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