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

  8: typedef struct {
  9:   PetscInt    id; /* node id */
 10:   PetscScalar inj; /* current injection (A) */
 11:   PetscBool   gr; /* grounded ? */
 12: } Node;

 14: typedef struct {
 15:   PetscInt    id;  /* branch id */
 16:   PetscScalar r;   /* resistance (ohms) */
 17:   PetscScalar bat; /* battery (V) */
 18: } Branch;

 20: typedef struct Edge {
 21:   PetscInt      n;
 22:   PetscInt      i;
 23:   PetscInt      j;
 24:   struct Edge   *next;
 25: } Edge;

 27: PetscReal findDistance(PetscReal x1, PetscReal x2, PetscReal y1, PetscReal y2)
 28: {
 29:   return PetscSqrtReal(PetscPowReal(x2-x1,2.0) + PetscPowReal(y2-y1,2.0));
 30: }

 32: /*
 33:   The algorithm for network formation is based on the paper:
 34:   Routing of Multipoint Connections, Bernard M. Waxman. 1988
 35: */

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

 50:   PetscRandomCreate(PETSC_COMM_SELF,&rnd);
 51:   PetscRandomSetFromOptions(rnd);

 53:   PetscRandomSetSeed(rnd, seed);
 54:   PetscRandomSeed(rnd);

 56:   /* These parameters might be modified for experimentation */
 57:   nbat  = (PetscInt)(0.1*nvertex);
 58:   ncurr = (PetscInt)(0.1*nvertex);
 59:   alpha = 0.6;
 60:   beta  = 0.2;

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

 64:   PetscRandomSetInterval(rnd,0.0,xmax);
 65:   for (i=0; i<nvertex; i++) {
 66:     PetscRandomGetValueReal(rnd,&x[i]);
 67:     PetscRandomGetValueReal(rnd,&y[i]);
 68:   }

 70:   /* find maximum distance */
 71:   for (i=0; i<nvertex; i++) {
 72:     for (j=0; j<nvertex; j++) {
 73:       dist = findDistance(x[i],x[j],y[i],y[j]);
 74:       if (dist >= maxdist) maxdist = dist;
 75:     }
 76:   }

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

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

109:   for (aux = head; aux; aux = aux->next) {
110:     edgelist[(aux->n)*2]     = aux->i;
111:     edgelist[(aux->n)*2 + 1] = aux->j;
112:   }

114:   aux = head;
115:   while (aux != NULL) {
116:     nnew = aux;
117:     aux = aux->next;
118:     PetscFree(nnew);
119:   }

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

123:   for (i = 0; i < nvertex; i++) {
124:     node[i].id  = i;
125:     node[i].inj = 0;
126:     node[i].gr = PETSC_FALSE;
127:   }

129:   for (i = 0; i < nedges; i++) {
130:     branch[i].id  = i;
131:     branch[i].r   = 1.0;
132:     branch[i].bat = 0;
133:   }

135:   /* Chose random node as ground voltage */
136:   PetscRandomSetInterval(rnd,0.0,nvertex);
137:   PetscRandomGetValueReal(rnd,&value);
138:   node[(int)value].gr = PETSC_TRUE;

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

150:   for (i=0; i<nbat; i++) {
151:     PetscRandomSetInterval(rnd,0.0,nedges);
152:     PetscRandomGetValueReal(rnd,&value);
153:     branch[(int)value].bat += 1.0;
154:   }

156:   PetscFree2(x,y);
157:   PetscRandomDestroy(&rnd);

159:   /* assign pointers */
160:   *pnbranch  = nedges;
161:   *pedgelist = edgelist;
162:   *pbranch   = branch;
163:   *pnode     = node;
164:   return 0;
165: }

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

178:   DMGetLocalVector(networkdm,&localb);
179:   VecSet(b,0.0);
180:   VecSet(localb,0.0);
181:   MatZeroEntries(A);

183:   VecGetArray(localb,&barr);

185:   /*
186:     We can define the current as a "edge characteristic" and the voltage
187:     and the voltage as a "vertex characteristic". With that, we can iterate
188:     the list of edges and vertices, query the associated voltages and currents
189:     and use them to write the Kirchoff equations.
190:   */

192:   /* Branch equations: i/r + uj - ui = battery */
193:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
194:   for (e = 0; e < eEnd; e++) {
195:     DMNetworkGetComponent(networkdm,e,0,NULL,(void**)&branch,NULL);
196:     DMNetworkGetLocalVecOffset(networkdm,e,ALL_COMPONENTS,&lofst);

198:     DMNetworkGetConnectedVertices(networkdm,e,&cone);
199:     DMNetworkGetLocalVecOffset(networkdm,cone[0],ALL_COMPONENTS,&lofst_fr);
200:     DMNetworkGetLocalVecOffset(networkdm,cone[1],ALL_COMPONENTS,&lofst_to);

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

204:     row[0] = lofst;
205:     col[0] = lofst;     val[0] =  1;
206:     col[1] = lofst_to;  val[1] =  1;
207:     col[2] = lofst_fr;  val[2] = -1;
208:     MatSetValuesLocal(A,1,row,3,col,val,ADD_VALUES);

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

213:     if (!node->gr) {
214:       row[0] = lofst_fr;
215:       col[0] = lofst;   val[0] =  1;
216:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
217:     }

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

222:     if (!node->gr) {
223:       row[0] = lofst_to;
224:       col[0] = lofst;   val[0] =  -1;
225:       MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
226:     }
227:   }

229:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
230:   for (v = vStart; v < vEnd; v++) {
231:     DMNetworkIsGhostVertex(networkdm,v,&ghost);
232:     if (!ghost) {
233:       DMNetworkGetComponent(networkdm,v,0,NULL,(void**)&node,NULL);
234:       DMNetworkGetLocalVecOffset(networkdm,v,ALL_COMPONENTS,&lofst);

236:       if (node->gr) {
237:         row[0] = lofst;
238:         col[0] = lofst;   val[0] =  1;
239:         MatSetValuesLocal(A,1,row,1,col,val,ADD_VALUES);
240:       } else {
241:         barr[lofst] -= node->inj;
242:       }
243:     }
244:   }

246:   VecRestoreArray(localb,&barr);

248:   DMLocalToGlobalBegin(networkdm,localb,ADD_VALUES,b);
249:   DMLocalToGlobalEnd(networkdm,localb,ADD_VALUES,b);
250:   DMRestoreLocalVector(networkdm,&localb);

252:   MatAssemblyBegin(A,MAT_FINAL_ASSEMBLY);
253:   MatAssemblyEnd(A,MAT_FINAL_ASSEMBLY);
254:   return 0;
255: }

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

274:   PetscInitialize(&argc,&argv,(char*)0,help);
275:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
276:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

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

280:   PetscLogStageRegister("Network Creation", &stage[0]);
281:   PetscLogStageRegister("DMNetwork data structures", &stage[1]);
282:   PetscLogStageRegister("KSP", &stage[2]);

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

293:   PetscLogStagePush(stage[1]);
294:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
295:   DMNetworkRegisterComponent(networkdm,"nstr",sizeof(Node),&componentkey[0]);
296:   DMNetworkRegisterComponent(networkdm,"bsrt",sizeof(Branch),&componentkey[1]);

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

302:   /* Set up the network layout */
303:   DMNetworkLayoutSetUp(networkdm);

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

312:     DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
313:     for (i = vStart; i < vEnd; i++) {
314:       DMNetworkAddComponent(networkdm,i,componentkey[0],&node[i-vStart],1);
315:     }
316:   }

318:   /* Network partitioning and distribution of data */
319:   DMSetUp(networkdm);
320:   DMNetworkDistribute(&networkdm,0);
321:   DMNetworkAssembleGraphStructures(networkdm);

323:   /* We don't use these data structures anymore since they have been copied to networkdm */
324:   if (rank == 0) {
325:     PetscFree(edgelist);
326:     PetscFree2(node,branch);
327:   }

329:   /* Create vectors and matrix */
330:   DMCreateGlobalVector(networkdm,&x);
331:   VecDuplicate(x,&b);
332:   DMCreateMatrix(networkdm,&A);

334:   PetscLogStagePop();

336:   PetscLogStagePush(stage[2]);
337:   /* Assembly system of equations */
338:   FormOperator(networkdm,A,b);

340:   /* Solve linear system: A x = b */
341:   KSPCreate(PETSC_COMM_WORLD, &ksp);
342:   KSPSetOperators(ksp, A, A);
343:   KSPSetFromOptions(ksp);
344:   KSPSolve(ksp, b, x);

346:   PetscLogStagePop();

348:   /* Free work space */
349:   VecDestroy(&x);
350:   VecDestroy(&b);
351:   MatDestroy(&A);
352:   KSPDestroy(&ksp);
353:   DMDestroy(&networkdm);
354:   PetscFinalize();
355:   return 0;
356: }

358: /*TEST

360:    build:
361:       requires: !single double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

363:    test:
364:       args: -ksp_converged_reason

366:    test:
367:       suffix: 2
368:       nsize: 2
369:       args: -petscpartitioner_type simple -pc_type asm -sub_pc_type ilu -ksp_converged_reason

371:    test:
372:       suffix: 3
373:       nsize: 4
374:       args: -petscpartitioner_type simple -pc_type asm -sub_pc_type lu -sub_pc_factor_shift_type nonzero -ksp_converged_reason

376:    test:
377:       suffix: graphindex
378:       args: -n 20 -vertex_global_section_view -edge_global_section_view

380:    test:
381:       suffix: graphindex_2
382:       nsize: 2
383:       args: -petscpartitioner_type simple -n 20 -vertex_global_section_view -edge_global_section_view

385: TEST*/