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: #include "power.h"
9: #include <petscdmnetwork.h>
11: PetscErrorCode FormFunction_Subnet(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
12: {
13: UserCtx_Power *User = (UserCtx_Power*)appctx;
14: PetscInt e,v,vfrom,vto;
15: const PetscScalar *xarr;
16: PetscScalar *farr;
17: PetscInt offsetfrom,offsetto,offset;
19: VecGetArrayRead(localX,&xarr);
20: VecGetArray(localF,&farr);
22: for (v=0; v<nv; v++) {
23: PetscInt i,j,key;
24: PetscScalar Vm;
25: PetscScalar Sbase = User->Sbase;
26: VERTEX_Power bus = NULL;
27: GEN gen;
28: LOAD load;
29: PetscBool ghostvtex;
30: PetscInt numComps;
31: void* component;
33: DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
34: DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
35: DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset);
36: for (j = 0; j < numComps; j++) {
37: DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component,NULL);
38: if (key == 1) {
39: PetscInt nconnedges;
40: const PetscInt *connedges;
42: bus = (VERTEX_Power)(component);
43: /* Handle reference bus constrained dofs */
44: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
45: farr[offset] = xarr[offset] - bus->va*PETSC_PI/180.0;
46: farr[offset+1] = xarr[offset+1] - bus->vm;
47: break;
48: }
50: if (!ghostvtex) {
51: Vm = xarr[offset+1];
53: /* Shunt injections */
54: farr[offset] += Vm*Vm*bus->gl/Sbase;
55: if (bus->ide != PV_BUS) farr[offset+1] += -Vm*Vm*bus->bl/Sbase;
56: }
58: DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
59: for (i=0; i < nconnedges; i++) {
60: EDGE_Power branch;
61: PetscInt keye;
62: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
63: const PetscInt *cone;
64: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
66: e = connedges[i];
67: DMNetworkGetComponent(networkdm,e,0,&keye,(void**)&branch,NULL);
68: if (!branch->status) continue;
69: Gff = branch->yff[0];
70: Bff = branch->yff[1];
71: Gft = branch->yft[0];
72: Bft = branch->yft[1];
73: Gtf = branch->ytf[0];
74: Btf = branch->ytf[1];
75: Gtt = branch->ytt[0];
76: Btt = branch->ytt[1];
78: DMNetworkGetConnectedVertices(networkdm,e,&cone);
79: vfrom = cone[0];
80: vto = cone[1];
82: DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
83: DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
85: thetaf = xarr[offsetfrom];
86: Vmf = xarr[offsetfrom+1];
87: thetat = xarr[offsetto];
88: Vmt = xarr[offsetto+1];
89: thetaft = thetaf - thetat;
90: thetatf = thetat - thetaf;
92: if (vfrom == vtx[v]) {
93: farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft));
94: farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
95: } else {
96: farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf));
97: farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
98: }
99: }
100: } else if (key == 2) {
101: if (!ghostvtex) {
102: gen = (GEN)(component);
103: if (!gen->status) continue;
104: farr[offset] += -gen->pg/Sbase;
105: farr[offset+1] += -gen->qg/Sbase;
106: }
107: } else if (key == 3) {
108: if (!ghostvtex) {
109: load = (LOAD)(component);
110: farr[offset] += load->pl/Sbase;
111: farr[offset+1] += load->ql/Sbase;
112: }
113: }
114: }
115: if (bus && bus->ide == PV_BUS) {
116: farr[offset+1] = xarr[offset+1] - bus->vm;
117: }
118: }
119: VecRestoreArrayRead(localX,&xarr);
120: VecRestoreArray(localF,&farr);
121: return 0;
122: }
124: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
125: {
126: DM networkdm;
127: Vec localX,localF;
128: PetscInt nv,ne;
129: const PetscInt *vtx,*edges;
131: SNESGetDM(snes,&networkdm);
132: DMGetLocalVector(networkdm,&localX);
133: DMGetLocalVector(networkdm,&localF);
134: VecSet(F,0.0);
136: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
137: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
139: DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
140: DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);
142: /* Form Function for first subnetwork */
143: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
144: FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);
146: /* Form Function for second subnetwork */
147: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
148: FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);
150: DMRestoreLocalVector(networkdm,&localX);
152: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
153: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
154: DMRestoreLocalVector(networkdm,&localF);
155: return 0;
156: }
158: PetscErrorCode FormJacobian_Subnet(DM networkdm,Vec localX, Mat J, Mat Jpre, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
159: {
160: UserCtx_Power *User=(UserCtx_Power*)appctx;
161: PetscInt e,v,vfrom,vto;
162: const PetscScalar *xarr;
163: PetscInt offsetfrom,offsetto,goffsetfrom,goffsetto;
164: PetscInt row[2],col[8];
165: PetscScalar values[8];
167: VecGetArrayRead(localX,&xarr);
169: for (v=0; v<nv; v++) {
170: PetscInt i,j,key;
171: PetscInt offset,goffset;
172: PetscScalar Vm;
173: PetscScalar Sbase=User->Sbase;
174: VERTEX_Power bus;
175: PetscBool ghostvtex;
176: PetscInt numComps;
177: void* component;
179: DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
180: DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
181: for (j = 0; j < numComps; j++) {
182: DMNetworkGetLocalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&offset);
183: DMNetworkGetGlobalVecOffset(networkdm,vtx[v],ALL_COMPONENTS,&goffset);
184: DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component,NULL);
185: if (key == 1) {
186: PetscInt nconnedges;
187: const PetscInt *connedges;
189: bus = (VERTEX_Power)(component);
190: if (!ghostvtex) {
191: /* Handle reference bus constrained dofs */
192: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
193: row[0] = goffset; row[1] = goffset+1;
194: col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
195: values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
196: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
197: break;
198: }
200: Vm = xarr[offset+1];
202: /* Shunt injections */
203: row[0] = goffset; row[1] = goffset+1;
204: col[0] = goffset; col[1] = goffset+1;
205: values[0] = values[1] = values[2] = values[3] = 0.0;
206: if (bus->ide != PV_BUS) {
207: values[1] = 2.0*Vm*bus->gl/Sbase;
208: values[3] = -2.0*Vm*bus->bl/Sbase;
209: }
210: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
211: }
213: DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
214: for (i=0; i < nconnedges; i++) {
215: EDGE_Power branch;
216: VERTEX_Power busf,bust;
217: PetscInt keyf,keyt;
218: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
219: const PetscInt *cone;
220: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
222: e = connedges[i];
223: DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch,NULL);
224: if (!branch->status) continue;
226: Gff = branch->yff[0];
227: Bff = branch->yff[1];
228: Gft = branch->yft[0];
229: Bft = branch->yft[1];
230: Gtf = branch->ytf[0];
231: Btf = branch->ytf[1];
232: Gtt = branch->ytt[0];
233: Btt = branch->ytt[1];
235: DMNetworkGetConnectedVertices(networkdm,e,&cone);
236: vfrom = cone[0];
237: vto = cone[1];
239: DMNetworkGetLocalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&offsetfrom);
240: DMNetworkGetLocalVecOffset(networkdm,vto,ALL_COMPONENTS,&offsetto);
241: DMNetworkGetGlobalVecOffset(networkdm,vfrom,ALL_COMPONENTS,&goffsetfrom);
242: DMNetworkGetGlobalVecOffset(networkdm,vto,ALL_COMPONENTS,&goffsetto);
244: if (goffsetto < 0) goffsetto = -goffsetto - 1;
246: thetaf = xarr[offsetfrom];
247: Vmf = xarr[offsetfrom+1];
248: thetat = xarr[offsetto];
249: Vmt = xarr[offsetto+1];
250: thetaft = thetaf - thetat;
251: thetatf = thetat - thetaf;
253: DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf,NULL);
254: DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust,NULL);
256: if (vfrom == vtx[v]) {
257: if (busf->ide != REF_BUS) {
258: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */
259: row[0] = goffsetfrom;
260: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
261: values[0] = Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
262: values[1] = 2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
263: values[2] = Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
264: values[3] = Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */
266: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
267: }
268: if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
269: row[0] = goffsetfrom+1;
270: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
271: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
272: values[0] = Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
273: values[1] = -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
274: values[2] = Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
275: values[3] = Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
277: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
278: }
279: } else {
280: if (bust->ide != REF_BUS) {
281: row[0] = goffsetto;
282: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
283: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
284: values[0] = Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
285: values[1] = 2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
286: values[2] = Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
287: values[3] = Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */
289: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
290: }
291: if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
292: row[0] = goffsetto+1;
293: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
294: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
295: values[0] = Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
296: values[1] = -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
297: values[2] = Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
298: values[3] = Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
300: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
301: }
302: }
303: }
304: if (!ghostvtex && bus->ide == PV_BUS) {
305: row[0] = goffset+1; col[0] = goffset+1;
306: values[0] = 1.0;
307: MatSetValues(J,1,row,1,col,values,ADD_VALUES);
308: }
309: }
310: }
311: }
312: VecRestoreArrayRead(localX,&xarr);
313: return 0;
314: }
316: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
317: {
318: DM networkdm;
319: Vec localX;
320: PetscInt ne,nv;
321: const PetscInt *vtx,*edges;
323: MatZeroEntries(J);
325: SNESGetDM(snes,&networkdm);
326: DMGetLocalVector(networkdm,&localX);
328: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
329: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
331: /* Form Jacobian for first subnetwork */
332: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
333: FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);
335: /* Form Jacobian for second subnetwork */
336: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
337: FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);
339: DMRestoreLocalVector(networkdm,&localX);
341: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
342: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
343: return 0;
344: }
346: PetscErrorCode SetInitialValues_Subnet(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx)
347: {
348: VERTEX_Power bus;
349: PetscInt i;
350: GEN gen;
351: PetscBool ghostvtex;
352: PetscScalar *xarr;
353: PetscInt key,numComps,j,offset;
354: void* component;
356: VecGetArray(localX,&xarr);
357: for (i = 0; i < nv; i++) {
358: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
359: if (ghostvtex) continue;
361: DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
362: DMNetworkGetNumComponents(networkdm,vtx[i],&numComps);
363: for (j=0; j < numComps; j++) {
364: DMNetworkGetComponent(networkdm,vtx[i],j,&key,&component,NULL);
365: if (key == 1) {
366: bus = (VERTEX_Power)(component);
367: xarr[offset] = bus->va*PETSC_PI/180.0;
368: xarr[offset+1] = bus->vm;
369: } else if (key == 2) {
370: gen = (GEN)(component);
371: if (!gen->status) continue;
372: xarr[offset+1] = gen->vs;
373: break;
374: }
375: }
376: }
377: VecRestoreArray(localX,&xarr);
378: return 0;
379: }
381: PetscErrorCode SetInitialValues(DM networkdm, Vec X,void* appctx)
382: {
383: PetscInt nv,ne;
384: const PetscInt *vtx,*edges;
385: Vec localX;
387: DMGetLocalVector(networkdm,&localX);
389: VecSet(X,0.0);
390: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
391: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
393: /* Set initial guess for first subnetwork */
394: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
395: SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);
397: /* Set initial guess for second subnetwork */
398: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
399: SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);
401: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
402: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
403: DMRestoreLocalVector(networkdm,&localX);
404: return 0;
405: }
407: int main(int argc,char ** argv)
408: {
409: char pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
410: PFDATA *pfdata1,*pfdata2;
411: PetscInt numEdges1=0,numEdges2=0;
412: PetscInt *edgelist1 = NULL,*edgelist2 = NULL,componentkey[4];
413: DM networkdm;
414: UserCtx_Power User;
415: #if defined(PETSC_USE_LOG)
416: PetscLogStage stage1,stage2;
417: #endif
418: PetscMPIInt rank;
419: PetscInt nsubnet = 2,nv,ne,i,j,genj,loadj;
420: const PetscInt *vtx,*edges;
421: Vec X,F;
422: Mat J;
423: SNES snes;
425: PetscInitialize(&argc,&argv,"poweroptions",help);
426: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
427: {
428: /* 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 */
429: /* this is an experiment to see how the analyzer reacts */
430: const PetscMPIInt crank = rank;
432: /* Create an empty network object */
433: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
435: /* Register the components in the network */
436: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0]);
437: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1]);
438: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
439: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);
441: PetscLogStageRegister("Read Data",&stage1);
442: PetscLogStagePush(stage1);
443: /* READ THE DATA */
444: if (!crank) {
445: /* Only rank 0 reads the data */
446: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);
447: /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */
449: /* READ DATA FOR THE FIRST SUBNETWORK */
450: PetscNew(&pfdata1);
451: PFReadMatPowerData(pfdata1,pfdata_file);
452: User.Sbase = pfdata1->sbase;
454: numEdges1 = pfdata1->nbranch;
455: PetscMalloc1(2*numEdges1,&edgelist1);
456: GetListofEdges_Power(pfdata1,edgelist1);
458: /* READ DATA FOR THE SECOND SUBNETWORK */
459: PetscNew(&pfdata2);
460: PFReadMatPowerData(pfdata2,pfdata_file);
461: User.Sbase = pfdata2->sbase;
463: numEdges2 = pfdata2->nbranch;
464: PetscMalloc1(2*numEdges2,&edgelist2);
465: GetListofEdges_Power(pfdata2,edgelist2);
466: }
468: PetscLogStagePop();
469: MPI_Barrier(PETSC_COMM_WORLD);
470: PetscLogStageRegister("Create network",&stage2);
471: PetscLogStagePush(stage2);
473: /* Set number of nodes/edges and edge connectivity */
474: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,nsubnet);
475: DMNetworkAddSubnetwork(networkdm,"",numEdges1,edgelist1,NULL);
476: DMNetworkAddSubnetwork(networkdm,"",numEdges2,edgelist2,NULL);
478: /* Set up the network layout */
479: DMNetworkLayoutSetUp(networkdm);
481: /* Add network components only process 0 has any data to add*/
482: if (!crank) {
483: genj=0; loadj=0;
485: /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
486: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
488: for (i = 0; i < ne; i++) {
489: DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata1->branch[i],0);
490: }
492: for (i = 0; i < nv; i++) {
493: DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata1->bus[i],2);
494: if (pfdata1->bus[i].ngen) {
495: for (j = 0; j < pfdata1->bus[i].ngen; j++) {
496: DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata1->gen[genj++],0);
497: }
498: }
499: if (pfdata1->bus[i].nload) {
500: for (j=0; j < pfdata1->bus[i].nload; j++) {
501: DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata1->load[loadj++],0);
502: }
503: }
504: }
506: genj=0; loadj=0;
508: /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
509: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
511: for (i = 0; i < ne; i++) {
512: DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata2->branch[i],0);
513: }
515: for (i = 0; i < nv; i++) {
516: DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata2->bus[i],2);
517: if (pfdata2->bus[i].ngen) {
518: for (j = 0; j < pfdata2->bus[i].ngen; j++) {
519: DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata2->gen[genj++],0);
520: }
521: }
522: if (pfdata2->bus[i].nload) {
523: for (j=0; j < pfdata2->bus[i].nload; j++) {
524: DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata2->load[loadj++],0);
525: }
526: }
527: }
528: }
530: /* Set up DM for use */
531: DMSetUp(networkdm);
533: if (!crank) {
534: PetscFree(edgelist1);
535: PetscFree(edgelist2);
536: }
538: if (!crank) {
539: PetscFree(pfdata1->bus);
540: PetscFree(pfdata1->gen);
541: PetscFree(pfdata1->branch);
542: PetscFree(pfdata1->load);
543: PetscFree(pfdata1);
545: PetscFree(pfdata2->bus);
546: PetscFree(pfdata2->gen);
547: PetscFree(pfdata2->branch);
548: PetscFree(pfdata2->load);
549: PetscFree(pfdata2);
550: }
552: /* Distribute networkdm to multiple processes */
553: DMNetworkDistribute(&networkdm,0);
555: PetscLogStagePop();
557: /* Broadcast Sbase to all processors */
558: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
560: DMCreateGlobalVector(networkdm,&X);
561: VecDuplicate(X,&F);
563: DMCreateMatrix(networkdm,&J);
564: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
566: SetInitialValues(networkdm,X,&User);
568: /* HOOK UP SOLVER */
569: SNESCreate(PETSC_COMM_WORLD,&snes);
570: SNESSetDM(snes,networkdm);
571: SNESSetFunction(snes,F,FormFunction,&User);
572: SNESSetJacobian(snes,J,J,FormJacobian,&User);
573: SNESSetFromOptions(snes);
575: SNESSolve(snes,NULL,X);
576: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
578: VecDestroy(&X);
579: VecDestroy(&F);
580: MatDestroy(&J);
582: SNESDestroy(&snes);
583: DMDestroy(&networkdm);
584: }
585: PetscFinalize();
586: return 0;
587: }
589: /*TEST
591: build:
592: depends: PFReadData.c pffunctions.c
593: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
595: test:
596: args: -snes_rtol 1.e-3
597: localrunfiles: poweroptions case9.m
598: output_file: output/power_1.out
600: test:
601: suffix: 2
602: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
603: nsize: 4
604: localrunfiles: poweroptions case9.m
605: output_file: output/power_1.out
607: TEST*/