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";
7
8 #include "power.h"
9 #include <petscdmnetwork.h>
10
FormFunction_Subnet(DM networkdm,Vec localX,Vec localF,PetscInt nv,PetscInt ne,const PetscInt * vtx,const PetscInt * edges,void * appctx)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;
18
19 PetscFunctionBegin;
20 PetscCall(VecGetArrayRead(localX, &xarr));
21 PetscCall(VecGetArray(localF, &farr));
22
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;
33
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;
42
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 }
50
51 if (!ghostvtex) {
52 Vm = xarr[offset + 1];
53
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 }
58
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;
66
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];
78
79 PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
80 vfrom = cone[0];
81 vto = cone[1];
82
83 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vfrom, ALL_COMPONENTS, &offsetfrom));
84 PetscCall(DMNetworkGetLocalVecOffset(networkdm, vto, ALL_COMPONENTS, &offsetto));
85
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;
92
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 }
122
FormFunction(SNES snes,Vec X,Vec F,void * appctx)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;
129
130 PetscFunctionBegin;
131 PetscCall(SNESGetDM(snes, &networkdm));
132 PetscCall(DMGetLocalVector(networkdm, &localX));
133 PetscCall(DMGetLocalVector(networkdm, &localF));
134 PetscCall(VecSet(F, 0.0));
135
136 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
137 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
138
139 PetscCall(DMGlobalToLocalBegin(networkdm, F, INSERT_VALUES, localF));
140 PetscCall(DMGlobalToLocalEnd(networkdm, F, INSERT_VALUES, localF));
141
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));
145
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));
149
150 PetscCall(DMRestoreLocalVector(networkdm, &localX));
151
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 }
157
FormJacobian_Subnet(DM networkdm,Vec localX,Mat J,Mat Jpre,PetscInt nv,PetscInt ne,const PetscInt * vtx,const PetscInt * edges,void * appctx)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];
166
167 PetscFunctionBegin;
168 PetscCall(VecGetArrayRead(localX, &xarr));
169
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;
179
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;
189
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 }
207
208 Vm = xarr[offset + 1];
209
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 }
222
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;
231
232 e = connedges[i];
233 PetscCall(DMNetworkGetComponent(networkdm, e, 0, &key, (void **)&branch, NULL));
234 if (!branch->status) continue;
235
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];
244
245 PetscCall(DMNetworkGetConnectedVertices(networkdm, e, &cone));
246 vfrom = cone[0];
247 vto = cone[1];
248
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));
253
254 if (goffsetto < 0) goffsetto = -goffsetto - 1;
255
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;
262
263 PetscCall(DMNetworkGetComponent(networkdm, vfrom, 0, &keyf, (void **)&busf, NULL));
264 PetscCall(DMNetworkGetComponent(networkdm, vto, 0, &keyt, (void **)&bust, NULL));
265
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 */
278
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));
292
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 */
307
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));
321
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 }
338
FormJacobian(SNES snes,Vec X,Mat J,Mat Jpre,void * appctx)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;
345
346 PetscFunctionBegin;
347 PetscCall(MatZeroEntries(J));
348
349 PetscCall(SNESGetDM(snes, &networkdm));
350 PetscCall(DMGetLocalVector(networkdm, &localX));
351
352 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
353 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
354
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));
358
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));
362
363 PetscCall(DMRestoreLocalVector(networkdm, &localX));
364
365 PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
366 PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
367 PetscFunctionReturn(PETSC_SUCCESS);
368 }
369
SetInitialValues_Subnet(DM networkdm,Vec localX,PetscInt nv,PetscInt ne,const PetscInt * vtx,const PetscInt * edges,void * appctx)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;
379
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;
385
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 }
405
SetInitialValues(DM networkdm,Vec X,void * appctx)406 PetscErrorCode SetInitialValues(DM networkdm, Vec X, void *appctx)
407 {
408 PetscInt nv, ne;
409 const PetscInt *vtx, *edges;
410 Vec localX;
411
412 PetscFunctionBegin;
413 PetscCall(DMGetLocalVector(networkdm, &localX));
414
415 PetscCall(VecSet(X, 0.0));
416 PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
417 PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
418
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));
422
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));
426
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 }
432
main(int argc,char ** argv)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;
448
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;
456
457 /* Create an empty network object */
458 PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
459
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]));
465
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 */
473
474 /* READ DATA FOR THE FIRST SUBNETWORK */
475 PetscCall(PetscNew(&pfdata1));
476 PetscCall(PFReadMatPowerData(pfdata1, pfdata_file));
477 User.Sbase = pfdata1->sbase;
478
479 numEdges1 = pfdata1->nbranch;
480 PetscCall(PetscMalloc1(2 * numEdges1, &edgelist1));
481 PetscCall(GetListofEdges_Power(pfdata1, edgelist1));
482
483 /* READ DATA FOR THE SECOND SUBNETWORK */
484 PetscCall(PetscNew(&pfdata2));
485 PetscCall(PFReadMatPowerData(pfdata2, pfdata_file));
486 User.Sbase = pfdata2->sbase;
487
488 numEdges2 = pfdata2->nbranch;
489 PetscCall(PetscMalloc1(2 * numEdges2, &edgelist2));
490 PetscCall(GetListofEdges_Power(pfdata2, edgelist2));
491 }
492
493 PetscCall(PetscLogStagePop());
494 PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
495 PetscCall(PetscLogStageRegister("Create network", &stage2));
496 PetscCall(PetscLogStagePush(stage2));
497
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));
502
503 /* Set up the network layout */
504 PetscCall(DMNetworkLayoutSetUp(networkdm));
505
506 /* Add network components only process 0 has any data to add*/
507 if (!crank) {
508 genj = 0;
509 loadj = 0;
510
511 /* ADD VARIABLES AND COMPONENTS FOR THE FIRST SUBNETWORK */
512 PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
513
514 for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata1->branch[i], 0));
515
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 }
525
526 genj = 0;
527 loadj = 0;
528
529 /* ADD VARIABLES AND COMPONENTS FOR THE SECOND SUBNETWORK */
530 PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
531
532 for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], componentkey[0], &pfdata2->branch[i], 0));
533
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 }
544
545 /* Set up DM for use */
546 PetscCall(DMSetUp(networkdm));
547
548 if (!crank) {
549 PetscCall(PetscFree(edgelist1));
550 PetscCall(PetscFree(edgelist2));
551 }
552
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));
559
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 }
566
567 /* Distribute networkdm to multiple processes */
568 PetscCall(DMNetworkDistribute(&networkdm, 0));
569
570 PetscCall(PetscLogStagePop());
571
572 /* Broadcast Sbase to all processors */
573 PetscCallMPI(MPI_Bcast(&User.Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD));
574
575 PetscCall(DMCreateGlobalVector(networkdm, &X));
576 PetscCall(VecDuplicate(X, &F));
577
578 PetscCall(DMCreateMatrix(networkdm, &J));
579 PetscCall(MatSetOption(J, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_FALSE));
580
581 PetscCall(SetInitialValues(networkdm, X, &User));
582
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));
589
590 PetscCall(SNESSolve(snes, NULL, X));
591 /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
592
593 PetscCall(VecDestroy(&X));
594 PetscCall(VecDestroy(&F));
595 PetscCall(MatDestroy(&J));
596
597 PetscCall(SNESDestroy(&snes));
598 PetscCall(DMDestroy(&networkdm));
599 }
600 PetscCall(PetscFinalize());
601 return 0;
602 }
603
604 /*TEST
605
606 build:
607 depends: PFReadData.c pffunctions.c
608 requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
609
610 test:
611 args: -snes_rtol 1.e-3
612 localrunfiles: poweroptions case9.m
613 output_file: output/power_1.out
614
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
621
622 TEST*/
623