Actual source code: power2.c
petsc-3.12.2 2019-11-22
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: {
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: DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);
43: for (j = 0; j < numComps; j++) {
44: DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);
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);
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: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
90: DMNetworkGetVariableOffset(networkdm,vto,&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: }
132: PetscErrorCode FormFunction(SNES snes,Vec X, Vec F,void *appctx)
133: {
135: DM networkdm;
136: Vec localX,localF;
137: PetscInt nv,ne;
138: const PetscInt *vtx,*edges;
141: SNESGetDM(snes,&networkdm);
142: DMGetLocalVector(networkdm,&localX);
143: DMGetLocalVector(networkdm,&localF);
144: VecSet(F,0.0);
146: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
147: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
149: DMGlobalToLocalBegin(networkdm,F,INSERT_VALUES,localF);
150: DMGlobalToLocalEnd(networkdm,F,INSERT_VALUES,localF);
152: /* Form Function for first subnetwork */
153: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
154: FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);
156: /* Form Function for second subnetwork */
157: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
158: FormFunction_Subnet(networkdm,localX,localF,nv,ne,vtx,edges,appctx);
160: DMRestoreLocalVector(networkdm,&localX);
162: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
163: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
164: DMRestoreLocalVector(networkdm,&localF);
165: return(0);
166: }
168: PetscErrorCode FormJacobian_Subnet(DM networkdm,Vec localX, Mat J, Mat Jpre, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
169: {
171: UserCtx_Power *User=(UserCtx_Power*)appctx;
172: PetscInt e,v,vfrom,vto;
173: const PetscScalar *xarr;
174: PetscInt offsetfrom,offsetto,goffsetfrom,goffsetto;
175: PetscInt row[2],col[8];
176: PetscScalar values[8];
179: VecGetArrayRead(localX,&xarr);
181: for (v=0; v<nv; v++) {
182: PetscInt i,j,key;
183: PetscInt offset,goffset;
184: PetscScalar Vm;
185: PetscScalar Sbase=User->Sbase;
186: VERTEX_Power bus;
187: PetscBool ghostvtex;
188: PetscInt numComps;
189: void* component;
191: DMNetworkIsGhostVertex(networkdm,vtx[v],&ghostvtex);
192: DMNetworkGetNumComponents(networkdm,vtx[v],&numComps);
193: for (j = 0; j < numComps; j++) {
194: DMNetworkGetVariableOffset(networkdm,vtx[v],&offset);
195: DMNetworkGetVariableGlobalOffset(networkdm,vtx[v],&goffset);
196: DMNetworkGetComponent(networkdm,vtx[v],j,&key,&component);
197: if (key == 1) {
198: PetscInt nconnedges;
199: const PetscInt *connedges;
201: bus = (VERTEX_Power)(component);
202: if (!ghostvtex) {
203: /* Handle reference bus constrained dofs */
204: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
205: row[0] = goffset; row[1] = goffset+1;
206: col[0] = goffset; col[1] = goffset+1; col[2] = goffset; col[3] = goffset+1;
207: values[0] = 1.0; values[1] = 0.0; values[2] = 0.0; values[3] = 1.0;
208: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
209: break;
210: }
212: Vm = xarr[offset+1];
214: /* Shunt injections */
215: row[0] = goffset; row[1] = goffset+1;
216: col[0] = goffset; col[1] = goffset+1;
217: values[0] = values[1] = values[2] = values[3] = 0.0;
218: if (bus->ide != PV_BUS) {
219: values[1] = 2.0*Vm*bus->gl/Sbase;
220: values[3] = -2.0*Vm*bus->bl/Sbase;
221: }
222: MatSetValues(J,2,row,2,col,values,ADD_VALUES);
223: }
225: DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
226: for (i=0; i < nconnedges; i++) {
227: EDGE_Power branch;
228: VERTEX_Power busf,bust;
229: PetscInt keyf,keyt;
230: PetscScalar Gff,Bff,Gft,Bft,Gtf,Btf,Gtt,Btt;
231: const PetscInt *cone;
232: PetscScalar Vmf,Vmt,thetaf,thetat,thetaft,thetatf;
234: e = connedges[i];
235: DMNetworkGetComponent(networkdm,e,0,&key,(void**)&branch);
236: if (!branch->status) continue;
237:
238: Gff = branch->yff[0];
239: Bff = branch->yff[1];
240: Gft = branch->yft[0];
241: Bft = branch->yft[1];
242: Gtf = branch->ytf[0];
243: Btf = branch->ytf[1];
244: Gtt = branch->ytt[0];
245: Btt = branch->ytt[1];
247: DMNetworkGetConnectedVertices(networkdm,e,&cone);
248: vfrom = cone[0];
249: vto = cone[1];
251: DMNetworkGetVariableOffset(networkdm,vfrom,&offsetfrom);
252: DMNetworkGetVariableOffset(networkdm,vto,&offsetto);
253: DMNetworkGetVariableGlobalOffset(networkdm,vfrom,&goffsetfrom);
254: DMNetworkGetVariableGlobalOffset(networkdm,vto,&goffsetto);
256: if (goffsetto < 0) goffsetto = -goffsetto - 1;
258: thetaf = xarr[offsetfrom];
259: Vmf = xarr[offsetfrom+1];
260: thetat = xarr[offsetto];
261: Vmt = xarr[offsetto+1];
262: thetaft = thetaf - thetat;
263: thetatf = thetat - thetaf;
265: DMNetworkGetComponent(networkdm,vfrom,0,&keyf,(void**)&busf);
266: DMNetworkGetComponent(networkdm,vto,0,&keyt,(void**)&bust);
268: if (vfrom == vtx[v]) {
269: if (busf->ide != REF_BUS) {
270: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */
271: row[0] = goffsetfrom;
272: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
273: values[0] = Vmf*Vmt*(Gft*-PetscSinScalar(thetaft) + Bft*PetscCosScalar(thetaft)); /* df_dthetaf */
274: values[1] = 2.0*Gff*Vmf + Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmf */
275: values[2] = Vmf*Vmt*(Gft*PetscSinScalar(thetaft) + Bft*-PetscCosScalar(thetaft)); /* df_dthetat */
276: values[3] = Vmf*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); /* df_dVmt */
277:
278: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
279: }
280: if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
281: row[0] = goffsetfrom+1;
282: col[0] = goffsetfrom; col[1] = goffsetfrom+1; col[2] = goffsetto; col[3] = goffsetto+1;
283: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
284: values[0] = Vmf*Vmt*(Bft*PetscSinScalar(thetaft) + Gft*PetscCosScalar(thetaft));
285: values[1] = -2.0*Bff*Vmf + Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
286: values[2] = Vmf*Vmt*(-Bft*PetscSinScalar(thetaft) + Gft*-PetscCosScalar(thetaft));
287: values[3] = Vmf*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft));
288:
289: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
290: }
291: } else {
292: if (bust->ide != REF_BUS) {
293: row[0] = goffsetto;
294: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
295: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
296: values[0] = Vmt*Vmf*(Gtf*-PetscSinScalar(thetatf) + Btf*PetscCosScalar(thetaft)); /* df_dthetat */
297: values[1] = 2.0*Gtt*Vmt + Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmt */
298: values[2] = Vmt*Vmf*(Gtf*PetscSinScalar(thetatf) + Btf*-PetscCosScalar(thetatf)); /* df_dthetaf */
299: values[3] = Vmt*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); /* df_dVmf */
300:
301: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
302: }
303: if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
304: row[0] = goffsetto+1;
305: col[0] = goffsetto; col[1] = goffsetto+1; col[2] = goffsetfrom; col[3] = goffsetfrom+1;
306: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
307: values[0] = Vmt*Vmf*(Btf*PetscSinScalar(thetatf) + Gtf*PetscCosScalar(thetatf));
308: values[1] = -2.0*Btt*Vmt + Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
309: values[2] = Vmt*Vmf*(-Btf*PetscSinScalar(thetatf) + Gtf*-PetscCosScalar(thetatf));
310: values[3] = Vmt*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf));
311:
312: MatSetValues(J,1,row,4,col,values,ADD_VALUES);
313: }
314: }
315: }
316: if (!ghostvtex && bus->ide == PV_BUS) {
317: row[0] = goffset+1; col[0] = goffset+1;
318: values[0] = 1.0;
319: MatSetValues(J,1,row,1,col,values,ADD_VALUES);
320: }
321: }
322: }
323: }
324: VecRestoreArrayRead(localX,&xarr);
325: return(0);
326: }
328: PetscErrorCode FormJacobian(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
329: {
331: DM networkdm;
332: Vec localX;
333: PetscInt ne,nv;
334: const PetscInt *vtx,*edges;
337: MatZeroEntries(J);
339: SNESGetDM(snes,&networkdm);
340: DMGetLocalVector(networkdm,&localX);
342: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
343: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
345: /* Form Jacobian for first subnetwork */
346: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
347: FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);
349: /* Form Jacobian for second subnetwork */
350: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
351: FormJacobian_Subnet(networkdm,localX,J,Jpre,nv,ne,vtx,edges,appctx);
353: DMRestoreLocalVector(networkdm,&localX);
355: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
356: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
357: return(0);
358: }
360: PetscErrorCode SetInitialValues_Subnet(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx)
361: {
363: VERTEX_Power bus;
364: PetscInt i;
365: GEN gen;
366: PetscBool ghostvtex;
367: PetscScalar *xarr;
368: PetscInt key,numComps,j,offset;
369: void* component;
372: VecGetArray(localX,&xarr);
373: for (i = 0; i < nv; i++) {
374: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
375: if (ghostvtex) continue;
377: DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);
378: DMNetworkGetNumComponents(networkdm,vtx[i],&numComps);
379: for (j=0; j < numComps; j++) {
380: DMNetworkGetComponent(networkdm,vtx[i],j,&key,&component);
381: if (key == 1) {
382: bus = (VERTEX_Power)(component);
383: xarr[offset] = bus->va*PETSC_PI/180.0;
384: xarr[offset+1] = bus->vm;
385: } else if(key == 2) {
386: gen = (GEN)(component);
387: if (!gen->status) continue;
388: xarr[offset+1] = gen->vs;
389: break;
390: }
391: }
392: }
393: VecRestoreArray(localX,&xarr);
394: return(0);
395: }
397: PetscErrorCode SetInitialValues(DM networkdm, Vec X,void* appctx)
398: {
400: PetscInt nv,ne;
401: const PetscInt *vtx,*edges;
402: Vec localX;
405: DMGetLocalVector(networkdm,&localX);
407: VecSet(X,0.0);
408: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
409: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
411: /* Set initial guess for first subnetwork */
412: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
413: SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);
415: /* Set initial guess for second subnetwork */
416: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
417: SetInitialValues_Subnet(networkdm,localX,nv,ne,vtx,edges,appctx);
419: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
420: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
421: DMRestoreLocalVector(networkdm,&localX);
422: return(0);
423: }
425: int main(int argc,char ** argv)
426: {
427: PetscErrorCode ierr;
428: char pfdata_file[PETSC_MAX_PATH_LEN]="case9.m";
429: PFDATA *pfdata1,*pfdata2;
430: PetscInt numEdges1=0,numVertices1=0,numEdges2=0,numVertices2=0;
431: PetscInt *edgelist1 = NULL,*edgelist2 = NULL;
432: DM networkdm;
433: PetscInt componentkey[4];
434: UserCtx_Power User;
435: PetscLogStage stage1,stage2;
436: PetscMPIInt rank;
437: PetscInt nsubnet = 2;
438: PetscInt numVertices[2],numEdges[2];
439: PetscInt *edgelist[2];
440: PetscInt nv,ne;
441: const PetscInt *vtx;
442: const PetscInt *edges;
443: PetscInt i,j;
444: PetscInt genj,loadj;
445: Vec X,F;
446: Mat J;
447: SNES snes;
450: PetscInitialize(&argc,&argv,"poweroptions",help);if (ierr) return ierr;
451: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
452: {
453: /* 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 */
454: /* this is an experiment to see how the analyzer reacts */
455: const PetscMPIInt crank = rank;
457: /* Create an empty network object */
458: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
460: /* Register the components in the network */
461: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&componentkey[0]);
462: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&componentkey[1]);
463: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&componentkey[2]);
464: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&componentkey[3]);
466: PetscLogStageRegister("Read Data",&stage1);
467: PetscLogStagePush(stage1);
468: /* READ THE DATA */
469: if (!crank) {
470: /* Only rank 0 reads the data */
471: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
472: /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */
474: /* READ DATA FOR THE FIRST SUBNETWORK */
475: PetscNew(&pfdata1);
476: PFReadMatPowerData(pfdata1,pfdata_file);
477: User.Sbase = pfdata1->sbase;
479: numEdges1 = pfdata1->nbranch;
480: numVertices1 = pfdata1->nbus;
482: PetscMalloc1(2*numEdges1,&edgelist1);
483: GetListofEdges_Power(pfdata1,edgelist1);
485: /* READ DATA FOR THE SECOND SUBNETWORK */
486: PetscNew(&pfdata2);
487: PFReadMatPowerData(pfdata2,pfdata_file);
488: User.Sbase = pfdata2->sbase;
490: numEdges2 = pfdata2->nbranch;
491: numVertices2 = pfdata2->nbus;
493: PetscMalloc1(2*numEdges2,&edgelist2);
494: GetListofEdges_Power(pfdata2,edgelist2);
495: }
497: PetscLogStagePop();
498: MPI_Barrier(PETSC_COMM_WORLD);
499: PetscLogStageRegister("Create network",&stage2);
500: PetscLogStagePush(stage2);
502: /* Set number of nodes/edges */
503: numVertices[0] = numVertices1; numVertices[1] = numVertices2;
504: numEdges[0] = numEdges1; numEdges[1] = numEdges2;
505: DMNetworkSetSizes(networkdm,nsubnet,numVertices,numEdges,0,NULL);
507: edgelist[0] = edgelist1; edgelist[1] = edgelist2;
509: /* Add edge connectivity */
510: DMNetworkSetEdgeList(networkdm,edgelist,NULL);
512: /* Set up the network layout */
513: DMNetworkLayoutSetUp(networkdm);
515: /* Add network components only process 0 has any data to add*/
516: if (!crank) {
517: genj=0; loadj=0;
519: /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
520: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
522: for (i = 0; i < ne; i++) {
523: DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata1->branch[i]);
524: }
526: for (i = 0; i < nv; i++) {
527: DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata1->bus[i]);
528: if (pfdata1->bus[i].ngen) {
529: for (j = 0; j < pfdata1->bus[i].ngen; j++) {
530: DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata1->gen[genj++]);
531: }
532: }
533: if (pfdata1->bus[i].nload) {
534: for (j=0; j < pfdata1->bus[i].nload; j++) {
535: DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata1->load[loadj++]);
536: }
537: }
538: /* Add number of variables */
539: DMNetworkAddNumVariables(networkdm,vtx[i],2);
540: }
542: genj=0; loadj=0;
544: /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
545: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
547: for (i = 0; i < ne; i++) {
548: DMNetworkAddComponent(networkdm,edges[i],componentkey[0],&pfdata2->branch[i]);
549: }
551: for (i = 0; i < nv; i++) {
552: DMNetworkAddComponent(networkdm,vtx[i],componentkey[1],&pfdata2->bus[i]);
553: if (pfdata2->bus[i].ngen) {
554: for (j = 0; j < pfdata2->bus[i].ngen; j++) {
555: DMNetworkAddComponent(networkdm,vtx[i],componentkey[2],&pfdata2->gen[genj++]);
556: }
557: }
558: if (pfdata2->bus[i].nload) {
559: for (j=0; j < pfdata2->bus[i].nload; j++) {
560: DMNetworkAddComponent(networkdm,vtx[i],componentkey[3],&pfdata2->load[loadj++]);
561: }
562: }
563: /* Add number of variables */
564: DMNetworkAddNumVariables(networkdm,vtx[i],2);
565: }
566: }
568: /* Set up DM for use */
569: DMSetUp(networkdm);
571: if (!crank) {
572: PetscFree(edgelist1);
573: PetscFree(edgelist2);
574: }
576: if (!crank) {
577: PetscFree(pfdata1->bus);
578: PetscFree(pfdata1->gen);
579: PetscFree(pfdata1->branch);
580: PetscFree(pfdata1->load);
581: PetscFree(pfdata1);
583: PetscFree(pfdata2->bus);
584: PetscFree(pfdata2->gen);
585: PetscFree(pfdata2->branch);
586: PetscFree(pfdata2->load);
587: PetscFree(pfdata2);
588: }
590: /* Distribute networkdm to multiple processes */
591: DMNetworkDistribute(&networkdm,0);
593: PetscLogStagePop();
595: /* Broadcast Sbase to all processors */
596: MPI_Bcast(&User.Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
598: DMCreateGlobalVector(networkdm,&X);
599: VecDuplicate(X,&F);
601: DMCreateMatrix(networkdm,&J);
602: MatSetOption(J,MAT_NEW_NONZERO_ALLOCATION_ERR,PETSC_FALSE);
604: SetInitialValues(networkdm,X,&User);
606: /* HOOK UP SOLVER */
607: SNESCreate(PETSC_COMM_WORLD,&snes);
608: SNESSetDM(snes,networkdm);
609: SNESSetFunction(snes,F,FormFunction,&User);
610: SNESSetJacobian(snes,J,J,FormJacobian,&User);
611: SNESSetFromOptions(snes);
613: SNESSolve(snes,NULL,X);
614: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
616: VecDestroy(&X);
617: VecDestroy(&F);
618: MatDestroy(&J);
620: SNESDestroy(&snes);
621: DMDestroy(&networkdm);
622: }
623: PetscFinalize();
624: return ierr;
625: }
627: /*TEST
629: build:
630: depends: PFReadData.c pffunctions.c
631: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
634: test:
635: args: -snes_rtol 1.e-3
636: localrunfiles: poweroptions case9.m
637: output_file: output/power2_1.out
638: requires: double !complex
640: test:
641: suffix: 2
642: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
643: nsize: 4
644: localrunfiles: poweroptions case9.m
645: output_file: output/power2_1.out
646: requires: double !complex
648: TEST*/