Actual source code: pf.c
petsc-3.5.4 2015-05-23
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] = 0.0;
91: farr[offset+1] = 0.0;
92: break;
93: }
95: if (!ghostvtex) {
96: Vm = xarr[offset+1];
98: /* Shunt injections */
99: farr[offset] += Vm*Vm*bus->gl/Sbase;
100: 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] = 0.0;
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; col[2] = goffset; col[3] = goffset+1;
240: values[0] = 2*Vm*bus->gl/Sbase;
241: values[1] = values[2] = 0.0;
242: values[3] = -2*Vm*bus->bl/Sbase;
243: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
244: }
246: PetscInt nconnedges;
247: const PetscInt *connedges;
249: DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
250: for (i=0; i < nconnedges; i++) {
251: EDGEDATA branch;
252: VERTEXDATA busf,bust;
253: PetscInt offsetfd,offsettd,keyf,keyt;
254: e = connedges[i];
255: DMNetworkGetComponentTypeOffset(networkdm,e,0,&key,&offsetd);
256: branch = (EDGEDATA)(arr+offsetd);
257: if (!branch->status) continue;
258: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
259: Gff = branch->yff[0];
260: Bff = branch->yff[1];
261: Gft = branch->yft[0];
262: Bft = branch->yft[1];
263: Gtf = branch->ytf[0];
264: Btf = branch->ytf[1];
265: Gtt = branch->ytt[0];
266: Btt = branch->ytt[1];
268: const PetscInt *cone;
269: DMNetworkGetConnectedNodes(networkdm,e,&cone);
270: vfrom = cone[0];
271: vto = cone[1];
273: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
274: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
275: DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);
276: DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);
278: if (goffsetfrom < 0) goffsetfrom = -goffsetfrom - 1; /* Convert to actual global offset for ghost nodes, global offset is -(gstart+1) */
279: if (goffsetto < 0) goffsetto = -goffsetto - 1;
280: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
282: thetaf = xarr[offsetfrom];
283: Vmf = xarr[offsetfrom+1];
284: thetat = xarr[offsetto];
285: Vmt = xarr[offsetto+1];
286: thetaft = thetaf - thetat;
287: thetatf = thetat - thetaf;
289: DMNetworkGetComponentTypeOffset(networkdm,vfrom,0,&keyf,&offsetfd);
290: DMNetworkGetComponentTypeOffset(networkdm,vto,0,&keyt,&offsettd);
291: busf = (VERTEXDATA)(arr+offsetfd);
292: bust = (VERTEXDATA)(arr+offsettd);
294: if (vfrom == v) {
295: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*cos(thetaft) + Bft*sin(thetaft)); */
296: row[0] = goffsetfrom;
297: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
298: values[0] = Vmf*Vmt*(Gft*-sin(thetaft) + Bft*cos(thetaft)); /* df_dthetaf */
299: values[1] = 2*Gff*Vmf + Vmt*(Gft*cos(thetaft) + Bft*sin(thetaft)); /* df_dVmf */
300: values[2] = Vmf*Vmt*(Gft*sin(thetaft) + Bft*-cos(thetaft)); /* df_dthetat */
301: values[3] = Vmf*(Gft*cos(thetaft) + Bft*sin(thetaft)); /* df_dVmt */
302:
303: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
304:
305: if (busf->ide != PV_BUS) {
306: row[0] = goffsetfrom+1;
307: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
308: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*cos(thetaft) + Gft*sin(thetaft)); */
309: values[0] = Vmf*Vmt*(Bft*sin(thetaft) + Gft*cos(thetaft));
310: values[1] = -2*Bff*Vmf + Vmt*(-Bft*cos(thetaft) + Gft*sin(thetaft));
311: values[2] = Vmf*Vmt*(-Bft*sin(thetaft) + Gft*-cos(thetaft));
312: values[3] = Vmf*(-Bft*cos(thetaft) + Gft*sin(thetaft));
313:
314: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
315: }
316: } else {
317: row[0] = goffsetto;
318: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
319: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*cos(thetatf) + Btf*sin(thetatf)); */
320: values[0] = Vmt*Vmf*(Gtf*-sin(thetatf) + Btf*cos(thetaft)); /* df_dthetat */
321: values[1] = 2*Gtt*Vmt + Vmf*(Gtf*cos(thetatf) + Btf*sin(thetatf)); /* df_dVmt */
322: values[2] = Vmt*Vmf*(Gtf*sin(thetatf) + Btf*-cos(thetatf)); /* df_dthetaf */
323: values[3] = Vmt*(Gtf*cos(thetatf) + Btf*sin(thetatf)); /* df_dVmf */
324:
325: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
326:
327: if (bust->ide != PV_BUS) {
328: row[0] = goffsetto+1;
329: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
330: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*cos(thetatf) + Gtf*sin(thetatf)); */
331: values[0] = Vmt*Vmf*(Btf*sin(thetatf) + Gtf*cos(thetatf));
332: values[1] = -2*Btt*Vmt + Vmf*(-Btf*cos(thetatf) + Gtf*sin(thetatf));
333: values[2] = Vmt*Vmf*(-Btf*sin(thetatf) + Gtf*-cos(thetatf));
334: values[3] = Vmt*(-Btf*cos(thetatf) + Gtf*sin(thetatf));
335:
336: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
337: }
338: }
339: }
340: if (!ghostvtex && bus->ide == PV_BUS) {
341: row[0] = goffset+1; col[0] = goffset+1;
342: values[0] = 1.0;
343: MatSetValues(J,1,row,1,col,values,ADD_VALUES);
344: }
345: }
346: }
347: }
348: VecRestoreArrayRead(localX,&xarr);
349: DMRestoreLocalVector(networkdm,&localX);
351: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
352: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
354: return(0);
355: }
359: PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx)
360: {
362: VERTEXDATA bus;
363: GEN gen;
364: PetscInt v, vStart, vEnd, offset;
365: PetscBool ghostvtex;
366: Vec localX;
367: PetscScalar *xarr;
368: PetscInt key;
369: DMNetworkComponentGenericDataType *arr;
370:
372: DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);
374: DMGetLocalVector(networkdm,&localX);
376: VecSet(X,0.0);
377: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
378: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
380: VecGetArray(localX,&xarr);
381: DMNetworkGetComponentDataArray(networkdm,&arr);
382: for (v = vStart; v < vEnd; v++) {
383: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
384: if (ghostvtex) continue;
385: PetscInt numComps;
386: PetscInt j;
387: PetscInt offsetd;
388: DMNetworkGetVariableOffset(networkdm,v,&offset);
389: DMNetworkGetNumComponents(networkdm,v,&numComps);
390: for (j=0; j < numComps; j++) {
391: DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);
392: if (key == 1) {
393: bus = (VERTEXDATA)(arr+offsetd);
394: xarr[offset] = bus->va*PETSC_PI/180.0;
395: xarr[offset+1] = bus->vm;
396: } else if(key == 2) {
397: gen = (GEN)(arr+offsetd);
398: if (!gen->status) continue;
399: xarr[offset+1] = gen->vs;
400: break;
401: }
402: }
403: }
404: VecRestoreArray(localX,&xarr);
405: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
406: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
407: DMRestoreLocalVector(networkdm,&localX);
408: return(0);
409: }
414: int main(int argc,char ** argv)
415: {
416: PetscErrorCode ierr;
417: char pfdata_file[PETSC_MAX_PATH_LEN]="datafiles/case9.m";
418: PFDATA pfdata;
419: PetscInt numEdges=0,numVertices=0;
420: int *edges = NULL;
421: PetscInt i;
422: DM networkdm;
423: PetscInt componentkey[4];
424: UserCtx User;
425: PetscLogStage stage1,stage2;
426: PetscInt size;
428: PetscInitialize(&argc,&argv,"pfoptions",help);
430: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
432: /* Create an empty network object */
433: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
434: /* Register the components in the network */
435: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGEDATA),&componentkey[0]);
436: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEXDATA),&componentkey[1]);
437: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
438: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);
440: PetscLogStageRegister("Read Data",&stage1);
441: PetscLogStagePush(stage1);
442: /* READ THE DATA */
443: if (!rank) {
444: /* READ DATA */
445: /* Only rank 0 reads the data */
446: PetscOptionsGetString(PETSC_NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
447: PFReadMatPowerData(&pfdata,pfdata_file);
448: User.Sbase = pfdata.sbase;
450: numEdges = pfdata.nbranch;
451: numVertices = pfdata.nbus;
453: PetscMalloc(2*numEdges*sizeof(int),&edges);
454: GetListofEdges(pfdata.nbranch,pfdata.branch,edges);
455: }
456: PetscLogStagePop();
457: MPI_Barrier(PETSC_COMM_WORLD);
458: PetscLogStageRegister("Create network",&stage2);
459: PetscLogStagePush(stage2);
460: /* Set number of nodes/edges */
461: DMNetworkSetSizes(networkdm,numVertices,numEdges,PETSC_DETERMINE,PETSC_DETERMINE);
462: /* Add edge connectivity */
463: DMNetworkSetEdgeList(networkdm,edges);
464: /* Set up the network layout */
465: DMNetworkLayoutSetUp(networkdm);
467: if (!rank) {
468: PetscFree(edges);
469: }
470: /* Add network components */
471: PetscInt eStart, eEnd, vStart, vEnd,j;
472: PetscInt genj=0,loadj=0;
473: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
474: for (i = eStart; i < eEnd; i++) {
475: DMNetworkAddComponent(networkdm,i,componentkey[0],&pfdata.branch[i-eStart]);
476: }
477: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
478: for (i = vStart; i < vEnd; i++) {
479: DMNetworkAddComponent(networkdm,i,componentkey[1],&pfdata.bus[i-vStart]);
480: if (pfdata.bus[i-vStart].ngen) {
481: for (j = 0; j < pfdata.bus[i-vStart].ngen; j++) {
482: DMNetworkAddComponent(networkdm,i,componentkey[2],&pfdata.gen[genj++]);
483: }
484: }
485: if (pfdata.bus[i-vStart].nload) {
486: for (j=0; j < pfdata.bus[i-vStart].nload; j++) {
487: DMNetworkAddComponent(networkdm,i,componentkey[3],&pfdata.load[loadj++]);
488: }
489: }
490: /* Add number of variables */
491: DMNetworkAddNumVariables(networkdm,i,2);
492: }
493: /* Set up DM for use */
494: DMSetUp(networkdm);
496: if (!rank) {
497: PetscFree(pfdata.bus);
498: PetscFree(pfdata.gen);
499: PetscFree(pfdata.branch);
500: PetscFree(pfdata.load);
501: }
504: MPI_Comm_size(PETSC_COMM_WORLD,&size);
505: if (size > 1) {
506: DM distnetworkdm;
507: /* Network partitioning and distribution of data */
508: DMNetworkDistribute(networkdm,"chaco",0,&distnetworkdm);
509: DMDestroy(&networkdm);
510: networkdm = distnetworkdm;
511: }
513: PetscLogStagePop();
514: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
515: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
516:
517: #if 0
518: PetscInt numComponents;
519: EDGEDATA edge;
520: PetscInt offset,key;
521: DMNetworkComponentGenericDataType *arr;
522: for (i = eStart; i < eEnd; i++) {
523: DMNetworkGetComponentDataArray(networkdm,&arr);
524: DMNetworkGetComponentTypeOffset(networkdm,i,0,&key,&offset);
525: edge = (EDGEDATA)(arr+offset);
526: DMNetworkGetNumComponents(networkdm,i,&numComponents);
527: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",rank,numComponents,edge->internal_i,edge->internal_j);
528: }
530: VERTEXDATA bus;
531: GEN gen;
532: LOAD load;
533: PetscInt kk;
534: for (i = vStart; i < vEnd; i++) {
535: DMNetworkGetComponentDataArray(networkdm,&arr);
536: DMNetworkGetNumComponents(networkdm,i,&numComponents);
537: for (kk=0; kk < numComponents; kk++) {
538: DMNetworkGetComponentTypeOffset(networkdm,i,kk,&key,&offset);
539: if (key == 1) {
540: bus = (VERTEXDATA)(arr+offset);
541: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",rank,numComponents,bus->internal_i);
542: } else if (key == 2) {
543: gen = (GEN)(arr+offset);
544: PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",rank,gen->pg,gen->qg);
545: } else if (key == 3) {
546: load = (LOAD)(arr+offset);
547: PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pd = %f qd = %f\n",rank,load->pl,load->ql);
548: }
549: }
550: }
551: #endif
552: /* Broadcast Sbase to all processors */
553: MPI_Bcast(&User.Sbase,1,MPI_DOUBLE,0,PETSC_COMM_WORLD);
555: Vec X,F;
556: DMCreateGlobalVector(networkdm,&X);
557: VecDuplicate(X,&F);
559: Mat J;
560: DMCreateMatrix(networkdm,&J);
561: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
563: SetInitialValues(networkdm,X,&User);
565: SNES snes;
566: /* HOOK UP SOLVER */
567: SNESCreate(PETSC_COMM_WORLD,&snes);
568: SNESSetDM(snes,networkdm);
569: SNESSetFunction(snes,F,FormFunction,&User);
570: SNESSetJacobian(snes,J,J,FormJacobian,&User);
571: SNESSetFromOptions(snes);
573: SNESSolve(snes,NULL,X);
575: VecDestroy(&X);
576: VecDestroy(&F);
577: MatDestroy(&J);
579: SNESDestroy(&snes);
580: DMDestroy(&networkdm);
582: PetscFinalize();
583: return 0;
584: }