Actual source code: ex1.c

  1: static char help[] = "This example demonstrates the use of DMNetwork interface with subnetworks for solving a coupled nonlinear \n\
  2:                       electric power grid and water pipe problem.\n\
  3:                       The available solver options are in the ex1options file \n\
  4:                       and the data files are in the datafiles of subdirectories.\n\
  5:                       This example shows the use of subnetwork feature in DMNetwork. \n\
  6:                       Run this program: mpiexec -n <n> ./ex1 \n\\n";

  8: #include "power/power.h"
  9: #include "water/water.h"

 11: typedef struct {
 12:   UserCtx_Power appctx_power;
 13:   AppCtx_Water  appctx_water;
 14:   PetscInt      subsnes_id; /* snes solver id */
 15:   PetscInt      it;         /* iteration number */
 16:   Vec           localXold;  /* store previous solution, used by FormFunction_Dummy() */
 17: } UserCtx;

 19: PetscErrorCode UserMonitor(SNES snes, PetscInt its, PetscReal fnorm, void *appctx)
 20: {
 21:   UserCtx    *user = (UserCtx *)appctx;
 22:   Vec         X, localXold = user->localXold;
 23:   DM          networkdm;
 24:   PetscMPIInt rank;
 25:   MPI_Comm    comm;

 27:   PetscFunctionBegin;
 28:   PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
 29:   PetscCallMPI(MPI_Comm_rank(comm, &rank));
 30: #if 0
 31:   if (rank == 0) {
 32:     PetscInt       subsnes_id = user->subsnes_id;
 33:     if (subsnes_id == 2) {
 34:       PetscCall(PetscPrintf(PETSC_COMM_SELF," it %" PetscInt_FMT ", subsnes_id %" PetscInt_FMT ", fnorm %g\n",user->it,user->subsnes_id,(double)fnorm));
 35:     } else {
 36:       PetscCall(PetscPrintf(PETSC_COMM_SELF,"       subsnes_id %" PetscInt_FMT ", fnorm %g\n",user->subsnes_id,(double)fnorm));
 37:     }
 38:   }
 39: #endif
 40:   PetscCall(SNESGetSolution(snes, &X));
 41:   PetscCall(SNESGetDM(snes, &networkdm));
 42:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localXold));
 43:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localXold));
 44:   PetscFunctionReturn(PETSC_SUCCESS);
 45: }

 47: PetscErrorCode FormJacobian_subPower(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx)
 48: {
 49:   DM              networkdm;
 50:   Vec             localX;
 51:   PetscInt        nv, ne, i, j, offset, nvar, row;
 52:   const PetscInt *vtx, *edges;
 53:   PetscBool       ghostvtex;
 54:   PetscScalar     one = 1.0;
 55:   PetscMPIInt     rank;
 56:   MPI_Comm        comm;

 58:   PetscFunctionBegin;
 59:   PetscCall(SNESGetDM(snes, &networkdm));
 60:   PetscCall(DMGetLocalVector(networkdm, &localX));

 62:   PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
 63:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

 65:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
 66:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

 68:   PetscCall(MatZeroEntries(J));

 70:   /* Power subnetwork: copied from power/FormJacobian_Power() */
 71:   PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
 72:   PetscCall(FormJacobian_Power_private(networkdm, localX, J, nv, ne, vtx, edges, appctx));

 74:   /* Water subnetwork: Identity */
 75:   PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
 76:   for (i = 0; i < nv; i++) {
 77:     PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex));
 78:     if (ghostvtex) continue;

 80:     PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset));
 81:     PetscCall(DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar));
 82:     for (j = 0; j < nvar; j++) {
 83:       row = offset + j;
 84:       PetscCall(MatSetValues(J, 1, &row, 1, &row, &one, ADD_VALUES));
 85:     }
 86:   }
 87:   PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
 88:   PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));

 90:   PetscCall(DMRestoreLocalVector(networkdm, &localX));
 91:   PetscFunctionReturn(PETSC_SUCCESS);
 92: }

 94: /* Dummy equation localF(X) = localX - localXold */
 95: PetscErrorCode FormFunction_Dummy(DM networkdm, Vec localX, Vec localF, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
 96: {
 97:   const PetscScalar *xarr, *xoldarr;
 98:   PetscScalar       *farr;
 99:   PetscInt           i, j, offset, nvar;
100:   PetscBool          ghostvtex;
101:   UserCtx           *user      = (UserCtx *)appctx;
102:   Vec                localXold = user->localXold;

104:   PetscFunctionBegin;
105:   PetscCall(VecGetArrayRead(localX, &xarr));
106:   PetscCall(VecGetArrayRead(localXold, &xoldarr));
107:   PetscCall(VecGetArray(localF, &farr));

109:   for (i = 0; i < nv; i++) {
110:     PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex));
111:     if (ghostvtex) continue;

113:     PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset));
114:     PetscCall(DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar));
115:     for (j = 0; j < nvar; j++) farr[offset + j] = xarr[offset + j] - xoldarr[offset + j];
116:   }

118:   PetscCall(VecRestoreArrayRead(localX, &xarr));
119:   PetscCall(VecRestoreArrayRead(localXold, &xoldarr));
120:   PetscCall(VecRestoreArray(localF, &farr));
121:   PetscFunctionReturn(PETSC_SUCCESS);
122: }

124: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *appctx)
125: {
126:   DM              networkdm;
127:   Vec             localX, localF;
128:   PetscInt        nv, ne, v;
129:   const PetscInt *vtx, *edges;
130:   PetscMPIInt     rank;
131:   MPI_Comm        comm;
132:   UserCtx        *user         = (UserCtx *)appctx;
133:   UserCtx_Power   appctx_power = (*user).appctx_power;
134:   AppCtx_Water    appctx_water = (*user).appctx_water;

136:   PetscFunctionBegin;
137:   PetscCall(SNESGetDM(snes, &networkdm));
138:   PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
139:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

141:   PetscCall(DMGetLocalVector(networkdm, &localX));
142:   PetscCall(DMGetLocalVector(networkdm, &localF));
143:   PetscCall(VecSet(F, 0.0));
144:   PetscCall(VecSet(localF, 0.0));

146:   PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
147:   PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));

149:   /* Form Function for power subnetwork */
150:   PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
151:   if (user->subsnes_id == 1) { /* snes_water only */
152:     PetscCall(FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user));
153:   } else {
154:     PetscCall(FormFunction_Power(networkdm, localX, localF, nv, ne, vtx, edges, &appctx_power));
155:   }

157:   /* Form Function for water subnetwork */
158:   PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
159:   if (user->subsnes_id == 0) { /* snes_power only */
160:     PetscCall(FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user));
161:   } else {
162:     PetscCall(FormFunction_Water(networkdm, localX, localF, nv, ne, vtx, edges, NULL));
163:   }

165:   /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */
166:   PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
167:   for (v = 0; v < nv; v++) {
168:     PetscInt        key, ncomp, nvar, nconnedges, k, e, keye, goffset[3];
169:     void           *component;
170:     const PetscInt *connedges;

172:     PetscCall(DMNetworkGetComponent(networkdm, vtx[v], ALL_COMPONENTS, NULL, NULL, &nvar));
173:     PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &ncomp));
174:     /* printf("  [%d] coupling vertex[%" PetscInt_FMT "]: v %" PetscInt_FMT ", ncomp %" PetscInt_FMT "; nvar %" PetscInt_FMT "\n",rank,v,vtx[v], ncomp,nvar); */

176:     for (k = 0; k < ncomp; k++) {
177:       PetscCall(DMNetworkGetComponent(networkdm, vtx[v], k, &key, &component, &nvar));
178:       PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[v], k, &goffset[k]));

180:       /* Verify the coupling vertex is a powernet load vertex or a water vertex */
181:       switch (k) {
182:       case 0:
183:         PetscCheck(key == appctx_power.compkey_bus && nvar == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "key %" PetscInt_FMT " not a power bus vertex or nvar %" PetscInt_FMT " != 2", key, nvar);
184:         break;
185:       case 1:
186:         PetscCheck(key == appctx_power.compkey_load && nvar == 0 && goffset[1] == goffset[0] + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a power load vertex");
187:         break;
188:       case 2:
189:         PetscCheck(key == appctx_water.compkey_vtx && nvar == 1 && goffset[2] == goffset[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a water vertex");
190:         break;
191:       default:
192:         SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %" PetscInt_FMT " is wrong", k);
193:       }
194:       /* printf("  [%d] coupling vertex[%" PetscInt_FMT "]: key %" PetscInt_FMT "; nvar %" PetscInt_FMT ", goffset %" PetscInt_FMT "\n",rank,v,key,nvar,goffset[k]); */
195:     }

197:     /* Get its supporting edges */
198:     PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges));
199:     /* printf("\n[%d] coupling vertex: nconnedges %" PetscInt_FMT "\n",rank,nconnedges); */
200:     for (k = 0; k < nconnedges; k++) {
201:       e = connedges[k];
202:       PetscCall(DMNetworkGetNumComponents(networkdm, e, &ncomp));
203:       /* printf("\n  [%d] connected edge[%" PetscInt_FMT "]=%" PetscInt_FMT " has ncomp %" PetscInt_FMT "\n",rank,k,e,ncomp); */
204:       PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, &component, NULL));
205:       if (keye != appctx_water.compkey_edge) PetscCheck(keye == appctx_power.compkey_branch, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a power branch");
206:     }
207:   }

209:   PetscCall(DMRestoreLocalVector(networkdm, &localX));

211:   PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
212:   PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
213:   PetscCall(DMRestoreLocalVector(networkdm, &localF));
214: #if 0
215:   if (rank == 0) printf("F:\n");
216:   PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD));
217: #endif
218:   PetscFunctionReturn(PETSC_SUCCESS);
219: }

221: PetscErrorCode SetInitialGuess(DM networkdm, Vec X, void *appctx)
222: {
223:   PetscInt        nv, ne, i, j, ncomp, offset, key;
224:   const PetscInt *vtx, *edges;
225:   UserCtx        *user         = (UserCtx *)appctx;
226:   Vec             localX       = user->localXold;
227:   UserCtx_Power   appctx_power = (*user).appctx_power;
228:   AppCtx_Water    appctx_water = (*user).appctx_water;
229:   PetscBool       ghost;
230:   PetscScalar    *xarr;
231:   VERTEX_Power    bus;
232:   VERTEX_Water    vertex;
233:   void           *component;
234:   GEN             gen;

236:   PetscFunctionBegin;
237:   PetscCall(VecSet(X, 0.0));
238:   PetscCall(VecSet(localX, 0.0));

240:   /* Set initial guess for power subnetwork */
241:   PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
242:   PetscCall(SetInitialGuess_Power(networkdm, localX, nv, ne, vtx, edges, &appctx_power));

244:   /* Set initial guess for water subnetwork */
245:   PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
246:   PetscCall(SetInitialGuess_Water(networkdm, localX, nv, ne, vtx, edges, NULL));

248:   /* Set initial guess at the coupling vertex */
249:   PetscCall(VecGetArray(localX, &xarr));
250:   PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
251:   for (i = 0; i < nv; i++) {
252:     PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghost));
253:     if (ghost) continue;

255:     PetscCall(DMNetworkGetNumComponents(networkdm, vtx[i], &ncomp));
256:     for (j = 0; j < ncomp; j++) {
257:       PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], j, &offset));
258:       PetscCall(DMNetworkGetComponent(networkdm, vtx[i], j, &key, (void **)&component, NULL));
259:       if (key == appctx_power.compkey_bus) {
260:         bus              = (VERTEX_Power)(component);
261:         xarr[offset]     = bus->va * PETSC_PI / 180.0;
262:         xarr[offset + 1] = bus->vm;
263:       } else if (key == appctx_power.compkey_gen) {
264:         gen = (GEN)(component);
265:         if (!gen->status) continue;
266:         xarr[offset + 1] = gen->vs;
267:       } else if (key == appctx_water.compkey_vtx) {
268:         vertex = (VERTEX_Water)(component);
269:         if (vertex->type == VERTEX_TYPE_JUNCTION) {
270:           xarr[offset] = 100;
271:         } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
272:           xarr[offset] = vertex->res.head;
273:         } else {
274:           xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
275:         }
276:       }
277:     }
278:   }
279:   PetscCall(VecRestoreArray(localX, &xarr));

281:   PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
282:   PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
283:   PetscFunctionReturn(PETSC_SUCCESS);
284: }

286: /* Set coordinates */
287: static PetscErrorCode CoordinateVecSetUp(DM dm, Vec coords)
288: {
289:   DM              dmclone;
290:   PetscInt        i, gidx, offset, v, nv, Nsubnet;
291:   const PetscInt *vtx;
292:   PetscScalar    *carray;

294:   PetscFunctionBeginUser;
295:   PetscCall(DMGetCoordinateDM(dm, &dmclone));
296:   PetscCall(VecGetArrayWrite(coords, &carray));
297:   PetscCall(DMNetworkGetNumSubNetworks(dm, NULL, &Nsubnet));
298:   for (i = 0; i < Nsubnet; i++) {
299:     PetscCall(DMNetworkGetSubnetwork(dm, i, &nv, NULL, &vtx, NULL));
300:     for (v = 0; v < nv; v++) {
301:       PetscCall(DMNetworkGetGlobalVertexIndex(dm, vtx[v], &gidx));
302:       PetscCall(DMNetworkGetLocalVecOffset(dmclone, vtx[v], 0, &offset));
303:       switch (gidx) {
304:       case 0:
305:         carray[offset]     = -1.0;
306:         carray[offset + 1] = -1.0;
307:         break;
308:       case 1:
309:         carray[offset]     = -2.0;
310:         carray[offset + 1] = 2.0;
311:         break;
312:       case 2:
313:         carray[offset]     = 0.0;
314:         carray[offset + 1] = 2.0;
315:         break;
316:       case 3:
317:         carray[offset]     = -1.0;
318:         carray[offset + 1] = 0.0;
319:         break;
320:       case 4:
321:         carray[offset]     = 0.0;
322:         carray[offset + 1] = 0.0;
323:         break;
324:       case 5:
325:         carray[offset]     = 0.0;
326:         carray[offset + 1] = 1.0;
327:         break;
328:       case 6:
329:         carray[offset]     = -1.0;
330:         carray[offset + 1] = 1.0;
331:         break;
332:       case 7:
333:         carray[offset]     = -2.0;
334:         carray[offset + 1] = 1.0;
335:         break;
336:       case 8:
337:         carray[offset]     = -2.0;
338:         carray[offset + 1] = 0.0;
339:         break;
340:       case 9:
341:         carray[offset]     = 1.0;
342:         carray[offset + 1] = 0.0;
343:         break;
344:       case 10:
345:         carray[offset]     = 1.0;
346:         carray[offset + 1] = -1.0;
347:         break;
348:       case 11:
349:         carray[offset]     = 2.0;
350:         carray[offset + 1] = -1.0;
351:         break;
352:       case 12:
353:         carray[offset]     = 2.0;
354:         carray[offset + 1] = 0.0;
355:         break;
356:       case 13:
357:         carray[offset]     = 0.0;
358:         carray[offset + 1] = -1.0;
359:         break;
360:       case 14:
361:         carray[offset]     = 2.0;
362:         carray[offset + 1] = 1.0;
363:         break;
364:       default:
365:         PetscCheck(gidx < 15 && gidx > -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "gidx %" PetscInt_FMT "must between 0 and 14", gidx);
366:       }
367:     }
368:   }
369:   PetscCall(VecRestoreArrayWrite(coords, &carray));
370:   PetscFunctionReturn(PETSC_SUCCESS);
371: }

373: static PetscErrorCode CoordinatePrint(DM dm)
374: {
375:   DM                 dmclone;
376:   PetscInt           cdim, v, off, vglobal, vStart, vEnd;
377:   const PetscScalar *carray;
378:   Vec                coords;
379:   MPI_Comm           comm;
380:   PetscMPIInt        rank;

382:   PetscFunctionBegin;
383:   PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
384:   PetscCallMPI(MPI_Comm_rank(comm, &rank));

386:   PetscCall(DMGetCoordinateDM(dm, &dmclone));
387:   PetscCall(DMNetworkGetVertexRange(dm, &vStart, &vEnd));
388:   PetscCall(DMGetCoordinatesLocal(dm, &coords));

390:   PetscCall(DMGetCoordinateDim(dm, &cdim));
391:   PetscCall(VecGetArrayRead(coords, &carray));

393:   PetscCall(PetscPrintf(MPI_COMM_WORLD, "\nCoordinatePrint, cdim %" PetscInt_FMT ":\n", cdim));
394:   PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%i]\n", rank));
395:   for (v = vStart; v < vEnd; v++) {
396:     PetscCall(DMNetworkGetLocalVecOffset(dmclone, v, 0, &off));
397:     PetscCall(DMNetworkGetGlobalVertexIndex(dmclone, v, &vglobal));
398:     switch (cdim) {
399:     case 2:
400:       PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "Vertex: %" PetscInt_FMT ", x =  %f y = %f \n", vglobal, (double)PetscRealPart(carray[off]), (double)PetscRealPart(carray[off + 1])));
401:       break;
402:     default:
403:       PetscCheck(cdim == 2, MPI_COMM_WORLD, PETSC_ERR_SUP, "Only supports Network embedding dimension of 2, not supplied  %" PetscInt_FMT, cdim);
404:       break;
405:     }
406:   }
407:   PetscCall(PetscSynchronizedFlush(MPI_COMM_WORLD, NULL));
408:   PetscCall(VecRestoreArrayRead(coords, &carray));
409:   PetscFunctionReturn(PETSC_SUCCESS);
410: }

412: int main(int argc, char **argv)
413: {
414:   DM                  networkdm, dmclone;
415:   PetscMPIInt         rank, size;
416:   PetscInt            Nsubnet = 2, numVertices[2], numEdges[2], i, j, nv, ne, it_max = 10;
417:   PetscInt            vStart, vEnd, compkey;
418:   const PetscInt     *vtx, *edges;
419:   Vec                 X, F, coords;
420:   SNES                snes, snes_power, snes_water;
421:   Mat                 Jac;
422:   PetscBool           ghost, viewJ = PETSC_FALSE, viewX = PETSC_FALSE, test = PETSC_FALSE, distribute = PETSC_TRUE, flg, printCoord = PETSC_FALSE, viewCSV = PETSC_FALSE;
423:   UserCtx             user;
424:   SNESConvergedReason reason;
425: #if defined(PETSC_USE_LOG)
426:   PetscLogStage stage[4];
427: #endif

429:   /* Power subnetwork */
430:   UserCtx_Power *appctx_power                    = &user.appctx_power;
431:   char           pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
432:   PFDATA        *pfdata                          = NULL;
433:   PetscInt       genj, loadj, *edgelist_power = NULL, power_netnum;
434:   PetscScalar    Sbase = 0.0;

436:   /* Water subnetwork */
437:   AppCtx_Water *appctx_water                       = &user.appctx_water;
438:   WATERDATA    *waterdata                          = NULL;
439:   char          waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
440:   PetscInt     *edgelist_water                     = NULL, water_netnum;

442:   /* Shared vertices between subnetworks */
443:   PetscInt power_svtx, water_svtx;

445:   PetscFunctionBeginUser;
446:   PetscCall(PetscInitialize(&argc, &argv, "ex1options", help));
447:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
448:   PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));

450:   /* (1) Read Data - Only rank 0 reads the data */
451:   PetscCall(PetscLogStageRegister("Read Data", &stage[0]));
452:   PetscCall(PetscLogStagePush(stage[0]));

454:   for (i = 0; i < Nsubnet; i++) {
455:     numVertices[i] = 0;
456:     numEdges[i]    = 0;
457:   }

459:   /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
460:   /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
461:   PetscCall(PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, PETSC_MAX_PATH_LEN - 1, NULL));
462:   PetscCall(PetscNew(&pfdata));
463:   PetscCall(PFReadMatPowerData(pfdata, pfdata_file));
464:   if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Power network: nb = %" PetscInt_FMT ", ngen = %" PetscInt_FMT ", nload = %" PetscInt_FMT ", nbranch = %" PetscInt_FMT "\n", pfdata->nbus, pfdata->ngen, pfdata->nload, pfdata->nbranch));
465:   Sbase = pfdata->sbase;
466:   if (rank == 0) { /* proc[0] will create Electric Power Grid */
467:     numEdges[0]    = pfdata->nbranch;
468:     numVertices[0] = pfdata->nbus;

470:     PetscCall(PetscMalloc1(2 * numEdges[0], &edgelist_power));
471:     PetscCall(GetListofEdges_Power(pfdata, edgelist_power));
472:   }
473:   /* Broadcast power Sbase to all processors */
474:   PetscCallMPI(MPI_Bcast(&Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD));
475:   appctx_power->Sbase     = Sbase;
476:   appctx_power->jac_error = PETSC_FALSE;
477:   /* If external option activated. Introduce error in jacobian */
478:   PetscCall(PetscOptionsHasName(NULL, NULL, "-jac_error", &appctx_power->jac_error));

480:   /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
481:   /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
482:   PetscCall(PetscNew(&waterdata));
483:   PetscCall(PetscOptionsGetString(NULL, NULL, "-waterdata", waterdata_file, PETSC_MAX_PATH_LEN - 1, NULL));
484:   PetscCall(WaterReadData(waterdata, waterdata_file));
485:   if (size == 1 || (size > 1 && rank == 1)) {
486:     PetscCall(PetscCalloc1(2 * waterdata->nedge, &edgelist_water));
487:     PetscCall(GetListofEdges_Water(waterdata, edgelist_water));
488:     numEdges[1]    = waterdata->nedge;
489:     numVertices[1] = waterdata->nvertex;
490:   }
491:   PetscCall(PetscLogStagePop());

493:   /* (2) Create a network consist of two subnetworks */
494:   PetscCall(PetscLogStageRegister("Net Setup", &stage[1]));
495:   PetscCall(PetscLogStagePush(stage[1]));

497:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewCSV", &viewCSV, NULL));
498:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-printCoord", &printCoord, NULL));
499:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-test", &test, NULL));
500:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-distribute", &distribute, NULL));

502:   /* Create an empty network object */
503:   PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));

505:   /* Register the components in the network */
506:   PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &appctx_power->compkey_branch));
507:   PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &appctx_power->compkey_bus));
508:   PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &appctx_power->compkey_gen));
509:   PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &appctx_power->compkey_load));

511:   PetscCall(DMNetworkRegisterComponent(networkdm, "edge_water", sizeof(struct _p_EDGE_Water), &appctx_water->compkey_edge));
512:   PetscCall(DMNetworkRegisterComponent(networkdm, "vertex_water", sizeof(struct _p_VERTEX_Water), &appctx_water->compkey_vtx));

514:   PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] Total local nvertices %" PetscInt_FMT " + %" PetscInt_FMT " = %" PetscInt_FMT ", nedges %" PetscInt_FMT " + %" PetscInt_FMT " = %" PetscInt_FMT "\n", rank, numVertices[0], numVertices[1], numVertices[0] + numVertices[1], numEdges[0], numEdges[1], numEdges[0] + numEdges[1]));
515:   PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));

517:   PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, Nsubnet));
518:   PetscCall(DMNetworkAddSubnetwork(networkdm, "power", numEdges[0], edgelist_power, &power_netnum));
519:   PetscCall(DMNetworkAddSubnetwork(networkdm, "water", numEdges[1], edgelist_water, &water_netnum));

521:   /* vertex subnet[0].4 shares with vertex subnet[1].0 */
522:   power_svtx = 4;
523:   water_svtx = 0;
524:   PetscCall(DMNetworkAddSharedVertices(networkdm, power_netnum, water_netnum, 1, &power_svtx, &water_svtx));

526:   /* Set up the network layout */
527:   PetscCall(DMNetworkLayoutSetUp(networkdm));

529:   /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
530:   genj  = 0;
531:   loadj = 0;
532:   PetscCall(DMNetworkGetSubnetwork(networkdm, power_netnum, &nv, &ne, &vtx, &edges));

534:   for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx_power->compkey_branch, &pfdata->branch[i], 0));

536:   for (i = 0; i < nv; i++) {
537:     PetscCall(DMNetworkIsSharedVertex(networkdm, vtx[i], &flg));
538:     if (flg) continue;

540:     PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[i], 2));
541:     if (pfdata->bus[i].ngen) {
542:       for (j = 0; j < pfdata->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_gen, &pfdata->gen[genj++], 0));
543:     }
544:     if (pfdata->bus[i].nload) {
545:       for (j = 0; j < pfdata->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[loadj++], 0));
546:     }
547:   }

549:   /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
550:   PetscCall(DMNetworkGetSubnetwork(networkdm, water_netnum, &nv, &ne, &vtx, &edges));
551:   for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx_water->compkey_edge, &waterdata->edge[i], 0));

553:   for (i = 0; i < nv; i++) {
554:     PetscCall(DMNetworkIsSharedVertex(networkdm, vtx[i], &flg));
555:     if (flg) continue;

557:     PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[i], 1));
558:   }

560:   /* ADD VARIABLES AND COMPONENTS AT THE SHARED VERTEX: net[0].4 coupls with net[1].0 -- owning and all ghost ranks of the vertex do this */
561:   PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
562:   for (i = 0; i < nv; i++) {
563:     /* power */
564:     PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[4], 2));
565:     /* bus[4] is a load, add its component */
566:     PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[0], 0));

568:     /* water */
569:     PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[0], 1));
570:   }

572:   /* Set coordinates for visualization */
573:   PetscCall(DMSetCoordinateDim(networkdm, 2));
574:   PetscCall(DMGetCoordinateDM(networkdm, &dmclone));
575:   PetscCall(DMNetworkGetVertexRange(dmclone, &vStart, &vEnd));
576:   PetscCall(DMNetworkRegisterComponent(dmclone, "coordinates", 0, &compkey));
577:   for (i = vStart; i < vEnd; i++) PetscCall(DMNetworkAddComponent(dmclone, i, compkey, NULL, 2));
578:   PetscCall(DMNetworkFinalizeComponents(dmclone));

580:   PetscCall(DMCreateLocalVector(dmclone, &coords));
581:   PetscCall(DMSetCoordinatesLocal(networkdm, coords)); /* set/get coords to/from networkdm */
582:   PetscCall(CoordinateVecSetUp(networkdm, coords));
583:   if (printCoord) PetscCall(CoordinatePrint(networkdm));

585:   /* Set up DM for use */
586:   PetscCall(DMSetUp(networkdm));

588:   /* Free user objects */
589:   PetscCall(PetscFree(edgelist_power));
590:   PetscCall(PetscFree(pfdata->bus));
591:   PetscCall(PetscFree(pfdata->gen));
592:   PetscCall(PetscFree(pfdata->branch));
593:   PetscCall(PetscFree(pfdata->load));
594:   PetscCall(PetscFree(pfdata));

596:   PetscCall(PetscFree(edgelist_water));
597:   PetscCall(PetscFree(waterdata->vertex));
598:   PetscCall(PetscFree(waterdata->edge));
599:   PetscCall(PetscFree(waterdata));

601:   /* Re-distribute networkdm to multiple processes for better job balance */
602:   if (distribute) {
603:     PetscCall(DMNetworkDistribute(&networkdm, 0));

605:     if (printCoord) PetscCall(CoordinatePrint(networkdm));
606:     if (viewCSV) { /* CSV View of network with coordinates */
607:       PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_CSV));
608:       PetscCall(DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD));
609:       PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
610:     }
611:   }
612:   PetscCall(VecDestroy(&coords));

614:   /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
615:   if (test) {
616:     PetscInt v, gidx;
617:     PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
618:     for (i = 0; i < Nsubnet; i++) {
619:       PetscCall(DMNetworkGetSubnetwork(networkdm, i, &nv, &ne, &vtx, &edges));
620:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, subnet[%" PetscInt_FMT "] ne %" PetscInt_FMT ", nv %" PetscInt_FMT "\n", rank, i, ne, nv));
621:       PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));

623:       for (v = 0; v < nv; v++) {
624:         PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghost));
625:         PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx));
626:         PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] subnet[%" PetscInt_FMT "] v %" PetscInt_FMT " %" PetscInt_FMT "; ghost %d\n", rank, i, vtx[v], gidx, ghost));
627:       }
628:       PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
629:     }
630:     PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));

632:     PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
633:     PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, num of shared vertices nsv = %" PetscInt_FMT "\n", rank, nv));
634:     for (v = 0; v < nv; v++) {
635:       PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx));
636:       PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] sv %" PetscInt_FMT ", gidx=%" PetscInt_FMT "\n", rank, vtx[v], gidx));
637:     }
638:     PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
639:   }

641:   /* Create solution vector X */
642:   PetscCall(DMCreateGlobalVector(networkdm, &X));
643:   PetscCall(VecDuplicate(X, &F));
644:   PetscCall(DMGetLocalVector(networkdm, &user.localXold));
645:   PetscCall(PetscLogStagePop());

647:   /* (3) Setup Solvers */
648:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewJ", &viewJ, NULL));
649:   PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewX", &viewX, NULL));

651:   PetscCall(PetscLogStageRegister("SNES Setup", &stage[2]));
652:   PetscCall(PetscLogStagePush(stage[2]));

654:   PetscCall(SetInitialGuess(networkdm, X, &user));

656:   /* Create coupled snes */
657:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_coupled setup ......\n"));
658:   user.subsnes_id = Nsubnet;
659:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
660:   PetscCall(SNESSetDM(snes, networkdm));
661:   PetscCall(SNESSetOptionsPrefix(snes, "coupled_"));
662:   PetscCall(SNESSetFunction(snes, F, FormFunction, &user));
663:   PetscCall(SNESMonitorSet(snes, UserMonitor, &user, NULL));
664:   PetscCall(SNESSetFromOptions(snes));

666:   if (viewJ) {
667:     /* View Jac structure */
668:     PetscCall(SNESGetJacobian(snes, &Jac, NULL, NULL, NULL));
669:     PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
670:   }

672:   if (viewX) {
673:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution:\n"));
674:     PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
675:   }

677:   if (viewJ) {
678:     /* View assembled Jac */
679:     PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
680:   }

682:   /* Create snes_power */
683:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_power setup ......\n"));

685:   user.subsnes_id = 0;
686:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_power));
687:   PetscCall(SNESSetDM(snes_power, networkdm));
688:   PetscCall(SNESSetOptionsPrefix(snes_power, "power_"));
689:   PetscCall(SNESSetFunction(snes_power, F, FormFunction, &user));
690:   PetscCall(SNESMonitorSet(snes_power, UserMonitor, &user, NULL));

692:   /* Use user-provide Jacobian */
693:   PetscCall(DMCreateMatrix(networkdm, &Jac));
694:   PetscCall(SNESSetJacobian(snes_power, Jac, Jac, FormJacobian_subPower, &user));
695:   PetscCall(SNESSetFromOptions(snes_power));

697:   if (viewX) {
698:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Power Solution:\n"));
699:     PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
700:   }
701:   if (viewJ) {
702:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Power Jac:\n"));
703:     PetscCall(SNESGetJacobian(snes_power, &Jac, NULL, NULL, NULL));
704:     PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
705:     /* PetscCall(MatView(Jac,PETSC_VIEWER_STDOUT_WORLD)); */
706:   }

708:   /* Create snes_water */
709:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_water setup......\n"));

711:   user.subsnes_id = 1;
712:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_water));
713:   PetscCall(SNESSetDM(snes_water, networkdm));
714:   PetscCall(SNESSetOptionsPrefix(snes_water, "water_"));
715:   PetscCall(SNESSetFunction(snes_water, F, FormFunction, &user));
716:   PetscCall(SNESMonitorSet(snes_water, UserMonitor, &user, NULL));
717:   PetscCall(SNESSetFromOptions(snes_water));

719:   if (viewX) {
720:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Water Solution:\n"));
721:     PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
722:   }
723:   PetscCall(PetscLogStagePop());

725:   /* (4) Solve */
726:   PetscCall(PetscLogStageRegister("SNES Solve", &stage[3]));
727:   PetscCall(PetscLogStagePush(stage[3]));
728:   user.it = 0;
729:   reason  = SNES_DIVERGED_DTOL;
730:   while (user.it < it_max && (PetscInt)reason < 0) {
731: #if 0
732:     user.subsnes_id = 0;
733:     PetscCall(SNESSolve(snes_power,NULL,X));

735:     user.subsnes_id = 1;
736:     PetscCall(SNESSolve(snes_water,NULL,X));
737: #endif
738:     user.subsnes_id = Nsubnet;
739:     PetscCall(SNESSolve(snes, NULL, X));

741:     PetscCall(SNESGetConvergedReason(snes, &reason));
742:     user.it++;
743:   }
744:   PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coupled_SNES converged in %" PetscInt_FMT " iterations\n", user.it));
745:   if (viewX) {
746:     PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final Solution:\n"));
747:     PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
748:   }
749:   PetscCall(PetscLogStagePop());

751:   /* Free objects */
752:   /* -------------*/
753:   PetscCall(VecDestroy(&X));
754:   PetscCall(VecDestroy(&F));
755:   PetscCall(DMRestoreLocalVector(networkdm, &user.localXold));

757:   PetscCall(SNESDestroy(&snes));
758:   PetscCall(MatDestroy(&Jac));
759:   PetscCall(SNESDestroy(&snes_power));
760:   PetscCall(SNESDestroy(&snes_water));

762:   PetscCall(DMDestroy(&networkdm));
763:   PetscCall(PetscFinalize());
764:   return 0;
765: }

767: /*TEST

769:    build:
770:      requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
771:      depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c

773:    test:
774:       args: -coupled_snes_converged_reason -options_left no -dmnetwork_view -fp_trap 0
775:       localrunfiles: ex1options power/case9.m water/sample1.inp
776:       output_file: output/ex1.out

778:    test:
779:       suffix: 2
780:       nsize: 3
781:       args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis -fp_trap 0
782:       localrunfiles: ex1options power/case9.m water/sample1.inp
783:       output_file: output/ex1_2.out
784:       requires: parmetis

786:    test:
787:       suffix: 3
788:       nsize: 3
789:       args: -coupled_snes_converged_reason -options_left no -distribute false -fp_trap 0
790:       localrunfiles: ex1options power/case9.m water/sample1.inp
791:       output_file: output/ex1_2.out

793:    test:
794:       suffix: 4
795:       nsize: 4
796:       args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -dmnetwork_view -dmnetwork_view_distributed -fp_trap 0
797:       localrunfiles: ex1options power/case9.m water/sample1.inp
798:       output_file: output/ex1_4.out

800:    test:
801:       suffix: 5
802:       args: -coupled_snes_converged_reason -options_left no -viewCSV -fp_trap 0
803:       localrunfiles: ex1options power/case9.m water/sample1.inp
804:       output_file: output/ex1_5.out

806:    test:
807:       suffix: 6
808:       nsize: 3
809:       args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis -dmnetwork_view_distributed draw:null -fp_trap 0
810:       localrunfiles: ex1options power/case9.m water/sample1.inp
811:       output_file: output/ex1_2.out
812:       requires: parmetis

814: TEST*/