Actual source code: pf.c
petsc-3.8.4 2018-03-24
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 \n";
7: /* T
8: Concepts: DMNetwork
9: Concepts: PETSc SNES solver
10: */
12: #include "pf.h"
13: #include <petscdmnetwork.h>
15: PetscErrorCode GetListofEdges(PetscInt nbranches, EDGEDATA branch,int edges[])
16: {
17: PetscInt i, fbus,tbus;
20: for (i=0; i < nbranches; i++) {
21: fbus = branch[i].internal_i;
22: tbus = branch[i].internal_j;
23: edges[2*i] = fbus;
24: edges[2*i+1] = tbus;
25: }
26: return(0);
27: }
29: typedef struct{
30: PetscScalar Sbase;
31: PetscBool jac_error; /* introduce error in the jacobian */
32: }UserCtx;
34: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
35: {
37: DM networkdm;
38: UserCtx *User=(UserCtx*)appctx;
39: Vec localX,localF;
40: PetscInt e;
41: PetscInt v,vStart,vEnd,vfrom,vto;
42: const PetscScalar *xarr;
43: PetscScalar *farr;
44: PetscInt offsetfrom,offsetto,offset;
45: DMNetworkComponentGenericDataType *arr;
48: SNESGetDM(snes,&networkdm);
49: DMGetLocalVector(networkdm,&localX);
50: DMGetLocalVector(networkdm,&localF);
51: VecSet(F,0.0);
53: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
54: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
56: DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
57: DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);
59: VecGetArrayRead(localX,&xarr);
60: VecGetArray(localF,&farr);
62: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
63: DMNetworkGetComponentDataArray(networkdm,&arr);
65: for (v=vStart; v < vEnd; v++) {
66: PetscInt i,j,offsetd,key;
67: PetscScalar Vm;
68: PetscScalar Sbase=User->Sbase;
69: VERTEXDATA bus=NULL;
70: GEN gen;
71: LOAD load;
72: PetscBool ghostvtex;
73: PetscInt numComps;
75: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
76: DMNetworkGetNumComponents(networkdm,v,&numComps);
77: DMNetworkGetVariableOffset(networkdm,v,&offset);
78: for (j = 0; j < numComps; j++) {
79: DMNetworkGetComponentKeyOffset(networkdm,v,j,&key,&offsetd);
80: if (key == 1) {
81: PetscInt nconnedges;
82: const PetscInt *connedges;
84: bus = (VERTEXDATA)(arr+offsetd);
85: /* Handle reference bus constrained dofs */
86: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
87: farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
88: farr[offset+1] = xarr[offset+1] - bus->vm;
89: break;
90: }
92: if (!ghostvtex) {
93: Vm = xarr[offset+1];
95: /* Shunt injections */
96: farr[offset] += Vm*Vm*bus->gl/Sbase;
97: if(bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
98: }
100: DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
101: for (i=0; i < nconnedges; i++) {
102: EDGEDATA branch;
103: PetscInt keye;
104: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
105: const PetscInt *cone;
106: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
108: e = connedges[i];
109: DMNetworkGetComponentKeyOffset(networkdm,e,0,&keye,&offsetd);
110: branch = (EDGEDATA)(arr+offsetd);
111: if (!branch->status) continue;
112: Gff = branch->yff[0];
113: Bff = branch->yff[1];
114: Gft = branch->yft[0];
115: Bft = branch->yft[1];
116: Gtf = branch->ytf[0];
117: Btf = branch->ytf[1];
118: Gtt = branch->ytt[0];
119: Btt = branch->ytt[1];
121: DMNetworkGetConnectedVertices(networkdm,e,&cone);
122: vfrom = cone[0];
123: vto = cone[1];
125: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
126: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
128: thetaf = xarr[offsetfrom];
129: Vmf = xarr[offsetfrom+1];
130: thetat = xarr[offsetto];
131: Vmt = xarr[offsetto+1];
132: thetaft = thetaf - thetat;
133: thetatf = thetat - thetaf;
135: if (vfrom == v) {
136: farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));
137: farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
138: } else {
139: farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf));
140: farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
141: }
142: }
143: } else if (key == 2) {
144: if (!ghostvtex) {
145: gen = (GEN)(arr+offsetd);
146: if (!gen->status) continue;
147: farr[offset] += -gen->pg/Sbase;
148: farr[offset+1] += -gen->qg/Sbase;
149: }
150: } else if (key == 3) {
151: if (!ghostvtex) {
152: load = (LOAD)(arr+offsetd);
153: farr[offset] += load->pl/Sbase;
154: farr[offset+1] += load->ql/Sbase;
155: }
156: }
157: }
158: if (bus && bus->ide == PV_BUS) {
159: farr[offset+1] = xarr[offset+1] - bus->vm;
160: }
161: }
162: VecRestoreArrayRead(localX,&xarr);
163: VecRestoreArray(localF,&farr);
164: DMRestoreLocalVector(networkdm,&localX);
166: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
167: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
168: DMRestoreLocalVector(networkdm,&localF);
169: return(0);
170: }
172: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
173: {
175: DM networkdm;
176: UserCtx *User=(UserCtx*)appctx;
177: Vec localX;
178: PetscInt e;
179: PetscInt v,vStart,vEnd,vfrom,vto;
180: const PetscScalar *xarr;
181: PetscInt offsetfrom,offsetto,goffsetfrom,goffsetto;
182: DMNetworkComponentGenericDataType *arr;
183: PetscInt row[2],col[8];
184: PetscScalar values[8];
187: MatZeroEntries(J);
189: SNESGetDM(snes,&networkdm);
190: DMGetLocalVector(networkdm,&localX);
192: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
193: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
195: VecGetArrayRead(localX,&xarr);
197: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
198: DMNetworkGetComponentDataArray(networkdm,&arr);
200: for (v=vStart; v < vEnd; v++) {
201: PetscInt i,j,offsetd,key;
202: PetscInt offset,goffset;
203: PetscScalar Vm;
204: PetscScalar Sbase=User->Sbase;
205: VERTEXDATA bus;
206: PetscBool ghostvtex;
207: PetscInt numComps;
209: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
210: DMNetworkGetNumComponents(networkdm,v,&numComps);
211: for (j = 0; j < numComps; j++) {
212: DMNetworkGetVariableOffset(networkdm,v,&offset);
213: DMNetworkGetVariableGlobalOffset(networkdm,v,&goffset);
214: DMNetworkGetComponentKeyOffset(networkdm,v,j,&key,&offsetd);
215: if (key == 1) {
216: PetscInt nconnedges;
217: const PetscInt *connedges;
219: bus = (VERTEXDATA)(arr+offsetd);
220: if (!ghostvtex) {
221: /* Handle reference bus constrained dofs */
222: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
223: row[0] = goffset; row[1] = goffset+1;
224: col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
225: values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
226: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
227: break;
228: }
229:
230: Vm = xarr[offset+1];
231:
232: /* Shunt injections */
233: row[0] = goffset; row[1] = goffset+1;
234: col[0] = goffset; col[1] = goffset+1;
235: values[0] = values[1] = values[2] = values[3] = 0.0;
236: if (bus->ide != PV_BUS) {
237: values[1] = 2.0*Vm*bus->gl/Sbase;
238: values[3] = -2.0*Vm*bus->bl/Sbase;
239: }
240: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
241: }
243: DMNetworkGetSupportingEdges(networkdm,v,&nconnedges,&connedges);
244: for (i=0; i < nconnedges; i++) {
245: EDGEDATA branch;
246: VERTEXDATA busf,bust;
247: PetscInt offsetfd,offsettd,keyf,keyt;
248: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
249: const PetscInt *cone;
250: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
252: e = connedges[i];
253: DMNetworkGetComponentKeyOffset(networkdm,e,0,&key,&offsetd);
254: branch = (EDGEDATA)(arr+offsetd);
255: if (!branch->status) continue;
256:
257: Gff = branch->yff[0];
258: Bff = branch->yff[1];
259: Gft = branch->yft[0];
260: Bft = branch->yft[1];
261: Gtf = branch->ytf[0];
262: Btf = branch->ytf[1];
263: Gtt = branch->ytt[0];
264: Btt = branch->ytt[1];
266: DMNetworkGetConnectedVertices(networkdm,e,&cone);
267: vfrom = cone[0];
268: vto = cone[1];
270: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
271: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
272: DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);
273: DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);
275: if (goffsetto < 0) goffsetto = -goffsetto - 1;
277: thetaf = xarr[offsetfrom];
278: Vmf = xarr[offsetfrom+1];
279: thetat = xarr[offsetto];
280: Vmt = xarr[offsetto+1];
281: thetaft = thetaf - thetat;
282: thetatf = thetat - thetaf;
284: DMNetworkGetComponentKeyOffset(networkdm,vfrom,0,&keyf,&offsetfd);
285: DMNetworkGetComponentKeyOffset(networkdm,vto,0,&keyt,&offsettd);
286: busf = (VERTEXDATA)(arr+offsetfd);
287: bust = (VERTEXDATA)(arr+offsettd);
289: if (vfrom == v) {
290: if (busf->ide != REF_BUS) {
291: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */
292: row[0] = goffsetfrom;
293: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
294: values[0] = Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
295: values[1] = 2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
296: values[2] = Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
297: values[3] = Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */
298:
299: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
300: }
301: if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
302: row[0] = goffsetfrom+1;
303: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
304: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
305: values[0] = Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
306: values[1] = -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
307: values[2] = Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
308: values[3] = Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
309:
310: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
311: }
312: } else {
313: if (bust->ide != REF_BUS) {
314: row[0] = goffsetto;
315: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
316: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
317: values[0] = Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
318: values[1] = 2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
319: values[2] = Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
320: values[3] = Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */
321:
322: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
323: }
324: if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
325: row[0] = goffsetto+1;
326: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
327: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
328: values[0] = Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
329: values[1] = -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
330: values[2] = Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
331: values[3] = Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
332:
333: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
334: }
335: }
336: }
337: if (!ghostvtex && bus->ide == PV_BUS) {
338: row[0] = goffset+1; col[0] = goffset+1;
339: values[0] = 1.0;
340: if (User->jac_error) values[0] = 50.0;
341: MatSetValues(J,1,row,1,col,values,ADD_VALUES);
342: }
343: }
344: }
345: }
346: VecRestoreArrayRead(localX,&xarr);
347: DMRestoreLocalVector(networkdm,&localX);
349: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
350: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
351: return(0);
352: }
354: PetscErrorCode SetInitialValues(DM networkdm,Vec X,void* appctx)
355: {
357: VERTEXDATA bus;
358: GEN gen;
359: PetscInt v, vStart, vEnd, offset;
360: PetscBool ghostvtex;
361: Vec localX;
362: PetscScalar *xarr;
363: PetscInt key,numComps,j,offsetd;
364: DMNetworkComponentGenericDataType *arr;
365:
367: DMNetworkGetVertexRange(networkdm,&vStart, &vEnd);
369: DMGetLocalVector(networkdm,&localX);
371: VecSet(X,0.0);
372: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
373: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
375: VecGetArray(localX,&xarr);
376: DMNetworkGetComponentDataArray(networkdm,&arr);
377: for (v = vStart; v < vEnd; v++) {
378: DMNetworkIsGhostVertex(networkdm,v,&ghostvtex);
379: if (ghostvtex) continue;
381: DMNetworkGetVariableOffset(networkdm,v,&offset);
382: DMNetworkGetNumComponents(networkdm,v,&numComps);
383: for (j=0; j < numComps; j++) {
384: DMNetworkGetComponentKeyOffset(networkdm,v,j,&key,&offsetd);
385: if (key == 1) {
386: bus = (VERTEXDATA)(arr+offsetd);
387: xarr[offset] = bus->va*PETSC_PI/180.0;
388: xarr[offset+1] = bus->vm;
389: } else if(key == 2) {
390: gen = (GEN)(arr+offsetd);
391: if (!gen->status) continue;
392: xarr[offset+1] = gen->vs;
393: break;
394: }
395: }
396: }
397: VecRestoreArray(localX,&xarr);
398: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
399: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
400: DMRestoreLocalVector(networkdm,&localX);
401: return(0);
402: }
405: int main(int argc,char ** argv)
406: {
407: PetscErrorCode ierr;
408: char pfdata_file[PETSC_MAX_PATH_LEN]="datafiles/case9.m";
409: PFDATA *pfdata;
410: PetscInt numEdges=0,numVertices=0;
411: int *edges = NULL;
412: PetscInt i;
413: DM networkdm;
414: PetscInt componentkey[4];
415: UserCtx User;
416: PetscLogStage stage1,stage2;
417: PetscMPIInt rank;
418: PetscInt eStart, eEnd, vStart, vEnd,j;
419: PetscInt genj,loadj;
420: Vec X,F;
421: Mat J;
422: SNES snes;
424: PetscInitialize(&argc,&argv,"pfoptions",help);
425: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
426: {
427: /* 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 */
428: /* this is an experiment to see how the analyzer reacts */
429: const PetscMPIInt crank = rank;
431: /* Create an empty network object */
432: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
433: /* Register the components in the network */
434: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGEDATA),&componentkey[0]);
435: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEXDATA),&componentkey[1]);
436: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
437: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);
439: PetscLogStageRegister("Read Data",&stage1);
440: PetscLogStagePush(stage1);
441: /* READ THE DATA */
442: if (!crank) {
443: /* READ DATA */
444: /* Only rank 0 reads the data */
445: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
446: PetscNew(&pfdata);
447: PFReadMatPowerData(pfdata,pfdata_file);
448: User.Sbase = pfdata->sbase;
450: numEdges = pfdata->nbranch;
451: numVertices = pfdata->nbus;
453: PetscMalloc1(2*numEdges,&edges);
454: GetListofEdges(pfdata->nbranch,pfdata->branch,edges);
455: }
457: /* If external option activated. Introduce error in jacobian */
458: PetscOptionsHasName(NULL,NULL, "-jac_error", &User.jac_error);
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 (!crank) {
472: PetscFree(edges);
473: }
475: /* Add network components only process 0 has any data to add*/
476: if (!crank) {
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: }
500: /* Set up DM for use */
501: DMSetUp(networkdm);
503: if (!crank) {
504: PetscFree(pfdata->bus);
505: PetscFree(pfdata->gen);
506: PetscFree(pfdata->branch);
507: PetscFree(pfdata->load);
508: PetscFree(pfdata);
509: }
511: /* Distribute networkdm to multiple processes */
512: DMNetworkDistribute(&networkdm,0);
514: PetscLogStagePop();
515: DMNetworkGetEdgeRange(networkdm,&eStart,&eEnd);
516: DMNetworkGetVertexRange(networkdm,&vStart,&vEnd);
518: #if 0
519: PetscInt numComponents;
520: EDGEDATA edge;
521: PetscInt offset,key,kk;
522: DMNetworkComponentGenericDataType *arr;
523: VERTEXDATA bus;
524: GEN gen;
525: LOAD load;
527: for (i = eStart; i < eEnd; i++) {
528: DMNetworkGetComponentDataArray(networkdm,&arr);
529: DMNetworkGetComponentKeyOffset(networkdm,i,0,&key,&offset);
530: edge = (EDGEDATA)(arr+offset);
531: DMNetworkGetNumComponents(networkdm,i,&numComponents);
532: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Line %d ---- %d\n",crank,numComponents,edge->internal_i,edge->internal_j);
533: }
535: for (i = vStart; i < vEnd; i++) {
536: DMNetworkGetComponentDataArray(networkdm,&arr);
537: DMNetworkGetNumComponents(networkdm,i,&numComponents);
538: for (kk=0; kk < numComponents; kk++) {
539: DMNetworkGetComponentKeyOffset(networkdm,i,kk,&key,&offset);
540: if (key == 1) {
541: bus = (VERTEXDATA)(arr+offset);
542: PetscPrintf(PETSC_COMM_SELF,"Rank %d ncomps = %d Bus %d\n",crank,numComponents,bus->internal_i);
543: } else if (key == 2) {
544: gen = (GEN)(arr+offset);
545: PetscPrintf(PETSC_COMM_SELF,"Rank %d Gen pg = %f qg = %f\n",crank,gen->pg,gen->qg);
546: } else if (key == 3) {
547: load = (LOAD)(arr+offset);
548: PetscPrintf(PETSC_COMM_SELF,"Rank %d Load pl = %f ql = %f\n",crank,load->pl,load->ql);
549: }
550: }
551: }
552: #endif
553: /* Broadcast Sbase to all processors */
554: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
556: DMCreateGlobalVector(networkdm,&X);
557: VecDuplicate(X,&F);
559: DMCreateMatrix(networkdm,&J);
560: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
561:
562: SetInitialValues(networkdm,X,&User);
563:
564: /* HOOK UP SOLVER */
565: SNESCreate(PETSC_COMM_WORLD,&snes);
566: SNESSetDM(snes,networkdm);
567: SNESSetFunction(snes,F,FormFunction,&User);
568: SNESSetJacobian(snes,J,J,FormJacobian,&User);
569: SNESSetFromOptions(snes);
570:
571: SNESSolve(snes,NULL,X);
572: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
573:
574: VecDestroy(&X);
575: VecDestroy(&F);
576: MatDestroy(&J);
577:
578: SNESDestroy(&snes);
579: DMDestroy(&networkdm);
580: }
581: PetscFinalize();
582: return ierr;
583: }