Actual source code: ex1.c
1: static char help[] = "This example demonstrates the use of DMNetwork with subnetworks for solving a coupled nonlinear \n\
2: electric power grid and water pipe problem.\n\
3: The available solver options are in the ex1options file \n\
4: and the data files are in the datafiles of subdirectories.\n\
5: Run this program: mpiexec -n <n> ./ex1 \n\\n";
6: /*
7: Example:
8: mpiexec -n 3 ./ex1 -petscpartitioner_type parmetis -dmnetwork_view draw -dmnetwork_view_distributed draw -dmnetwork_view_rank_range 0,1,2
9: mpiexec -n 3 ./ex1 -petscpartitioner_type simple -dmnetwork_view_distributed draw -dmnetwork_view_zoomin_vertices 0 -dmnetwork_view_zoomin_vertices_padding 2 -dmnetwork_view_rank_range 0
10: mpiexec -n <n> ./ex1 -monitorIteration -monitorColor -power_snes_max_it 0 -water_snes_max_it 0 -coupled_snes_max_it 10 -draw_pause 5.0
11: */
13: #include "power/power.h"
14: #include "water/water.h"
16: typedef struct {
17: UserCtx_Power appctx_power;
18: AppCtx_Water appctx_water;
19: PetscInt subsnes_id; /* snes solver id */
20: PetscInt it; /* iteration number */
21: Vec localXold; /* store previous solution, used by FormFunction_Dummy() */
22: PetscBool monitorColor;
23: } UserCtx;
25: /*
26: UserMonitor -- called at the end of every SNES iteration via option `-monitorIteration' or `-monitorColor'
27: */
28: PetscErrorCode UserMonitor(SNES snes, PetscInt its, PetscReal fnorm, void *appctx)
29: {
30: UserCtx *user = (UserCtx *)appctx;
31: PetscMPIInt rank;
32: MPI_Comm comm;
33: PetscInt it;
35: PetscFunctionBegin;
36: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
37: PetscCallMPI(MPI_Comm_rank(comm, &rank));
39: PetscCall(SNESGetIterationNumber(snes, &it));
40: if (rank == 0) {
41: PetscCall(SNESGetIterationNumber(snes, &it));
42: if (user->subsnes_id == 0 || user->subsnes_id == 1) {
43: PetscCall(PetscPrintf(PETSC_COMM_SELF, " subsnes_id %" PetscInt_FMT ", it %" PetscInt_FMT ", fnorm %g\n", user->subsnes_id, it, (double)fnorm));
44: } else {
45: PetscCall(PetscPrintf(PETSC_COMM_SELF, " coupled_snes_it %" PetscInt_FMT ", total_snes_it %" PetscInt_FMT ", fnorm %g\n", it, user->it, (double)fnorm));
46: }
47: }
49: if (user->monitorColor) {
50: DM networkdm, dmcoords;
51: Vec F;
52: PetscInt v, vStart, vEnd, offset, gidx, rstart;
53: PetscReal *color;
54: PetscScalar *farr;
55: PetscBool ghost;
57: PetscCall(SNESGetDM(snes, &networkdm));
58: PetscCall(DMGetCoordinateDM(networkdm, &dmcoords));
60: PetscCall(SNESGetFunction(snes, &F, NULL, NULL));
61: PetscCall(VecGetOwnershipRange(F, &rstart, NULL));
63: PetscCall(VecGetArray(F, &farr));
64: PetscCall(DMNetworkGetVertexRange(dmcoords, &vStart, &vEnd));
66: PetscCall(PetscPrintf(MPI_COMM_WORLD, "\nColorPrint:\n"));
67: for (v = vStart; v < vEnd; v++) {
68: PetscCall(DMNetworkIsGhostVertex(networkdm, v, &ghost));
69: PetscCall(DMNetworkGetComponent(dmcoords, v, 0, NULL, (void **)&color, NULL));
70: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, v, &gidx));
71: if (!ghost) {
72: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, v, 0, &offset));
73: *color = (PetscRealPart(farr[offset - rstart]));
74: }
75: PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%d] v %" PetscInt_FMT ": color[%" PetscInt_FMT "] = %g\n", rank, gidx, offset - rstart, *color));
76: }
77: PetscCall(PetscSynchronizedFlush(MPI_COMM_WORLD, NULL));
78: PetscCall(VecRestoreArray(F, &farr));
80: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_CSV));
81: PetscCall(DMView(networkdm, PETSC_VIEWER_DRAW_WORLD));
82: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
83: }
84: PetscFunctionReturn(PETSC_SUCCESS);
85: }
87: PetscErrorCode FormJacobian_subPower(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx)
88: {
89: DM networkdm;
90: Vec localX;
91: PetscInt nv, ne, i, j, offset, nvar, row;
92: const PetscInt *vtx, *edges;
93: PetscBool ghostvtex;
94: PetscScalar one = 1.0;
95: PetscMPIInt rank;
96: MPI_Comm comm;
98: PetscFunctionBegin;
99: PetscCall(SNESGetDM(snes, &networkdm));
100: PetscCall(DMGetLocalVector(networkdm, &localX));
102: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
103: PetscCallMPI(MPI_Comm_rank(comm, &rank));
105: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
106: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
108: PetscCall(MatZeroEntries(J));
110: /* Power subnetwork: copied from power/FormJacobian_Power() */
111: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
112: PetscCall(FormJacobian_Power_private(networkdm, localX, J, nv, ne, vtx, edges, appctx));
114: /* Water subnetwork: Identity */
115: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
116: for (i = 0; i < nv; i++) {
117: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex));
118: if (ghostvtex) continue;
120: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset));
121: PetscCall(DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar));
122: for (j = 0; j < nvar; j++) {
123: row = offset + j;
124: PetscCall(MatSetValues(J, 1, &row, 1, &row, &one, ADD_VALUES));
125: }
126: }
127: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
128: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
130: PetscCall(DMRestoreLocalVector(networkdm, &localX));
131: PetscFunctionReturn(PETSC_SUCCESS);
132: }
134: /* Dummy equation localF(X) = localX - localXold */
135: PetscErrorCode FormFunction_Dummy(DM networkdm, Vec localX, Vec localF, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
136: {
137: const PetscScalar *xarr, *xoldarr;
138: PetscScalar *farr;
139: PetscInt i, j, offset, nvar;
140: PetscBool ghostvtex;
141: UserCtx *user = (UserCtx *)appctx;
142: Vec localXold = user->localXold;
144: PetscFunctionBegin;
145: PetscCall(VecGetArrayRead(localX, &xarr));
146: PetscCall(VecGetArrayRead(localXold, &xoldarr));
147: PetscCall(VecGetArray(localF, &farr));
149: for (i = 0; i < nv; i++) {
150: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex));
151: if (ghostvtex) continue;
153: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset));
154: PetscCall(DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar));
155: for (j = 0; j < nvar; j++) farr[offset + j] = xarr[offset + j] - xoldarr[offset + j];
156: }
158: PetscCall(VecRestoreArrayRead(localX, &xarr));
159: PetscCall(VecRestoreArrayRead(localXold, &xoldarr));
160: PetscCall(VecRestoreArray(localF, &farr));
161: PetscFunctionReturn(PETSC_SUCCESS);
162: }
164: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *appctx)
165: {
166: DM networkdm;
167: Vec localX, localF;
168: PetscInt nv, ne, v;
169: const PetscInt *vtx, *edges;
170: PetscMPIInt rank;
171: MPI_Comm comm;
172: UserCtx *user = (UserCtx *)appctx;
173: UserCtx_Power appctx_power = (*user).appctx_power;
174: AppCtx_Water appctx_water = (*user).appctx_water;
176: PetscFunctionBegin;
177: PetscCall(SNESGetDM(snes, &networkdm));
178: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
179: PetscCallMPI(MPI_Comm_rank(comm, &rank));
181: PetscCall(DMGetLocalVector(networkdm, &localX));
182: PetscCall(DMGetLocalVector(networkdm, &localF));
183: PetscCall(VecSet(F, 0.0));
184: PetscCall(VecSet(localF, 0.0));
186: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
187: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
189: /* Form Function for power subnetwork */
190: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
191: if (user->subsnes_id == 1) { /* snes_water only */
192: PetscCall(FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user));
193: } else {
194: PetscCall(FormFunction_Power(networkdm, localX, localF, nv, ne, vtx, edges, &appctx_power));
195: }
197: /* Form Function for water subnetwork */
198: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
199: if (user->subsnes_id == 0) { /* snes_power only */
200: PetscCall(FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user));
201: } else {
202: PetscCall(FormFunction_Water(networkdm, localX, localF, nv, ne, vtx, edges, NULL));
203: }
205: /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */
206: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
207: for (v = 0; v < nv; v++) {
208: PetscInt key, ncomp, nvar, nconnedges, k, e, keye, goffset[3];
209: void *component;
210: const PetscInt *connedges;
212: PetscCall(DMNetworkGetComponent(networkdm, vtx[v], ALL_COMPONENTS, NULL, NULL, &nvar));
213: PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &ncomp));
214: /* printf(" [%d] coupling vertex[%" PetscInt_FMT "]: v %" PetscInt_FMT ", ncomp %" PetscInt_FMT "; nvar %" PetscInt_FMT "\n",rank,v,vtx[v], ncomp,nvar); */
216: for (k = 0; k < ncomp; k++) {
217: PetscCall(DMNetworkGetComponent(networkdm, vtx[v], k, &key, &component, &nvar));
218: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[v], k, &goffset[k]));
220: /* Verify the coupling vertex is a powernet load vertex or a water vertex */
221: switch (k) {
222: case 0:
223: PetscCheck(key == appctx_power.compkey_bus && nvar == 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "key %" PetscInt_FMT " not a power bus vertex or nvar %" PetscInt_FMT " != 2", key, nvar);
224: break;
225: case 1:
226: PetscCheck(key == appctx_power.compkey_load && nvar == 0 && goffset[1] == goffset[0] + 2, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a power load vertex");
227: break;
228: case 2:
229: PetscCheck(key == appctx_water.compkey_vtx && nvar == 1 && goffset[2] == goffset[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a water vertex");
230: break;
231: default:
232: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %" PetscInt_FMT " is wrong", k);
233: }
234: /* printf(" [%d] coupling vertex[%" PetscInt_FMT "]: key %" PetscInt_FMT "; nvar %" PetscInt_FMT ", goffset %" PetscInt_FMT "\n",rank,v,key,nvar,goffset[k]); */
235: }
237: /* Get its supporting edges */
238: PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges));
239: /* printf("\n[%d] coupling vertex: nconnedges %" PetscInt_FMT "\n",rank,nconnedges); */
240: for (k = 0; k < nconnedges; k++) {
241: e = connedges[k];
242: PetscCall(DMNetworkGetNumComponents(networkdm, e, &ncomp));
243: /* printf("\n [%d] connected edge[%" PetscInt_FMT "]=%" PetscInt_FMT " has ncomp %" PetscInt_FMT "\n",rank,k,e,ncomp); */
244: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, &component, NULL));
245: if (keye != appctx_water.compkey_edge) PetscCheck(keye == appctx_power.compkey_branch, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a power branch");
246: }
247: }
249: PetscCall(DMRestoreLocalVector(networkdm, &localX));
251: PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
252: PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
253: PetscCall(DMRestoreLocalVector(networkdm, &localF));
254: #if 0
255: if (rank == 0) printf("F:\n");
256: PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD));
257: #endif
258: PetscFunctionReturn(PETSC_SUCCESS);
259: }
261: PetscErrorCode SetInitialGuess(DM networkdm, Vec X, void *appctx)
262: {
263: PetscInt nv, ne, i, j, ncomp, offset, key;
264: const PetscInt *vtx, *edges;
265: UserCtx *user = (UserCtx *)appctx;
266: Vec localX = user->localXold;
267: UserCtx_Power appctx_power = (*user).appctx_power;
268: AppCtx_Water appctx_water = (*user).appctx_water;
269: PetscBool ghost;
270: PetscScalar *xarr;
271: VERTEX_Power bus;
272: VERTEX_Water vertex;
273: void *component;
274: GEN gen;
276: PetscFunctionBegin;
277: PetscCall(VecSet(X, 0.0));
278: PetscCall(VecSet(localX, 0.0));
280: /* Set initial guess for power subnetwork */
281: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
282: PetscCall(SetInitialGuess_Power(networkdm, localX, nv, ne, vtx, edges, &appctx_power));
284: /* Set initial guess for water subnetwork */
285: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
286: PetscCall(SetInitialGuess_Water(networkdm, localX, nv, ne, vtx, edges, NULL));
288: /* Set initial guess at the coupling vertex */
289: PetscCall(VecGetArray(localX, &xarr));
290: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
291: for (i = 0; i < nv; i++) {
292: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghost));
293: if (ghost) continue;
295: PetscCall(DMNetworkGetNumComponents(networkdm, vtx[i], &ncomp));
296: for (j = 0; j < ncomp; j++) {
297: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], j, &offset));
298: PetscCall(DMNetworkGetComponent(networkdm, vtx[i], j, &key, (void **)&component, NULL));
299: if (key == appctx_power.compkey_bus) {
300: bus = (VERTEX_Power)(component);
301: xarr[offset] = bus->va * PETSC_PI / 180.0;
302: xarr[offset + 1] = bus->vm;
303: } else if (key == appctx_power.compkey_gen) {
304: gen = (GEN)(component);
305: if (!gen->status) continue;
306: xarr[offset + 1] = gen->vs;
307: } else if (key == appctx_water.compkey_vtx) {
308: vertex = (VERTEX_Water)(component);
309: if (vertex->type == VERTEX_TYPE_JUNCTION) {
310: xarr[offset] = 100;
311: } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
312: xarr[offset] = vertex->res.head;
313: } else {
314: xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
315: }
316: }
317: }
318: }
319: PetscCall(VecRestoreArray(localX, &xarr));
321: PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
322: PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
323: PetscFunctionReturn(PETSC_SUCCESS);
324: }
326: /* Set coordinates */
327: static PetscErrorCode CoordinateVecSetUp(DM dmcoords, Vec coords)
328: {
329: PetscInt i, gidx, offset, v, nv, Nsubnet;
330: const PetscInt *vtx;
331: PetscScalar *carray;
332: PetscReal *color;
334: PetscFunctionBeginUser;
335: PetscCall(VecGetArrayWrite(coords, &carray));
336: PetscCall(DMNetworkGetNumSubNetworks(dmcoords, NULL, &Nsubnet));
337: for (i = 0; i < Nsubnet; i++) {
338: PetscCall(DMNetworkGetSubnetwork(dmcoords, i, &nv, NULL, &vtx, NULL));
339: for (v = 0; v < nv; v++) {
340: PetscCall(DMNetworkGetGlobalVertexIndex(dmcoords, vtx[v], &gidx));
341: PetscCall(DMNetworkGetLocalVecOffset(dmcoords, vtx[v], 0, &offset));
342: PetscCall(DMNetworkGetComponent(dmcoords, vtx[v], 0, NULL, (void **)&color, NULL));
343: *color = 0.0;
344: switch (gidx) {
345: case 0:
346: carray[offset] = -1.0;
347: carray[offset + 1] = -1.0;
348: break;
349: case 1:
350: carray[offset] = -2.0;
351: carray[offset + 1] = 2.0;
352: break;
353: case 2:
354: carray[offset] = 0.0;
355: carray[offset + 1] = 2.0;
356: break;
357: case 3:
358: carray[offset] = -1.0;
359: carray[offset + 1] = 0.0;
360: break;
361: case 4:
362: carray[offset] = 0.0;
363: carray[offset + 1] = 0.0;
364: break;
365: case 5:
366: carray[offset] = 0.0;
367: carray[offset + 1] = 1.0;
368: break;
369: case 6:
370: carray[offset] = -1.0;
371: carray[offset + 1] = 1.0;
372: break;
373: case 7:
374: carray[offset] = -2.0;
375: carray[offset + 1] = 1.0;
376: break;
377: case 8:
378: carray[offset] = -2.0;
379: carray[offset + 1] = 0.0;
380: break;
381: case 9:
382: carray[offset] = 1.0;
383: carray[offset + 1] = 0.0;
384: break;
385: case 10:
386: carray[offset] = 1.0;
387: carray[offset + 1] = -1.0;
388: break;
389: case 11:
390: carray[offset] = 2.0;
391: carray[offset + 1] = -1.0;
392: break;
393: case 12:
394: carray[offset] = 2.0;
395: carray[offset + 1] = 0.0;
396: break;
397: case 13:
398: carray[offset] = 0.0;
399: carray[offset + 1] = -1.0;
400: break;
401: case 14:
402: carray[offset] = 2.0;
403: carray[offset + 1] = 1.0;
404: break;
405: default:
406: PetscCheck(gidx < 15 && gidx > -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "gidx %" PetscInt_FMT "must between 0 and 14", gidx);
407: }
408: }
409: }
410: PetscCall(VecRestoreArrayWrite(coords, &carray));
411: PetscFunctionReturn(PETSC_SUCCESS);
412: }
414: static PetscErrorCode CoordinatePrint(DM dm)
415: {
416: DM dmcoords;
417: PetscInt cdim, v, off, vglobal, vStart, vEnd;
418: const PetscScalar *carray;
419: Vec coords;
420: MPI_Comm comm;
421: PetscMPIInt rank;
423: PetscFunctionBegin;
424: /* get info from dm */
425: PetscCall(DMGetCoordinateDim(dm, &cdim));
426: PetscCall(DMGetCoordinatesLocal(dm, &coords));
428: PetscCall(DMGetCoordinateDM(dm, &dmcoords));
429: PetscCall(PetscObjectGetComm((PetscObject)dmcoords, &comm));
430: PetscCallMPI(MPI_Comm_rank(comm, &rank));
432: /* print coordinates from dmcoords */
433: PetscCall(PetscPrintf(MPI_COMM_WORLD, "\nCoordinatePrint, cdim %" PetscInt_FMT ":\n", cdim));
434: PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%d]\n", rank));
436: PetscCall(DMNetworkGetVertexRange(dmcoords, &vStart, &vEnd));
437: PetscCall(VecGetArrayRead(coords, &carray));
438: for (v = vStart; v < vEnd; v++) {
439: PetscCall(DMNetworkGetLocalVecOffset(dmcoords, v, 0, &off));
440: PetscCall(DMNetworkGetGlobalVertexIndex(dmcoords, v, &vglobal));
441: switch (cdim) {
442: case 2:
443: PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "Vertex: %" PetscInt_FMT ", x = %f y = %f \n", vglobal, (double)PetscRealPart(carray[off]), (double)PetscRealPart(carray[off + 1])));
444: break;
445: default:
446: PetscCheck(cdim == 2, MPI_COMM_WORLD, PETSC_ERR_SUP, "Only supports Network embedding dimension of 2, not supplied %" PetscInt_FMT, cdim);
447: break;
448: }
449: }
450: PetscCall(PetscSynchronizedFlush(MPI_COMM_WORLD, NULL));
451: PetscCall(VecRestoreArrayRead(coords, &carray));
452: PetscFunctionReturn(PETSC_SUCCESS);
453: }
455: int main(int argc, char **argv)
456: {
457: DM networkdm, dmcoords;
458: PetscMPIInt rank, size;
459: PetscInt Nsubnet = 2, numVertices[2], numEdges[2], i, j, nv, ne, it_max = 10;
460: PetscInt vStart, vEnd, compkey;
461: const PetscInt *vtx, *edges;
462: PetscReal *color;
463: Vec X, F, coords;
464: SNES snes, snes_power, snes_water;
465: Mat Jac;
466: PetscBool ghost, viewJ = PETSC_FALSE, viewX = PETSC_FALSE, test = PETSC_FALSE, distribute = PETSC_TRUE, flg, printCoord = PETSC_FALSE, viewCSV = PETSC_FALSE, monitorIt = PETSC_FALSE;
467: UserCtx user;
468: SNESConvergedReason reason;
469: PetscLogStage stage[4];
471: /* Power subnetwork */
472: UserCtx_Power *appctx_power = &user.appctx_power;
473: char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
474: PFDATA *pfdata = NULL;
475: PetscInt genj, loadj, *edgelist_power = NULL, power_netnum;
476: PetscScalar Sbase = 0.0;
478: /* Water subnetwork */
479: AppCtx_Water *appctx_water = &user.appctx_water;
480: WATERDATA *waterdata = NULL;
481: char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
482: PetscInt *edgelist_water = NULL, water_netnum;
484: /* Shared vertices between subnetworks */
485: PetscInt power_svtx, water_svtx;
487: PetscFunctionBeginUser;
488: PetscCall(PetscInitialize(&argc, &argv, "ex1options", help));
489: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
490: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
492: /* (1) Read Data - Only rank 0 reads the data */
493: PetscCall(PetscLogStageRegister("Read Data", &stage[0]));
494: PetscCall(PetscLogStagePush(stage[0]));
496: for (i = 0; i < Nsubnet; i++) {
497: numVertices[i] = 0;
498: numEdges[i] = 0;
499: }
501: /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
502: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
503: PetscCall(PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, PETSC_MAX_PATH_LEN - 1, NULL));
504: PetscCall(PetscNew(&pfdata));
505: PetscCall(PFReadMatPowerData(pfdata, pfdata_file));
506: if (rank == 0) PetscCall(PetscPrintf(PETSC_COMM_SELF, "Power network: nb = %" PetscInt_FMT ", ngen = %" PetscInt_FMT ", nload = %" PetscInt_FMT ", nbranch = %" PetscInt_FMT "\n", pfdata->nbus, pfdata->ngen, pfdata->nload, pfdata->nbranch));
507: Sbase = pfdata->sbase;
508: if (rank == 0) { /* proc[0] will create Electric Power Grid */
509: numEdges[0] = pfdata->nbranch;
510: numVertices[0] = pfdata->nbus;
512: PetscCall(PetscMalloc1(2 * numEdges[0], &edgelist_power));
513: PetscCall(GetListofEdges_Power(pfdata, edgelist_power));
514: }
515: /* Broadcast power Sbase to all processors */
516: PetscCallMPI(MPI_Bcast(&Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD));
517: appctx_power->Sbase = Sbase;
518: appctx_power->jac_error = PETSC_FALSE;
519: /* If external option activated. Introduce error in jacobian */
520: PetscCall(PetscOptionsHasName(NULL, NULL, "-jac_error", &appctx_power->jac_error));
522: /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
523: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
524: PetscCall(PetscNew(&waterdata));
525: PetscCall(PetscOptionsGetString(NULL, NULL, "-waterdata", waterdata_file, PETSC_MAX_PATH_LEN - 1, NULL));
526: PetscCall(WaterReadData(waterdata, waterdata_file));
527: if (size == 1 || (size > 1 && rank == 1)) {
528: PetscCall(PetscCalloc1(2 * waterdata->nedge, &edgelist_water));
529: PetscCall(GetListofEdges_Water(waterdata, edgelist_water));
530: numEdges[1] = waterdata->nedge;
531: numVertices[1] = waterdata->nvertex;
532: }
533: PetscCall(PetscLogStagePop());
535: /* (2) Create a network consist of two subnetworks */
536: PetscCall(PetscLogStageRegister("Net Setup", &stage[1]));
537: PetscCall(PetscLogStagePush(stage[1]));
539: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewCSV", &viewCSV, NULL));
540: PetscCall(PetscOptionsGetBool(NULL, NULL, "-printCoord", &printCoord, NULL));
541: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test", &test, NULL));
542: PetscCall(PetscOptionsGetBool(NULL, NULL, "-distribute", &distribute, NULL));
544: /* Create an empty network object */
545: PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
547: /* Register the components in the network */
548: PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &appctx_power->compkey_branch));
549: PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &appctx_power->compkey_bus));
550: PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &appctx_power->compkey_gen));
551: PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &appctx_power->compkey_load));
553: PetscCall(DMNetworkRegisterComponent(networkdm, "edge_water", sizeof(struct _p_EDGE_Water), &appctx_water->compkey_edge));
554: PetscCall(DMNetworkRegisterComponent(networkdm, "vertex_water", sizeof(struct _p_VERTEX_Water), &appctx_water->compkey_vtx));
556: PetscCall(PetscSynchronizedPrintf(PETSC_COMM_WORLD, "[%d] Total local nvertices %" PetscInt_FMT " + %" PetscInt_FMT " = %" PetscInt_FMT ", nedges %" PetscInt_FMT " + %" PetscInt_FMT " = %" PetscInt_FMT "\n", rank, numVertices[0], numVertices[1], numVertices[0] + numVertices[1], numEdges[0], numEdges[1], numEdges[0] + numEdges[1]));
557: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
559: PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, Nsubnet));
560: PetscCall(DMNetworkAddSubnetwork(networkdm, "power", numEdges[0], edgelist_power, &power_netnum));
561: PetscCall(DMNetworkAddSubnetwork(networkdm, "water", numEdges[1], edgelist_water, &water_netnum));
563: /* vertex subnet[0].4 shares with vertex subnet[1].0 */
564: power_svtx = 4;
565: water_svtx = 0;
566: PetscCall(DMNetworkAddSharedVertices(networkdm, power_netnum, water_netnum, 1, &power_svtx, &water_svtx));
568: /* Set up the network layout */
569: PetscCall(DMNetworkLayoutSetUp(networkdm));
571: /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
572: genj = 0;
573: loadj = 0;
574: PetscCall(DMNetworkGetSubnetwork(networkdm, power_netnum, &nv, &ne, &vtx, &edges));
576: for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx_power->compkey_branch, &pfdata->branch[i], 0));
578: for (i = 0; i < nv; i++) {
579: PetscCall(DMNetworkIsSharedVertex(networkdm, vtx[i], &flg));
580: if (flg) continue;
582: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[i], 2));
583: if (pfdata->bus[i].ngen) {
584: for (j = 0; j < pfdata->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_gen, &pfdata->gen[genj++], 0));
585: }
586: if (pfdata->bus[i].nload) {
587: for (j = 0; j < pfdata->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[loadj++], 0));
588: }
589: }
591: /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
592: PetscCall(DMNetworkGetSubnetwork(networkdm, water_netnum, &nv, &ne, &vtx, &edges));
593: for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx_water->compkey_edge, &waterdata->edge[i], 0));
595: for (i = 0; i < nv; i++) {
596: PetscCall(DMNetworkIsSharedVertex(networkdm, vtx[i], &flg));
597: if (flg) continue;
599: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[i], 1));
600: }
602: /* ADD VARIABLES AND COMPONENTS AT THE SHARED VERTEX: net[0].4 coupls with net[1].0 -- owning and all ghost ranks of the vertex do this */
603: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
604: for (i = 0; i < nv; i++) {
605: /* power */
606: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[4], 2));
607: /* bus[4] is a load, add its component */
608: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[0], 0));
610: /* water */
611: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[0], 1));
612: }
614: /* Set coordinates for visualization */
615: PetscCall(DMSetCoordinateDim(networkdm, 2));
616: PetscCall(DMGetCoordinateDM(networkdm, &dmcoords));
617: PetscCall(DMNetworkGetVertexRange(dmcoords, &vStart, &vEnd));
619: PetscCall(PetscCalloc1(vEnd - vStart, &color));
620: PetscCall(DMNetworkRegisterComponent(dmcoords, "coordinate&color", sizeof(PetscReal), &compkey));
621: for (i = vStart; i < vEnd; i++) PetscCall(DMNetworkAddComponent(dmcoords, i, compkey, &color[i - vStart], 2));
622: PetscCall(DMNetworkFinalizeComponents(dmcoords));
624: PetscCall(DMCreateLocalVector(dmcoords, &coords));
625: PetscCall(DMSetCoordinatesLocal(networkdm, coords)); /* set/get coords to/from networkdm */
626: PetscCall(CoordinateVecSetUp(dmcoords, coords));
627: if (printCoord) PetscCall(CoordinatePrint(networkdm));
629: /* Set up DM for use */
630: PetscCall(DMSetUp(networkdm));
632: /* Free user objects */
633: PetscCall(PetscFree(edgelist_power));
634: PetscCall(PetscFree(pfdata->bus));
635: PetscCall(PetscFree(pfdata->gen));
636: PetscCall(PetscFree(pfdata->branch));
637: PetscCall(PetscFree(pfdata->load));
638: PetscCall(PetscFree(pfdata));
640: PetscCall(PetscFree(edgelist_water));
641: PetscCall(PetscFree(waterdata->vertex));
642: PetscCall(PetscFree(waterdata->edge));
643: PetscCall(PetscFree(waterdata));
645: /* Re-distribute networkdm to multiple processes for better job balance */
646: if (distribute) {
647: PetscCall(DMNetworkDistribute(&networkdm, 0));
649: if (printCoord) PetscCall(CoordinatePrint(networkdm));
650: if (viewCSV) { /* CSV View of network with coordinates */
651: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_CSV));
652: PetscCall(DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD));
653: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
654: }
655: }
656: PetscCall(VecDestroy(&coords));
658: /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
659: if (test) {
660: PetscInt v, gidx;
661: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
662: for (i = 0; i < Nsubnet; i++) {
663: PetscCall(DMNetworkGetSubnetwork(networkdm, i, &nv, &ne, &vtx, &edges));
664: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, subnet[%" PetscInt_FMT "] ne %" PetscInt_FMT ", nv %" PetscInt_FMT "\n", rank, i, ne, nv));
665: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
667: for (v = 0; v < nv; v++) {
668: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghost));
669: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx));
670: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] subnet[%" PetscInt_FMT "] v %" PetscInt_FMT " %" PetscInt_FMT "; ghost %d\n", rank, i, vtx[v], gidx, ghost));
671: }
672: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
673: }
674: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
676: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
677: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, num of shared vertices nsv = %" PetscInt_FMT "\n", rank, nv));
678: for (v = 0; v < nv; v++) {
679: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx));
680: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] sv %" PetscInt_FMT ", gidx=%" PetscInt_FMT "\n", rank, vtx[v], gidx));
681: }
682: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
683: }
685: /* Create solution vector X */
686: PetscCall(DMCreateGlobalVector(networkdm, &X));
687: PetscCall(VecDuplicate(X, &F));
688: PetscCall(DMGetLocalVector(networkdm, &user.localXold));
689: PetscCall(PetscLogStagePop());
691: /* (3) Setup Solvers */
692: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewJ", &viewJ, NULL));
693: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewX", &viewX, NULL));
695: PetscCall(PetscLogStageRegister("SNES Setup", &stage[2]));
696: PetscCall(PetscLogStagePush(stage[2]));
698: PetscCall(SetInitialGuess(networkdm, X, &user));
700: /* Create coupled snes */
701: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_coupled setup ......\n"));
702: user.subsnes_id = Nsubnet;
703: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
704: PetscCall(SNESSetDM(snes, networkdm));
705: PetscCall(SNESSetOptionsPrefix(snes, "coupled_"));
706: PetscCall(SNESSetFunction(snes, F, FormFunction, &user));
707: /* set maxit=1 which can be changed via option '-coupled_snes_max_it <>', see ex1options */
708: PetscCall(SNESSetTolerances(snes, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1, PETSC_DEFAULT));
709: PetscCall(SNESSetFromOptions(snes));
711: if (viewJ) {
712: /* View Jac structure */
713: PetscCall(SNESGetJacobian(snes, &Jac, NULL, NULL, NULL));
714: PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
715: }
717: if (viewX) {
718: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution:\n"));
719: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
720: }
722: if (viewJ) {
723: /* View assembled Jac */
724: PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
725: }
727: /* Create snes_power */
728: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_power setup ......\n"));
729: user.subsnes_id = 0;
730: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_power));
731: PetscCall(SNESSetDM(snes_power, networkdm));
732: PetscCall(SNESSetOptionsPrefix(snes_power, "power_"));
733: PetscCall(SNESSetFunction(snes_power, F, FormFunction, &user));
734: /* set maxit=1 which can be changed via option '-power_snes_max_it <>', see ex1options */
735: PetscCall(SNESSetTolerances(snes_power, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1, PETSC_DEFAULT));
737: /* Use user-provide Jacobian */
738: PetscCall(DMCreateMatrix(networkdm, &Jac));
739: PetscCall(SNESSetJacobian(snes_power, Jac, Jac, FormJacobian_subPower, &user));
740: PetscCall(SNESSetFromOptions(snes_power));
742: if (viewX) {
743: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Power Solution:\n"));
744: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
745: }
746: if (viewJ) {
747: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Power Jac:\n"));
748: PetscCall(SNESGetJacobian(snes_power, &Jac, NULL, NULL, NULL));
749: PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
750: }
752: /* Create snes_water */
753: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_water setup......\n"));
755: user.subsnes_id = 1;
756: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_water));
757: PetscCall(SNESSetDM(snes_water, networkdm));
758: PetscCall(SNESSetOptionsPrefix(snes_water, "water_"));
759: PetscCall(SNESSetFunction(snes_water, F, FormFunction, &user));
760: /* set maxit=1 which can be changed via option '-water_snes_max_it <>', see ex1options */
761: PetscCall(SNESSetTolerances(snes_water, PETSC_DEFAULT, PETSC_DEFAULT, PETSC_DEFAULT, 1, PETSC_DEFAULT));
762: PetscCall(SNESSetFromOptions(snes_water));
764: if (viewX) {
765: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Water Solution:\n"));
766: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
767: }
769: /* Monitor snes, snes_power and snes_water iterations */
770: PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitorIteration", &monitorIt, NULL));
771: user.monitorColor = PETSC_FALSE;
772: PetscCall(PetscOptionsGetBool(NULL, NULL, "-monitorColor", &user.monitorColor, NULL));
773: if (user.monitorColor) monitorIt = PETSC_TRUE; /* require installation of pandas and matplotlib */
774: if (monitorIt) {
775: PetscCall(SNESMonitorSet(snes_power, UserMonitor, &user, NULL));
776: PetscCall(SNESMonitorSet(snes_water, UserMonitor, &user, NULL));
777: PetscCall(SNESMonitorSet(snes, UserMonitor, &user, NULL));
778: }
779: PetscCall(PetscLogStagePop());
781: /* (4) Solve: we must update user.localXold after each call of SNESSolve().
782: See "PETSc DMNetwork: A Library for Scalable Network PDE-Based Multiphysics Simulations",
783: https://dl.acm.org/doi/10.1145/3344587
784: */
785: PetscCall(PetscLogStageRegister("SNES Solve", &stage[3]));
786: PetscCall(PetscLogStagePush(stage[3]));
787: user.it = 0; /* total_snes_it */
788: reason = SNES_DIVERGED_DTOL;
789: while (user.it < it_max && (PetscInt)reason < 0) {
790: user.subsnes_id = 0;
791: PetscCall(SNESSolve(snes_power, NULL, X));
792: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, user.localXold));
793: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, user.localXold));
795: user.subsnes_id = 1;
796: PetscCall(SNESSolve(snes_water, NULL, X));
797: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, user.localXold));
798: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, user.localXold));
800: user.subsnes_id = Nsubnet;
801: PetscCall(SNESSolve(snes, NULL, X));
802: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, user.localXold));
803: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, user.localXold));
805: PetscCall(SNESGetConvergedReason(snes, &reason));
806: user.it++;
807: }
808: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coupled_SNES converged in %" PetscInt_FMT " iterations\n", user.it));
809: if (viewX) {
810: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final Solution:\n"));
811: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
812: }
813: PetscCall(PetscLogStagePop());
815: /* Free objects */
816: /* -------------*/
817: PetscCall(PetscFree(color));
818: PetscCall(VecDestroy(&X));
819: PetscCall(VecDestroy(&F));
820: PetscCall(DMRestoreLocalVector(networkdm, &user.localXold));
822: PetscCall(SNESDestroy(&snes));
823: PetscCall(MatDestroy(&Jac));
824: PetscCall(SNESDestroy(&snes_power));
825: PetscCall(SNESDestroy(&snes_water));
827: PetscCall(DMDestroy(&networkdm));
828: PetscCall(PetscFinalize());
829: return 0;
830: }
832: /*TEST
834: build:
835: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
836: depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c
838: test:
839: args: -options_left no -dmnetwork_view -fp_trap 0
840: localrunfiles: ex1options power/case9.m water/sample1.inp
841: output_file: output/ex1.out
843: test:
844: suffix: 2
845: nsize: 3
846: args: -options_left no -petscpartitioner_type parmetis -fp_trap 0
847: localrunfiles: ex1options power/case9.m water/sample1.inp
848: output_file: output/ex1_2.out
849: requires: parmetis
851: test:
852: suffix: 3
853: nsize: 3
854: args: -options_left no -distribute false -fp_trap 0
855: localrunfiles: ex1options power/case9.m water/sample1.inp
856: output_file: output/ex1_2.out
858: test:
859: suffix: 4
860: nsize: 4
861: args: -options_left no -petscpartitioner_type simple -dmnetwork_view -dmnetwork_view_distributed -fp_trap 0
862: localrunfiles: ex1options power/case9.m water/sample1.inp
863: output_file: output/ex1_4.out
865: test:
866: suffix: 5
867: args: -options_left no -viewCSV -fp_trap 0
868: localrunfiles: ex1options power/case9.m water/sample1.inp
869: output_file: output/ex1_5.out
871: test:
872: suffix: 6
873: nsize: 3
874: args: -options_left no -petscpartitioner_type parmetis -dmnetwork_view_distributed draw:null -fp_trap 0
875: localrunfiles: ex1options power/case9.m water/sample1.inp
876: output_file: output/ex1_2.out
877: requires: parmetis
879: TEST*/