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) farr[offset + 1] = xarr[offset + 1] - bus->vm;
116: }
117: VecRestoreArrayRead(localX, &xarr);
118: VecRestoreArray(localF, &farr);
119: return 0;
120: }
122: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *appctx)
123: {
124: DM networkdm;
125: Vec localX, localF;
126: PetscInt nv, ne;
127: const PetscInt *vtx, *edges;
129: SNESGetDM(snes, &networkdm);
130: DMGetLocalVector(networkdm, &localX);
131: DMGetLocalVector(networkdm, &localF);
132: VecSet(F, 0.0);
134: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
135: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
137: DMGlobalToLocalBegin(networkdm, F, INSERT_VALUES, localF);
138: DMGlobalToLocalEnd(networkdm, F, INSERT_VALUES, localF);
140: /* Form Function for first subnetwork */
141: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
142: FormFunction_Subnet(networkdm, localX, localF, nv, ne, vtx, edges, appctx);
144: /* Form Function for second subnetwork */
145: DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges);
146: FormFunction_Subnet(networkdm, localX, localF, nv, ne, vtx, edges, appctx);
148: DMRestoreLocalVector(networkdm, &localX);
150: DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F);
151: DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F);
152: DMRestoreLocalVector(networkdm, &localF);
153: return 0;
154: }
156: PetscErrorCode FormJacobian_Subnet(DM networkdm, Vec localX, Mat J, Mat Jpre, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
157: {
158: UserCtx_Power *User = (UserCtx_Power *)appctx;
159: PetscInt e, v, vfrom, vto;
160: const PetscScalar *xarr;
161: PetscInt offsetfrom, offsetto, goffsetfrom, goffsetto;
162: PetscInt row[2], col[8];
163: PetscScalar values[8];
165: VecGetArrayRead(localX, &xarr);
167: for (v = 0; v < nv; v++) {
168: PetscInt i, j, key;
169: PetscInt offset, goffset;
170: PetscScalar Vm;
171: PetscScalar Sbase = User->Sbase;
172: VERTEX_Power bus;
173: PetscBool ghostvtex;
174: PetscInt numComps;
175: void *component;
177: DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex);
178: DMNetworkGetNumComponents(networkdm, vtx[v], &numComps);
179: for (j = 0; j < numComps; j++) {
180: DMNetworkGetLocalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &offset);
181: DMNetworkGetGlobalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &goffset);
182: DMNetworkGetComponent(networkdm, vtx[v], j, &key, &component, NULL);
183: if (key == 1) {
184: PetscInt nconnedges;
185: const PetscInt *connedges;
187: bus = (VERTEX_Power)(component);
188: if (!ghostvtex) {
189: /* Handle reference bus constrained dofs */
190: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
191: row[0] = goffset;
192: row[1] = goffset + 1;
193: col[0] = goffset;
194: col[1] = goffset + 1;
195: col[2] = goffset;
196: col[3] = goffset + 1;
197: values[0] = 1.0;
198: values[1] = 0.0;
199: values[2] = 0.0;
200: 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;
209: row[1] = goffset + 1;
210: col[0] = goffset;
211: col[1] = goffset + 1;
212: values[0] = values[1] = values[2] = values[3] = 0.0;
213: if (bus->ide != PV_BUS) {
214: values[1] = 2.0 * Vm * bus->gl / Sbase;
215: values[3] = -2.0 * Vm * bus->bl / Sbase;
216: }
217: MatSetValues(J, 2, row, 2, col, values, ADD_VALUES);
218: }
220: DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges);
221: for (i = 0; i < nconnedges; i++) {
222: EDGE_Power branch;
223: VERTEX_Power busf, bust;
224: PetscInt keyf, keyt;
225: PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt;
226: const PetscInt *cone;
227: PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf;
229: e = connedges[i];
230: DMNetworkGetComponent(networkdm, e, 0, &key, (void **)&branch, NULL);
231: if (!branch->status) continue;
233: Gff = branch->yff[0];
234: Bff = branch->yff[1];
235: Gft = branch->yft[0];
236: Bft = branch->yft[1];
237: Gtf = branch->ytf[0];
238: Btf = branch->ytf[1];
239: Gtt = branch->ytt[0];
240: Btt = branch->ytt[1];
242: DMNetworkGetConnectedVertices(networkdm, e, &cone);
243: vfrom = cone[0];
244: vto = cone[1];
246: DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom);
247: DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto);
248: DMNetworkGetGlobalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &goffsetfrom);
249: DMNetworkGetGlobalVecOffset(networkdm, vto, ALL_COMPONENTS, &goffsetto);
251: if (goffsetto < 0) goffsetto = -goffsetto - 1;
253: thetaf = xarr[offsetfrom];
254: Vmf = xarr[offsetfrom + 1];
255: thetat = xarr[offsetto];
256: Vmt = xarr[offsetto + 1];
257: thetaft = thetaf - thetat;
258: thetatf = thetat - thetaf;
260: DMNetworkGetComponent(networkdm, vfrom, 0, &keyf, (void **)&busf, NULL);
261: DMNetworkGetComponent(networkdm, vto, 0, &keyt, (void **)&bust, NULL);
263: if (vfrom == vtx[v]) {
264: if (busf->ide != REF_BUS) {
265: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */
266: row[0] = goffsetfrom;
267: col[0] = goffsetfrom;
268: col[1] = goffsetfrom + 1;
269: col[2] = goffsetto;
270: col[3] = goffsetto + 1;
271: values[0] = Vmf * Vmt * (Gft * -PetscSinScalar(thetaft) + Bft * PetscCosScalar(thetaft)); /* df_dthetaf */
272: values[1] = 2.0 * Gff * Vmf + Vmt * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmf */
273: values[2] = Vmf * Vmt * (Gft * PetscSinScalar(thetaft) + Bft * -PetscCosScalar(thetaft)); /* df_dthetat */
274: values[3] = Vmf * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmt */
276: MatSetValues(J, 1, row, 4, col, values, ADD_VALUES);
277: }
278: if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
279: row[0] = goffsetfrom + 1;
280: col[0] = goffsetfrom;
281: col[1] = goffsetfrom + 1;
282: col[2] = goffsetto;
283: col[3] = goffsetto + 1;
284: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
285: values[0] = Vmf * Vmt * (Bft * PetscSinScalar(thetaft) + Gft * PetscCosScalar(thetaft));
286: values[1] = -2.0 * Bff * Vmf + Vmt * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft));
287: values[2] = Vmf * Vmt * (-Bft * PetscSinScalar(thetaft) + Gft * -PetscCosScalar(thetaft));
288: values[3] = Vmf * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft));
290: MatSetValues(J, 1, row, 4, col, values, ADD_VALUES);
291: }
292: } else {
293: if (bust->ide != REF_BUS) {
294: row[0] = goffsetto;
295: col[0] = goffsetto;
296: col[1] = goffsetto + 1;
297: col[2] = goffsetfrom;
298: col[3] = goffsetfrom + 1;
299: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
300: values[0] = Vmt * Vmf * (Gtf * -PetscSinScalar(thetatf) + Btf * PetscCosScalar(thetaft)); /* df_dthetat */
301: values[1] = 2.0 * Gtt * Vmt + Vmf * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmt */
302: values[2] = Vmt * Vmf * (Gtf * PetscSinScalar(thetatf) + Btf * -PetscCosScalar(thetatf)); /* df_dthetaf */
303: values[3] = Vmt * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmf */
305: MatSetValues(J, 1, row, 4, col, values, ADD_VALUES);
306: }
307: if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
308: row[0] = goffsetto + 1;
309: col[0] = goffsetto;
310: col[1] = goffsetto + 1;
311: col[2] = goffsetfrom;
312: col[3] = goffsetfrom + 1;
313: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
314: values[0] = Vmt * Vmf * (Btf * PetscSinScalar(thetatf) + Gtf * PetscCosScalar(thetatf));
315: values[1] = -2.0 * Btt * Vmt + Vmf * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf));
316: values[2] = Vmt * Vmf * (-Btf * PetscSinScalar(thetatf) + Gtf * -PetscCosScalar(thetatf));
317: values[3] = Vmt * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf));
319: MatSetValues(J, 1, row, 4, col, values, ADD_VALUES);
320: }
321: }
322: }
323: if (!ghostvtex && bus->ide == PV_BUS) {
324: row[0] = goffset + 1;
325: col[0] = goffset + 1;
326: values[0] = 1.0;
327: MatSetValues(J, 1, row, 1, col, values, ADD_VALUES);
328: }
329: }
330: }
331: }
332: VecRestoreArrayRead(localX, &xarr);
333: return 0;
334: }
336: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx)
337: {
338: DM networkdm;
339: Vec localX;
340: PetscInt ne, nv;
341: const PetscInt *vtx, *edges;
343: MatZeroEntries(J);
345: SNESGetDM(snes, &networkdm);
346: DMGetLocalVector(networkdm, &localX);
348: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
349: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
351: /* Form Jacobian for first subnetwork */
352: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
353: FormJacobian_Subnet(networkdm, localX, J, Jpre, nv, ne, vtx, edges, appctx);
355: /* Form Jacobian for second subnetwork */
356: DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges);
357: FormJacobian_Subnet(networkdm, localX, J, Jpre, nv, ne, vtx, edges, appctx);
359: DMRestoreLocalVector(networkdm, &localX);
361: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
362: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
363: return 0;
364: }
366: PetscErrorCode SetInitialValues_Subnet(DM networkdm, Vec localX, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
367: {
368: VERTEX_Power bus;
369: PetscInt i;
370: GEN gen;
371: PetscBool ghostvtex;
372: PetscScalar *xarr;
373: PetscInt key, numComps, j, offset;
374: void *component;
376: VecGetArray(localX, &xarr);
377: for (i = 0; i < nv; i++) {
378: DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex);
379: if (ghostvtex) continue;
381: DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset);
382: DMNetworkGetNumComponents(networkdm, vtx[i], &numComps);
383: for (j = 0; j < numComps; j++) {
384: DMNetworkGetComponent(networkdm, vtx[i], j, &key, &component, NULL);
385: if (key == 1) {
386: bus = (VERTEX_Power)(component);
387: xarr[offset] = bus->va * PETSC_PI / 180.0;
388: xarr[offset + 1] = bus->vm;
389: } else if (key == 2) {
390: gen = (GEN)(component);
391: if (!gen->status) continue;
392: xarr[offset + 1] = gen->vs;
393: break;
394: }
395: }
396: }
397: VecRestoreArray(localX, &xarr);
398: return 0;
399: }
401: PetscErrorCode SetInitialValues(DM networkdm, Vec X, void *appctx)
402: {
403: PetscInt nv, ne;
404: const PetscInt *vtx, *edges;
405: Vec localX;
407: DMGetLocalVector(networkdm, &localX);
409: VecSet(X, 0.0);
410: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
411: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
413: /* Set initial guess for first subnetwork */
414: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
415: SetInitialValues_Subnet(networkdm, localX, nv, ne, vtx, edges, appctx);
417: /* Set initial guess for second subnetwork */
418: DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges);
419: SetInitialValues_Subnet(networkdm, localX, nv, ne, vtx, edges, appctx);
421: DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X);
422: DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X);
423: DMRestoreLocalVector(networkdm, &localX);
424: return 0;
425: }
427: int main(int argc, char **argv)
428: {
429: char pfdata_file[PETSC_MAX_PATH_LEN] = "case9.m";
430: PFDATA *pfdata1, *pfdata2;
431: PetscInt numEdges1 = 0, numEdges2 = 0;
432: PetscInt *edgelist1 = NULL, *edgelist2 = NULL, componentkey[4];
433: DM networkdm;
434: UserCtx_Power User;
435: #if defined(PETSC_USE_LOG)
436: PetscLogStage stage1, stage2;
437: #endif
438: PetscMPIInt rank;
439: PetscInt nsubnet = 2, nv, ne, i, j, genj, loadj;
440: const PetscInt *vtx, *edges;
441: Vec X, F;
442: Mat J;
443: SNES snes;
446: PetscInitialize(&argc, &argv, "poweroptions", help);
447: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
448: {
449: /* 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 */
450: /* this is an experiment to see how the analyzer reacts */
451: const PetscMPIInt crank = rank;
453: /* Create an empty network object */
454: DMNetworkCreate(PETSC_COMM_WORLD, &networkdm);
456: /* Register the components in the network */
457: DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &componentkey[0]);
458: DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &componentkey[1]);
459: DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &componentkey[2]);
460: DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &componentkey[3]);
462: PetscLogStageRegister("Read Data", &stage1);
463: PetscLogStagePush(stage1);
464: /* READ THE DATA */
465: if (!crank) {
466: /* Only rank 0 reads the data */
467: PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, sizeof(pfdata_file), NULL);
468: /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */
470: /* READ DATA FOR THE FIRST SUBNETWORK */
471: PetscNew(&pfdata1);
472: PFReadMatPowerData(pfdata1, pfdata_file);
473: User.Sbase = pfdata1->sbase;
475: numEdges1 = pfdata1->nbranch;
476: PetscMalloc1(2 * numEdges1, &edgelist1);
477: GetListofEdges_Power(pfdata1, edgelist1);
479: /* READ DATA FOR THE SECOND SUBNETWORK */
480: PetscNew(&pfdata2);
481: PFReadMatPowerData(pfdata2, pfdata_file);
482: User.Sbase = pfdata2->sbase;
484: numEdges2 = pfdata2->nbranch;
485: PetscMalloc1(2 * numEdges2, &edgelist2);
486: GetListofEdges_Power(pfdata2, edgelist2);
487: }
489: PetscLogStagePop();
490: MPI_Barrier(PETSC_COMM_WORLD);
491: PetscLogStageRegister("Create network", &stage2);
492: PetscLogStagePush(stage2);
494: /* Set number of nodes/edges and edge connectivity */
495: DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, nsubnet);
496: DMNetworkAddSubnetwork(networkdm, "", numEdges1, edgelist1, NULL);
497: DMNetworkAddSubnetwork(networkdm, "", numEdges2, edgelist2, NULL);
499: /* Set up the network layout */
500: DMNetworkLayoutSetUp(networkdm);
502: /* Add network components only process 0 has any data to add*/
503: if (!crank) {
504: genj = 0;
505: loadj = 0;
507: /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
508: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
510: for (i = 0; i < ne; i++) DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata1->branch[i], 0);
512: for (i = 0; i < nv; i++) {
513: DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata1->bus[i], 2);
514: if (pfdata1->bus[i].ngen) {
515: for (j = 0; j < pfdata1->bus[i].ngen; j++) DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata1->gen[genj++], 0);
516: }
517: if (pfdata1->bus[i].nload) {
518: for (j = 0; j < pfdata1->bus[i].nload; j++) DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata1->load[loadj++], 0);
519: }
520: }
522: genj = 0;
523: loadj = 0;
525: /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
526: DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges);
528: for (i = 0; i < ne; i++) DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata2->branch[i], 0);
530: for (i = 0; i < nv; i++) {
531: DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata2->bus[i], 2);
532: if (pfdata2->bus[i].ngen) {
533: for (j = 0; j < pfdata2->bus[i].ngen; j++) DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata2->gen[genj++], 0);
534: }
535: if (pfdata2->bus[i].nload) {
536: for (j = 0; j < pfdata2->bus[i].nload; j++) DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata2->load[loadj++], 0);
537: }
538: }
539: }
541: /* Set up DM for use */
542: DMSetUp(networkdm);
544: if (!crank) {
545: PetscFree(edgelist1);
546: PetscFree(edgelist2);
547: }
549: if (!crank) {
550: PetscFree(pfdata1->bus);
551: PetscFree(pfdata1->gen);
552: PetscFree(pfdata1->branch);
553: PetscFree(pfdata1->load);
554: PetscFree(pfdata1);
556: PetscFree(pfdata2->bus);
557: PetscFree(pfdata2->gen);
558: PetscFree(pfdata2->branch);
559: PetscFree(pfdata2->load);
560: PetscFree(pfdata2);
561: }
563: /* Distribute networkdm to multiple processes */
564: DMNetworkDistribute(&networkdm, 0);
566: PetscLogStagePop();
568: /* Broadcast Sbase to all processors */
569: MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD);
571: DMCreateGlobalVector(networkdm, &X);
572: VecDuplicate(X, &F);
574: DMCreateMatrix(networkdm, &J);
575: MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE);
577: SetInitialValues(networkdm, X, &User);
579: /* HOOK UP SOLVER */
580: SNESCreate(PETSC_COMM_WORLD, &snes);
581: SNESSetDM(snes, networkdm);
582: SNESSetFunction(snes, F, FormFunction, &User);
583: SNESSetJacobian(snes, J, J, FormJacobian, &User);
584: SNESSetFromOptions(snes);
586: SNESSolve(snes, NULL, X);
587: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
589: VecDestroy(&X);
590: VecDestroy(&F);
591: MatDestroy(&J);
593: SNESDestroy(&snes);
594: DMDestroy(&networkdm);
595: }
596: PetscFinalize();
597: return 0;
598: }
600: /*TEST
602: build:
603: depends: PFReadData.c pffunctions.c
604: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
606: test:
607: args: -snes_rtol 1.e-3
608: localrunfiles: poweroptions case9.m
609: output_file: output/power_1.out
611: test:
612: suffix: 2
613: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
614: nsize: 4
615: localrunfiles: poweroptions case9.m
616: output_file: output/power_1.out
618: TEST*/