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: {
 26:   UserCtx        *user = (UserCtx*)appctx;
 27:   Vec            X,localXold = user->localXold;
 28:   DM             networkdm;
 29:   PetscMPIInt    rank;
 30:   MPI_Comm       comm;

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

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

 62:   SNESGetDM(snes,&networkdm);
 63:   DMGetLocalVector(networkdm,&localX);

 65:   PetscObjectGetComm((PetscObject)networkdm,&comm);
 66:   MPI_Comm_rank(comm,&rank);

 68:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 69:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 71:   MatZeroEntries(J);

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

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

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

 93:   DMRestoreLocalVector(networkdm,&localX);
 94:   return 0;
 95: }

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

107:   VecGetArrayRead(localX,&xarr);
108:   VecGetArrayRead(localXold,&xoldarr);
109:   VecGetArray(localF,&farr);

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

115:     DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
116:     DMNetworkGetComponent(networkdm,vtx[i],ALL_COMPONENTS,NULL,NULL,&nvar);
117:     for (j=0; j<nvar; j++) {
118:       farr[offset+j] = xarr[offset+j] - xoldarr[offset+j];
119:     }
120:   }

122:   VecRestoreArrayRead(localX,&xarr);
123:   VecRestoreArrayRead(localXold,&xoldarr);
124:   VecRestoreArray(localF,&farr);
125:   return 0;
126: }

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

140:   SNESGetDM(snes,&networkdm);
141:   PetscObjectGetComm((PetscObject)networkdm,&comm);
142:   MPI_Comm_rank(comm,&rank);

144:   DMGetLocalVector(networkdm,&localX);
145:   DMGetLocalVector(networkdm,&localF);
146:   VecSet(F,0.0);
147:   VecSet(localF,0.0);

149:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
150:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

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

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

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

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

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

183:       /* Verify the coupling vertex is a powernet load vertex or a water vertex */
184:       switch (k) {
185:       case 0:
187:         break;
188:       case 1:
190:         break;
191:       case 2:
193:         break;
194:       default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %D is wrong",k);
195:       }
196:       /* printf("  [%d] coupling vertex[%D]: key %D; nvar %D, goffset %D\n",rank,v,key,nvar,goffset[k]); */
197:     }

199:     /* Get its supporting edges */
200:     DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
201:     /* printf("\n[%d] coupling vertex: nconnedges %D\n",rank,nconnedges); */
202:     for (k=0; k<nconnedges; k++) {
203:       e = connedges[k];
204:       DMNetworkGetNumComponents(networkdm,e,&ncomp);
205:       /* printf("\n  [%d] connected edge[%D]=%D has ncomp %D\n",rank,k,e,ncomp); */
206:       DMNetworkGetComponent(networkdm,e,0,&keye,&component,NULL);
207:       if (keye == appctx_water.compkey_edge) { /* water_compkey_edge */
208:         EDGE_Water        edge=(EDGE_Water)component;
209:         if (edge->type == EDGE_TYPE_PUMP) {
210:           /* printf("  connected edge[%D]=%D has keye=%D, is appctx_water.compkey_edge with EDGE_TYPE_PUMP\n",k,e,keye); */
211:         }
212:       } else { /* ower->compkey_branch */
214:       }
215:     }
216:   }

218:   DMRestoreLocalVector(networkdm,&localX);

220:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
221:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
222:   DMRestoreLocalVector(networkdm,&localF);
223: #if 0
224:   if (rank == 0) printf("F:\n");
225:   VecView(F,PETSC_VIEWER_STDOUT_WORLD);
226: #endif
227:   return 0;
228: }

230: PetscErrorCode SetInitialGuess(DM networkdm,Vec X,void* appctx)
231: {
232:   PetscInt       nv,ne,i,j,ncomp,offset,key;
233:   const PetscInt *vtx,*edges;
234:   UserCtx        *user = (UserCtx*)appctx;
235:   Vec            localX = user->localXold;
236:   UserCtx_Power  appctx_power = (*user).appctx_power;
237:   AppCtx_Water   appctx_water = (*user).appctx_water;
238:   PetscBool      ghost;
239:   PetscScalar    *xarr;
240:   VERTEX_Power   bus;
241:   VERTEX_Water   vertex;
242:   void*          component;
243:   GEN            gen;

245:   VecSet(X,0.0);
246:   VecSet(localX,0.0);

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

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

256:   /* Set initial guess at the coupling vertex */
257:   VecGetArray(localX,&xarr);
258:   DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
259:   for (i=0; i<nv; i++) {
260:     DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);
261:     if (ghost) continue;

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

289:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
290:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
291:   return 0;
292: }

294: int main(int argc,char **argv)
295: {
296:   DM               networkdm;
297:   PetscLogStage    stage[4];
298:   PetscMPIInt      rank,size;
299:   PetscInt         Nsubnet=2,numVertices[2],numEdges[2],i,j,nv,ne,it_max=10;
300:   const PetscInt   *vtx,*edges;
301:   Vec              X,F;
302:   SNES             snes,snes_power,snes_water;
303:   Mat              Jac;
304:   PetscBool        ghost,viewJ=PETSC_FALSE,viewX=PETSC_FALSE,viewDM=PETSC_FALSE,test=PETSC_FALSE,distribute=PETSC_TRUE,flg;
305:   UserCtx          user;
306:   SNESConvergedReason reason;

308:   /* Power subnetwork */
309:   UserCtx_Power       *appctx_power  = &user.appctx_power;
310:   char                pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
311:   PFDATA              *pfdata = NULL;
312:   PetscInt            genj,loadj,*edgelist_power = NULL,power_netnum;
313:   PetscScalar         Sbase = 0.0;

315:   /* Water subnetwork */
316:   AppCtx_Water        *appctx_water = &user.appctx_water;
317:   WATERDATA           *waterdata = NULL;
318:   char                waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
319:   PetscInt            *edgelist_water = NULL,water_netnum;

321:   /* Shared vertices between subnetworks */
322:   PetscInt           power_svtx,water_svtx;

324:   PetscInitialize(&argc,&argv,"ex1options",help);
325:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
326:   MPI_Comm_size(PETSC_COMM_WORLD,&size);

328:   /* (1) Read Data - Only rank 0 reads the data */
329:   /*--------------------------------------------*/
330:   PetscLogStageRegister("Read Data",&stage[0]);
331:   PetscLogStagePush(stage[0]);

333:   for (i=0; i<Nsubnet; i++) {
334:     numVertices[i] = 0;
335:     numEdges[i]    = 0;
336:   }

338:   /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
339:   /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
340:   PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
341:   PetscNew(&pfdata);
342:   PFReadMatPowerData(pfdata,pfdata_file);
343:   if (rank == 0) {
344:     PetscPrintf(PETSC_COMM_SELF,"Power network: nb = %D, ngen = %D, nload = %D, nbranch = %D\n",pfdata->nbus,pfdata->ngen,pfdata->nload,pfdata->nbranch);
345:   }
346:   Sbase = pfdata->sbase;
347:   if (rank == 0) { /* proc[0] will create Electric Power Grid */
348:     numEdges[0]    = pfdata->nbranch;
349:     numVertices[0] = pfdata->nbus;

351:     PetscMalloc1(2*numEdges[0],&edgelist_power);
352:     GetListofEdges_Power(pfdata,edgelist_power);
353:   }
354:   /* Broadcast power Sbase to all processors */
355:   MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
356:   appctx_power->Sbase     = Sbase;
357:   appctx_power->jac_error = PETSC_FALSE;
358:   /* If external option activated. Introduce error in jacobian */
359:   PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);

361:   /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
362:   /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
363:   PetscNew(&waterdata);
364:   PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,PETSC_MAX_PATH_LEN-1,NULL);
365:   WaterReadData(waterdata,waterdata_file);
366:   if (size == 1 || (size > 1 && rank == 1)) {
367:     PetscCalloc1(2*waterdata->nedge,&edgelist_water);
368:     GetListofEdges_Water(waterdata,edgelist_water);
369:     numEdges[1]    = waterdata->nedge;
370:     numVertices[1] = waterdata->nvertex;
371:   }
372:   PetscLogStagePop();

374:   /* (2) Create a network consist of two subnetworks */
375:   /*-------------------------------------------------*/
376:   PetscLogStageRegister("Net Setup",&stage[1]);
377:   PetscLogStagePush(stage[1]);

379:   PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);
380:   PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);
381:   PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);

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

386:   /* Register the components in the network */
387:   DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);
388:   DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);
389:   DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);
390:   DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);

392:   DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);
393:   DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);
394: #if 0
395:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_branch %d\n",appctx_power->compkey_branch);
396:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_bus    %d\n",appctx_power->compkey_bus);
397:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_gen    %d\n",appctx_power->compkey_gen);
398:   PetscPrintf(PETSC_COMM_WORLD,"power->compkey_load   %d\n",appctx_power->compkey_load);
399:   PetscPrintf(PETSC_COMM_WORLD,"water->compkey_edge   %d\n",appctx_water->compkey_edge);
400:   PetscPrintf(PETSC_COMM_WORLD,"water->compkey_vtx    %d\n",appctx_water->compkey_vtx);
401: #endif
402:   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]);
403:   PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);

405:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,Nsubnet);
406:   DMNetworkAddSubnetwork(networkdm,"power",numEdges[0],edgelist_power,&power_netnum);
407:   DMNetworkAddSubnetwork(networkdm,"water",numEdges[1],edgelist_water,&water_netnum);

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

413:   /* Set up the network layout */
414:   DMNetworkLayoutSetUp(networkdm);

416:   /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
417:   /*-------------------------------------------------------*/
418:   genj = 0; loadj = 0;
419:   DMNetworkGetSubnetwork(networkdm,power_netnum,&nv,&ne,&vtx,&edges);

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

425:   for (i = 0; i < nv; i++) {
426:     DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
427:     if (flg) continue;

429:     DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i],2);
430:     if (pfdata->bus[i].ngen) {
431:       for (j = 0; j < pfdata->bus[i].ngen; j++) {
432:         DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++],0);
433:       }
434:     }
435:     if (pfdata->bus[i].nload) {
436:       for (j=0; j < pfdata->bus[i].nload; j++) {
437:         DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++],0);
438:       }
439:     }
440:   }

442:   /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
443:   /*-------------------------------------------------------*/
444:   DMNetworkGetSubnetwork(networkdm,water_netnum,&nv,&ne,&vtx,&edges);
445:   for (i = 0; i < ne; i++) {
446:     DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i],0);
447:   }

449:   for (i = 0; i < nv; i++) {
450:     DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
451:     if (flg) continue;

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

456:   /* 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 */
457:   /*----------------------------------------------------------------------------------------------------------------------------*/
458:   DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
459:   for (i = 0; i < nv; i++) {
460:     /* power */
461:     DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[4],2);
462:     /* bus[4] is a load, add its component */
463:     DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[0],0);

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

469:   /* Set up DM for use */
470:   DMSetUp(networkdm);
471:   if (viewDM) {
472:     PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
473:     DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
474:   }

476:   /* Free user objects */
477:   PetscFree(edgelist_power);
478:   PetscFree(pfdata->bus);
479:   PetscFree(pfdata->gen);
480:   PetscFree(pfdata->branch);
481:   PetscFree(pfdata->load);
482:   PetscFree(pfdata);

484:   PetscFree(edgelist_water);
485:   PetscFree(waterdata->vertex);
486:   PetscFree(waterdata->edge);
487:   PetscFree(waterdata);

489:   /* Re-distribute networkdm to multiple processes for better job balance */
490:   if (size >1 && distribute) {
491:     DMNetworkDistribute(&networkdm,0);
492:     if (viewDM) {
493:       PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
494:       DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
495:     }
496:   }

498:   /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
499:   if (test) {
500:     PetscInt  v,gidx;
501:     MPI_Barrier(PETSC_COMM_WORLD);
502:     for (i=0; i<Nsubnet; i++) {
503:       DMNetworkGetSubnetwork(networkdm,i,&nv,&ne,&vtx,&edges);
504:       PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, subnet[%d] ne %d, nv %d\n",rank,i,ne,nv);
505:       MPI_Barrier(PETSC_COMM_WORLD);

507:       for (v=0; v<nv; v++) {
508:         DMNetworkIsGhostVertex(networkdm,vtx[v],&ghost);
509:         DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
510:         PetscPrintf(PETSC_COMM_SELF,"[%d] subnet[%d] v %d %d; ghost %d\n",rank,i,vtx[v],gidx,ghost);
511:       }
512:       MPI_Barrier(PETSC_COMM_WORLD);
513:     }
514:     MPI_Barrier(PETSC_COMM_WORLD);

516:     DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
517:     PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, num of shared vertices nsv = %d\n",rank,nv);
518:     for (v=0; v<nv; v++) {
519:       DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
520:       PetscPrintf(PETSC_COMM_SELF,"[%d] sv %d, gidx=%d\n",rank,vtx[v],gidx);
521:     }
522:     MPI_Barrier(PETSC_COMM_WORLD);
523:   }

525:   /* Create solution vector X */
526:   DMCreateGlobalVector(networkdm,&X);
527:   VecDuplicate(X,&F);
528:   DMGetLocalVector(networkdm,&user.localXold);
529:   PetscLogStagePop();

531:   /* (3) Setup Solvers */
532:   /*-------------------*/
533:   PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);
534:   PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);

536:   PetscLogStageRegister("SNES Setup",&stage[2]);
537:   PetscLogStagePush(stage[2]);

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

541:   /* Create coupled snes */
542:   /*-------------------- */
543:   PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");
544:   user.subsnes_id = Nsubnet;
545:   SNESCreate(PETSC_COMM_WORLD,&snes);
546:   SNESSetDM(snes,networkdm);
547:   SNESSetOptionsPrefix(snes,"coupled_");
548:   SNESSetFunction(snes,F,FormFunction,&user);
549:   SNESMonitorSet(snes,UserMonitor,&user,NULL);
550:   SNESSetFromOptions(snes);

552:   if (viewJ) {
553:     /* View Jac structure */
554:     SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
555:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
556:   }

558:   if (viewX) {
559:     PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");
560:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
561:   }

563:   if (viewJ) {
564:     /* View assembled Jac */
565:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
566:   }

568:   /* Create snes_power */
569:   /*-------------------*/
570:   PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");

572:   user.subsnes_id = 0;
573:   SNESCreate(PETSC_COMM_WORLD,&snes_power);
574:   SNESSetDM(snes_power,networkdm);
575:   SNESSetOptionsPrefix(snes_power,"power_");
576:   SNESSetFunction(snes_power,F,FormFunction,&user);
577:   SNESMonitorSet(snes_power,UserMonitor,&user,NULL);

579:   /* Use user-provide Jacobian */
580:   DMCreateMatrix(networkdm,&Jac);
581:   SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);
582:   SNESSetFromOptions(snes_power);

584:   if (viewX) {
585:     PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");
586:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
587:   }
588:   if (viewJ) {
589:     PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");
590:     SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);
591:     MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
592:     /* MatView(Jac,PETSC_VIEWER_STDOUT_WORLD); */
593:   }

595:   /* Create snes_water */
596:   /*-------------------*/
597:   PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");

599:   user.subsnes_id = 1;
600:   SNESCreate(PETSC_COMM_WORLD,&snes_water);
601:   SNESSetDM(snes_water,networkdm);
602:   SNESSetOptionsPrefix(snes_water,"water_");
603:   SNESSetFunction(snes_water,F,FormFunction,&user);
604:   SNESMonitorSet(snes_water,UserMonitor,&user,NULL);
605:   SNESSetFromOptions(snes_water);

607:   if (viewX) {
608:     PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");
609:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
610:   }
611:   PetscLogStagePop();

613:   /* (4) Solve */
614:   /*-----------*/
615:   PetscLogStageRegister("SNES Solve",&stage[3]);
616:   PetscLogStagePush(stage[3]);
617:   user.it = 0;
618:   reason  = SNES_DIVERGED_DTOL;
619:   while (user.it < it_max && (PetscInt)reason<0) {
620: #if 0
621:     user.subsnes_id = 0;
622:     SNESSolve(snes_power,NULL,X);

624:     user.subsnes_id = 1;
625:     SNESSolve(snes_water,NULL,X);
626: #endif
627:     user.subsnes_id = Nsubnet;
628:     SNESSolve(snes,NULL,X);

630:     SNESGetConvergedReason(snes,&reason);
631:     user.it++;
632:   }
633:   PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);
634:   if (viewX) {
635:     PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");
636:     VecView(X,PETSC_VIEWER_STDOUT_WORLD);
637:   }
638:   PetscLogStagePop();

640:   /* Free objects */
641:   /* -------------*/
642:   VecDestroy(&X);
643:   VecDestroy(&F);
644:   DMRestoreLocalVector(networkdm,&user.localXold);

646:   SNESDestroy(&snes);
647:   MatDestroy(&Jac);
648:   SNESDestroy(&snes_power);
649:   SNESDestroy(&snes_water);

651:   DMDestroy(&networkdm);
652:   PetscFinalize();
653:   return 0;
654: }

656: /*TEST

658:    build:
659:      requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
660:      depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c

662:    test:
663:       args: -coupled_snes_converged_reason -options_left no -viewDM
664:       localrunfiles: ex1options power/case9.m water/sample1.inp
665:       output_file: output/ex1.out

667:    test:
668:       suffix: 2
669:       nsize: 3
670:       args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis
671:       localrunfiles: ex1options power/case9.m water/sample1.inp
672:       output_file: output/ex1_2.out
673:       requires: parmetis

675: #   test:
676: #      suffix: 3
677: #      nsize: 3
678: #      args: -coupled_snes_converged_reason -options_left no -distribute false
679: #      localrunfiles: ex1options power/case9.m water/sample1.inp
680: #      output_file: output/ex1_2.out

682:    test:
683:       suffix: 4
684:       nsize: 4
685:       args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -viewDM
686:       localrunfiles: ex1options power/case9.m water/sample1.inp
687:       output_file: output/ex1_4.out

689: TEST*/