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: #include "power.h"
  9: #include <petscdmnetwork.h>

 11: PetscErrorCode FormFunction_Subnet(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
 12: {
 13:   UserCtx_Power     *User = (UserCtx_Power*)appctx;
 14:   PetscInt          e,v,vfrom,vto;
 15:   const PetscScalar *xarr;
 16:   PetscScalar       *farr;
 17:   PetscInt          offsetfrom,offsetto,offset;

 19:   VecGetArrayRead(localX,&xarr);
 20:   VecGetArray(localF,&farr);

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

 33:     DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
 34:     DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
 35:     DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset);
 36:     for (j = 0; j < numComps; j++) {
 37:       DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component,NULL);
 38:       if (key == 1) {
 39:         PetscInt       nconnedges;
 40:         const PetscInt *connedges;

 42:         bus = (VERTEX_Power)(component);
 43:         /* Handle reference bus constrained dofs */
 44:         if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
 45:           farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
 46:           farr[offset+1] = xarr[offset+1] - bus->vm;
 47:           break;
 48:         }

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

 53:           /* Shunt injections */
 54:           farr[offset] += Vm*Vm*bus->gl/Sbase;
 55:           if (bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
 56:         }

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

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

 78:           DMNetworkGetConnectedVertices(networkdm,e,&cone);
 79:           vfrom = cone[0];
 80:           vto   = cone[1];

 82:           DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
 83:           DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);

 85:           thetaf = xarr[offsetfrom];
 86:           Vmf     = xarr[offsetfrom+1];
 87:           thetat  = xarr[offsetto];
 88:           Vmt     = xarr[offsetto+1];
 89:           thetaft = thetaf - thetat;
 90:           thetatf = thetat - thetaf;

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

124: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
125: {
126:   DM             networkdm;
127:   Vec            localX,localF;
128:   PetscInt       nv,ne;
129:   const PetscInt *vtx,*edges;

131:   SNESGetDM(snes,&networkdm);
132:   DMGetLocalVector(networkdm,&localX);
133:   DMGetLocalVector(networkdm,&localF);
134:   VecSet(F,0.0);

136:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
137:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

139:   DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
140:   DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);

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

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

150:   DMRestoreLocalVector(networkdm,&localX);

152:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
153:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
154:   DMRestoreLocalVector(networkdm,&localF);
155:   return 0;
156: }

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

167:   VecGetArrayRead(localX,&xarr);

169:   for (v=0; v<nv; v++) {
170:     PetscInt    i,j,key;
171:     PetscInt    offset,goffset;
172:     PetscScalar Vm;
173:     PetscScalar Sbase=User->Sbase;
174:     VERTEX_Power bus;
175:     PetscBool   ghostvtex;
176:     PetscInt    numComps;
177:     void*       component;

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

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

200:           Vm = xarr[offset+1];

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

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

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

226:           Gff = branch->yff[0];
227:           Bff = branch->yff[1];
228:           Gft = branch->yft[0];
229:           Bft = branch->yft[1];
230:           Gtf = branch->ytf[0];
231:           Btf = branch->ytf[1];
232:           Gtt = branch->ytt[0];
233:           Btt = branch->ytt[1];

235:           DMNetworkGetConnectedVertices(networkdm,e,&cone);
236:           vfrom = cone[0];
237:           vto   = cone[1];

239:           DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
240:           DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
241:           DMNetworkGetGlobalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&goffsetfrom);
242:           DMNetworkGetGlobalVecOffset(networkdm,vto,ALL_COMPONENTS,&goffsetto);

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

246:           thetaf = xarr[offsetfrom];
247:           Vmf     = xarr[offsetfrom+1];
248:           thetat = xarr[offsetto];
249:           Vmt     = xarr[offsetto+1];
250:           thetaft = thetaf - thetat;
251:           thetatf = thetat - thetaf;

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

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

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

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

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

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

316: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
317: {
318:   DM             networkdm;
319:   Vec            localX;
320:   PetscInt       ne,nv;
321:   const PetscInt *vtx,*edges;

323:   MatZeroEntries(J);

325:   SNESGetDM(snes,&networkdm);
326:   DMGetLocalVector(networkdm,&localX);

328:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
329:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

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

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

339:   DMRestoreLocalVector(networkdm,&localX);

341:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
342:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
343:   return 0;
344: }

346: PetscErrorCode SetInitialValues_Subnet(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx)
347: {
348:   VERTEX_Power   bus;
349:   PetscInt       i;
350:   GEN            gen;
351:   PetscBool      ghostvtex;
352:   PetscScalar    *xarr;
353:   PetscInt       key,numComps,j,offset;
354:   void*          component;

356:   VecGetArray(localX,&xarr);
357:   for (i = 0; i < nv; i++) {
358:     DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
359:     if (ghostvtex) continue;

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

381: PetscErrorCode SetInitialValues(DM networkdm, Vec X,void* appctx)
382: {
383:   PetscInt       nv,ne;
384:   const PetscInt *vtx,*edges;
385:   Vec            localX;

387:   DMGetLocalVector(networkdm,&localX);

389:   VecSet(X,0.0);
390:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
391:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

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

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

401:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
402:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
403:   DMRestoreLocalVector(networkdm,&localX);
404:   return 0;
405: }

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

425:   PetscInitialize(&argc,&argv,"poweroptions",help);
426:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
427:   {
428:     /* 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 */
429:     /* this is an experiment to see how the analyzer reacts */
430:     const PetscMPIInt crank = rank;

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

435:     /* Register the components in the network */
436:     DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0]);
437:     DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1]);
438:     DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
439:     DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);

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

449:       /*    READ DATA FOR THE FIRST SUBNETWORK */
450:       PetscNew(&pfdata1);
451:       PFReadMatPowerData(pfdata1,pfdata_file);
452:       User.Sbase = pfdata1->sbase;

454:       numEdges1 = pfdata1->nbranch;
455:       PetscMalloc1(2*numEdges1,&edgelist1);
456:       GetListofEdges_Power(pfdata1,edgelist1);

458:       /*    READ DATA FOR THE SECOND SUBNETWORK */
459:       PetscNew(&pfdata2);
460:       PFReadMatPowerData(pfdata2,pfdata_file);
461:       User.Sbase = pfdata2->sbase;

463:       numEdges2 = pfdata2->nbranch;
464:       PetscMalloc1(2*numEdges2,&edgelist2);
465:       GetListofEdges_Power(pfdata2,edgelist2);
466:     }

468:     PetscLogStagePop();
469:     MPI_Barrier(PETSC_COMM_WORLD);
470:     PetscLogStageRegister("Create network",&stage2);
471:     PetscLogStagePush(stage2);

473:     /* Set number of nodes/edges and edge connectivity */
474:     DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,nsubnet);
475:     DMNetworkAddSubnetwork(networkdm,"",numEdges1,edgelist1,NULL);
476:     DMNetworkAddSubnetwork(networkdm,"",numEdges2,edgelist2,NULL);

478:     /* Set up the network layout */
479:     DMNetworkLayoutSetUp(networkdm);

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

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

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

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

506:       genj=0; loadj=0;

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

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

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

530:     /* Set up DM for use */
531:     DMSetUp(networkdm);

533:     if (!crank) {
534:       PetscFree(edgelist1);
535:       PetscFree(edgelist2);
536:     }

538:     if (!crank) {
539:       PetscFree(pfdata1->bus);
540:       PetscFree(pfdata1->gen);
541:       PetscFree(pfdata1->branch);
542:       PetscFree(pfdata1->load);
543:       PetscFree(pfdata1);

545:       PetscFree(pfdata2->bus);
546:       PetscFree(pfdata2->gen);
547:       PetscFree(pfdata2->branch);
548:       PetscFree(pfdata2->load);
549:       PetscFree(pfdata2);
550:     }

552:     /* Distribute networkdm to multiple processes */
553:     DMNetworkDistribute(&networkdm,0);

555:     PetscLogStagePop();

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

560:     DMCreateGlobalVector(networkdm,&X);
561:     VecDuplicate(X,&F);

563:     DMCreateMatrix(networkdm,&J);
564:     MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);

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

568:     /* HOOK UP SOLVER */
569:     SNESCreate(PETSC_COMM_WORLD,&snes);
570:     SNESSetDM(snes,networkdm);
571:     SNESSetFunction(snes,F,FormFunction,&User);
572:     SNESSetJacobian(snes,J,J,FormJacobian,&User);
573:     SNESSetFromOptions(snes);

575:     SNESSolve(snes,NULL,X);
576:     /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

578:     VecDestroy(&X);
579:     VecDestroy(&F);
580:     MatDestroy(&J);

582:     SNESDestroy(&snes);
583:     DMDestroy(&networkdm);
584:   }
585:   PetscFinalize();
586:   return 0;
587: }

589: /*TEST

591:    build:
592:      depends: PFReadData.c pffunctions.c
593:      requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

595:    test:
596:      args: -snes_rtol 1.e-3
597:      localrunfiles: poweroptions case9.m
598:      output_file: output/power_1.out

600:    test:
601:      suffix: 2
602:      args: -snes_rtol 1.e-3 -petscpartitioner_type simple
603:      nsize: 4
604:      localrunfiles: poweroptions case9.m
605:      output_file: output/power_1.out

607: TEST*/