Actual source code: ex1_nest.c

  1: static char help[] = "This example is based on ex1 using MATNEST. \n";

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

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

 11: /* The topology looks like:

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

 28:   Nodes:          (1), ... (4)
 29:   Branches:       b1, ... b6
 30:   Resistances:    R
 31:   Voltage source: V

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

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

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

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

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

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

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

 69:   PetscCalloc1(nnode,&node);
 70:   PetscCalloc1(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 would have that:
 97:       edgelist[2*i]     = from node
 98:       edgelist[2*i + 1] = to node
 99:   */

101:   PetscCalloc1(2*nbranch, &edgelist);

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

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

143: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
144: {
145:   Vec               localb;
146:   Branch            *branch;
147:   Node              *node;
148:   PetscInt          e;
149:   PetscInt          v,vStart,vEnd;
150:   PetscInt          eStart, eEnd;
151:   PetscBool         ghost;
152:   const PetscInt    *cone;
153:   PetscScalar       *barr;
154:   PetscInt          lofst, lofst_to, lofst_fr;
155:   PetscInt          key;
156:   PetscInt          row[2],col[6];
157:   PetscScalar       val[6];
158:   Mat               e11, c12, c21, v22;
159:   Mat               **mats;

161:   DMGetLocalVector(networkdm,&localb);
162:   VecSet(b,0.0);
163:   VecSet(localb,0.0);

165:   VecGetArray(localb,&barr);

167:   /* Get nested submatrices */
168:   MatNestGetSubMats(A,NULL,NULL,&mats);
169:   e11 = mats[0][0];  /* edges */
170:   c12 = mats[0][1];  /* couplings */
171:   c21 = mats[1][0];  /* couplings */
172:   v22 = mats[1][1];  /* vertices */

174:   /* Get vertices and edge range */
175:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
176:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

178:   for (e = 0; e < eEnd; e++) {
179:     DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
180:     DMNetworkGetEdgeOffset(networkdm,e,&lofst);

182:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
183:     DMNetworkGetVertexOffset(networkdm,cone[0],&lofst_fr);
184:     DMNetworkGetVertexOffset(networkdm,cone[1],&lofst_to);

186:     /* These are edge-edge and go to e11 */
187:     row[0] = lofst;
188:     col[0] = lofst;     val[0] =  1;
189:     MatSetValuesLocal(e11,1,row,1,col,val,INSERT_VALUES);

191:     /* These are edge-vertex and go to c12 */
192:     col[0] = lofst_to;  val[0] =  1;
193:     col[1] = lofst_fr;  val[1] = -1;
194:     MatSetValuesLocal(c12,1,row,2,col,val,INSERT_VALUES);

196:     /* These are edge-vertex and go to c12 */
197:     /* from node */
198:     DMNetworkGetComponent(networkdm,cone[0],0,&key,(void**)&node,NULL);

200:     if (!node->gr) {
201:       row[0] = lofst_fr;
202:       col[0] = lofst;   val[0] =  1;
203:       MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
204:     }

206:     /* to node */
207:     DMNetworkGetComponent(networkdm,cone[1],0,&key,(void**)&node,NULL);

209:     if (!node->gr) {
210:       row[0] = lofst_to;
211:       col[0] = lofst;   val[0] =  -1;
212:       MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
213:     }

215:     /* TODO: this is not a nested vector. Need to implement nested vector */
216:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);
217:     barr[lofst] = branch->bat;
218:   }

220:   for (v = vStart; v < vEnd; v++) {
221:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
222:     if (!ghost) {
223:       DMNetworkGetComponent(networkdm,v,0,&key,(void**)&node,NULL);
224:       DMNetworkGetVertexOffset(networkdm,v,&lofst);

226:       if (node->gr) {
227:         row[0] = lofst;
228:         col[0] = lofst;   val[0] =  1;
229:         MatSetValuesLocal(v22,1,row,1,col,val,INSERT_VALUES);
230:       } else {
231:         /* TODO: this is not a nested vector. Need to implement nested vector */
232:         DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);
233:         barr[lofst] -= node->inj;
234:       }
235:     }
236:   }

238:   VecRestoreArray(localb,&barr);

240:   DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
241:   DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
242:   DMRestoreLocalVector(networkdm,&localb);

244:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
245:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
246:   return 0;
247: }

249: int main(int argc,char ** argv)
250: {
251:   PetscInt          i, nnode = 0, nbranch = 0;
252:   PetscInt          eStart, eEnd, vStart, vEnd;
253:   PetscMPIInt       size, rank;
254:   DM                networkdm;
255:   Vec               x, b;
256:   Mat               A;
257:   KSP               ksp;
258:   PetscInt          *edgelist = NULL;
259:   PetscInt          componentkey[2];
260:   Node              *node;
261:   Branch            *branch;

263:   PetscInitialize(&argc,&argv,(char*)0,help);
264:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
265:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

267:   /* "read" data only for processor 0 */
268:   if (rank == 0) {
269:     read_data(&nnode, &nbranch, &node, &branch, &edgelist);
270:   }

272:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
273:   DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
274:   DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);

276:   /* Set number of nodes/edges, add edge connectivity */
277:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
278:   DMNetworkAddSubnetwork(networkdm,"",nbranch,edgelist,NULL);

280:   /* Set up the network layout */
281:   DMNetworkLayoutSetUp(networkdm);

283:   /* Add network components (physical parameters of nodes and branches) and num of variables */
284:   if (rank == 0) {
285:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
286:     for (i = eStart; i < eEnd; i++) {
287:       DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
288:     }

290:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
291:     for (i = vStart; i < vEnd; i++) {
292:       DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
293:     }
294:   }

296:   /* Network partitioning and distribution of data */
297:   DMSetUp(networkdm);
298:   DMNetworkDistribute(&networkdm,0);

300:   DMNetworkAssembleGraphStructures(networkdm);

302:   /* We don't use these data structures anymore since they have been copied to networkdm */
303:   if (rank == 0) {
304:     PetscFree(edgelist);
305:     PetscFree(node);
306:     PetscFree(branch);
307:   }

309:   DMCreateGlobalVector(networkdm,&x);
310:   VecDuplicate(x,&b);

312:   DMSetMatType(networkdm, MATNEST);
313:   DMCreateMatrix(networkdm,&A);

315:   /* Assembly system of equations */
316:   FormOperator(networkdm,A,b);

318:   KSPCreate(PETSC_COMM_WORLD, &ksp);
319:   KSPSetOperators(ksp, A, A);
320:   KSPSetFromOptions(ksp);
321:   KSPSolve(ksp, b, x);
322:   VecView(x, 0);

324:   VecDestroy(&x);
325:   VecDestroy(&b);
326:   MatDestroy(&A);
327:   KSPDestroy(&ksp);
328:   DMDestroy(&networkdm);
329:   PetscFinalize();
330:   return 0;
331: }

333: /*TEST

335:    build:
336:       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

338:    test:
339:       args: -ksp_converged_reason

341:    test:
342:       suffix: 2
343:       nsize: 2
344:       args: -petscpartitioner_type simple -ksp_converged_reason

346: TEST*/