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;
206:     val[0] = 1;
207:     col[1] = lofst_to;
208:     val[1] = 1;
209:     col[2] = lofst_fr;
210:     val[2] = -1;
211:     MatSetValuesLocal(A, 1, row, 3, col, val, ADD_VALUES);

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

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

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

226:     if (!node->gr) {
227:       row[0] = lofst_to;
228:       col[0] = lofst;
229:       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;
244:         val[0] = 1;
245:         MatSetValuesLocal(A, 1, row, 1, col, val, ADD_VALUES);
246:       } else {
247:         barr[lofst] -= node->inj;
248:       }
249:     }
250:   }

252:   VecRestoreArray(localb, &barr);

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

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

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

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

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

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

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

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

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

309:   /* Set up the network layout */
310:   DMNetworkLayoutSetUp(networkdm);

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

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

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

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

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

337:   PetscLogStagePop();

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

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

349:   PetscLogStagePop();

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

361: /*TEST

363:    build:
364:       requires: !single double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

366:    test:
367:       args: -ksp_converged_reason

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

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

379:    test:
380:       suffix: graphindex
381:       args: -n 20 -vertex_global_section_view -edge_global_section_view

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

388: TEST*/