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: #include <petscdmnetwork.h>
  5: #include <petscksp.h>

  7: /* The topology looks like:

  9:             (0)
 10:             /|\
 11:            / | \
 12:           /  |  \
 13:          R   R   V
 14:         /    |b3  \
 15:     b0 /    (3)    \ b1
 16:       /    /   \    R
 17:      /   R       R   \
 18:     /  /           \  \
 19:    / / b4        b5  \ \
 20:   //                   \\
 21: (1)--------- R -------- (2)
 22:              b2

 24:   Nodes:          (0), ... (3)
 25:   Branches:       b0, ... b5
 26:   Resistances:    R
 27:   Voltage source: V

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

 32: /*
 33:   Structures containing physical data of circuit.
 34:   Note that no topology is defined
 35: */

 37: typedef struct {
 38:   PetscInt    id;  /* node id */
 39:   PetscScalar inj; /* current injection (A) */
 40:   PetscBool   gr;  /* boundary node */
 41: } Node;

 43: typedef struct {
 44:   PetscInt    id;  /* branch id */
 45:   PetscScalar r;   /* resistance (ohms) */
 46:   PetscScalar bat; /* battery (V) */
 47: } Branch;

 49: /*
 50:   read_data: this routine fills data structures with problem data.
 51:   This can be substituted by an external parser.
 52: */

 54: PetscErrorCode read_data(PetscInt *pnnode,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist)
 55: {
 56:   PetscInt          nnode, nbranch, i;
 57:   Branch            *branch;
 58:   Node              *node;
 59:   PetscInt          *edgelist;

 62:   nnode   = 4;
 63:   nbranch = 6;

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

 67:   for (i = 0; i < nnode; i++) {
 68:     node[i].id  = i;
 69:     node[i].inj = 0;
 70:     node[i].gr = PETSC_FALSE;
 71:   }

 73:   for (i = 0; i < nbranch; i++) {
 74:     branch[i].id  = i;
 75:     branch[i].r   = 1.0;
 76:     branch[i].bat = 0;
 77:   }

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

 89:   /*
 90:     edgelist is an array of len = 2*nbranch that defines the
 91:     topology of the network. For branch i we have:
 92:       edgelist[2*i]     = from node
 93:       edgelist[2*i + 1] = to node.
 94:   */
 95:   PetscCalloc1(2*nbranch, &edgelist);

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

128:   /* assign pointers */
129:   *pnnode    = nnode;
130:   *pnbranch  = nbranch;
131:   *pedgelist = edgelist;
132:   *pbranch   = branch;
133:   *pnode     = node;
134:   return 0;
135: }

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

147:   MatZeroEntries(A);

149:   VecSet(b,0.0);
150:   VecGetArray(b,&barr);

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

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

165:     DMNetworkGetConnectedVertices(dmnetwork,e,&cone);
166:     DMNetworkGetLocalVecOffset(dmnetwork,cone[0],ALL_COMPONENTS,&lofst_fr);
167:     DMNetworkGetLocalVecOffset(dmnetwork,cone[1],ALL_COMPONENTS,&lofst_to);

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

172:     /* set Branch equation */
173:     row[0] = lofst;
174:     col[0] = lofst;     val[0] =  1./branch->r;
175:     col[1] = lofst_to;  val[1] =  1;
176:     col[2] = lofst_fr;  val[2] = -1;
177:     MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);

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

182:     /* from node */
183:     if (!node->gr) { /* not a boundary node */
184:       row[0] = lofst_fr;
185:       col[0] = lofst;   val[0] = -1;
186:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
187:     }

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

192:     if (!node->gr) { /* not a boundary node */
193:       row[0] = lofst_to;
194:       col[0] = lofst;   val[0] = 1;
195:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
196:     }
197:   }

199:   /* set rhs b for Node equation */
200:   DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
201:   for (v = vStart; v < vEnd; v++) {
202:     DMNetworkIsGhostVertex(dmnetwork,v,&ghost);
203:     if (!ghost) {
204:       DMNetworkGetComponent(dmnetwork,v,0,NULL,(void**)&node,NULL);
205:       DMNetworkGetLocalVecOffset(dmnetwork,v,ALL_COMPONENTS,&lofst);

207:       if (node->gr) { /* a boundary node */
208:         row[0] = lofst;
209:         col[0] = lofst;   val[0] = 1;
210:         MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
211:       } else {       /* not a boundary node */
212:         barr[lofst] += node->inj;
213:       }
214:     }
215:   }

217:   VecRestoreArray(b,&barr);

219:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
220:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
221:   return 0;
222: }

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

238:   PetscInitialize(&argc,&argv,(char*)0,help);
239:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
240:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

242:   /* "Read" data only for processor 0 */
243:   if (rank == 0) {
244:     read_data(&nnode, &nbranch, &node, &branch, &edgelist);
245:   }

247:   DMNetworkCreate(PETSC_COMM_WORLD,&dmnetwork);
248:   DMNetworkRegisterComponent(dmnetwork,"nstr",sizeof(Node),&componentkey[0]);
249:   DMNetworkRegisterComponent(dmnetwork,"bsrt",sizeof(Branch),&componentkey[1]);

251:   /* Set local number of nodes/edges, add edge connectivity */
252:   nE[0] = nbranch;
253:   DMNetworkSetNumSubNetworks(dmnetwork,PETSC_DECIDE,1);
254:   DMNetworkAddSubnetwork(dmnetwork,"",nE[0],edgelist,NULL);

256:   /* Set up the network layout */
257:   DMNetworkLayoutSetUp(dmnetwork);

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

266:     DMNetworkGetVertexRange(dmnetwork,&vStart,&vEnd);
267:     for (i = vStart; i < vEnd; i++) {
268:       DMNetworkAddComponent(dmnetwork,i,componentkey[0],&node[i-vStart],1);
269:     }
270:   }

272:   /* Network partitioning and distribution of data */
273:   DMSetUp(dmnetwork);
274:   DMNetworkDistribute(&dmnetwork,0);

276:   /* We do not use these data structures anymore since they have been copied to dmnetwork */
277:   if (rank == 0) {
278:     PetscFree(edgelist);
279:     PetscFree2(node,branch);
280:   }

282:   /* Create vectors and matrix */
283:   DMCreateGlobalVector(dmnetwork,&x);
284:   VecDuplicate(x,&b);
285:   DMCreateMatrix(dmnetwork,&A);

287:   /* Assembly system of equations */
288:   FormOperator(dmnetwork,A,b);

290:   /* Solve linear system: A x = b */
291:   KSPCreate(PETSC_COMM_WORLD, &ksp);
292:   KSPSetOperators(ksp, A, A);
293:   KSPSetFromOptions(ksp);
294:   KSPSolve(ksp, b, x);
295:   VecView(x, PETSC_VIEWER_STDOUT_WORLD);

297:   /* Free work space */
298:   VecDestroy(&x);
299:   VecDestroy(&b);
300:   MatDestroy(&A);
301:   KSPDestroy(&ksp);
302:   DMDestroy(&dmnetwork);
303:   PetscFinalize();
304:   return 0;
305: }

307: /*TEST

309:    build:
310:       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

312:    test:
313:       args: -ksp_monitor_short

315:    test:
316:       suffix: 2
317:       nsize: 2
318:       args: -petscpartitioner_type simple -ksp_converged_reason

320: TEST*/