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:   PetscErrorCode    ierr;
 61:   PetscInt          nnode, nbranch, i;
 62:   Branch            *branch;
 63:   Node              *node;
 64:   PetscInt          *edgelist;

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

 70:   PetscCalloc1(nnode,&node);
 71:   PetscCalloc1(nbranch,&branch);

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

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

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

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

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

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

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

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

164:   DMGetLocalVector(networkdm,&localb);
165:   VecSet(b,0.0);
166:   VecSet(localb,0.0);

168:   VecGetArray(localb,&barr);

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

177:   /* Get vertices and edge range */
178:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
179:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);

181:   for (e = 0; e < eEnd; e++) {
182:     DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
183:     DMNetworkGetEdgeOffset(networkdm,e,&lofst);

185:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
186:     DMNetworkGetVertexOffset(networkdm,cone[0],&lofst_fr);
187:     DMNetworkGetVertexOffset(networkdm,cone[1],&lofst_to);

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

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

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

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

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

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

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

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

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

241:   VecRestoreArray(localb,&barr);

243:   DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
244:   DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
245:   DMRestoreLocalVector(networkdm,&localb);

247:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
248:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
249:   return(0);
250: }

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

267:   PetscInitialize(&argc,&argv,(char*)0,help);if (ierr) return ierr;
268:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
269:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

271:   /* "read" data only for processor 0 */
272:   if (rank == 0) {
273:     read_data(&nnode, &nbranch, &node, &branch, &edgelist);
274:   }

276:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
277:   DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
278:   DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);

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

284:   /* Set up the network layout */
285:   DMNetworkLayoutSetUp(networkdm);

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

294:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
295:     for (i = vStart; i < vEnd; i++) {
296:       DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
297:     }
298:   }

300:   /* Network partitioning and distribution of data */
301:   DMSetUp(networkdm);
302:   DMNetworkDistribute(&networkdm,0);

304:   DMNetworkAssembleGraphStructures(networkdm);

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

313:   DMCreateGlobalVector(networkdm,&x);
314:   VecDuplicate(x,&b);

316:   DMSetMatType(networkdm, MATNEST);
317:   DMCreateMatrix(networkdm,&A);

319:   /* Assembly system of equations */
320:   FormOperator(networkdm,A,b);

322:   KSPCreate(PETSC_COMM_WORLD, &ksp);
323:   KSPSetOperators(ksp, A, A);
324:   KSPSetFromOptions(ksp);
325:   KSPSolve(ksp, b, x);
326:   VecView(x, 0);

328:   VecDestroy(&x);
329:   VecDestroy(&b);
330:   MatDestroy(&A);
331:   KSPDestroy(&ksp);
332:   DMDestroy(&networkdm);
333:   PetscFinalize();
334:   return ierr;
335: }

337: /*TEST

339:    build:
340:       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

342:    test:
343:       args: -ksp_converged_reason

345:    test:
346:       suffix: 2
347:       nsize: 2
348:       args: -petscpartitioner_type simple -ksp_converged_reason

350: TEST*/