Actual source code: ex2.c

  1: static char help[] = "This example is based on ex1.c, but generates a random network of chosen sizes and parameters. \n\
  2:   Usage: -n determines number of nodes. The nonnegative seed can be specified with the flag -seed, otherwise the program generates a random seed.\n\n";

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

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

 13: typedef struct {
 14:   PetscInt    id; /* node id */
 15:   PetscScalar inj; /* current injection (A) */
 16:   PetscBool   gr; /* grounded ? */
 17: } Node;

 19: typedef struct {
 20:   PetscInt    id;  /* branch id */
 21:   PetscScalar r;   /* resistance (ohms) */
 22:   PetscScalar bat; /* battery (V) */
 23: } Branch;

 25: typedef struct Edge {
 26:   PetscInt      n;
 27:   PetscInt      i;
 28:   PetscInt      j;
 29:   struct Edge   *next;
 30: } Edge;

 32: PetscReal findDistance(PetscReal x1, PetscReal x2, PetscReal y1, PetscReal y2)
 33: {
 34:   return PetscSqrtReal(PetscPowReal(x2-x1,2.0) + PetscPowReal(y2-y1,2.0));
 35: }

 37: /*
 38:   The algorithm for network formation is based on the paper:
 39:   Routing of Multipoint Connections, Bernard M. Waxman. 1988
 40: */

 42: PetscErrorCode random_network(PetscInt nvertex,PetscInt *pnbranch,Node **pnode,Branch **pbranch,PetscInt **pedgelist,PetscInt seed)
 43: {
 44:   PetscInt       i, j, nedges = 0;
 45:   PetscInt       *edgelist;
 46:   PetscInt       nbat, ncurr, fr, to;
 47:   PetscReal      *x, *y, value, xmax = 10.0; /* generate points in square */
 48:   PetscReal      maxdist = 0.0, dist, alpha, beta, prob;
 49:   PetscRandom    rnd;
 50:   Branch         *branch;
 51:   Node           *node;
 52:   Edge           *head = NULL, *nnew= NULL, *aux= NULL;

 55:   PetscRandomCreate(PETSC_COMM_SELF,&rnd);
 56:   PetscRandomSetFromOptions(rnd);

 58:   PetscRandomSetSeed(rnd, seed);
 59:   PetscRandomSeed(rnd);

 61:   /* These parameters might be modified for experimentation */
 62:   nbat  = (PetscInt)(0.1*nvertex);
 63:   ncurr = (PetscInt)(0.1*nvertex);
 64:   alpha = 0.6;
 65:   beta  = 0.2;

 67:   PetscMalloc2(nvertex,&x,nvertex,&y);

 69:   PetscRandomSetInterval(rnd,0.0,xmax);
 70:   for (i=0; i<nvertex; i++) {
 71:     PetscRandomGetValueReal(rnd,&x[i]);
 72:     PetscRandomGetValueReal(rnd,&y[i]);
 73:   }

 75:   /* find maximum distance */
 76:   for (i=0; i<nvertex; i++) {
 77:     for (j=0; j<nvertex; j++) {
 78:       dist = findDistance(x[i],x[j],y[i],y[j]);
 79:       if (dist >= maxdist) maxdist = dist;
 80:     }
 81:   }

 83:   PetscRandomSetInterval(rnd,0.0,1.0);
 84:   for (i=0; i<nvertex; i++) {
 85:     for (j=0; j<nvertex; j++) {
 86:       if (j != i) {
 87:         dist = findDistance(x[i],x[j],y[i],y[j]);
 88:         prob = beta*PetscExpScalar(-dist/(maxdist*alpha));
 89:         PetscRandomGetValueReal(rnd,&value);
 90:         if (value <= prob) {
 91:           PetscMalloc1(1,&nnew);
 92:           if (head == NULL) {
 93:             head       = nnew;
 94:             head->next = NULL;
 95:             head->n    = nedges;
 96:             head->i    = i;
 97:             head->j    = j;
 98:           } else {
 99:             aux = head;
100:             head = nnew;
101:             head->n    = nedges;
102:             head->next = aux;
103:             head->i    = i;
104:             head->j    = j;
105:           }
106:           nedges += 1;
107:         }
108:       }
109:     }
110:   }

112:   PetscMalloc1(2*nedges,&edgelist);

114:   for (aux = head; aux; aux = aux->next) {
115:     edgelist[(aux->n)*2]     = aux->i;
116:     edgelist[(aux->n)*2 + 1] = aux->j;
117:   }

119:   aux = head;
120:   while (aux != NULL) {
121:     nnew = aux;
122:     aux = aux->next;
123:     PetscFree(nnew);
124:   }

126:   PetscCalloc2(nvertex,&node,nedges,&branch);

128:   for (i = 0; i < nvertex; i++) {
129:     node[i].id  = i;
130:     node[i].inj = 0;
131:     node[i].gr = PETSC_FALSE;
132:   }

134:   for (i = 0; i < nedges; i++) {
135:     branch[i].id  = i;
136:     branch[i].r   = 1.0;
137:     branch[i].bat = 0;
138:   }

140:   /* Chose random node as ground voltage */
141:   PetscRandomSetInterval(rnd,0.0,nvertex);
142:   PetscRandomGetValueReal(rnd,&value);
143:   node[(int)value].gr = PETSC_TRUE;

145:   /* Create random current and battery injectionsa */
146:   for (i=0; i<ncurr; i++) {
147:     PetscRandomSetInterval(rnd,0.0,nvertex);
148:     PetscRandomGetValueReal(rnd,&value);
149:     fr   = edgelist[(int)value*2];
150:     to   = edgelist[(int)value*2 + 1];
151:     node[fr].inj += 1.0;
152:     node[to].inj -= 1.0;
153:   }

155:   for (i=0; i<nbat; i++) {
156:     PetscRandomSetInterval(rnd,0.0,nedges);
157:     PetscRandomGetValueReal(rnd,&value);
158:     branch[(int)value].bat += 1.0;
159:   }

161:   PetscFree2(x,y);
162:   PetscRandomDestroy(&rnd);

164:   /* assign pointers */
165:   *pnbranch  = nedges;
166:   *pedgelist = edgelist;
167:   *pbranch   = branch;
168:   *pnode     = node;
169:   return 0;
170: }

172: PetscErrorCode FormOperator(DM networkdm,Mat A,Vec b)
173: {
174:   Vec               localb;
175:   Branch            *branch;
176:   Node              *node;
177:   PetscInt          e,v,vStart,vEnd,eStart, eEnd;
178:   PetscInt          lofst,lofst_to,lofst_fr,row[2],col[6];
179:   PetscBool         ghost;
180:   const PetscInt    *cone;
181:   PetscScalar       *barr,val[6];

183:   DMGetLocalVector(networkdm,&localb);
184:   VecSet(b,0.0);
185:   VecSet(localb,0.0);
186:   MatZeroEntries(A);

188:   VecGetArray(localb,&barr);

190:   /*
191:     We can define the current as a "edge characteristic" and the voltage
192:     and the voltage as a "vertex characteristic". With that, we can iterate
193:     the list of edges and vertices, query the associated voltages and currents
194:     and use them to write the Kirchoff equations.
195:   */

197:   /* Branch equations: i/r + uj - ui = battery */
198:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
199:   for (e = 0; e < eEnd; e++) {
200:     DMNetworkGetComponent(networkdm,e,0,NULL,(void**)&branch,NULL);
201:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);

203:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
204:     DMNetworkGetLocalVecOffset(networkdm,cone[0],ALL_COMPONENTS,&lofst_fr);
205:     DMNetworkGetLocalVecOffset(networkdm,cone[1],ALL_COMPONENTS,&lofst_to);

207:     barr[lofst] = branch->bat;

209:     row[0] = lofst;
210:     col[0] = lofst;     val[0] =  1;
211:     col[1] = lofst_to;  val[1] =  1;
212:     col[2] = lofst_fr;  val[2] = -1;
213:     MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);

215:     /* from node */
216:     DMNetworkGetComponent(networkdm,cone[0],0,NULL,(void**)&node,NULL);

218:     if (!node->gr) {
219:       row[0] = lofst_fr;
220:       col[0] = lofst;   val[0] =  1;
221:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
222:     }

224:     /* to node */
225:     DMNetworkGetComponent(networkdm,cone[1],0,NULL,(void**)&node,NULL);

227:     if (!node->gr) {
228:       row[0] = lofst_to;
229:       col[0] = lofst;   val[0] =  -1;
230:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
231:     }
232:   }

234:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
235:   for (v = vStart; v < vEnd; v++) {
236:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
237:     if (!ghost) {
238:       DMNetworkGetComponent(networkdm,v,0,NULL,(void**)&node,NULL);
239:       DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);

241:       if (node->gr) {
242:         row[0] = lofst;
243:         col[0] = lofst;   val[0] =  1;
244:         MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
245:       } else {
246:         barr[lofst] -= node->inj;
247:       }
248:     }
249:   }

251:   VecRestoreArray(localb,&barr);

253:   DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
254:   DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
255:   DMRestoreLocalVector(networkdm,&localb);

257:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
258:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
259:   return 0;
260: }

262: int main(int argc,char ** argv)
263: {
264:   PetscInt          i, nbranch = 0, eStart, eEnd, vStart, vEnd;
265:   PetscInt          seed = 0, nnode = 0;
266:   PetscMPIInt       size, rank;
267:   DM                networkdm;
268:   Vec               x, b;
269:   Mat               A;
270:   KSP               ksp;
271:   PetscInt          *edgelist = NULL;
272:   PetscInt          componentkey[2];
273:   Node              *node;
274:   Branch            *branch;
275: #if defined(PETSC_USE_LOG)
276:   PetscLogStage stage[3];
277: #endif

279:   PetscInitialize(&argc,&argv,(char*)0,help);
280:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
281:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

283:   PetscOptionsGetInt(NULL,NULL,"-seed",&seed,NULL);

285:   PetscLogStageRegister("Network Creation", &stage[0]);
286:   PetscLogStageRegister("DMNetwork data structures", &stage[1]);
287:   PetscLogStageRegister("KSP", &stage[2]);

289:   PetscLogStagePush(stage[0]);
290:   /* "read" data only for processor 0 */
291:   if (rank == 0) {
292:     nnode = 100;
293:     PetscOptionsGetInt(NULL,NULL,"-n",&nnode,NULL);
294:     random_network(nnode, &nbranch, &node, &branch, &edgelist, seed);
295:   }
296:   PetscLogStagePop();

298:   PetscLogStagePush(stage[1]);
299:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
300:   DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
301:   DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);

303:   /* Set number of nodes/edges and edge connectivity */
304:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
305:   DMNetworkAddSubnetwork(networkdm,"",nbranch,edgelist,NULL);

307:   /* Set up the network layout */
308:   DMNetworkLayoutSetUp(networkdm);

310:   /* Add network components (physical parameters of nodes and branches) and num of variables */
311:   if (rank == 0) {
312:     DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
313:     for (i = eStart; i < eEnd; i++) {
314:       DMNetworkAddComponent(networkdm,i,componentkey[1],&branch[i-eStart],1);
315:     }

317:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
318:     for (i = vStart; i < vEnd; i++) {
319:       DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
320:     }
321:   }

323:   /* Network partitioning and distribution of data */
324:   DMSetUp(networkdm);
325:   DMNetworkDistribute(&networkdm,0);
326:   DMNetworkAssembleGraphStructures(networkdm);

328:   /* We don't use these data structures anymore since they have been copied to networkdm */
329:   if (rank == 0) {
330:     PetscFree(edgelist);
331:     PetscFree2(node,branch);
332:   }

334:   /* Create vectors and matrix */
335:   DMCreateGlobalVector(networkdm,&x);
336:   VecDuplicate(x,&b);
337:   DMCreateMatrix(networkdm,&A);

339:   PetscLogStagePop();

341:   PetscLogStagePush(stage[2]);
342:   /* Assembly system of equations */
343:   FormOperator(networkdm,A,b);

345:   /* Solve linear system: A x = b */
346:   KSPCreate(PETSC_COMM_WORLD, &ksp);
347:   KSPSetOperators(ksp, A, A);
348:   KSPSetFromOptions(ksp);
349:   KSPSolve(ksp, b, x);

351:   PetscLogStagePop();

353:   /* Free work space */
354:   VecDestroy(&x);
355:   VecDestroy(&b);
356:   MatDestroy(&A);
357:   KSPDestroy(&ksp);
358:   DMDestroy(&networkdm);
359:   PetscFinalize();
360:   return 0;
361: }

363: /*TEST

365:    build:
366:       requires: !single double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

368:    test:
369:       args: -ksp_converged_reason

371:    test:
372:       suffix: 2
373:       nsize: 2
374:       args: -petscpartitioner_type simple -pc_type asm -sub_pc_type ilu -ksp_converged_reason

376:    test:
377:       suffix: 3
378:       nsize: 4
379:       args: -petscpartitioner_type simple -pc_type asm -sub_pc_type lu -sub_pc_factor_shift_type nonzero -ksp_converged_reason

381:    test:
382:       suffix: graphindex
383:       args: -n 20 -vertex_global_section_view -edge_global_section_view

385:    test:
386:       suffix: graphindex_2
387:       nsize: 2
388:       args: -petscpartitioner_type simple -n 20 -vertex_global_section_view -edge_global_section_view

390: TEST*/