Actual source code: pf.c
petsc-3.6.4 2016-04-12
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: }