Actual source code: power2.c
1: static char help[] = "This example demonstrates the use of DMNetwork interface with subnetworks for solving a nonlinear electric power grid problem.\n\
2: The available solver options are in the poweroptions 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: This example shows the use of subnetwork feature in DMNetwork. It creates duplicates of the same network which are treated as subnetworks.\n\
5: Run this program: mpiexec -n <n> ./pf2\n\
6: mpiexec -n <n> ./pf2 \n";
8: /* T
9: Concepts: DMNetwork
10: Concepts: PETSc SNES solver
11: */
13: #include "power.h"
14: #include <petscdmnetwork.h>
16: PetscErrorCode FormFunction_Subnet(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
17: {
18: PetscErrorCode ierr;
19: UserCtx_Power *User = (UserCtx_Power*)appctx;
20: PetscInt e,v,vfrom,vto;
21: const PetscScalar *xarr;
22: PetscScalar *farr;
23: PetscInt offsetfrom,offsetto,offset;
26: VecGetArrayRead(localX,&xarr);
27: VecGetArray(localF,&farr);
29: for (v=0; v<nv; v++) {
30: PetscInt i,j,key;
31: PetscScalar Vm;
32: PetscScalar Sbase = User->Sbase;
33: VERTEX_Power bus = NULL;
34: GEN gen;
35: LOAD load;
36: PetscBool ghostvtex;
37: PetscInt numComps;
38: void* component;
40: DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
41: DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
42: DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset);
43: for (j = 0; j < numComps; j++) {
44: DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component,NULL);
45: if (key == 1) {
46: PetscInt nconnedges;
47: const PetscInt *connedges;
49: bus = (VERTEX_Power)(component);
50: /* Handle reference bus constrained dofs */
51: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
52: farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
53: farr[offset+1] = xarr[offset+1] - bus->vm;
54: break;
55: }
57: if (!ghostvtex) {
58: Vm = xarr[offset+1];
60: /* Shunt injections */
61: farr[offset] += Vm*Vm*bus->gl/Sbase;
62: if (bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
63: }
65: DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
66: for (i=0; i < nconnedges; i++) {
67: EDGE_Power branch;
68: PetscInt keye;
69: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
70: const PetscInt *cone;
71: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
73: e = connedges[i];
74: DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);
75: if (!branch->status) continue;
76: Gff = branch->yff[0];
77: Bff = branch->yff[1];
78: Gft = branch->yft[0];
79: Bft = branch->yft[1];
80: Gtf = branch->ytf[0];
81: Btf = branch->ytf[1];
82: Gtt = branch->ytt[0];
83: Btt = branch->ytt[1];
85: DMNetworkGetConnectedVertices(networkdm,e,&cone);
86: vfrom = cone[0];
87: vto = cone[1];
89: DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
90: DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
92: thetaf = xarr[offsetfrom];
93: Vmf = xarr[offsetfrom+1];
94: thetat = xarr[offsetto];
95: Vmt = xarr[offsetto+1];
96: thetaft = thetaf - thetat;
97: thetatf = thetat - thetaf;
99: if (vfrom == vtx[v]) {
100: farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));
101: farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
102: } else {
103: farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf));
104: farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
105: }
106: }
107: } else if (key == 2) {
108: if (!ghostvtex) {
109: gen = (GEN)(component);
110: if (!gen->status) continue;
111: farr[offset] += -gen->pg/Sbase;
112: farr[offset+1] += -gen->qg/Sbase;
113: }
114: } else if (key == 3) {
115: if (!ghostvtex) {
116: load = (LOAD)(component);
117: farr[offset] += load->pl/Sbase;
118: farr[offset+1] += load->ql/Sbase;
119: }
120: }
121: }
122: if (bus && bus->ide == PV_BUS) {
123: farr[offset+1] = xarr[offset+1] - bus->vm;
124: }
125: }
126: VecRestoreArrayRead(localX,&xarr);
127: VecRestoreArray(localF,&farr);
128: return(0);
129: }
131: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
132: {
134: DM networkdm;
135: Vec localX,localF;
136: PetscInt nv,ne;
137: const PetscInt *vtx,*edges;
140: SNESGetDM(snes,&networkdm);
141: DMGetLocalVector(networkdm,&localX);
142: DMGetLocalVector(networkdm,&localF);
143: VecSet(F,0.0);
145: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
146: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
148: DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
149: DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);
151: /* Form Function for first subnetwork */
152: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
153: FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);
155: /* Form Function for second subnetwork */
156: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
157: FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);
159: DMRestoreLocalVector(networkdm,&localX);
161: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
162: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
163: DMRestoreLocalVector(networkdm,&localF);
164: return(0);
165: }
167: PetscErrorCode FormJacobian_Subnet(DM networkdm,Vec localX, Mat J, Mat Jpre, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
168: {
169: PetscErrorCode ierr;
170: UserCtx_Power *User=(UserCtx_Power*)appctx;
171: PetscInt e,v,vfrom,vto;
172: const PetscScalar *xarr;
173: PetscInt offsetfrom,offsetto,goffsetfrom,goffsetto;
174: PetscInt row[2],col[8];
175: PetscScalar values[8];
178: VecGetArrayRead(localX,&xarr);
180: for (v=0; v<nv; v++) {
181: PetscInt i,j,key;
182: PetscInt offset,goffset;
183: PetscScalar Vm;
184: PetscScalar Sbase=User->Sbase;
185: VERTEX_Power bus;
186: PetscBool ghostvtex;
187: PetscInt numComps;
188: void* component;
190: DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
191: DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
192: for (j = 0; j < numComps; j++) {
193: DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset);
194: DMNetworkGetGlobalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&goffset);
195: DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component,NULL);
196: if (key == 1) {
197: PetscInt nconnedges;
198: const PetscInt *connedges;
200: bus = (VERTEX_Power)(component);
201: if (!ghostvtex) {
202: /* Handle reference bus constrained dofs */
203: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
204: row[0] = goffset; row[1] = goffset+1;
205: col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
206: values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
207: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
208: break;
209: }
211: Vm = xarr[offset+1];
213: /* Shunt injections */
214: row[0] = goffset; row[1] = goffset+1;
215: col[0] = goffset; col[1] = goffset+1;
216: values[0] = values[1] = values[2] = values[3] = 0.0;
217: if (bus->ide != PV_BUS) {
218: values[1] = 2.0*Vm*bus->gl/Sbase;
219: values[3] = -2.0*Vm*bus->bl/Sbase;
220: }
221: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
222: }
224: DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
225: for (i=0; i < nconnedges; i++) {
226: EDGE_Power branch;
227: VERTEX_Power busf,bust;
228: PetscInt keyf,keyt;
229: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
230: const PetscInt *cone;
231: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
233: e = connedges[i];
234: DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
235: if (!branch->status) continue;
237: Gff = branch->yff[0];
238: Bff = branch->yff[1];
239: Gft = branch->yft[0];
240: Bft = branch->yft[1];
241: Gtf = branch->ytf[0];
242: Btf = branch->ytf[1];
243: Gtt = branch->ytt[0];
244: Btt = branch->ytt[1];
246: DMNetworkGetConnectedVertices(networkdm,e,&cone);
247: vfrom = cone[0];
248: vto = cone[1];
250: DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
251: DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
252: DMNetworkGetGlobalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&goffsetfrom);
253: DMNetworkGetGlobalVecOffset(networkdm,vto,ALL_COMPONENTS,&goffsetto);
255: if (goffsetto < 0) goffsetto = -goffsetto - 1;
257: thetaf = xarr[offsetfrom];
258: Vmf = xarr[offsetfrom+1];
259: thetat = xarr[offsetto];
260: Vmt = xarr[offsetto+1];
261: thetaft = thetaf - thetat;
262: thetatf = thetat - thetaf;
264: DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf,NULL);
265: DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust,NULL);
267: if (vfrom == vtx[v]) {
268: if (busf->ide != REF_BUS) {
269: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */
270: row[0] = goffsetfrom;
271: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
272: values[0] = Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
273: values[1] = 2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
274: values[2] = Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
275: values[3] = Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */
277: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
278: }
279: if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
280: row[0] = goffsetfrom+1;
281: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
282: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
283: values[0] = Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
284: values[1] = -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
285: values[2] = Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
286: values[3] = Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
288: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
289: }
290: } else {
291: if (bust->ide != REF_BUS) {
292: row[0] = goffsetto;
293: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
294: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
295: values[0] = Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
296: values[1] = 2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
297: values[2] = Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
298: values[3] = Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */
300: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
301: }
302: if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
303: row[0] = goffsetto+1;
304: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
305: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
306: values[0] = Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
307: values[1] = -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
308: values[2] = Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
309: values[3] = Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
311: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
312: }
313: }
314: }
315: if (!ghostvtex && bus->ide == PV_BUS) {
316: row[0] = goffset+1; col[0] = goffset+1;
317: values[0] = 1.0;
318: MatSetValues(J,1,row,1,col,values,ADD_VALUES);
319: }
320: }
321: }
322: }
323: VecRestoreArrayRead(localX,&xarr);
324: return(0);
325: }
327: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
328: {
330: DM networkdm;
331: Vec localX;
332: PetscInt ne,nv;
333: const PetscInt *vtx,*edges;
336: MatZeroEntries(J);
338: SNESGetDM(snes,&networkdm);
339: DMGetLocalVector(networkdm,&localX);
341: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
342: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
344: /* Form Jacobian for first subnetwork */
345: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
346: FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);
348: /* Form Jacobian for second subnetwork */
349: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
350: FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);
352: DMRestoreLocalVector(networkdm,&localX);
354: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
355: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
356: return(0);
357: }
359: PetscErrorCode SetInitialValues_Subnet(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx)
360: {
362: VERTEX_Power bus;
363: PetscInt i;
364: GEN gen;
365: PetscBool ghostvtex;
366: PetscScalar *xarr;
367: PetscInt key,numComps,j,offset;
368: void* component;
371: VecGetArray(localX,&xarr);
372: for (i = 0; i < nv; i++) {
373: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
374: if (ghostvtex) continue;
376: DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
377: DMNetworkGetNumComponents(networkdm,vtx[i],&numComps);
378: for (j=0; j < numComps; j++) {
379: DMNetworkGetComponent(networkdm,vtx[i],j,&key,&component,NULL);
380: if (key == 1) {
381: bus = (VERTEX_Power)(component);
382: xarr[offset] = bus->va*PETSC_PI/180.0;
383: xarr[offset+1] = bus->vm;
384: } else if (key == 2) {
385: gen = (GEN)(component);
386: if (!gen->status) continue;
387: xarr[offset+1] = gen->vs;
388: break;
389: }
390: }
391: }
392: VecRestoreArray(localX,&xarr);
393: return(0);
394: }
396: PetscErrorCode SetInitialValues(DM networkdm, Vec X,void* appctx)
397: {
399: PetscInt nv,ne;
400: const PetscInt *vtx,*edges;
401: Vec localX;
404: DMGetLocalVector(networkdm,&localX);
406: VecSet(X,0.0);
407: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
408: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
410: /* Set initial guess for first subnetwork */
411: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
412: SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);
414: /* Set initial guess for second subnetwork */
415: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
416: SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);
418: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
419: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
420: DMRestoreLocalVector(networkdm,&localX);
421: return(0);
422: }
424: int main(int argc,char ** argv)
425: {
426: PetscErrorCode ierr;
427: char pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
428: PFDATA *pfdata1,*pfdata2;
429: PetscInt numEdges1=0,numEdges2=0;
430: PetscInt *edgelist1 = NULL,*edgelist2 = NULL,componentkey[4];
431: DM networkdm;
432: UserCtx_Power User;
433: #if defined(PETSC_USE_LOG)
434: PetscLogStage stage1,stage2;
435: #endif
436: PetscMPIInt rank;
437: PetscInt nsubnet = 2,nv,ne,i,j,genj,loadj;
438: const PetscInt *vtx,*edges;
439: Vec X,F;
440: Mat J;
441: SNES snes;
443: PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr;
444: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
445: {
446: /* 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 */
447: /* this is an experiment to see how the analyzer reacts */
448: const PetscMPIInt crank = rank;
450: /* Create an empty network object */
451: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
453: /* Register the components in the network */
454: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0]);
455: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1]);
456: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
457: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);
459: PetscLogStageRegister("Read Data",&stage1);
460: PetscLogStagePush(stage1);
461: /* READ THE DATA */
462: if (!crank) {
463: /* Only rank 0 reads the data */
464: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);
465: /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */
467: /* READ DATA FOR THE FIRST SUBNETWORK */
468: PetscNew(&pfdata1);
469: PFReadMatPowerData(pfdata1,pfdata_file);
470: User.Sbase = pfdata1->sbase;
472: numEdges1 = pfdata1->nbranch;
473: PetscMalloc1(2*numEdges1,&edgelist1);
474: GetListofEdges_Power(pfdata1,edgelist1);
476: /* READ DATA FOR THE SECOND SUBNETWORK */
477: PetscNew(&pfdata2);
478: PFReadMatPowerData(pfdata2,pfdata_file);
479: User.Sbase = pfdata2->sbase;
481: numEdges2 = pfdata2->nbranch;
482: PetscMalloc1(2*numEdges2,&edgelist2);
483: GetListofEdges_Power(pfdata2,edgelist2);
484: }
486: PetscLogStagePop();
487: MPI_Barrier(PETSC_COMM_WORLD);
488: PetscLogStageRegister("Create network",&stage2);
489: PetscLogStagePush(stage2);
491: /* Set number of nodes/edges and edge connectivity */
492: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,nsubnet);
493: DMNetworkAddSubnetwork(networkdm,"",numEdges1,edgelist1,NULL);
494: DMNetworkAddSubnetwork(networkdm,"",numEdges2,edgelist2,NULL);
496: /* Set up the network layout */
497: DMNetworkLayoutSetUp(networkdm);
499: /* Add network components only process 0 has any data to add*/
500: if (!crank) {
501: genj=0; loadj=0;
503: /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
504: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
506: for (i = 0; i < ne; i++) {
507: DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata1->branch[i],0);
508: }
510: for (i = 0; i < nv; i++) {
511: DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata1->bus[i],2);
512: if (pfdata1->bus[i].ngen) {
513: for (j = 0; j < pfdata1->bus[i].ngen; j++) {
514: DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata1->gen[genj++],0);
515: }
516: }
517: if (pfdata1->bus[i].nload) {
518: for (j=0; j < pfdata1->bus[i].nload; j++) {
519: DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata1->load[loadj++],0);
520: }
521: }
522: }
524: genj=0; loadj=0;
526: /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
527: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
529: for (i = 0; i < ne; i++) {
530: DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata2->branch[i],0);
531: }
533: for (i = 0; i < nv; i++) {
534: DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata2->bus[i],2);
535: if (pfdata2->bus[i].ngen) {
536: for (j = 0; j < pfdata2->bus[i].ngen; j++) {
537: DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata2->gen[genj++],0);
538: }
539: }
540: if (pfdata2->bus[i].nload) {
541: for (j=0; j < pfdata2->bus[i].nload; j++) {
542: DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata2->load[loadj++],0);
543: }
544: }
545: }
546: }
548: /* Set up DM for use */
549: DMSetUp(networkdm);
551: if (!crank) {
552: PetscFree(edgelist1);
553: PetscFree(edgelist2);
554: }
556: if (!crank) {
557: PetscFree(pfdata1->bus);
558: PetscFree(pfdata1->gen);
559: PetscFree(pfdata1->branch);
560: PetscFree(pfdata1->load);
561: PetscFree(pfdata1);
563: PetscFree(pfdata2->bus);
564: PetscFree(pfdata2->gen);
565: PetscFree(pfdata2->branch);
566: PetscFree(pfdata2->load);
567: PetscFree(pfdata2);
568: }
570: /* Distribute networkdm to multiple processes */
571: DMNetworkDistribute(&networkdm,0);
573: PetscLogStagePop();
575: /* Broadcast Sbase to all processors */
576: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
578: DMCreateGlobalVector(networkdm,&X);
579: VecDuplicate(X,&F);
581: DMCreateMatrix(networkdm,&J);
582: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
584: SetInitialValues(networkdm,X,&User);
586: /* HOOK UP SOLVER */
587: SNESCreate(PETSC_COMM_WORLD,&snes);
588: SNESSetDM(snes,networkdm);
589: SNESSetFunction(snes,F,FormFunction,&User);
590: SNESSetJacobian(snes,J,J,FormJacobian,&User);
591: SNESSetFromOptions(snes);
593: SNESSolve(snes,NULL,X);
594: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
596: VecDestroy(&X);
597: VecDestroy(&F);
598: MatDestroy(&J);
600: SNESDestroy(&snes);
601: DMDestroy(&networkdm);
602: }
603: PetscFinalize();
604: return ierr;
605: }
607: /*TEST
609: build:
610: depends: PFReadData.c pffunctions.c
611: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
613: test:
614: args: -snes_rtol 1.e-3
615: localrunfiles: poweroptions case9.m
616: output_file: output/power_1.out
618: test:
619: suffix: 2
620: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
621: nsize: 4
622: localrunfiles: poweroptions case9.m
623: output_file: output/power_1.out
625: TEST*/