Actual source code: pf.c

petsc-3.6.1 2015-08-06
Report Typos and Errors
  1: static char help[] = "This example demonstrates the use of DMNetwork interface for solving a nonlinear electric power grid problem.\n\
  2:                       The available solver options are in the pfoptions 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:                       Run this program: mpiexec -n <n> ./PF\n\
  5:                       mpiexec -n <n> ./PF -pfdata <filename>\n";

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

 12:  #include pf.h
 13: #include <petscdmnetwork.h>

 15: PetscMPIInt rank;

 19: PetscErrorCode GetListofEdges(PetscInt nbranches, EDGEDATA branch,int edges[])
 20: {
 21:   PetscInt       i, fbus,tbus;

 24: 
 25:   for (i=0; i < nbranches; i++) {
 26:     fbus = branch[i].internal_i;
 27:     tbus = branch[i].internal_j;
 28:     edges[2*i]   = fbus;
 29:     edges[2*i+1] = tbus;
 30:   }
 31:   return(0);
 32: }

 34: typedef struct{
 35:   PetscScalar  Sbase;
 36: }UserCtx;

 40: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
 41: {
 43:   DM             networkdm;
 44:   UserCtx       *User=(UserCtx*)appctx;
 45:   Vec           localX,localF;
 46:   PetscInt      e;
 47:   PetscInt      v,vStart,vEnd,vfrom,vto;
 48:   const PetscScalar *xarr;
 49:   PetscScalar   *farr;
 50:   PetscInt      offsetfrom,offsetto,offset;
 51:   DMNetworkComponentGenericDataType *arr;

 54:   SNESGetDM(snes,&networkdm);
 55:   DMGetLocalVector(networkdm,&localX);
 56:   DMGetLocalVector(networkdm,&localF);
 57:   VecSet(F,0.0);

 59:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
 60:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

 62:   DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
 63:   DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);

 65:   VecGetArrayRead(localX,&xarr);
 66:   VecGetArray(localF,&farr);

 68:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
 69:   DMNetworkGetComponentDataArray(networkdm,&arr);

 71:   for (v=vStart; v < vEnd; v++) {
 72:     PetscInt    i,j,offsetd,key;
 73:     PetscScalar Vm;
 74:     PetscScalar Sbase=User->Sbase;
 75:     VERTEXDATA  bus;
 76:     GEN         gen;
 77:     LOAD        load;
 78:     PetscBool   ghostvtex;
 79:     PetscInt    numComps;

 81:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
 82:     DMNetworkGetNumComponents(networkdm,v,&numComps);
 83:     DMNetworkGetVariableOffset(networkdm,v,&offset);
 84:     for (j = 0; j < numComps; j++) {
 85:       DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);
 86:       if (key == 1) {
 87:         bus = (VERTEXDATA)(arr+offsetd);
 88:         /* Handle reference bus constrained dofs */
 89:         if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
 90:           farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
 91:           farr[offset+1] = xarr[offset+1] - bus->vm;
 92:           break;
 93:         }

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

 98:           /* Shunt injections */
 99:           farr[offset] += Vm*Vm*bus->gl/Sbase;
100:           if(bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
101:         }
102:         PetscInt nconnedges;
103:         const PetscInt *connedges;

105:         DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
106:         for (i=0; i < nconnedges; i++) {
107:           EDGEDATA branch;
108:           PetscInt keye;
109:           e = connedges[i];
110:           DMNetworkGetComponentTypeOffset(networkdm,e,0,&keye,&offsetd);
111:           branch = (EDGEDATA)(arr+offsetd);
112:           if (!branch->status) continue;
113:           PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
114:           Gff = branch->yff[0];
115:           Bff = branch->yff[1];
116:           Gft = branch->yft[0];
117:           Bft = branch->yft[1];
118:           Gtf = branch->ytf[0];
119:           Btf = branch->ytf[1];
120:           Gtt = branch->ytt[0];
121:           Btt = branch->ytt[1];

123:           const PetscInt *cone;
124:           DMNetworkGetConnectedNodes(networkdm,e,&cone);
125:           vfrom = cone[0];
126:           vto   = cone[1];

128:           DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
129:           DMNetworkGetVariableOffset(networkdm,vto,&offsetto);

131:           PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;

133:           thetaf = xarr[offsetfrom];
134:           Vmf     = xarr[offsetfrom+1];
135:           thetat = xarr[offsetto];
136:           Vmt     = xarr[offsetto+1];
137:           thetaft = thetaf - thetat;
138:           thetatf = thetat - thetaf;

140:           if (vfrom == v) {
141:             farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*cos(thetaft) + Bft*sin(thetaft));
142:             farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*cos(thetaft) + Gft*sin(thetaft));
143:           } else {
144:             farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*cos(thetatf) + Btf*sin(thetatf));
145:             farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*cos(thetatf) + Gtf*sin(thetatf));
146:           }
147:         }
148:       } else if (key == 2) {
149:         if (!ghostvtex) {
150:           gen = (GEN)(arr+offsetd);
151:           if (!gen->status) continue;
152:           farr[offset] += -gen->pg/Sbase;
153:           farr[offset+1] += -gen->qg/Sbase;
154:         }
155:       } else if (key == 3) {
156:         if (!ghostvtex) {
157:           load = (LOAD)(arr+offsetd);
158:           farr[offset] += load->pl/Sbase;
159:           farr[offset+1] += load->ql/Sbase;
160:         }
161:       }
162:     }
163:     if (bus->ide == PV_BUS) {
164:       farr[offset+1] = xarr[offset+1] - bus->vm;
165:     }
166:   }
167:   VecRestoreArrayRead(localX,&xarr);
168:   VecRestoreArray(localF,&farr);
169:   DMRestoreLocalVector(networkdm,&localX);

171:   DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
172:   DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
173:   DMRestoreLocalVector(networkdm,&localF);

175:   return(0);
176: }

180: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
181: {
183:   DM            networkdm;
184:   UserCtx       *User=(UserCtx*)appctx;
185:   Vec           localX;
186:   PetscInt      e;
187:   PetscInt      v,vStart,vEnd,vfrom,vto;
188:   const PetscScalar *xarr;
189:   PetscInt      offsetfrom,offsetto,goffsetfrom,goffsetto;
190:   DMNetworkComponentGenericDataType *arr;
191:   PetscInt      row[2],col[8];
192:   PetscScalar   values[8];

195:   MatZeroEntries(J);

197:   SNESGetDM(snes,&networkdm);
198:   DMGetLocalVector(networkdm,&localX);

200:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
201:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

203:   VecGetArrayRead(localX,&xarr);

205:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
206:   DMNetworkGetComponentDataArray(networkdm,&arr);

208:   for (v=vStart; v < vEnd; v++) {
209:     PetscInt    i,j,offsetd,key;
210:     PetscInt    offset,goffset;
211:     PetscScalar Vm;
212:     PetscScalar Sbase=User->Sbase;
213:     VERTEXDATA  bus;
214:     PetscBool   ghostvtex;
215:     PetscInt    numComps;

217:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
218:     DMNetworkGetNumComponents(networkdm,v,&numComps);
219:     for (j = 0; j < numComps; j++) {
220:       DMNetworkGetVariableOffset(networkdm,v,&offset);
221:       DMNetworkGetVariableGlobalOffset(networkdm,v,&goffset);
222:       DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);
223:       if (key == 1) {
224:         bus = (VERTEXDATA)(arr+offsetd);
225:         if (!ghostvtex) {
226:           /* Handle reference bus constrained dofs */
227:           if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
228:             row[0] = goffset; row[1] = goffset+1;
229:             col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
230:             values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
231:             MatSetValues(J,2,row,2,col,values,ADD_VALUES);
232:             break;
233:           }
234: 
235:           Vm = xarr[offset+1];
236: 
237:           /* Shunt injections */
238:           row[0] = goffset; row[1] = goffset+1;
239:           col[0] = goffset; col[1] = goffset+1;
240:           values[0] = values[1] = values[2] = values[3] = 0.0;
241:           if(bus->ide != PV_BUS) {
242:             values[1] = 2*Vm*bus->gl/Sbase;
243:             values[3] = -2*Vm*bus->bl/Sbase;
244:           }
245:           MatSetValues(J,2,row,2,col,values,ADD_VALUES);
246:         }

248:         PetscInt nconnedges;
249:         const PetscInt *connedges;

251:         DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
252:         for (i=0; i < nconnedges; i++) {
253:           EDGEDATA branch;
254:           VERTEXDATA busf,bust;
255:           PetscInt   offsetfd,offsettd,keyf,keyt;
256:           e = connedges[i];
257:           DMNetworkGetComponentTypeOffset(networkdm,e,0,&key,&offsetd);
258:           branch = (EDGEDATA)(arr+offsetd);
259:           if (!branch->status) continue;
260:           PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
261:           Gff = branch->yff[0];
262:           Bff = branch->yff[1];
263:           Gft = branch->yft[0];
264:           Bft = branch->yft[1];
265:           Gtf = branch->ytf[0];
266:           Btf = branch->ytf[1];
267:           Gtt = branch->ytt[0];
268:           Btt = branch->ytt[1];

270:           const PetscInt *cone;
271:           DMNetworkGetConnectedNodes(networkdm,e,&cone);
272:           vfrom = cone[0];
273:           vto   = cone[1];

275:           DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
276:           DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
277:           DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);
278:           DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);

280:           if (goffsetfrom < 0) goffsetfrom = -goffsetfrom - 1; /* Convert to actual global offset for ghost nodes, global offset is -(gstart+1) */
281:           if (goffsetto < 0) goffsetto = -goffsetto - 1;
282:           PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;

284:           thetaf = xarr[offsetfrom];
285:           Vmf     = xarr[offsetfrom+1];
286:           thetat = xarr[offsetto];
287:           Vmt     = xarr[offsetto+1];
288:           thetaft = thetaf - thetat;
289:           thetatf = thetat - thetaf;

291:           DMNetworkGetComponentTypeOffset(networkdm,vfrom,0,&keyf,&offsetfd);
292:           DMNetworkGetComponentTypeOffset(networkdm,vto,0,&keyt,&offsettd);
293:           busf = (VERTEXDATA)(arr+offsetfd);
294:           bust = (VERTEXDATA)(arr+offsettd);

296:           if (vfrom == v) {
297:             if (busf->ide != REF_BUS) {
298:               /*    farr[offsetfrom]   += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*cos(thetaft) + Bft*sin(thetaft));  */
299:               row[0]  = goffsetfrom;
300:               col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
301:               values[0] =  Vmf*Vmt*(Gft*-sin(thetaft) + Bft*cos(thetaft)); /* df_dthetaf */
302:               values[1] =  2*Gff*Vmf + Vmt*(Gft*cos(thetaft) + Bft*sin(thetaft)); /* df_dVmf */
303:               values[2] =  Vmf*Vmt*(Gft*sin(thetaft) + Bft*-cos(thetaft)); /* df_dthetat */
304:               values[3] =  Vmf*(Gft*cos(thetaft) + Bft*sin(thetaft)); /* df_dVmt */
305: 
306:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
307:             }
308:             if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
309:               row[0] = goffsetfrom+1;
310:               col[0]  = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
311:               /*    farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*cos(thetaft) + Gft*sin(thetaft)); */
312:               values[0] =  Vmf*Vmt*(Bft*sin(thetaft) + Gft*cos(thetaft));
313:               values[1] =  -2*Bff*Vmf + Vmt*(-Bft*cos(thetaft) + Gft*sin(thetaft));
314:               values[2] =  Vmf*Vmt*(-Bft*sin(thetaft) + Gft*-cos(thetaft));
315:               values[3] =  Vmf*(-Bft*cos(thetaft) + Gft*sin(thetaft));
316: 
317:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
318:             }
319:           } else {
320:             if(bust->ide != REF_BUS) {
321:               row[0] = goffsetto;
322:               col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
323:               /*    farr[offsetto]   += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*cos(thetatf) + Btf*sin(thetatf)); */
324:               values[0] =  Vmt*Vmf*(Gtf*-sin(thetatf) + Btf*cos(thetaft)); /* df_dthetat */
325:               values[1] =  2*Gtt*Vmt + Vmf*(Gtf*cos(thetatf) + Btf*sin(thetatf)); /* df_dVmt */
326:               values[2] =  Vmt*Vmf*(Gtf*sin(thetatf) + Btf*-cos(thetatf)); /* df_dthetaf */
327:               values[3] =  Vmt*(Gtf*cos(thetatf) + Btf*sin(thetatf)); /* df_dVmf */
328: 
329:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
330:             }
331:             if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
332:               row[0] = goffsetto+1;
333:               col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
334:               /*    farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*cos(thetatf) + Gtf*sin(thetatf)); */
335:               values[0] =  Vmt*Vmf*(Btf*sin(thetatf) + Gtf*cos(thetatf));
336:               values[1] =  -2*Btt*Vmt + Vmf*(-Btf*cos(thetatf) + Gtf*sin(thetatf));
337:               values[2] =  Vmt*Vmf*(-Btf*sin(thetatf) + Gtf*-cos(thetatf));
338:               values[3] =  Vmt*(-Btf*cos(thetatf) + Gtf*sin(thetatf));
339: 
340:               MatSetValues(J,1,row,4,col,values,ADD_VALUES);
341:             }
342:           }
343:         }
344:         if (!ghostvtex && bus->ide == PV_BUS) {
345:           row[0] = goffset+1; col[0] = goffset+1;
346:           values[0]  = 1.0;
347:           MatSetValues(J,1,row,1,col,values,ADD_VALUES);
348:         }
349:       }
350:     }
351:   }
352:   VecRestoreArrayRead(localX,&xarr);
353:   DMRestoreLocalVector(networkdm,&localX);

355:   MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
356:   MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);

358:   return(0);
359: }

363: PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx)
364: {
366:   VERTEXDATA     bus;
367:   GEN            gen;
368:   PetscInt       v, vStart, vEnd, offset;
369:   PetscBool      ghostvtex;
370:   Vec            localX;
371:   PetscScalar    *xarr;
372:   PetscInt       key;
373:   DMNetworkComponentGenericDataType *arr;
374: 
376:   DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);

378:   DMGetLocalVector(networkdm,&localX);

380:   VecSet(X,0.0);
381:   DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
382:   DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);

384:   VecGetArray(localX,&xarr);
385:   DMNetworkGetComponentDataArray(networkdm,&arr);
386:   for (v = vStart; v < vEnd; v++) {
387:     DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
388:     if (ghostvtex) continue;
389:     PetscInt numComps;
390:     PetscInt j;
391:     PetscInt offsetd;
392:     DMNetworkGetVariableOffset(networkdm,v,&offset);
393:     DMNetworkGetNumComponents(networkdm,v,&numComps);
394:     for (j=0; j < numComps; j++) {
395:       DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);
396:       if (key == 1) {
397:         bus = (VERTEXDATA)(arr+offsetd);
398:         xarr[offset] = bus->va*PETSC_PI/180.0;
399:         xarr[offset+1] = bus->vm;
400:       } else if(key == 2) {
401:         gen = (GEN)(arr+offsetd);
402:         if (!gen->status) continue;
403:         xarr[offset+1] = gen->vs;
404:         break;
405:       }
406:     }
407:   }
408:   VecRestoreArray(localX,&xarr);
409:   DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
410:   DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
411:   DMRestoreLocalVector(networkdm,&localX);
412:   return(0);
413: }


418: int main(int argc,char ** argv)
419: {
420:   PetscErrorCode       ierr;
421:   char                 pfdata_file[PETSC_MAX_PATH_LEN]="datafiles/case9.m";
422:   PFDATA               pfdata;
423:   PetscInt             numEdges=0,numVertices=0;
424:   int                  *edges = NULL;
425:   PetscInt             i;
426:   DM                   networkdm;
427:   PetscInt             componentkey[4];
428:   UserCtx              User;
429:   PetscLogStage        stage1,stage2;
430:   PetscInt             size;

432:   PetscInitialize(&argc,&argv,"pfoptions",help);

434:   MPI_Comm_rank(PETSC_COMM_WORLD,&rank);

436:   /* Create an empty network object */
437:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
438:   /* Register the components in the network */
439:   DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGEDATA),&componentkey[0]);
440:   DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEXDATA),&componentkey[1]);
441:   DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
442:   DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);

444:   PetscLogStageRegister("Read Data",&stage1);
445:   PetscLogStagePush(stage1);
446:   /* READ THE DATA */
447:   if (!rank) {
448:     /*    READ DATA */
449:     /* Only rank 0 reads the data */
450:     PetscOptionsGetString(PETSC_NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
451:     PFReadMatPowerData(&pfdata,pfdata_file);
452:     User.Sbase = pfdata.sbase;

454:     numEdges = pfdata.nbranch;
455:     numVertices = pfdata.nbus;

457:     PetscMalloc(2*numEdges*sizeof(int),&edges);
458:     GetListofEdges(pfdata.nbranch,pfdata.branch,edges);
459:   }
460:   PetscLogStagePop();
461:   MPI_Barrier(PETSC_COMM_WORLD);
462:   PetscLogStageRegister("Create network",&stage2);
463:   PetscLogStagePush(stage2);
464:   /* Set number of nodes/edges */
465:   DMNetworkSetSizes(networkdm,numVertices,numEdges,PETSC_DETERMINE,PETSC_DETERMINE);
466:   /* Add edge connectivity */
467:   DMNetworkSetEdgeList(networkdm,edges);
468:   /* Set up the network layout */
469:   DMNetworkLayoutSetUp(networkdm);

471:   if (!rank) {
472:     PetscFree(edges);
473:   }
474:   /* Add network components */
475:   PetscInt eStart, eEnd, vStart, vEnd,j;
476:   PetscInt genj=0,loadj=0;
477:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
478:   for (i = eStart; i < eEnd; i++) {
479:     DMNetworkAddComponent(networkdm,i,componentkey[0],&pfdata.branch[i-eStart]);
480:   }
481:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
482:   for (i = vStart; i < vEnd; i++) {
483:     DMNetworkAddComponent(networkdm,i,componentkey[1],&pfdata.bus[i-vStart]);
484:     if (pfdata.bus[i-vStart].ngen) {
485:       for (j = 0; j < pfdata.bus[i-vStart].ngen; j++) {
486:         DMNetworkAddComponent(networkdm,i,componentkey[2],&pfdata.gen[genj++]);
487:       }
488:     }
489:     if (pfdata.bus[i-vStart].nload) {
490:       for (j=0; j < pfdata.bus[i-vStart].nload; j++) {
491:         DMNetworkAddComponent(networkdm,i,componentkey[3],&pfdata.load[loadj++]);
492:       }
493:     }
494:     /* Add number of variables */
495:     DMNetworkAddNumVariables(networkdm,i,2);
496:   }
497:   /* Set up DM for use */
498:   DMSetUp(networkdm);

500:   if (!rank) {
501:     PetscFree(pfdata.bus);
502:     PetscFree(pfdata.gen);
503:     PetscFree(pfdata.branch);
504:     PetscFree(pfdata.load);
505:   }


508:   MPI_Comm_size(PETSC_COMM_WORLD,&size);
509:   if (size > 1) {
510:     DM distnetworkdm;
511:     /* Network partitioning and distribution of data */
512:     DMNetworkDistribute(networkdm,0,&distnetworkdm);
513:     DMDestroy(&networkdm);
514:     networkdm = distnetworkdm;
515:   }

517:   PetscLogStagePop();
518:   DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
519:   DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
520: 
521: #if 0
522:   PetscInt numComponents;
523:   EDGEDATA edge;
524:   PetscInt offset,key;
525:   DMNetworkComponentGenericDataType *arr;
526:   for (i = eStart; i < eEnd; i++) {
527:     DMNetworkGetComponentDataArray(networkdm,&arr);
528:     DMNetworkGetComponentTypeOffset(networkdm,i,0,&key,&offset);
529:     edge = (EDGEDATA)(arr+offset);
530:     DMNetworkGetNumComponents(networkdm,i,&numComponents);
531:     PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",rank,numComponents,edge->internal_i,edge->internal_j);
532:   }

534:   VERTEXDATA bus;
535:   GEN        gen;
536:   LOAD       load;
537:   PetscInt   kk;
538:   for (i = vStart; i < vEnd; i++) {
539:     DMNetworkGetComponentDataArray(networkdm,&arr);
540:     DMNetworkGetNumComponents(networkdm,i,&numComponents);
541:     for (kk=0; kk < numComponents; kk++) {
542:       DMNetworkGetComponentTypeOffset(networkdm,i,kk,&key,&offset);
543:       if (key == 1) {
544:         bus = (VERTEXDATA)(arr+offset);
545:         PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",rank,numComponents,bus->internal_i);
546:       } else if (key == 2) {
547:         gen = (GEN)(arr+offset);
548:         PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",rank,gen->pg,gen->qg);
549:       } else if (key == 3) {
550:         load = (LOAD)(arr+offset);
551:         PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pd = %f qd = %f\n",rank,load->pl,load->ql);
552:       }
553:     }
554:   }
555: #endif  
556:   /* Broadcast Sbase to all processors */
557:   MPI_Bcast(&User.Sbase,1,MPI_DOUBLE,0,PETSC_COMM_WORLD);

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

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

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

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

577:   SNESSolve(snes,NULL,X);

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

583:   SNESDestroy(&snes);
584:   DMDestroy(&networkdm);

586:   PetscFinalize();
587:   return 0;
588: }