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: PetscFunctionBegin;
20: PetscCall(VecGetArrayRead(localX, &xarr));
21: PetscCall(VecGetArray(localF, &farr));
23: for (v = 0; v < nv; v++) {
24: PetscInt i, j, key;
25: PetscScalar Vm;
26: PetscScalar Sbase = User->Sbase;
27: VERTEX_Power bus = NULL;
28: GEN gen;
29: LOAD load;
30: PetscBool ghostvtex;
31: PetscInt numComps;
32: void *component;
34: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex));
35: PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &numComps));
36: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &offset));
37: for (j = 0; j < numComps; j++) {
38: PetscCall(DMNetworkGetComponent(networkdm, vtx[v], j, &key, &component, NULL));
39: if (key == 1) {
40: PetscInt nconnedges;
41: const PetscInt *connedges;
43: bus = (VERTEX_Power)(component);
44: /* Handle reference bus constrained dofs */
45: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
46: farr[offset] = xarr[offset] - bus->va * PETSC_PI / 180.0;
47: farr[offset + 1] = xarr[offset + 1] - bus->vm;
48: break;
49: }
51: if (!ghostvtex) {
52: Vm = xarr[offset + 1];
54: /* Shunt injections */
55: farr[offset] += Vm * Vm * bus->gl / Sbase;
56: if (bus->ide != PV_BUS) farr[offset + 1] += -Vm * Vm * bus->bl / Sbase;
57: }
59: PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges));
60: for (i = 0; i < nconnedges; i++) {
61: EDGE_Power branch;
62: PetscInt keye;
63: PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt;
64: const PetscInt *cone;
65: PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf;
67: e = connedges[i];
68: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, (void **)&branch, NULL));
69: if (!branch->status) continue;
70: Gff = branch->yff[0];
71: Bff = branch->yff[1];
72: Gft = branch->yft[0];
73: Bft = branch->yft[1];
74: Gtf = branch->ytf[0];
75: Btf = branch->ytf[1];
76: Gtt = branch->ytt[0];
77: Btt = branch->ytt[1];
79: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
80: vfrom = cone[0];
81: vto = cone[1];
83: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
84: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));
86: thetaf = xarr[offsetfrom];
87: Vmf = xarr[offsetfrom + 1];
88: thetat = xarr[offsetto];
89: Vmt = xarr[offsetto + 1];
90: thetaft = thetaf - thetat;
91: thetatf = thetat - thetaf;
93: if (vfrom == vtx[v]) {
94: farr[offsetfrom] += Gff * Vmf * Vmf + Vmf * Vmt * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft));
95: farr[offsetfrom + 1] += -Bff * Vmf * Vmf + Vmf * Vmt * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft));
96: } else {
97: farr[offsetto] += Gtt * Vmt * Vmt + Vmt * Vmf * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf));
98: farr[offsetto + 1] += -Btt * Vmt * Vmt + Vmt * Vmf * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf));
99: }
100: }
101: } else if (key == 2) {
102: if (!ghostvtex) {
103: gen = (GEN)(component);
104: if (!gen->status) continue;
105: farr[offset] += -gen->pg / Sbase;
106: farr[offset + 1] += -gen->qg / Sbase;
107: }
108: } else if (key == 3) {
109: if (!ghostvtex) {
110: load = (LOAD)(component);
111: farr[offset] += load->pl / Sbase;
112: farr[offset + 1] += load->ql / Sbase;
113: }
114: }
115: }
116: if (bus && bus->ide == PV_BUS) farr[offset + 1] = xarr[offset + 1] - bus->vm;
117: }
118: PetscCall(VecRestoreArrayRead(localX, &xarr));
119: PetscCall(VecRestoreArray(localF, &farr));
120: PetscFunctionReturn(PETSC_SUCCESS);
121: }
123: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *appctx)
124: {
125: DM networkdm;
126: Vec localX, localF;
127: PetscInt nv, ne;
128: const PetscInt *vtx, *edges;
130: PetscFunctionBegin;
131: PetscCall(SNESGetDM(snes, &networkdm));
132: PetscCall(DMGetLocalVector(networkdm, &localX));
133: PetscCall(DMGetLocalVector(networkdm, &localF));
134: PetscCall(VecSet(F, 0.0));
136: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
137: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
139: PetscCall(DMGlobalToLocalBegin(networkdm, F, INSERT_VALUES, localF));
140: PetscCall(DMGlobalToLocalEnd(networkdm, F, INSERT_VALUES, localF));
142: /* Form Function for first subnetwork */
143: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
144: PetscCall(FormFunction_Subnet(networkdm, localX, localF, nv, ne, vtx, edges, appctx));
146: /* Form Function for second subnetwork */
147: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
148: PetscCall(FormFunction_Subnet(networkdm, localX, localF, nv, ne, vtx, edges, appctx));
150: PetscCall(DMRestoreLocalVector(networkdm, &localX));
152: PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
153: PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
154: PetscCall(DMRestoreLocalVector(networkdm, &localF));
155: PetscFunctionReturn(PETSC_SUCCESS);
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: PetscFunctionBegin;
168: PetscCall(VecGetArrayRead(localX, &xarr));
170: for (v = 0; v < nv; v++) {
171: PetscInt i, j, key;
172: PetscInt offset, goffset;
173: PetscScalar Vm;
174: PetscScalar Sbase = User->Sbase;
175: VERTEX_Power bus;
176: PetscBool ghostvtex;
177: PetscInt numComps;
178: void *component;
180: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghostvtex));
181: PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &numComps));
182: for (j = 0; j < numComps; j++) {
183: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &offset));
184: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[v], ALL_COMPONENTS, &goffset));
185: PetscCall(DMNetworkGetComponent(networkdm, vtx[v], j, &key, &component, NULL));
186: if (key == 1) {
187: PetscInt nconnedges;
188: const PetscInt *connedges;
190: bus = (VERTEX_Power)(component);
191: if (!ghostvtex) {
192: /* Handle reference bus constrained dofs */
193: if (bus->ide == REF_BUS || bus->ide == ISOLATED_BUS) {
194: row[0] = goffset;
195: row[1] = goffset + 1;
196: col[0] = goffset;
197: col[1] = goffset + 1;
198: col[2] = goffset;
199: col[3] = goffset + 1;
200: values[0] = 1.0;
201: values[1] = 0.0;
202: values[2] = 0.0;
203: values[3] = 1.0;
204: PetscCall(MatSetValues(J, 2, row, 2, col, values, ADD_VALUES));
205: break;
206: }
208: Vm = xarr[offset + 1];
210: /* Shunt injections */
211: row[0] = goffset;
212: row[1] = goffset + 1;
213: col[0] = goffset;
214: col[1] = goffset + 1;
215: values[0] = values[1] = values[2] = values[3] = 0.0;
216: if (bus->ide != PV_BUS) {
217: values[1] = 2.0 * Vm * bus->gl / Sbase;
218: values[3] = -2.0 * Vm * bus->bl / Sbase;
219: }
220: PetscCall(MatSetValues(J, 2, row, 2, col, values, ADD_VALUES));
221: }
223: PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges));
224: for (i = 0; i < nconnedges; i++) {
225: EDGE_Power branch;
226: VERTEX_Power busf, bust;
227: PetscInt keyf, keyt;
228: PetscScalar Gff, Bff, Gft, Bft, Gtf, Btf, Gtt, Btt;
229: const PetscInt *cone;
230: PetscScalar Vmf, Vmt, thetaf, thetat, thetaft, thetatf;
232: e = connedges[i];
233: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &key, (void **)&branch, NULL));
234: if (!branch->status) continue;
236: Gff = branch->yff[0];
237: Bff = branch->yff[1];
238: Gft = branch->yft[0];
239: Bft = branch->yft[1];
240: Gtf = branch->ytf[0];
241: Btf = branch->ytf[1];
242: Gtt = branch->ytt[0];
243: Btt = branch->ytt[1];
245: PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
246: vfrom = cone[0];
247: vto = cone[1];
249: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
250: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));
251: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &goffsetfrom));
252: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vto, ALL_COMPONENTS, &goffsetto));
254: if (goffsetto < 0) goffsetto = -goffsetto - 1;
256: thetaf = xarr[offsetfrom];
257: Vmf = xarr[offsetfrom + 1];
258: thetat = xarr[offsetto];
259: Vmt = xarr[offsetto + 1];
260: thetaft = thetaf - thetat;
261: thetatf = thetat - thetaf;
263: PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &keyf, (void **)&busf, NULL));
264: PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &keyt, (void **)&bust, NULL));
266: if (vfrom == vtx[v]) {
267: if (busf->ide != REF_BUS) {
268: /* farr[offsetfrom] += Gff*Vmf*Vmf + Vmf*Vmt*(Gft*PetscCosScalar(thetaft) + Bft*PetscSinScalar(thetaft)); */
269: row[0] = goffsetfrom;
270: col[0] = goffsetfrom;
271: col[1] = goffsetfrom + 1;
272: col[2] = goffsetto;
273: col[3] = goffsetto + 1;
274: values[0] = Vmf * Vmt * (Gft * -PetscSinScalar(thetaft) + Bft * PetscCosScalar(thetaft)); /* df_dthetaf */
275: values[1] = 2.0 * Gff * Vmf + Vmt * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmf */
276: values[2] = Vmf * Vmt * (Gft * PetscSinScalar(thetaft) + Bft * -PetscCosScalar(thetaft)); /* df_dthetat */
277: values[3] = Vmf * (Gft * PetscCosScalar(thetaft) + Bft * PetscSinScalar(thetaft)); /* df_dVmt */
279: PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES));
280: }
281: if (busf->ide != PV_BUS && busf->ide != REF_BUS) {
282: row[0] = goffsetfrom + 1;
283: col[0] = goffsetfrom;
284: col[1] = goffsetfrom + 1;
285: col[2] = goffsetto;
286: col[3] = goffsetto + 1;
287: /* farr[offsetfrom+1] += -Bff*Vmf*Vmf + Vmf*Vmt*(-Bft*PetscCosScalar(thetaft) + Gft*PetscSinScalar(thetaft)); */
288: values[0] = Vmf * Vmt * (Bft * PetscSinScalar(thetaft) + Gft * PetscCosScalar(thetaft));
289: values[1] = -2.0 * Bff * Vmf + Vmt * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft));
290: values[2] = Vmf * Vmt * (-Bft * PetscSinScalar(thetaft) + Gft * -PetscCosScalar(thetaft));
291: values[3] = Vmf * (-Bft * PetscCosScalar(thetaft) + Gft * PetscSinScalar(thetaft));
293: PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES));
294: }
295: } else {
296: if (bust->ide != REF_BUS) {
297: row[0] = goffsetto;
298: col[0] = goffsetto;
299: col[1] = goffsetto + 1;
300: col[2] = goffsetfrom;
301: col[3] = goffsetfrom + 1;
302: /* farr[offsetto] += Gtt*Vmt*Vmt + Vmt*Vmf*(Gtf*PetscCosScalar(thetatf) + Btf*PetscSinScalar(thetatf)); */
303: values[0] = Vmt * Vmf * (Gtf * -PetscSinScalar(thetatf) + Btf * PetscCosScalar(thetaft)); /* df_dthetat */
304: values[1] = 2.0 * Gtt * Vmt + Vmf * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmt */
305: values[2] = Vmt * Vmf * (Gtf * PetscSinScalar(thetatf) + Btf * -PetscCosScalar(thetatf)); /* df_dthetaf */
306: values[3] = Vmt * (Gtf * PetscCosScalar(thetatf) + Btf * PetscSinScalar(thetatf)); /* df_dVmf */
308: PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES));
309: }
310: if (bust->ide != PV_BUS && bust->ide != REF_BUS) {
311: row[0] = goffsetto + 1;
312: col[0] = goffsetto;
313: col[1] = goffsetto + 1;
314: col[2] = goffsetfrom;
315: col[3] = goffsetfrom + 1;
316: /* farr[offsetto+1] += -Btt*Vmt*Vmt + Vmt*Vmf*(-Btf*PetscCosScalar(thetatf) + Gtf*PetscSinScalar(thetatf)); */
317: values[0] = Vmt * Vmf * (Btf * PetscSinScalar(thetatf) + Gtf * PetscCosScalar(thetatf));
318: values[1] = -2.0 * Btt * Vmt + Vmf * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf));
319: values[2] = Vmt * Vmf * (-Btf * PetscSinScalar(thetatf) + Gtf * -PetscCosScalar(thetatf));
320: values[3] = Vmt * (-Btf * PetscCosScalar(thetatf) + Gtf * PetscSinScalar(thetatf));
322: PetscCall(MatSetValues(J, 1, row, 4, col, values, ADD_VALUES));
323: }
324: }
325: }
326: if (!ghostvtex && bus->ide == PV_BUS) {
327: row[0] = goffset + 1;
328: col[0] = goffset + 1;
329: values[0] = 1.0;
330: PetscCall(MatSetValues(J, 1, row, 1, col, values, ADD_VALUES));
331: }
332: }
333: }
334: }
335: PetscCall(VecRestoreArrayRead(localX, &xarr));
336: PetscFunctionReturn(PETSC_SUCCESS);
337: }
339: PetscErrorCode FormJacobian(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx)
340: {
341: DM networkdm;
342: Vec localX;
343: PetscInt ne, nv;
344: const PetscInt *vtx, *edges;
346: PetscFunctionBegin;
347: PetscCall(MatZeroEntries(J));
349: PetscCall(SNESGetDM(snes, &networkdm));
350: PetscCall(DMGetLocalVector(networkdm, &localX));
352: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
353: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
355: /* Form Jacobian for first subnetwork */
356: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
357: PetscCall(FormJacobian_Subnet(networkdm, localX, J, Jpre, nv, ne, vtx, edges, appctx));
359: /* Form Jacobian for second subnetwork */
360: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
361: PetscCall(FormJacobian_Subnet(networkdm, localX, J, Jpre, nv, ne, vtx, edges, appctx));
363: PetscCall(DMRestoreLocalVector(networkdm, &localX));
365: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
366: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
367: PetscFunctionReturn(PETSC_SUCCESS);
368: }
370: PetscErrorCode SetInitialValues_Subnet(DM networkdm, Vec localX, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
371: {
372: VERTEX_Power bus;
373: PetscInt i;
374: GEN gen;
375: PetscBool ghostvtex;
376: PetscScalar *xarr;
377: PetscInt key, numComps, j, offset;
378: void *component;
380: PetscFunctionBegin;
381: PetscCall(VecGetArray(localX, &xarr));
382: for (i = 0; i < nv; i++) {
383: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex));
384: if (ghostvtex) continue;
386: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset));
387: PetscCall(DMNetworkGetNumComponents(networkdm, vtx[i], &numComps));
388: for (j = 0; j < numComps; j++) {
389: PetscCall(DMNetworkGetComponent(networkdm, vtx[i], j, &key, &component, NULL));
390: if (key == 1) {
391: bus = (VERTEX_Power)(component);
392: xarr[offset] = bus->va * PETSC_PI / 180.0;
393: xarr[offset + 1] = bus->vm;
394: } else if (key == 2) {
395: gen = (GEN)(component);
396: if (!gen->status) continue;
397: xarr[offset + 1] = gen->vs;
398: break;
399: }
400: }
401: }
402: PetscCall(VecRestoreArray(localX, &xarr));
403: PetscFunctionReturn(PETSC_SUCCESS);
404: }
406: PetscErrorCode SetInitialValues(DM networkdm, Vec X, void *appctx)
407: {
408: PetscInt nv, ne;
409: const PetscInt *vtx, *edges;
410: Vec localX;
412: PetscFunctionBegin;
413: PetscCall(DMGetLocalVector(networkdm, &localX));
415: PetscCall(VecSet(X, 0.0));
416: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
417: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
419: /* Set initial guess for first subnetwork */
420: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
421: PetscCall(SetInitialValues_Subnet(networkdm, localX, nv, ne, vtx, edges, appctx));
423: /* Set initial guess for second subnetwork */
424: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
425: PetscCall(SetInitialValues_Subnet(networkdm, localX, nv, ne, vtx, edges, appctx));
427: PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
428: PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
429: PetscCall(DMRestoreLocalVector(networkdm, &localX));
430: PetscFunctionReturn(PETSC_SUCCESS);
431: }
433: int main(int argc, char **argv)
434: {
435: char pfdata_file[PETSC_MAX_PATH_LEN] = "case9.m";
436: PFDATA *pfdata1, *pfdata2;
437: PetscInt numEdges1 = 0, numEdges2 = 0;
438: PetscInt *edgelist1 = NULL, *edgelist2 = NULL, componentkey[4];
439: DM networkdm;
440: UserCtx_Power User;
441: PetscLogStage stage1, stage2;
442: PetscMPIInt rank;
443: PetscInt nsubnet = 2, nv, ne, i, j, genj, loadj;
444: const PetscInt *vtx, *edges;
445: Vec X, F;
446: Mat J;
447: SNES snes;
449: PetscFunctionBeginUser;
450: PetscCall(PetscInitialize(&argc, &argv, "poweroptions", help));
451: PetscCallMPI(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: PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
460: /* Register the components in the network */
461: PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &componentkey[0]));
462: PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &componentkey[1]));
463: PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &componentkey[2]));
464: PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &componentkey[3]));
466: PetscCall(PetscLogStageRegister("Read Data", &stage1));
467: PetscCall(PetscLogStagePush(stage1));
468: /* READ THE DATA */
469: if (!crank) {
470: /* Only rank 0 reads the data */
471: PetscCall(PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, sizeof(pfdata_file), NULL));
472: /* HERE WE CREATE COPIES OF THE SAME NETWORK THAT WILL BE TREATED AS SUBNETWORKS */
474: /* READ DATA FOR THE FIRST SUBNETWORK */
475: PetscCall(PetscNew(&pfdata1));
476: PetscCall(PFReadMatPowerData(pfdata1, pfdata_file));
477: User.Sbase = pfdata1->sbase;
479: numEdges1 = pfdata1->nbranch;
480: PetscCall(PetscMalloc1(2 * numEdges1, &edgelist1));
481: PetscCall(GetListofEdges_Power(pfdata1, edgelist1));
483: /* READ DATA FOR THE SECOND SUBNETWORK */
484: PetscCall(PetscNew(&pfdata2));
485: PetscCall(PFReadMatPowerData(pfdata2, pfdata_file));
486: User.Sbase = pfdata2->sbase;
488: numEdges2 = pfdata2->nbranch;
489: PetscCall(PetscMalloc1(2 * numEdges2, &edgelist2));
490: PetscCall(GetListofEdges_Power(pfdata2, edgelist2));
491: }
493: PetscCall(PetscLogStagePop());
494: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
495: PetscCall(PetscLogStageRegister("Create network", &stage2));
496: PetscCall(PetscLogStagePush(stage2));
498: /* Set number of nodes/edges and edge connectivity */
499: PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, nsubnet));
500: PetscCall(DMNetworkAddSubnetwork(networkdm, "", numEdges1, edgelist1, NULL));
501: PetscCall(DMNetworkAddSubnetwork(networkdm, "", numEdges2, edgelist2, NULL));
503: /* Set up the network layout */
504: PetscCall(DMNetworkLayoutSetUp(networkdm));
506: /* Add network components only process 0 has any data to add*/
507: if (!crank) {
508: genj = 0;
509: loadj = 0;
511: /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
512: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
514: for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata1->branch[i], 0));
516: for (i = 0; i < nv; i++) {
517: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata1->bus[i], 2));
518: if (pfdata1->bus[i].ngen) {
519: for (j = 0; j < pfdata1->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata1->gen[genj++], 0));
520: }
521: if (pfdata1->bus[i].nload) {
522: for (j = 0; j < pfdata1->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata1->load[loadj++], 0));
523: }
524: }
526: genj = 0;
527: loadj = 0;
529: /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
530: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
532: for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata2->branch[i], 0));
534: for (i = 0; i < nv; i++) {
535: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[1], &pfdata2->bus[i], 2));
536: if (pfdata2->bus[i].ngen) {
537: for (j = 0; j < pfdata2->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[2], &pfdata2->gen[genj++], 0));
538: }
539: if (pfdata2->bus[i].nload) {
540: for (j = 0; j < pfdata2->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], componentkey[3], &pfdata2->load[loadj++], 0));
541: }
542: }
543: }
545: /* Set up DM for use */
546: PetscCall(DMSetUp(networkdm));
548: if (!crank) {
549: PetscCall(PetscFree(edgelist1));
550: PetscCall(PetscFree(edgelist2));
551: }
553: if (!crank) {
554: PetscCall(PetscFree(pfdata1->bus));
555: PetscCall(PetscFree(pfdata1->gen));
556: PetscCall(PetscFree(pfdata1->branch));
557: PetscCall(PetscFree(pfdata1->load));
558: PetscCall(PetscFree(pfdata1));
560: PetscCall(PetscFree(pfdata2->bus));
561: PetscCall(PetscFree(pfdata2->gen));
562: PetscCall(PetscFree(pfdata2->branch));
563: PetscCall(PetscFree(pfdata2->load));
564: PetscCall(PetscFree(pfdata2));
565: }
567: /* Distribute networkdm to multiple processes */
568: PetscCall(DMNetworkDistribute(&networkdm, 0));
570: PetscCall(PetscLogStagePop());
572: /* Broadcast Sbase to all processors */
573: PetscCallMPI(MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD));
575: PetscCall(DMCreateGlobalVector(networkdm, &X));
576: PetscCall(VecDuplicate(X, &F));
578: PetscCall(DMCreateMatrix(networkdm, &J));
579: PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
581: PetscCall(SetInitialValues(networkdm, X, &User));
583: /* HOOK UP SOLVER */
584: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
585: PetscCall(SNESSetDM(snes, networkdm));
586: PetscCall(SNESSetFunction(snes, F, FormFunction, &User));
587: PetscCall(SNESSetJacobian(snes, J, J, FormJacobian, &User));
588: PetscCall(SNESSetFromOptions(snes));
590: PetscCall(SNESSolve(snes, NULL, X));
591: /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
593: PetscCall(VecDestroy(&X));
594: PetscCall(VecDestroy(&F));
595: PetscCall(MatDestroy(&J));
597: PetscCall(SNESDestroy(&snes));
598: PetscCall(DMDestroy(&networkdm));
599: }
600: PetscCall(PetscFinalize());
601: return 0;
602: }
604: /*TEST
606: build:
607: depends: PFReadData.c pffunctions.c
608: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
610: test:
611: args: -snes_rtol 1.e-3
612: localrunfiles: poweroptions case9.m
613: output_file: output/power_1.out
615: test:
616: suffix: 2
617: args: -snes_rtol 1.e-3 -petscpartitioner_type simple
618: nsize: 4
619: localrunfiles: poweroptions case9.m
620: output_file: output/power_1.out
622: TEST*/