Actual source code: ex1_nest.c

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

  3: #include <petscdmnetwork.h>
  4: #include <petscksp.h>

  6: /* The topology looks like:

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

 23:   Nodes:          (1), ... (4)
 24:   Branches:       b1, ... b6
 25:   Resistances:    R
 26:   Voltage source: V

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

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

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

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

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

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

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

 64:   PetscCalloc1(nnode,&node);
 65:   PetscCalloc1(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 would have that:
 92:       edgelist[2*i]     = from node
 93:       edgelist[2*i + 1] = to node
 94:   */

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

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

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

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

156:   DMGetLocalVector(networkdm,&localb);
157:   VecSet(b,0.0);
158:   VecSet(localb,0.0);

160:   VecGetArray(localb,&barr);

162:   /* Get nested submatrices */
163:   MatNestGetSubMats(A,NULL,NULL,&mats);
164:   e11 = mats[0][0];  /* edges */
165:   c12 = mats[0][1];  /* couplings */
166:   c21 = mats[1][0];  /* couplings */
167:   v22 = mats[1][1];  /* vertices */

169:   /* Get vertices and edge range */
170:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
171:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

173:   for (e = 0; e < eEnd; e++) {
174:     DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
175:     DMNetworkGetEdgeOffset(networkdm,e,&lofst);

177:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
178:     DMNetworkGetVertexOffset(networkdm,cone[0],&lofst_fr);
179:     DMNetworkGetVertexOffset(networkdm,cone[1],&lofst_to);

181:     /* These are edge-edge and go to e11 */
182:     row[0] = lofst;
183:     col[0] = lofst;     val[0] =  1;
184:     MatSetValuesLocal(e11,1,row,1,col,val,INSERT_VALUES);

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

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

195:     if (!node->gr) {
196:       row[0] = lofst_fr;
197:       col[0] = lofst;   val[0] =  1;
198:       MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
199:     }

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

204:     if (!node->gr) {
205:       row[0] = lofst_to;
206:       col[0] = lofst;   val[0] =  -1;
207:       MatSetValuesLocal(c21,1,row,1,col,val,INSERT_VALUES);
208:     }

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

215:   for (v = vStart; v < vEnd; v++) {
216:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
217:     if (!ghost) {
218:       DMNetworkGetComponent(networkdm,v,0,&key,(void**)&node,NULL);
219:       DMNetworkGetVertexOffset(networkdm,v,&lofst);

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

233:   VecRestoreArray(localb,&barr);

235:   DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
236:   DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
237:   DMRestoreLocalVector(networkdm,&localb);

239:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
240:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
241:   return 0;
242: }

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

258:   PetscInitialize(&argc,&argv,(char*)0,help);
259:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
260:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

262:   /* "read" data only for processor 0 */
263:   if (rank == 0) {
264:     read_data(&nnode, &nbranch, &node, &branch, &edgelist);
265:   }

267:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
268:   DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
269:   DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);

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

275:   /* Set up the network layout */
276:   DMNetworkLayoutSetUp(networkdm);

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

285:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
286:     for (i = vStart; i < vEnd; i++) {
287:       DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
288:     }
289:   }

291:   /* Network partitioning and distribution of data */
292:   DMSetUp(networkdm);
293:   DMNetworkDistribute(&networkdm,0);

295:   DMNetworkAssembleGraphStructures(networkdm);

297:   /* We don't use these data structures anymore since they have been copied to networkdm */
298:   if (rank == 0) {
299:     PetscFree(edgelist);
300:     PetscFree(node);
301:     PetscFree(branch);
302:   }

304:   DMCreateGlobalVector(networkdm,&x);
305:   VecDuplicate(x,&b);

307:   DMSetMatType(networkdm, MATNEST);
308:   DMCreateMatrix(networkdm,&A);

310:   /* Assembly system of equations */
311:   FormOperator(networkdm,A,b);

313:   KSPCreate(PETSC_COMM_WORLD, &ksp);
314:   KSPSetOperators(ksp, A, A);
315:   KSPSetFromOptions(ksp);
316:   KSPSolve(ksp, b, x);
317:   VecView(x, 0);

319:   VecDestroy(&x);
320:   VecDestroy(&b);
321:   MatDestroy(&A);
322:   KSPDestroy(&ksp);
323:   DMDestroy(&networkdm);
324:   PetscFinalize();
325:   return 0;
326: }

328: /*TEST

330:    build:
331:       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

333:    test:
334:       args: -ksp_converged_reason

336:    test:
337:       suffix: 2
338:       nsize: 2
339:       args: -petscpartitioner_type simple -ksp_converged_reason

341: TEST*/