Actual source code: ex1.c

  1: static char help[] = "This example demonstrates the use of DMNetwork interface for solving a simple electric circuit. \n\
  2:                       The example can be found in p.150 of 'Strang, Gilbert. Computational Science and Engineering. Wellesley, MA'.\n\n";

  4: /* T
  5:   Concepts: DMNetwork
  6:   Concepts: KSP
  7: */

  9: #include <petscdmnetwork.h>
 10: #include <petscksp.h>

 12: /* The topology looks like:

 14:             (0)
 15:             /|\
 16:            / | \
 17:           /  |  \
 18:          R   R   V
 19:         /    |b3  \
 20:     b0 /    (3)    \ b1
 21:       /    /   \    R
 22:      /   R       R   \
 23:     /  /           \  \
 24:    / / b4        b5  \ \
 25:   //                   \\
 26: (1)--------- R -------- (2)
 27:              b2

 29:   Nodes:          (0), ... (3)
 30:   Branches:       b0, ... b5
 31:   Resistances:    R
 32:   Voltage source: V

 34:   Additionally, there is a current source from (1) to (0).
 35: */

 37: /*
 38:   Structures containing physical data of circuit.
 39:   Note that no topology is defined
 40: */

 42: typedef struct {
 43:   PetscInt    id;  /* node id */
 44:   PetscScalar inj; /* current injection (A) */
 45:   PetscBool   gr;  /* boundary node */
 46: } Node;

 48: typedef struct {
 49:   PetscInt    id;  /* branch id */
 50:   PetscScalar r;   /* resistance (ohms) */
 51:   PetscScalar bat; /* battery (V) */
 52: } Branch;

 54: /*
 55:   read_data: this routine fills data structures with problem data.
 56:   This can be substituted by an external parser.
 57: */

 59: PetscErrorCode read_data(PetscInt *pnnode,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist)
 60: {
 61:   PetscInt          nnode, nbranch, i;
 62:   Branch            *branch;
 63:   Node              *node;
 64:   PetscInt          *edgelist;

 67:   nnode   = 4;
 68:   nbranch = 6;

 70:   PetscCalloc2(nnode,&node,nbranch,&branch);

 72:   for (i = 0; i < nnode; i++) {
 73:     node[i].id  = i;
 74:     node[i].inj = 0;
 75:     node[i].gr = PETSC_FALSE;
 76:   }

 78:   for (i = 0; i < nbranch; i++) {
 79:     branch[i].id  = i;
 80:     branch[i].r   = 1.0;
 81:     branch[i].bat = 0;
 82:   }

 84:   /*
 85:     Branch 1 contains a voltage source of 12.0 V
 86:     From node 0 to 1 there exists a current source of 4.0 A
 87:     Node 3 is grounded, hence the voltage is 0.
 88:   */
 89:   branch[1].bat = 12.0;
 90:   node[0].inj   = -4.0;
 91:   node[1].inj   =  4.0;
 92:   node[3].gr    =  PETSC_TRUE;

 94:   /*
 95:     edgelist is an array of len = 2*nbranch that defines the
 96:     topology of the network. For branch i we have:
 97:       edgelist[2*i]     = from node
 98:       edgelist[2*i + 1] = to node.
 99:   */
100:   PetscCalloc1(2*nbranch, &edgelist);

102:   for (i = 0; i < nbranch; i++) {
103:     switch (i) {
104:       case 0:
105:         edgelist[2*i]     = 0;
106:         edgelist[2*i + 1] = 1;
107:         break;
108:       case 1:
109:         edgelist[2*i]     = 0;
110:         edgelist[2*i + 1] = 2;
111:         break;
112:       case 2:
113:         edgelist[2*i]     = 1;
114:         edgelist[2*i + 1] = 2;
115:         break;
116:       case 3:
117:         edgelist[2*i]     = 0;
118:         edgelist[2*i + 1] = 3;
119:         break;
120:       case 4:
121:         edgelist[2*i]     = 1;
122:         edgelist[2*i + 1] = 3;
123:         break;
124:       case 5:
125:         edgelist[2*i]     = 2;
126:         edgelist[2*i + 1] = 3;
127:         break;
128:       default:
129:         break;
130:     }
131:   }

133:   /* assign pointers */
134:   *pnnode    = nnode;
135:   *pnbranch  = nbranch;
136:   *pedgelist = edgelist;
137:   *pbranch   = branch;
138:   *pnode     = node;
139:   return 0;
140: }

142: PetscErrorCode FormOperator(DM dmnetwork,Mat A,Vec b)
143: {
144:   Branch            *branch;
145:   Node              *node;
146:   PetscInt          e,v,vStart,vEnd,eStart, eEnd;
147:   PetscInt          lofst,lofst_to,lofst_fr,row[2],col[6];
148:   PetscBool         ghost;
149:   const PetscInt    *cone;
150:   PetscScalar       *barr,val[6];

152:   MatZeroEntries(A);

154:   VecSet(b,0.0);
155:   VecGetArray(b,&barr);

157:   /*
158:     We define the current i as an "edge characteristic" and the voltage v as a "vertex characteristic".
159:     We iterate the list of edges and vertices, query the associated voltages and currents
160:     and use them to write the Kirchoff equations:

162:     Branch equations: i/r + v_to - v_from     = v_source (battery)
163:     Node equations:   sum(i_to) - sum(i_from) = i_source
164:    */
165:   DMNetworkGetEdgeRange(dmnetwork,&eStart,&eEnd);
166:   for (e = 0; e < eEnd; e++) {
167:     DMNetworkGetComponent(dmnetwork,e,0,NULL,(void**)&branch,NULL);
168:     DMNetworkGetLocalVecOffset(dmnetwork,e,ALL_COMPONENTS,&lofst);

170:     DMNetworkGetConnectedVertices(dmnetwork,e,&cone);
171:     DMNetworkGetLocalVecOffset(dmnetwork,cone[0],ALL_COMPONENTS,&lofst_fr);
172:     DMNetworkGetLocalVecOffset(dmnetwork,cone[1],ALL_COMPONENTS,&lofst_to);

174:     /* set rhs b for Branch equation */
175:     barr[lofst] = branch->bat;

177:     /* set Branch equation */
178:     row[0] = lofst;
179:     col[0] = lofst;     val[0] =  1./branch->r;
180:     col[1] = lofst_to;  val[1] =  1;
181:     col[2] = lofst_fr;  val[2] = -1;
182:     MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);

184:     /* set Node equation */
185:     DMNetworkGetComponent(dmnetwork,cone[0],0,NULL,(void**)&node,NULL);

187:     /* from node */
188:     if (!node->gr) { /* not a boundary node */
189:       row[0] = lofst_fr;
190:       col[0] = lofst;   val[0] = -1;
191:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
192:     }

194:     /* to node */
195:     DMNetworkGetComponent(dmnetwork,cone[1],0,NULL,(void**)&node,NULL);

197:     if (!node->gr) { /* not a boundary node */
198:       row[0] = lofst_to;
199:       col[0] = lofst;   val[0] = 1;
200:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
201:     }
202:   }

204:   /* set rhs b for Node equation */
205:   DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
206:   for (v = vStart; v < vEnd; v++) {
207:     DMNetworkIsGhostVertex(dmnetwork,v,&ghost);
208:     if (!ghost) {
209:       DMNetworkGetComponent(dmnetwork,v,0,NULL,(void**)&node,NULL);
210:       DMNetworkGetLocalVecOffset(dmnetwork,v,ALL_COMPONENTS,&lofst);

212:       if (node->gr) { /* a boundary node */
213:         row[0] = lofst;
214:         col[0] = lofst;   val[0] = 1;
215:         MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
216:       } else {       /* not a boundary node */
217:         barr[lofst] += node->inj;
218:       }
219:     }
220:   }

222:   VecRestoreArray(b,&barr);

224:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
225:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
226:   return 0;
227: }

229: int main(int argc,char ** argv)
230: {
231:   PetscInt          i, nnode = 0, nbranch = 0, eStart, eEnd, vStart, vEnd;
232:   PetscMPIInt       size, rank;
233:   DM                dmnetwork;
234:   Vec               x, b;
235:   Mat               A;
236:   KSP               ksp;
237:   PetscInt          *edgelist = NULL;
238:   PetscInt          componentkey[2];
239:   Node              *node;
240:   Branch            *branch;
241:   PetscInt          nE[1];

243:   PetscInitialize(&argc,&argv,(char*)0,help);
244:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
245:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

247:   /* "Read" data only for processor 0 */
248:   if (rank == 0) {
249:     read_data(&nnode, &nbranch, &node, &branch, &edgelist);
250:   }

252:   DMNetworkCreate(PETSC_COMM_WORLD,&dmnetwork);
253:   DMNetworkRegisterComponent(dmnetwork,"nstr",sizeof(Node),&componentkey[0]);
254:   DMNetworkRegisterComponent(dmnetwork,"bsrt",sizeof(Branch),&componentkey[1]);

256:   /* Set local number of nodes/edges, add edge connectivity */
257:   nE[0] = nbranch;
258:   DMNetworkSetNumSubNetworks(dmnetwork,PETSC_DECIDE,1);
259:   DMNetworkAddSubnetwork(dmnetwork,"",nE[0],edgelist,NULL);

261:   /* Set up the network layout */
262:   DMNetworkLayoutSetUp(dmnetwork);

264:   /* Add network components (physical parameters of nodes and branches) and num of variables */
265:   if (rank == 0) {
266:     DMNetworkGetEdgeRange(dmnetwork,&eStart,&eEnd);
267:     for (i = eStart; i < eEnd; i++) {
268:       DMNetworkAddComponent(dmnetwork,i,componentkey[1],&branch[i-eStart],1);
269:     }

271:     DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
272:     for (i = vStart; i < vEnd; i++) {
273:       DMNetworkAddComponent(dmnetwork,i,componentkey[0],&node[i-vStart],1);
274:     }
275:   }

277:   /* Network partitioning and distribution of data */
278:   DMSetUp(dmnetwork);
279:   DMNetworkDistribute(&dmnetwork,0);

281:   /* We do not use these data structures anymore since they have been copied to dmnetwork */
282:   if (rank == 0) {
283:     PetscFree(edgelist);
284:     PetscFree2(node,branch);
285:   }

287:   /* Create vectors and matrix */
288:   DMCreateGlobalVector(dmnetwork,&x);
289:   VecDuplicate(x,&b);
290:   DMCreateMatrix(dmnetwork,&A);

292:   /* Assembly system of equations */
293:   FormOperator(dmnetwork,A,b);

295:   /* Solve linear system: A x = b */
296:   KSPCreate(PETSC_COMM_WORLD, &ksp);
297:   KSPSetOperators(ksp, A, A);
298:   KSPSetFromOptions(ksp);
299:   KSPSolve(ksp, b, x);
300:   VecView(x, PETSC_VIEWER_STDOUT_WORLD);

302:   /* Free work space */
303:   VecDestroy(&x);
304:   VecDestroy(&b);
305:   MatDestroy(&A);
306:   KSPDestroy(&ksp);
307:   DMDestroy(&dmnetwork);
308:   PetscFinalize();
309:   return 0;
310: }

312: /*TEST

314:    build:
315:       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

317:    test:
318:       args: -ksp_monitor_short

320:    test:
321:       suffix: 2
322:       nsize: 2
323:       args: -petscpartitioner_type simple -ksp_converged_reason

325: TEST*/