Actual source code: power2.c

  1: static char help[] = "This example demonstrates the use of DMNetwork interface with subnetworks for solving a nonlinear electric power grid problem.\n\
  2:                       The available solver options are in the poweroptions file and the data files are in the datafiles directory.\n\
  3:                       The data file format used is from the MatPower package (http://www.pserc.cornell.edu//matpower/).\n\
  4:                       This example shows the use of subnetwork feature in DMNetwork. It creates duplicates of the same network which are treated as subnetworks.\n\
  5:                       Run this program: mpiexec -n <n> ./pf2\n\
  6:                       mpiexec -n <n> ./pf2 \n";

  8: /* T
  9:    Concepts: DMNetwork
 10:    Concepts: PETSc SNES solver
 11: */

 13: #include "power.h"
 14: #include <petscdmnetwork.h>

 16: PetscErrorCode FormFunction_Subnet(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
 17: {
 18:   UserCtx_Power     *User = (UserCtx_Power*)appctx;
 19:   PetscInt          e,v,vfrom,vto;
 20:   const PetscScalar *xarr;
 21:   PetscScalar       *farr;
 22:   PetscInt          offsetfrom,offsetto,offset;

 24:   VecGetArrayRead(localX,&xarr);
 25:   VecGetArray(localF,&farr);

 27:   for (v=0; v<nv; v++) {
 28:     PetscInt      i,j,key;
 29:     PetscScalar   Vm;
 30:     PetscScalar   Sbase = User->Sbase;
 31:     VERTEX_Power  bus = NULL;
 32:     GEN           gen;
 33:     LOAD          load;
 34:     PetscBool     ghostvtex;
 35:     PetscInt      numComps;
 36:     void*         component;

 38:     DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
 39:     DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
 40:     DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset);
 41:     for (j = 0; j < numComps; j++) {
 42:       DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component,NULL);
 43:       if (key == 1) {
 44:         PetscInt       nconnedges;
 45:         const PetscInt *connedges;

 47:         bus = (VERTEX_Power)(component);
 48:         /* Handle reference bus constrained dofs */
 49:         if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
 50:           farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
 51:           farr[offset+1] = xarr[offset+1] - bus->vm;
 52:           break;
 53:         }

 55:         if (!ghostvtex) {
 56:           Vm = xarr[offset+1];

 58:           /* Shunt injections */
 59:           farr[offset] += Vm*Vm*bus->gl/Sbase;
 60:           if (bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
 61:         }

 63:         DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
 64:         for (i=0; i < nconnedges; i++) {
 65:           EDGE_Power     branch;
 66:           PetscInt       keye;
 67:           PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
 68:           const PetscInt *cone;
 69:           PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;

 71:           e = connedges[i];
 72:           DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);
 73:           if (!branch->status) continue;
 74:           Gff = branch->yff[0];
 75:           Bff = branch->yff[1];
 76:           Gft = branch->yft[0];
 77:           Bft = branch->yft[1];
 78:           Gtf = branch->ytf[0];
 79:           Btf = branch->ytf[1];
 80:           Gtt = branch->ytt[0];
 81:           Btt = branch->ytt[1];

 83:           DMNetworkGetConnectedVertices(networkdm,e,&cone);
 84:           vfrom = cone[0];
 85:           vto   = cone[1];

 87:           DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
 88:           DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);

 90:           thetaf = xarr[offsetfrom];
 91:           Vmf     = xarr[offsetfrom+1];
 92:           thetat  = xarr[offsetto];
 93:           Vmt     = xarr[offsetto+1];
 94:           thetaft = thetaf - thetat;
 95:           thetatf = thetat - thetaf;

 97:           if (vfrom == vtx[v]) {
 98:             farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));
 99:             farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
100:           } else {
101:             farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf));
102:             farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
103:           }
104:         }
105:       } else if (key == 2) {
106:         if (!ghostvtex) {
107:           gen = (GEN)(component);
108:           if (!gen->status) continue;
109:           farr[offset] += -gen->pg/Sbase;
110:           farr[offset+1] += -gen->qg/Sbase;
111:         }
112:       } else if (key == 3) {
113:         if (!ghostvtex) {
114:           load = (LOAD)(component);
115:           farr[offset] += load->pl/Sbase;
116:           farr[offset+1] += load->ql/Sbase;
117:         }
118:       }
119:     }
120:     if (bus && bus->ide == PV_BUS) {
121:       farr[offset+1] = xarr[offset+1] - bus->vm;
122:     }
123:   }
124:   VecRestoreArrayRead(localX,&xarr);
125:   VecRestoreArray(localF,&farr);
126:   return 0;
127: }

129: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
130: {
131:   DM             networkdm;
132:   Vec            localX,localF;
133:   PetscInt       nv,ne;
134:   const PetscInt *vtx,*edges;

136:   SNESGetDM(snes,&networkdm);
137:   DMGetLocalVector(networkdm,&localX);
138:   DMGetLocalVector(networkdm,&localF);
139:   VecSet(F,0.0);

141:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
142:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

144:   DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
145:   DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);

147:   /* Form Function for first subnetwork */
148:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
149:   FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);

151:   /* Form Function for second subnetwork */
152:   DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
153:   FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);

155:   DMRestoreLocalVector(networkdm,&localX);

157:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
158:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
159:   DMRestoreLocalVector(networkdm,&localF);
160:   return 0;
161: }

163: PetscErrorCode FormJacobian_Subnet(DM networkdm,Vec localX, Mat J, Mat Jpre, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
164: {
165:   UserCtx_Power     *User=(UserCtx_Power*)appctx;
166:   PetscInt          e,v,vfrom,vto;
167:   const PetscScalar *xarr;
168:   PetscInt          offsetfrom,offsetto,goffsetfrom,goffsetto;
169:   PetscInt          row[2],col[8];
170:   PetscScalar       values[8];

172:   VecGetArrayRead(localX,&xarr);

174:   for (v=0; v<nv; v++) {
175:     PetscInt    i,j,key;
176:     PetscInt    offset,goffset;
177:     PetscScalar Vm;
178:     PetscScalar Sbase=User->Sbase;
179:     VERTEX_Power bus;
180:     PetscBool   ghostvtex;
181:     PetscInt    numComps;
182:     void*       component;

184:     DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
185:     DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
186:     for (j = 0; j < numComps; j++) {
187:       DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset);
188:       DMNetworkGetGlobalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&goffset);
189:       DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component,NULL);
190:       if (key == 1) {
191:         PetscInt       nconnedges;
192:         const PetscInt *connedges;

194:         bus = (VERTEX_Power)(component);
195:         if (!ghostvtex) {
196:           /* Handle reference bus constrained dofs */
197:           if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
198:             row[0] = goffset; row[1] = goffset+1;
199:             col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
200:             values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
201:             MatSetValues(J,2,row,2,col,values,ADD_VALUES);
202:             break;
203:           }

205:           Vm = xarr[offset+1];

207:           /* Shunt injections */
208:           row[0] = goffset; row[1] = goffset+1;
209:           col[0] = goffset; col[1] = goffset+1;
210:           values[0] = values[1] = values[2] = values[3] = 0.0;
211:           if (bus->ide != PV_BUS) {
212:             values[1] = 2.0*Vm*bus->gl/Sbase;
213:             values[3] = -2.0*Vm*bus->bl/Sbase;
214:           }
215:           MatSetValues(J,2,row,2,col,values,ADD_VALUES);
216:         }

218:         DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
219:         for (i=0; i < nconnedges; i++) {
220:           EDGE_Power       branch;
221:           VERTEX_Power     busf,bust;
222:           PetscInt       keyf,keyt;
223:           PetscScalar    Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
224:           const PetscInt *cone;
225:           PetscScalar    Vmf,Vmt,thetaf,thetat,thetaft,thetatf;

227:           e = connedges[i];
228:           DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
229:           if (!branch->status) continue;

231:           Gff = branch->yff[0];
232:           Bff = branch->yff[1];
233:           Gft = branch->yft[0];
234:           Bft = branch->yft[1];
235:           Gtf = branch->ytf[0];
236:           Btf = branch->ytf[1];
237:           Gtt = branch->ytt[0];
238:           Btt = branch->ytt[1];

240:           DMNetworkGetConnectedVertices(networkdm,e,&cone);
241:           vfrom = cone[0];
242:           vto   = cone[1];

244:           DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
245:           DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
246:           DMNetworkGetGlobalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&goffsetfrom);
247:           DMNetworkGetGlobalVecOffset(networkdm,vto,ALL_COMPONENTS,&goffsetto);

249:           if (goffsetto < 0) goffsetto = -goffsetto - 1;

251:           thetaf = xarr[offsetfrom];
252:           Vmf     = xarr[offsetfrom+1];
253:           thetat = xarr[offsetto];
254:           Vmt     = xarr[offsetto+1];
255:           thetaft = thetaf - thetat;
256:           thetatf = thetat - thetaf;

258:           DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf,NULL);
259:           DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust,NULL);

261:           if (vfrom == vtx[v]) {
262:             if (busf->ide != REF_BUS) {
263:               /*    farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));  */
264:               row[0]  = goffsetfrom;
265:               col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
266:               values[0] =  Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
267:               values[1] =  2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
268:               values[2] =  Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
269:               values[3] =  Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */

271:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
272:             }
273:             if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
274:               row[0] = goffsetfrom+1;
275:               col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
276:               /*    farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
277:               values[0] =  Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
278:               values[1] =  -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
279:               values[2] =  Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
280:               values[3] =  Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));

282:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
283:             }
284:           } else {
285:             if (bust->ide != REF_BUS) {
286:               row[0] = goffsetto;
287:               col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
288:               /*    farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
289:               values[0] =  Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
290:               values[1] =  2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
291:               values[2] =  Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
292:               values[3] =  Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */

294:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
295:             }
296:             if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
297:               row[0] = goffsetto+1;
298:               col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
299:               /*    farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
300:               values[0] =  Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
301:               values[1] =  -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
302:               values[2] =  Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
303:               values[3] =  Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));

305:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
306:             }
307:           }
308:         }
309:         if (!ghostvtex && bus->ide == PV_BUS) {
310:           row[0] = goffset+1; col[0] = goffset+1;
311:           values[0]  = 1.0;
312:           MatSetValues(J,1,row,1,col,values,ADD_VALUES);
313:         }
314:       }
315:     }
316:   }
317:   VecRestoreArrayRead(localX,&xarr);
318:   return 0;
319: }

321: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
322: {
323:   DM             networkdm;
324:   Vec            localX;
325:   PetscInt       ne,nv;
326:   const PetscInt *vtx,*edges;

328:   MatZeroEntries(J);

330:   SNESGetDM(snes,&networkdm);
331:   DMGetLocalVector(networkdm,&localX);

333:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
334:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

336:   /* Form Jacobian for first subnetwork */
337:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
338:   FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);

340:   /* Form Jacobian for second subnetwork */
341:   DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
342:   FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);

344:   DMRestoreLocalVector(networkdm,&localX);

346:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
347:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
348:   return 0;
349: }

351: PetscErrorCode SetInitialValues_Subnet(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx)
352: {
353:   VERTEX_Power   bus;
354:   PetscInt       i;
355:   GEN            gen;
356:   PetscBool      ghostvtex;
357:   PetscScalar    *xarr;
358:   PetscInt       key,numComps,j,offset;
359:   void*          component;

361:   VecGetArray(localX,&xarr);
362:   for (i = 0; i < nv; i++) {
363:     DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
364:     if (ghostvtex) continue;

366:     DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
367:     DMNetworkGetNumComponents(networkdm,vtx[i],&numComps);
368:     for (j=0; j < numComps; j++) {
369:       DMNetworkGetComponent(networkdm,vtx[i],j,&key,&component,NULL);
370:       if (key == 1) {
371:         bus = (VERTEX_Power)(component);
372:         xarr[offset] = bus->va*PETSC_PI/180.0;
373:         xarr[offset+1] = bus->vm;
374:       } else if (key == 2) {
375:         gen = (GEN)(component);
376:         if (!gen->status) continue;
377:         xarr[offset+1] = gen->vs;
378:         break;
379:       }
380:     }
381:   }
382:   VecRestoreArray(localX,&xarr);
383:   return 0;
384: }

386: PetscErrorCode SetInitialValues(DM networkdm, Vec X,void* appctx)
387: {
388:   PetscInt       nv,ne;
389:   const PetscInt *vtx,*edges;
390:   Vec            localX;

392:   DMGetLocalVector(networkdm,&localX);

394:   VecSet(X,0.0);
395:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
396:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

398:   /* Set initial guess for first subnetwork */
399:   DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
400:   SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);

402:   /* Set initial guess for second subnetwork */
403:   DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
404:   SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);

406:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
407:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
408:   DMRestoreLocalVector(networkdm,&localX);
409:   return 0;
410: }

412: int main(int argc,char ** argv)
413: {
414:   char             pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
415:   PFDATA           *pfdata1,*pfdata2;
416:   PetscInt         numEdges1=0,numEdges2=0;
417:   PetscInt         *edgelist1 = NULL,*edgelist2 = NULL,componentkey[4];
418:   DM               networkdm;
419:   UserCtx_Power    User;
420: #if defined(PETSC_USE_LOG)
421:   PetscLogStage    stage1,stage2;
422: #endif
423:   PetscMPIInt      rank;
424:   PetscInt         nsubnet = 2,nv,ne,i,j,genj,loadj;
425:   const PetscInt   *vtx,*edges;
426:   Vec              X,F;
427:   Mat              J;
428:   SNES             snes;

430:   PetscInitialize(&argc,&argv,"poweroptions",help);
431:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
432:   {
433:     /* introduce the const crank so the clang static analyzer realizes that if it enters any of the if (crank) then it must have entered the first */
434:     /* this is an experiment to see how the analyzer reacts */
435:     const PetscMPIInt crank = rank;

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

440:     /* Register the components in the network */
441:     DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0]);
442:     DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1]);
443:     DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
444:     DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);

446:     PetscLogStageRegister("Read Data",&stage1);
447:     PetscLogStagePush(stage1);
448:     /* READ THE DATA */
449:     if (!crank) {
450:       /* Only rank 0 reads the data */
451:       PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);
452:       /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */

454:       /*    READ DATA FOR THE FIRST SUBNETWORK */
455:       PetscNew(&pfdata1);
456:       PFReadMatPowerData(pfdata1,pfdata_file);
457:       User.Sbase = pfdata1->sbase;

459:       numEdges1 = pfdata1->nbranch;
460:       PetscMalloc1(2*numEdges1,&edgelist1);
461:       GetListofEdges_Power(pfdata1,edgelist1);

463:       /*    READ DATA FOR THE SECOND SUBNETWORK */
464:       PetscNew(&pfdata2);
465:       PFReadMatPowerData(pfdata2,pfdata_file);
466:       User.Sbase = pfdata2->sbase;

468:       numEdges2 = pfdata2->nbranch;
469:       PetscMalloc1(2*numEdges2,&edgelist2);
470:       GetListofEdges_Power(pfdata2,edgelist2);
471:     }

473:     PetscLogStagePop();
474:     MPI_Barrier(PETSC_COMM_WORLD);
475:     PetscLogStageRegister("Create network",&stage2);
476:     PetscLogStagePush(stage2);

478:     /* Set number of nodes/edges and edge connectivity */
479:     DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,nsubnet);
480:     DMNetworkAddSubnetwork(networkdm,"",numEdges1,edgelist1,NULL);
481:     DMNetworkAddSubnetwork(networkdm,"",numEdges2,edgelist2,NULL);

483:     /* Set up the network layout */
484:     DMNetworkLayoutSetUp(networkdm);

486:     /* Add network components only process 0 has any data to add*/
487:     if (!crank) {
488:       genj=0; loadj=0;

490:       /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
491:       DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);

493:       for (i = 0; i < ne; i++) {
494:         DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata1->branch[i],0);
495:       }

497:       for (i = 0; i < nv; i++) {
498:         DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata1->bus[i],2);
499:         if (pfdata1->bus[i].ngen) {
500:           for (j = 0; j < pfdata1->bus[i].ngen; j++) {
501:             DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata1->gen[genj++],0);
502:           }
503:         }
504:         if (pfdata1->bus[i].nload) {
505:           for (j=0; j < pfdata1->bus[i].nload; j++) {
506:             DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata1->load[loadj++],0);
507:           }
508:         }
509:       }

511:       genj=0; loadj=0;

513:       /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
514:       DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);

516:       for (i = 0; i < ne; i++) {
517:         DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata2->branch[i],0);
518:       }

520:       for (i = 0; i < nv; i++) {
521:         DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata2->bus[i],2);
522:         if (pfdata2->bus[i].ngen) {
523:           for (j = 0; j < pfdata2->bus[i].ngen; j++) {
524:             DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata2->gen[genj++],0);
525:           }
526:         }
527:         if (pfdata2->bus[i].nload) {
528:           for (j=0; j < pfdata2->bus[i].nload; j++) {
529:             DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata2->load[loadj++],0);
530:           }
531:         }
532:       }
533:     }

535:     /* Set up DM for use */
536:     DMSetUp(networkdm);

538:     if (!crank) {
539:       PetscFree(edgelist1);
540:       PetscFree(edgelist2);
541:     }

543:     if (!crank) {
544:       PetscFree(pfdata1->bus);
545:       PetscFree(pfdata1->gen);
546:       PetscFree(pfdata1->branch);
547:       PetscFree(pfdata1->load);
548:       PetscFree(pfdata1);

550:       PetscFree(pfdata2->bus);
551:       PetscFree(pfdata2->gen);
552:       PetscFree(pfdata2->branch);
553:       PetscFree(pfdata2->load);
554:       PetscFree(pfdata2);
555:     }

557:     /* Distribute networkdm to multiple processes */
558:     DMNetworkDistribute(&networkdm,0);

560:     PetscLogStagePop();

562:     /* Broadcast Sbase to all processors */
563:     MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);

565:     DMCreateGlobalVector(networkdm,&X);
566:     VecDuplicate(X,&F);

568:     DMCreateMatrix(networkdm,&J);
569:     MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

571:     SetInitialValues(networkdm,X,&User);

573:     /* HOOK UP SOLVER */
574:     SNESCreate(PETSC_COMM_WORLD,&snes);
575:     SNESSetDM(snes,networkdm);
576:     SNESSetFunction(snes,F,FormFunction,&User);
577:     SNESSetJacobian(snes,J,J,FormJacobian,&User);
578:     SNESSetFromOptions(snes);

580:     SNESSolve(snes,NULL,X);
581:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

583:     VecDestroy(&X);
584:     VecDestroy(&F);
585:     MatDestroy(&J);

587:     SNESDestroy(&snes);
588:     DMDestroy(&networkdm);
589:   }
590:   PetscFinalize();
591:   return 0;
592: }

594: /*TEST

596:    build:
597:      depends: PFReadData.c pffunctions.c
598:      requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

600:    test:
601:      args: -snes_rtol 1.e-3
602:      localrunfiles: poweroptions case9.m
603:      output_file: output/power_1.out

605:    test:
606:      suffix: 2
607:      args: -snes_rtol 1.e-3 -petscpartitioner_type simple
608:      nsize: 4
609:      localrunfiles: poweroptions case9.m
610:      output_file: output/power_1.out

612: TEST*/