Actual source code: pf.c
petsc-3.7.3 2016-08-01
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
14: PetscMPIInt rank;
18: PetscErrorCode GetListofEdges(PetscInt nbranches, EDGEDATA branch,int edges[])
19: {
20: PetscInt i, fbus,tbus;
23: for (i=0; i < nbranches; i++) {
24: fbus = branch[i].internal_i;
25: tbus = branch[i].internal_j;
26: edges[2*i] = fbus;
27: edges[2*i+1] = tbus;
28: }
29: return(0);
30: }
32: typedef struct{
33: PetscScalar Sbase;
34: }UserCtx;
38: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
39: {
41: DM networkdm;
42: UserCtx *User=(UserCtx*)appctx;
43: Vec localX,localF;
44: PetscInt e;
45: PetscInt v,vStart,vEnd,vfrom,vto;
46: const PetscScalar *xarr;
47: PetscScalar *farr;
48: PetscInt offsetfrom,offsetto,offset;
49: DMNetworkComponentGenericDataType *arr;
52: SNESGetDM(snes,&networkdm);
53: DMGetLocalVector(networkdm,&localX);
54: DMGetLocalVector(networkdm,&localF);
55: VecSet(F,0.0);
57: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
58: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
60: DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
61: DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);
63: VecGetArrayRead(localX,&xarr);
64: VecGetArray(localF,&farr);
66: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
67: DMNetworkGetComponentDataArray(networkdm,&arr);
69: for (v=vStart; v < vEnd; v++) {
70: PetscInt i,j,offsetd,key;
71: PetscScalar Vm;
72: PetscScalar Sbase=User->Sbase;
73: VERTEXDATA bus=NULL;
74: GEN gen;
75: LOAD load;
76: PetscBool ghostvtex;
77: PetscInt numComps;
79: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
80: DMNetworkGetNumComponents(networkdm,v,&numComps);
81: DMNetworkGetVariableOffset(networkdm,v,&offset);
82: for (j = 0; j < numComps; j++) {
83: DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);
84: if (key == 1) {
85: PetscInt nconnedges;
86: const PetscInt *connedges;
88: bus = (VERTEXDATA)(arr+offsetd);
89: /* Handle reference bus constrained dofs */
90: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
91: farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
92: farr[offset+1] = xarr[offset+1] - bus->vm;
93: break;
94: }
96: if (!ghostvtex) {
97: Vm = xarr[offset+1];
99: /* Shunt injections */
100: farr[offset] += Vm*Vm*bus->gl/Sbase;
101: if(bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
102: }
104: DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
105: for (i=0; i < nconnedges; i++) {
106: EDGEDATA branch;
107: PetscInt keye;
108: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
109: const PetscInt *cone;
110: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
112: e = connedges[i];
113: DMNetworkGetComponentTypeOffset(networkdm,e,0,&keye,&offsetd);
114: branch = (EDGEDATA)(arr+offsetd);
115: if (!branch->status) continue;
116: Gff = branch->yff[0];
117: Bff = branch->yff[1];
118: Gft = branch->yft[0];
119: Bft = branch->yft[1];
120: Gtf = branch->ytf[0];
121: Btf = branch->ytf[1];
122: Gtt = branch->ytt[0];
123: Btt = branch->ytt[1];
125: DMNetworkGetConnectedNodes(networkdm,e,&cone);
126: vfrom = cone[0];
127: vto = cone[1];
129: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
130: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
132: thetaf = xarr[offsetfrom];
133: Vmf = xarr[offsetfrom+1];
134: thetat = xarr[offsetto];
135: Vmt = xarr[offsetto+1];
136: thetaft = thetaf - thetat;
137: thetatf = thetat - thetaf;
139: if (vfrom == v) {
140: farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));
141: farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
142: } else {
143: farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf));
144: farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
145: }
146: }
147: } else if (key == 2) {
148: if (!ghostvtex) {
149: gen = (GEN)(arr+offsetd);
150: if (!gen->status) continue;
151: farr[offset] += -gen->pg/Sbase;
152: farr[offset+1] += -gen->qg/Sbase;
153: }
154: } else if (key == 3) {
155: if (!ghostvtex) {
156: load = (LOAD)(arr+offsetd);
157: farr[offset] += load->pl/Sbase;
158: farr[offset+1] += load->ql/Sbase;
159: }
160: }
161: }
162: if (bus && bus->ide == PV_BUS) {
163: farr[offset+1] = xarr[offset+1] - bus->vm;
164: }
165: }
166: VecRestoreArrayRead(localX,&xarr);
167: VecRestoreArray(localF,&farr);
168: DMRestoreLocalVector(networkdm,&localX);
170: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
171: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
172: DMRestoreLocalVector(networkdm,&localF);
173: return(0);
174: }
178: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
179: {
181: DM networkdm;
182: UserCtx *User=(UserCtx*)appctx;
183: Vec localX;
184: PetscInt e;
185: PetscInt v,vStart,vEnd,vfrom,vto;
186: const PetscScalar *xarr;
187: PetscInt offsetfrom,offsetto,goffsetfrom,goffsetto;
188: DMNetworkComponentGenericDataType *arr;
189: PetscInt row[2],col[8];
190: PetscScalar values[8];
193: MatZeroEntries(J);
195: SNESGetDM(snes,&networkdm);
196: DMGetLocalVector(networkdm,&localX);
198: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
199: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
201: VecGetArrayRead(localX,&xarr);
203: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
204: DMNetworkGetComponentDataArray(networkdm,&arr);
206: for (v=vStart; v < vEnd; v++) {
207: PetscInt i,j,offsetd,key;
208: PetscInt offset,goffset;
209: PetscScalar Vm;
210: PetscScalar Sbase=User->Sbase;
211: VERTEXDATA bus;
212: PetscBool ghostvtex;
213: PetscInt numComps;
215: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
216: DMNetworkGetNumComponents(networkdm,v,&numComps);
217: for (j = 0; j < numComps; j++) {
218: DMNetworkGetVariableOffset(networkdm,v,&offset);
219: DMNetworkGetVariableGlobalOffset(networkdm,v,&goffset);
220: DMNetworkGetComponentTypeOffset(networkdm,v,j,&key,&offsetd);
221: if (key == 1) {
222: PetscInt nconnedges;
223: const PetscInt *connedges;
225: bus = (VERTEXDATA)(arr+offsetd);
226: if (!ghostvtex) {
227: /* Handle reference bus constrained dofs */
228: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
229: row[0] = goffset; row[1] = goffset+1;
230: col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
231: values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
232: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
233: break;
234: }
235:
236: Vm = xarr[offset+1];
237:
238: /* Shunt injections */
239: row[0] = goffset; row[1] = goffset+1;
240: col[0] = goffset; col[1] = goffset+1;
241: values[0] = values[1] = values[2] = values[3] = 0.0;
242: if (bus->ide != PV_BUS) {
243: values[1] = 2.0*Vm*bus->gl/Sbase;
244: values[3] = -2.0*Vm*bus->bl/Sbase;
245: }
246: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
247: }
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: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
255: const PetscInt *cone;
256: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
258: e = connedges[i];
259: DMNetworkGetComponentTypeOffset(networkdm,e,0,&key,&offsetd);
260: branch = (EDGEDATA)(arr+offsetd);
261: if (!branch->status) continue;
262:
263: Gff = branch->yff[0];
264: Bff = branch->yff[1];
265: Gft = branch->yft[0];
266: Bft = branch->yft[1];
267: Gtf = branch->ytf[0];
268: Btf = branch->ytf[1];
269: Gtt = branch->ytt[0];
270: Btt = branch->ytt[1];
272: DMNetworkGetConnectedNodes(networkdm,e,&cone);
273: vfrom = cone[0];
274: vto = cone[1];
276: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
277: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
278: DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);
279: DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);
281: if (goffsetto < 0) goffsetto = -goffsetto - 1;
283: thetaf = xarr[offsetfrom];
284: Vmf = xarr[offsetfrom+1];
285: thetat = xarr[offsetto];
286: Vmt = xarr[offsetto+1];
287: thetaft = thetaf - thetat;
288: thetatf = thetat - thetaf;
290: DMNetworkGetComponentTypeOffset(networkdm,vfrom,0,&keyf,&offsetfd);
291: DMNetworkGetComponentTypeOffset(networkdm,vto,0,&keyt,&offsettd);
292: busf = (VERTEXDATA)(arr+offsetfd);
293: bust = (VERTEXDATA)(arr+offsettd);
295: if (vfrom == v) {
296: if (busf->ide != REF_BUS) {
297: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */
298: row[0] = goffsetfrom;
299: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
300: values[0] = Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
301: values[1] = 2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
302: values[2] = Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
303: values[3] = Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */
304:
305: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
306: }
307: if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
308: row[0] = goffsetfrom+1;
309: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
310: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
311: values[0] = Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
312: values[1] = -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
313: values[2] = Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
314: values[3] = Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
315:
316: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
317: }
318: } else {
319: if (bust->ide != REF_BUS) {
320: row[0] = goffsetto;
321: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
322: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
323: values[0] = Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
324: values[1] = 2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
325: values[2] = Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
326: values[3] = Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */
327:
328: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
329: }
330: if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
331: row[0] = goffsetto+1;
332: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
333: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
334: values[0] = Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
335: values[1] = -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
336: values[2] = Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
337: values[3] = Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
338:
339: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
340: }
341: }
342: }
343: if (!ghostvtex && bus->ide == PV_BUS) {
344: row[0] = goffset+1; col[0] = goffset+1;
345: values[0] = 1.0;
346: MatSetValues(J,1,row,1,col,values,ADD_VALUES);
347: }
348: }
349: }
350: }
351: VecRestoreArrayRead(localX,&xarr);
352: DMRestoreLocalVector(networkdm,&localX);
354: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
355: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
356: return(0);
357: }
361: PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx)
362: {
364: VERTEXDATA bus;
365: GEN gen;
366: PetscInt v, vStart, vEnd, offset;
367: PetscBool ghostvtex;
368: Vec localX;
369: PetscScalar *xarr;
370: PetscInt key,numComps,j,offsetd;
371: DMNetworkComponentGenericDataType *arr;
372:
374: DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);
376: DMGetLocalVector(networkdm,&localX);
378: VecSet(X,0.0);
379: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
380: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
382: VecGetArray(localX,&xarr);
383: DMNetworkGetComponentDataArray(networkdm,&arr);
384: for (v = vStart; v < vEnd; v++) {
385: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
386: if (ghostvtex) continue;
387:
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: {
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: PetscMPIInt size;
427: PetscInt eStart, eEnd, vStart, vEnd,j;
428: PetscInt genj,loadj;
429: Vec X,F;
430: Mat J;
431: SNES snes;
433: 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(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
451: PetscNew(&pfdata);
452: PFReadMatPowerData(pfdata,pfdata_file);
453: User.Sbase = pfdata->sbase;
455: numEdges = pfdata->nbranch;
456: numVertices = pfdata->nbus;
458: PetscMalloc(2*numEdges*sizeof(int),&edges);
459: GetListofEdges(pfdata->nbranch,pfdata->branch,edges);
460: }
461: PetscLogStagePop();
462: MPI_Barrier(PETSC_COMM_WORLD);
463: PetscLogStageRegister("Create network",&stage2);
464: PetscLogStagePush(stage2);
465: /* Set number of nodes/edges */
466: DMNetworkSetSizes(networkdm,numVertices,numEdges,PETSC_DETERMINE,PETSC_DETERMINE);
467: /* Add edge connectivity */
468: DMNetworkSetEdgeList(networkdm,edges);
469: /* Set up the network layout */
470: DMNetworkLayoutSetUp(networkdm);
472: if (!rank) {
473: PetscFree(edges);
474: }
475: /* Add network components */
476:
477: genj=0; loadj=0;
478: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
479: for (i = eStart; i < eEnd; i++) {
480: DMNetworkAddComponent(networkdm,i,componentkey[0],&pfdata->branch[i-eStart]);
481: }
482: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
483: for (i = vStart; i < vEnd; i++) {
484: DMNetworkAddComponent(networkdm,i,componentkey[1],&pfdata->bus[i-vStart]);
485: if (pfdata->bus[i-vStart].ngen) {
486: for (j = 0; j < pfdata->bus[i-vStart].ngen; j++) {
487: DMNetworkAddComponent(networkdm,i,componentkey[2],&pfdata->gen[genj++]);
488: }
489: }
490: if (pfdata->bus[i-vStart].nload) {
491: for (j=0; j < pfdata->bus[i-vStart].nload; j++) {
492: DMNetworkAddComponent(networkdm,i,componentkey[3],&pfdata->load[loadj++]);
493: }
494: }
495: /* Add number of variables */
496: DMNetworkAddNumVariables(networkdm,i,2);
497: }
498: /* Set up DM for use */
499: DMSetUp(networkdm);
501: if (!rank) {
502: PetscFree(pfdata->bus);
503: PetscFree(pfdata->gen);
504: PetscFree(pfdata->branch);
505: PetscFree(pfdata->load);
506: PetscFree(pfdata);
507: }
510: MPI_Comm_size(PETSC_COMM_WORLD,&size);
511: if (size > 1) {
512: DM distnetworkdm;
513: /* Network partitioning and distribution of data */
514: DMNetworkDistribute(networkdm,0,&distnetworkdm);
515: DMDestroy(&networkdm);
516: networkdm = distnetworkdm;
517: }
519: PetscLogStagePop();
520: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
521: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
522:
523: #if 0
524: PetscInt numComponents;
525: EDGEDATA edge;
526: PetscInt offset,key,kk;
527: DMNetworkComponentGenericDataType *arr;
528: VERTEXDATA bus;
529: GEN gen;
530: LOAD load;
531:
532: for (i = eStart; i < eEnd; i++) {
533: DMNetworkGetComponentDataArray(networkdm,&arr);
534: DMNetworkGetComponentTypeOffset(networkdm,i,0,&key,&offset);
535: edge = (EDGEDATA)(arr+offset);
536: DMNetworkGetNumComponents(networkdm,i,&numComponents);
537: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",rank,numComponents,edge->internal_i,edge->internal_j);
538: }
540: for (i = vStart; i < vEnd; i++) {
541: DMNetworkGetComponentDataArray(networkdm,&arr);
542: DMNetworkGetNumComponents(networkdm,i,&numComponents);
543: for (kk=0; kk < numComponents; kk++) {
544: DMNetworkGetComponentTypeOffset(networkdm,i,kk,&key,&offset);
545: if (key == 1) {
546: bus = (VERTEXDATA)(arr+offset);
547: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",rank,numComponents,bus->internal_i);
548: } else if (key == 2) {
549: gen = (GEN)(arr+offset);
550: PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",rank,gen->pg,gen->qg);
551: } else if (key == 3) {
552: load = (LOAD)(arr+offset);
553: PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",rank,load->pl,load->ql);
554: }
555: }
556: }
557: #endif
558: /* Broadcast Sbase to all processors */
559: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
561: DMCreateGlobalVector(networkdm,&X);
562: VecDuplicate(X,&F);
564: DMCreateMatrix(networkdm,&J);
565: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
567: SetInitialValues(networkdm,X,&User);
569: /* HOOK UP SOLVER */
570: SNESCreate(PETSC_COMM_WORLD,&snes);
571: SNESSetDM(snes,networkdm);
572: SNESSetFunction(snes,F,FormFunction,&User);
573: SNESSetJacobian(snes,J,J,FormJacobian,&User);
574: SNESSetFromOptions(snes);
576: SNESSolve(snes,NULL,X);
578: VecDestroy(&X);
579: VecDestroy(&F);
580: MatDestroy(&J);
582: SNESDestroy(&snes);
583: DMDestroy(&networkdm);
585: PetscFinalize();
586: return 0;
587: }