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: /* T
  9:    Concepts: DMNetwork
 10:    Concepts: PETSc SNES solver
 11: */

 13: #include "power/power.h"
 14: #include "water/water.h"

 16: typedef struct{
 17:   UserCtx_Power appctx_power;
 18:   AppCtx_Water  appctx_water;
 19:   PetscInt      subsnes_id; /* snes solver id */
 20:   PetscInt      it;         /* iteration number */
 21:   Vec           localXold;  /* store previous solution, used by FormFunction_Dummy() */
 22: } UserCtx;

 24: PetscErrorCode UserMonitor(SNES snes,PetscInt its,PetscReal fnorm ,void *appctx)
 25: {
 27:   UserCtx        *user = (UserCtx*)appctx;
 28:   Vec            X,localXold = user->localXold;
 29:   DM             networkdm;
 30:   PetscMPIInt    rank;
 31:   MPI_Comm       comm;

 34:   PetscObjectGetComm((PetscObject)snes,&comm);
 35:   MPI_Comm_rank(comm,&rank);
 36: #if 0
 37:   if (!rank) {
 38:     PetscInt       subsnes_id = user->subsnes_id;
 39:     if (subsnes_id == 2) {
 40:       PetscPrintf(PETSC_COMM_SELF," it %D, subsnes_id %D, fnorm %g\n",user->it,user->subsnes_id,(double)fnorm);
 41:     } else {
 42:       PetscPrintf(PETSC_COMM_SELF,"       subsnes_id %D, fnorm %g\n",user->subsnes_id,(double)fnorm);
 43:     }
 44:   }
 45: #endif
 46:   SNESGetSolution(snes,&X);
 47:   SNESGetDM(snes,&networkdm);
 48:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localXold);
 49:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localXold);
 50:   return(0);
 51: }

 53: PetscErrorCode FormJacobian_subPower(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
 54: {
 56:   DM             networkdm;
 57:   Vec            localX;
 58:   PetscInt       nv,ne,i,j,offset,nvar,row;
 59:   const PetscInt *vtx,*edges;
 60:   PetscBool      ghostvtex;
 61:   PetscScalar    one = 1.0;
 62:   PetscMPIInt    rank;
 63:   MPI_Comm       comm;

 66:   SNESGetDM(snes,&networkdm);
 67:   DMGetLocalVector(networkdm,&localX);

 69:   PetscObjectGetComm((PetscObject)networkdm,&comm);
 70:   MPI_Comm_rank(comm,&rank);

 72:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 73:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 75:   MatZeroEntries(J);

 77:   /* Power subnetwork: copied from power/FormJacobian_Power() */
 78:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
 79:   FormJacobian_Power_private(networkdm,localX,J,nv,ne,vtx,edges,appctx);

 81:   /* Water subnetwork: Identity */
 82:   DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
 83:   for (i=0; i<nv; i++) {
 84:     DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
 85:     if (ghostvtex) continue;

 87:     DMNetworkGetGlobalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
 88:     DMNetworkGetComponent(networkdm,vtx[i],ALL_COMPONENTS,NULL,NULL,&nvar);
 89:     for (j=0; j<nvar; j++) {
 90:       row = offset + j;
 91:       MatSetValues(J,1,&row,1,&row,&one,ADD_VALUES);
 92:     }
 93:   }
 94:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
 95:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

 97:   DMRestoreLocalVector(networkdm,&localX);
 98:   return(0);
 99: }

101: /* Dummy equation localF(X) = localX - localXold */
102: PetscErrorCode FormFunction_Dummy(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
103: {
104:   PetscErrorCode    ierr;
105:   const PetscScalar *xarr,*xoldarr;
106:   PetscScalar       *farr;
107:   PetscInt          i,j,offset,nvar;
108:   PetscBool         ghostvtex;
109:   UserCtx           *user = (UserCtx*)appctx;
110:   Vec               localXold = user->localXold;

113:   VecGetArrayRead(localX,&xarr);
114:   VecGetArrayRead(localXold,&xoldarr);
115:   VecGetArray(localF,&farr);

117:   for (i=0; i<nv; i++) {
118:     DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
119:     if (ghostvtex) continue;

121:     DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
122:     DMNetworkGetComponent(networkdm,vtx[i],ALL_COMPONENTS,NULL,NULL,&nvar);
123:     for (j=0; j<nvar; j++) {
124:       farr[offset+j] = xarr[offset+j] - xoldarr[offset+j];
125:     }
126:   }

128:   VecRestoreArrayRead(localX,&xarr);
129:   VecRestoreArrayRead(localXold,&xoldarr);
130:   VecRestoreArray(localF,&farr);
131:   return(0);
132: }

134: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *appctx)
135: {
137:   DM             networkdm;
138:   Vec            localX,localF;
139:   PetscInt       nv,ne,v;
140:   const PetscInt *vtx,*edges;
141:   PetscMPIInt    rank;
142:   MPI_Comm       comm;
143:   UserCtx        *user = (UserCtx*)appctx;
144:   UserCtx_Power  appctx_power = (*user).appctx_power;
145:   AppCtx_Water   appctx_water = (*user).appctx_water;

148:   SNESGetDM(snes,&networkdm);
149:   PetscObjectGetComm((PetscObject)networkdm,&comm);
150:   MPI_Comm_rank(comm,&rank);

152:   DMGetLocalVector(networkdm,&localX);
153:   DMGetLocalVector(networkdm,&localF);
154:   VecSet(F,0.0);
155:   VecSet(localF,0.0);

157:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
158:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

160:   /* Form Function for power subnetwork */
161:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
162:   if (user->subsnes_id == 1) { /* snes_water only */
163:     FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
164:   } else {
165:     FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,&appctx_power);
166:   }

168:   /* Form Function for water subnetwork */
169:   DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
170:   if (user->subsnes_id == 0) { /* snes_power only */
171:     FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
172:   } else {
173:     FormFunction_Water(networkdm,localX,localF,nv,ne,vtx,edges,NULL);
174:   }

176:   /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */
177:   DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
178:   for (v=0; v<nv; v++) {
179:     PetscInt       key,ncomp,nvar,nconnedges,k,e,keye,goffset[3];
180:     void*          component;
181:     const PetscInt *connedges;

183:     DMNetworkGetComponent(networkdm,vtx[v],ALL_COMPONENTS,NULL,NULL,&nvar);
184:     DMNetworkGetNumComponents(networkdm,vtx[v],&ncomp);
185:     /* printf("  [%d] coupling vertex[%D]: v %D, ncomp %D; nvar %D\n",rank,v,vtx[v], ncomp,nvar); */

187:     for (k=0; k<ncomp; k++) {
188:       DMNetworkGetComponent(networkdm,vtx[v],k,&key,&component,&nvar);
189:       DMNetworkGetGlobalVecOffset(networkdm,vtx[v],k,&goffset[k]);

191:       /* Verify the coupling vertex is a powernet load vertex or a water vertex */
192:       switch (k) {
193:       case 0:
194:         if (key != appctx_power.compkey_bus || nvar != 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"key %D not a power bus vertex or nvar %D != 2",key,nvar);
195:         break;
196:       case 1:
197:         if (key != appctx_power.compkey_load || nvar != 0 || goffset[1] != goffset[0]+2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power load vertex");
198:         break;
199:       case 2:
200:         if (key != appctx_water.compkey_vtx || nvar != 1 || goffset[2] != goffset[1]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex");
201:         break;
202:       default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %D is wrong",k);
203:       }
204:       /* printf("  [%d] coupling vertex[%D]: key %D; nvar %D, goffset %D\n",rank,v,key,nvar,goffset[k]); */
205:     }

207:     /* Get its supporting edges */
208:     DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
209:     /* printf("\n[%d] coupling vertex: nconnedges %D\n",rank,nconnedges); */
210:     for (k=0; k<nconnedges; k++) {
211:       e = connedges[k];
212:       DMNetworkGetNumComponents(networkdm,e,&ncomp);
213:       /* printf("\n  [%d] connected edge[%D]=%D has ncomp %D\n",rank,k,e,ncomp); */
214:       DMNetworkGetComponent(networkdm,e,0,&keye,&component,NULL);
215:       if (keye == appctx_water.compkey_edge) { /* water_compkey_edge */
216:         EDGE_Water        edge=(EDGE_Water)component;
217:         if (edge->type == EDGE_TYPE_PUMP) {
218:           /* printf("  connected edge[%D]=%D has keye=%D, is appctx_water.compkey_edge with EDGE_TYPE_PUMP\n",k,e,keye); */
219:         }
220:       } else { /* ower->compkey_branch */
221:         if (keye != appctx_power.compkey_branch) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power branch");
222:       }
223:     }
224:   }

226:   DMRestoreLocalVector(networkdm,&localX);

228:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
229:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
230:   DMRestoreLocalVector(networkdm,&localF);
231: #if 0
232:   if (!rank) printf("F:\n");
233:   VecView(F,PETSC_VIEWER_STDOUT_WORLD);
234: #endif
235:   return(0);
236: }

238: PetscErrorCode SetInitialGuess(DM networkdm,Vec X,void* appctx)
239: {
241:   PetscInt       nv,ne,i,j,ncomp,offset,key;
242:   const PetscInt *vtx,*edges;
243:   UserCtx        *user = (UserCtx*)appctx;
244:   Vec            localX = user->localXold;
245:   UserCtx_Power  appctx_power = (*user).appctx_power;
246:   AppCtx_Water   appctx_water = (*user).appctx_water;
247:   PetscBool      ghost;
248:   PetscScalar    *xarr;
249:   VERTEX_Power   bus;
250:   VERTEX_Water   vertex;
251:   void*          component;
252:   GEN            gen;

255:   VecSet(X,0.0);
256:   VecSet(localX,0.0);

258:   /* Set initial guess for power subnetwork */
259:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
260:   SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,&appctx_power);

262:   /* Set initial guess for water subnetwork */
263:   DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
264:   SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);

266:   /* Set initial guess at the coupling vertex */
267:   VecGetArray(localX,&xarr);
268:   DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
269:   for (i=0; i<nv; i++) {
270:     DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);
271:     if (ghost) continue;

273:     DMNetworkGetNumComponents(networkdm,vtx[i],&ncomp);
274:     for (j=0; j<ncomp; j++) {
275:       DMNetworkGetLocalVecOffset(networkdm,vtx[i],j,&offset);
276:       DMNetworkGetComponent(networkdm,vtx[i],j,&key,(void**)&component,NULL);
277:       if (key == appctx_power.compkey_bus) {
278:         bus = (VERTEX_Power)(component);
279:         xarr[offset]   = bus->va*PETSC_PI/180.0;
280:         xarr[offset+1] = bus->vm;
281:       } else if (key == appctx_power.compkey_gen) {
282:         gen = (GEN)(component);
283:         if (!gen->status) continue;
284:         xarr[offset+1] = gen->vs;
285:       } else if (key == appctx_water.compkey_vtx) {
286:         vertex = (VERTEX_Water)(component);
287:         if (vertex->type == VERTEX_TYPE_JUNCTION) {
288:           xarr[offset] = 100;
289:         } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
290:           xarr[offset] = vertex->res.head;
291:         } else {
292:           xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
293:         }
294:       }
295:     }
296:   }
297:   VecRestoreArray(localX,&xarr);

299:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
300:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
301:   return(0);
302: }

304: int main(int argc,char **argv)
305: {
306:   PetscErrorCode   ierr;
307:   DM               networkdm;
308:   PetscLogStage    stage[4];
309:   PetscMPIInt      rank,size;
310:   PetscInt         Nsubnet=2,numVertices[2],numEdges[2],i,j,nv,ne,it_max=10;
311:   const PetscInt   *vtx,*edges;
312:   Vec              X,F;
313:   SNES             snes,snes_power,snes_water;
314:   Mat              Jac;
315:   PetscBool        ghost,viewJ=PETSC_FALSE,viewX=PETSC_FALSE,viewDM=PETSC_FALSE,test=PETSC_FALSE,distribute=PETSC_TRUE,flg;
316:   UserCtx          user;
317:   SNESConvergedReason reason;

319:   /* Power subnetwork */
320:   UserCtx_Power       *appctx_power  = &user.appctx_power;
321:   char                pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
322:   PFDATA              *pfdata = NULL;
323:   PetscInt            genj,loadj,*edgelist_power = NULL,power_netnum;
324:   PetscScalar         Sbase = 0.0;

326:   /* Water subnetwork */
327:   AppCtx_Water        *appctx_water = &user.appctx_water;
328:   WATERDATA           *waterdata = NULL;
329:   char                waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
330:   PetscInt            *edgelist_water = NULL,water_netnum;

332:   /* Shared vertices between subnetworks */
333:   PetscInt           power_svtx,water_svtx;

335:   PetscInitialize(&argc,&argv,"ex1options",help);if (ierr) return ierr;
336:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
337:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

339:   /* (1) Read Data - Only rank 0 reads the data */
340:   /*--------------------------------------------*/
341:   PetscLogStageRegister("Read Data",&stage[0]);
342:   PetscLogStagePush(stage[0]);

344:   for (i=0; i<Nsubnet; i++) {
345:     numVertices[i] = 0;
346:     numEdges[i]    = 0;
347:   }

349:   /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
350:   /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
351:   PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
352:   PetscNew(&pfdata);
353:   PFReadMatPowerData(pfdata,pfdata_file);
354:   if (!rank) {
355:     PetscPrintf(PETSC_COMM_SELF,"Power network: nb = %D, ngen = %D, nload = %D, nbranch = %D\n",pfdata->nbus,pfdata->ngen,pfdata->nload,pfdata->nbranch);
356:   }
357:   Sbase = pfdata->sbase;
358:   if (rank == 0) { /* proc[0] will create Electric Power Grid */
359:     numEdges[0]    = pfdata->nbranch;
360:     numVertices[0] = pfdata->nbus;

362:     PetscMalloc1(2*numEdges[0],&edgelist_power);
363:     GetListofEdges_Power(pfdata,edgelist_power);
364:   }
365:   /* Broadcast power Sbase to all processors */
366:   MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
367:   appctx_power->Sbase     = Sbase;
368:   appctx_power->jac_error = PETSC_FALSE;
369:   /* If external option activated. Introduce error in jacobian */
370:   PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);

372:   /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
373:   /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
374:   PetscNew(&waterdata);
375:   PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,PETSC_MAX_PATH_LEN-1,NULL);
376:   WaterReadData(waterdata,waterdata_file);
377:   if (size == 1 || (size > 1 && rank == 1)) {
378:     PetscCalloc1(2*waterdata->nedge,&edgelist_water);
379:     GetListofEdges_Water(waterdata,edgelist_water);
380:     numEdges[1]    = waterdata->nedge;
381:     numVertices[1] = waterdata->nvertex;
382:   }
383:   PetscLogStagePop();

385:   /* (2) Create a network consist of two subnetworks */
386:   /*-------------------------------------------------*/
387:   PetscLogStageRegister("Net Setup",&stage[1]);
388:   PetscLogStagePush(stage[1]);

390:   PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);
391:   PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);
392:   PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);

394:   /* Create an empty network object */
395:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);

397:   /* Register the components in the network */
398:   DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);
399:   DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);
400:   DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);
401:   DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);

403:   DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);
404:   DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);
405: #if 0
406:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_branch %d\n",appctx_power->compkey_branch);
407:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_bus    %d\n",appctx_power->compkey_bus);
408:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_gen    %d\n",appctx_power->compkey_gen);
409:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_load   %d\n",appctx_power->compkey_load);
410:   PetscPrintf(PETSC_COMM_WORLD,"water->compkey_edge   %d\n",appctx_water->compkey_edge);
411:   PetscPrintf(PETSC_COMM_WORLD,"water->compkey_vtx    %d\n",appctx_water->compkey_vtx);
412: #endif
413:   PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Total local nvertices %D + %D = %D, nedges %D + %D = %D\n",rank,numVertices[0],numVertices[1],numVertices[0]+numVertices[1],numEdges[0],numEdges[1],numEdges[0]+numEdges[1]);
414:   PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);

416:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,Nsubnet);
417:   DMNetworkAddSubnetwork(networkdm,"power",numVertices[0],numEdges[0],edgelist_power,&power_netnum);
418:   DMNetworkAddSubnetwork(networkdm,"water",numVertices[1],numEdges[1],edgelist_water,&water_netnum);

420:   /* vertex subnet[0].4 shares with vertex subnet[1].0 */
421:   power_svtx = 4; water_svtx = 0;
422:   DMNetworkAddSharedVertices(networkdm,power_netnum,water_netnum,1,&power_svtx,&water_svtx);

424:   /* Set up the network layout */
425:   DMNetworkLayoutSetUp(networkdm);

427:   /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
428:   /*-------------------------------------------------------*/
429:   genj = 0; loadj = 0;
430:   DMNetworkGetSubnetwork(networkdm,power_netnum,&nv,&ne,&vtx,&edges);

432:   for (i = 0; i < ne; i++) {
433:     DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i],0);
434:   }

436:   for (i = 0; i < nv; i++) {
437:     DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
438:     if (flg) continue;

440:     DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i],2);
441:     if (pfdata->bus[i].ngen) {
442:       for (j = 0; j < pfdata->bus[i].ngen; j++) {
443:         DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++],0);
444:       }
445:     }
446:     if (pfdata->bus[i].nload) {
447:       for (j=0; j < pfdata->bus[i].nload; j++) {
448:         DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++],0);
449:       }
450:     }
451:   }

453:   /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
454:   /*-------------------------------------------------------*/
455:   DMNetworkGetSubnetwork(networkdm,water_netnum,&nv,&ne,&vtx,&edges);
456:   for (i = 0; i < ne; i++) {
457:     DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i],0);
458:   }

460:   for (i = 0; i < nv; i++) {
461:     DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
462:     if (flg) continue;

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

467:   /* ADD VARIABLES AND COMPONENTS AT THE SHARED VERTEX: net[0].4 coupls with net[1].0 -- only the owner of the vertex does this */
468:   /*----------------------------------------------------------------------------------------------------------------------------*/
469:   DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
470:   for (i = 0; i < nv; i++) {
471:     DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);
472:     /* printf("[%d] coupling info: nv %d; sv[0] %d; ghost %d\n",rank,nv,vtx[0],ghost); */
473:     if (ghost) continue;

475:     /* power */
476:     DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[4],2);
477:     /* bus[4] is a load, add its component */
478:     DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[0],0);

480:     /* water */
481:     DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[0],1);
482:   }

484:   /* Set up DM for use */
485:   DMSetUp(networkdm);
486:   if (viewDM) {
487:     PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
488:     DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
489:   }

491:   /* Free user objects */
492:   PetscFree(edgelist_power);
493:   PetscFree(pfdata->bus);
494:   PetscFree(pfdata->gen);
495:   PetscFree(pfdata->branch);
496:   PetscFree(pfdata->load);
497:   PetscFree(pfdata);

499:   PetscFree(edgelist_water);
500:   PetscFree(waterdata->vertex);
501:   PetscFree(waterdata->edge);
502:   PetscFree(waterdata);

504:   /* Re-distribute networkdm to multiple processes for better job balance */
505:   if (size >1 && distribute) {
506:     DMNetworkDistribute(&networkdm,0);
507:     if (viewDM) {
508:       PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
509:       DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
510:     }
511:   }

513:   /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
514:   if (test) {
515:     PetscInt  v,gidx;
516:     MPI_Barrier(PETSC_COMM_WORLD);
517:     for (i=0; i<Nsubnet; i++) {
518:       DMNetworkGetSubnetwork(networkdm,i,&nv,&ne,&vtx,&edges);
519:       PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, subnet[%d] ne %d, nv %d\n",rank,i,ne,nv);
520:       MPI_Barrier(PETSC_COMM_WORLD);

522:       for (v=0; v<nv; v++) {
523:         DMNetworkIsGhostVertex(networkdm,vtx[v],&ghost);
524:         DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
525:         PetscPrintf(PETSC_COMM_SELF,"[%d] subnet[%d] v %d %d; ghost %d\n",rank,i,vtx[v],gidx,ghost);
526:       }
527:       MPI_Barrier(PETSC_COMM_WORLD);
528:     }
529:     MPI_Barrier(PETSC_COMM_WORLD);

531:     DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
532:     PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, num of shared vertices nsv = %d\n",rank,nv);
533:     for (v=0; v<nv; v++) {
534:       DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
535:       PetscPrintf(PETSC_COMM_SELF,"[%d] sv %d, gidx=%d\n",rank,vtx[v],gidx);
536:     }
537:     MPI_Barrier(PETSC_COMM_WORLD);
538:   }

540:   /* Create solution vector X */
541:   DMCreateGlobalVector(networkdm,&X);
542:   VecDuplicate(X,&F);
543:   DMGetLocalVector(networkdm,&user.localXold);
544:   PetscLogStagePop();

546:   /* (3) Setup Solvers */
547:   /*-------------------*/
548:   PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);
549:   PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);

551:   PetscLogStageRegister("SNES Setup",&stage[2]);
552:   PetscLogStagePush(stage[2]);

554:   SetInitialGuess(networkdm,X,&user);

556:   /* Create coupled snes */
557:   /*-------------------- */
558:   PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");
559:   user.subsnes_id = Nsubnet;
560:   SNESCreate(PETSC_COMM_WORLD,&snes);
561:   SNESSetDM(snes,networkdm);
562:   SNESSetOptionsPrefix(snes,"coupled_");
563:   SNESSetFunction(snes,F,FormFunction,&user);
564:   SNESMonitorSet(snes,UserMonitor,&user,NULL);
565:   SNESSetFromOptions(snes);

567:   if (viewJ) {
568:     /* View Jac structure */
569:     SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
570:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
571:   }

573:   if (viewX) {
574:     PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");
575:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
576:   }

578:   if (viewJ) {
579:     /* View assembled Jac */
580:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
581:   }

583:   /* Create snes_power */
584:   /*-------------------*/
585:   PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");

587:   user.subsnes_id = 0;
588:   SNESCreate(PETSC_COMM_WORLD,&snes_power);
589:   SNESSetDM(snes_power,networkdm);
590:   SNESSetOptionsPrefix(snes_power,"power_");
591:   SNESSetFunction(snes_power,F,FormFunction,&user);
592:   SNESMonitorSet(snes_power,UserMonitor,&user,NULL);

594:   /* Use user-provide Jacobian */
595:   DMCreateMatrix(networkdm,&Jac);
596:   SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);
597:   SNESSetFromOptions(snes_power);

599:   if (viewX) {
600:     PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");
601:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
602:   }
603:   if (viewJ) {
604:     PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");
605:     SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);
606:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
607:     /* MatView(Jac,PETSC_VIEWER_STDOUT_WORLD); */
608:   }

610:   /* Create snes_water */
611:   /*-------------------*/
612:   PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");

614:   user.subsnes_id = 1;
615:   SNESCreate(PETSC_COMM_WORLD,&snes_water);
616:   SNESSetDM(snes_water,networkdm);
617:   SNESSetOptionsPrefix(snes_water,"water_");
618:   SNESSetFunction(snes_water,F,FormFunction,&user);
619:   SNESMonitorSet(snes_water,UserMonitor,&user,NULL);
620:   SNESSetFromOptions(snes_water);

622:   if (viewX) {
623:     PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");
624:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
625:   }
626:   PetscLogStagePop();

628:   /* (4) Solve */
629:   /*-----------*/
630:   PetscLogStageRegister("SNES Solve",&stage[3]);
631:   PetscLogStagePush(stage[3]);
632:   user.it = 0;
633:   reason  = SNES_DIVERGED_DTOL;
634:   while (user.it < it_max && (PetscInt)reason<0) {
635: #if 0
636:     user.subsnes_id = 0;
637:     SNESSolve(snes_power,NULL,X);

639:     user.subsnes_id = 1;
640:     SNESSolve(snes_water,NULL,X);
641: #endif
642:     user.subsnes_id = Nsubnet;
643:     SNESSolve(snes,NULL,X);

645:     SNESGetConvergedReason(snes,&reason);
646:     user.it++;
647:   }
648:   PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);
649:   if (viewX) {
650:     PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");
651:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
652:   }
653:   PetscLogStagePop();

655:   /* Free objects */
656:   /* -------------*/
657:   VecDestroy(&X);
658:   VecDestroy(&F);
659:   DMRestoreLocalVector(networkdm,&user.localXold);

661:   SNESDestroy(&snes);
662:   MatDestroy(&Jac);
663:   SNESDestroy(&snes_power);
664:   SNESDestroy(&snes_water);

666:   DMDestroy(&networkdm);
667:   PetscFinalize();
668:   return ierr;
669: }

671: /*TEST

673:    build:
674:      requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
675:      depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c

677:    test:
678:       args: -coupled_snes_converged_reason -options_left no -viewDM
679:       localrunfiles: ex1options power/case9.m water/sample1.inp
680:       output_file: output/ex1.out

682:    test:
683:       suffix: 2
684:       nsize: 3
685:       args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis
686:       localrunfiles: ex1options power/case9.m water/sample1.inp
687:       output_file: output/ex1_2.out
688:       requires: parmetis

690: #   test:
691: #      suffix: 3
692: #      nsize: 3
693: #      args: -coupled_snes_converged_reason -options_left no -distribute false
694: #      localrunfiles: ex1options power/case9.m water/sample1.inp
695: #      output_file: output/ex1_2.out

697:    test:
698:       suffix: 4
699:       nsize: 4
700:       args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -viewDM
701:       localrunfiles: ex1options power/case9.m water/sample1.inp
702:       output_file: output/ex1_4.out

704: TEST*/