Actual source code: ex1.c
1: static char help[] = "This example demonstrates the use of DMNetwork interface 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: This example shows the use of subnetwork feature in DMNetwork. \n\
6: Run this program: mpiexec -n <n> ./ex1 \n\\n";
8: #include "power/power.h"
9: #include "water/water.h"
11: typedef struct {
12: UserCtx_Power appctx_power;
13: AppCtx_Water appctx_water;
14: PetscInt subsnes_id; /* snes solver id */
15: PetscInt it; /* iteration number */
16: Vec localXold; /* store previous solution, used by FormFunction_Dummy() */
17: } UserCtx;
19: PetscErrorCode UserMonitor(SNES snes, PetscInt its, PetscReal fnorm, void *appctx)
20: {
21: UserCtx *user = (UserCtx *)appctx;
22: Vec X, localXold = user->localXold;
23: DM networkdm;
24: PetscMPIInt rank;
25: MPI_Comm comm;
27: PetscFunctionBegin;
28: PetscCall(PetscObjectGetComm((PetscObject)snes, &comm));
29: PetscCallMPI(MPI_Comm_rank(comm, &rank));
30: #if 0
31: if (rank == 0) {
32: PetscInt subsnes_id = user->subsnes_id;
33: if (subsnes_id == 2) {
34: PetscCall(PetscPrintf(PETSC_COMM_SELF," it %" PetscInt_FMT ", subsnes_id %" PetscInt_FMT ", fnorm %g\n",user->it,user->subsnes_id,(double)fnorm));
35: } else {
36: PetscCall(PetscPrintf(PETSC_COMM_SELF," subsnes_id %" PetscInt_FMT ", fnorm %g\n",user->subsnes_id,(double)fnorm));
37: }
38: }
39: #endif
40: PetscCall(SNESGetSolution(snes, &X));
41: PetscCall(SNESGetDM(snes, &networkdm));
42: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localXold));
43: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localXold));
44: PetscFunctionReturn(PETSC_SUCCESS);
45: }
47: PetscErrorCode FormJacobian_subPower(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx)
48: {
49: DM networkdm;
50: Vec localX;
51: PetscInt nv, ne, i, j, offset, nvar, row;
52: const PetscInt *vtx, *edges;
53: PetscBool ghostvtex;
54: PetscScalar one = 1.0;
55: PetscMPIInt rank;
56: MPI_Comm comm;
58: PetscFunctionBegin;
59: PetscCall(SNESGetDM(snes, &networkdm));
60: PetscCall(DMGetLocalVector(networkdm, &localX));
62: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
63: PetscCallMPI(MPI_Comm_rank(comm, &rank));
65: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
66: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
68: PetscCall(MatZeroEntries(J));
70: /* Power subnetwork: copied from power/FormJacobian_Power() */
71: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
72: PetscCall(FormJacobian_Power_private(networkdm, localX, J, nv, ne, vtx, edges, appctx));
74: /* Water subnetwork: Identity */
75: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
76: for (i = 0; i < nv; i++) {
77: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex));
78: if (ghostvtex) continue;
80: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset));
81: PetscCall(DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar));
82: for (j = 0; j < nvar; j++) {
83: row = offset + j;
84: PetscCall(MatSetValues(J, 1, &row, 1, &row, &one, ADD_VALUES));
85: }
86: }
87: PetscCall(MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY));
88: PetscCall(MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY));
90: PetscCall(DMRestoreLocalVector(networkdm, &localX));
91: PetscFunctionReturn(PETSC_SUCCESS);
92: }
94: /* Dummy equation localF(X) = localX - localXold */
95: PetscErrorCode FormFunction_Dummy(DM networkdm, Vec localX, Vec localF, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
96: {
97: const PetscScalar *xarr, *xoldarr;
98: PetscScalar *farr;
99: PetscInt i, j, offset, nvar;
100: PetscBool ghostvtex;
101: UserCtx *user = (UserCtx *)appctx;
102: Vec localXold = user->localXold;
104: PetscFunctionBegin;
105: PetscCall(VecGetArrayRead(localX, &xarr));
106: PetscCall(VecGetArrayRead(localXold, &xoldarr));
107: PetscCall(VecGetArray(localF, &farr));
109: for (i = 0; i < nv; i++) {
110: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex));
111: if (ghostvtex) continue;
113: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset));
114: PetscCall(DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar));
115: for (j = 0; j < nvar; j++) farr[offset + j] = xarr[offset + j] - xoldarr[offset + j];
116: }
118: PetscCall(VecRestoreArrayRead(localX, &xarr));
119: PetscCall(VecRestoreArrayRead(localXold, &xoldarr));
120: PetscCall(VecRestoreArray(localF, &farr));
121: PetscFunctionReturn(PETSC_SUCCESS);
122: }
124: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *appctx)
125: {
126: DM networkdm;
127: Vec localX, localF;
128: PetscInt nv, ne, v;
129: const PetscInt *vtx, *edges;
130: PetscMPIInt rank;
131: MPI_Comm comm;
132: UserCtx *user = (UserCtx *)appctx;
133: UserCtx_Power appctx_power = (*user).appctx_power;
134: AppCtx_Water appctx_water = (*user).appctx_water;
136: PetscFunctionBegin;
137: PetscCall(SNESGetDM(snes, &networkdm));
138: PetscCall(PetscObjectGetComm((PetscObject)networkdm, &comm));
139: PetscCallMPI(MPI_Comm_rank(comm, &rank));
141: PetscCall(DMGetLocalVector(networkdm, &localX));
142: PetscCall(DMGetLocalVector(networkdm, &localF));
143: PetscCall(VecSet(F, 0.0));
144: PetscCall(VecSet(localF, 0.0));
146: PetscCall(DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX));
147: PetscCall(DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX));
149: /* Form Function for power subnetwork */
150: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
151: if (user->subsnes_id == 1) { /* snes_water only */
152: PetscCall(FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user));
153: } else {
154: PetscCall(FormFunction_Power(networkdm, localX, localF, nv, ne, vtx, edges, &appctx_power));
155: }
157: /* Form Function for water subnetwork */
158: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
159: if (user->subsnes_id == 0) { /* snes_power only */
160: PetscCall(FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user));
161: } else {
162: PetscCall(FormFunction_Water(networkdm, localX, localF, nv, ne, vtx, edges, NULL));
163: }
165: /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */
166: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
167: for (v = 0; v < nv; v++) {
168: PetscInt key, ncomp, nvar, nconnedges, k, e, keye, goffset[3];
169: void *component;
170: const PetscInt *connedges;
172: PetscCall(DMNetworkGetComponent(networkdm, vtx[v], ALL_COMPONENTS, NULL, NULL, &nvar));
173: PetscCall(DMNetworkGetNumComponents(networkdm, vtx[v], &ncomp));
174: /* printf(" [%d] coupling vertex[%" PetscInt_FMT "]: v %" PetscInt_FMT ", ncomp %" PetscInt_FMT "; nvar %" PetscInt_FMT "\n",rank,v,vtx[v], ncomp,nvar); */
176: for (k = 0; k < ncomp; k++) {
177: PetscCall(DMNetworkGetComponent(networkdm, vtx[v], k, &key, &component, &nvar));
178: PetscCall(DMNetworkGetGlobalVecOffset(networkdm, vtx[v], k, &goffset[k]));
180: /* Verify the coupling vertex is a powernet load vertex or a water vertex */
181: switch (k) {
182: case 0:
183: 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);
184: break;
185: case 1:
186: 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");
187: break;
188: case 2:
189: PetscCheck(key == appctx_water.compkey_vtx && nvar == 1 && goffset[2] == goffset[1], PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a water vertex");
190: break;
191: default:
192: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %" PetscInt_FMT " is wrong", k);
193: }
194: /* printf(" [%d] coupling vertex[%" PetscInt_FMT "]: key %" PetscInt_FMT "; nvar %" PetscInt_FMT ", goffset %" PetscInt_FMT "\n",rank,v,key,nvar,goffset[k]); */
195: }
197: /* Get its supporting edges */
198: PetscCall(DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges));
199: /* printf("\n[%d] coupling vertex: nconnedges %" PetscInt_FMT "\n",rank,nconnedges); */
200: for (k = 0; k < nconnedges; k++) {
201: e = connedges[k];
202: PetscCall(DMNetworkGetNumComponents(networkdm, e, &ncomp));
203: /* printf("\n [%d] connected edge[%" PetscInt_FMT "]=%" PetscInt_FMT " has ncomp %" PetscInt_FMT "\n",rank,k,e,ncomp); */
204: PetscCall(DMNetworkGetComponent(networkdm, e, 0, &keye, &component, NULL));
205: if (keye != appctx_water.compkey_edge) PetscCheck(keye == appctx_power.compkey_branch, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "Not a power branch");
206: }
207: }
209: PetscCall(DMRestoreLocalVector(networkdm, &localX));
211: PetscCall(DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F));
212: PetscCall(DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F));
213: PetscCall(DMRestoreLocalVector(networkdm, &localF));
214: #if 0
215: if (rank == 0) printf("F:\n");
216: PetscCall(VecView(F,PETSC_VIEWER_STDOUT_WORLD));
217: #endif
218: PetscFunctionReturn(PETSC_SUCCESS);
219: }
221: PetscErrorCode SetInitialGuess(DM networkdm, Vec X, void *appctx)
222: {
223: PetscInt nv, ne, i, j, ncomp, offset, key;
224: const PetscInt *vtx, *edges;
225: UserCtx *user = (UserCtx *)appctx;
226: Vec localX = user->localXold;
227: UserCtx_Power appctx_power = (*user).appctx_power;
228: AppCtx_Water appctx_water = (*user).appctx_water;
229: PetscBool ghost;
230: PetscScalar *xarr;
231: VERTEX_Power bus;
232: VERTEX_Water vertex;
233: void *component;
234: GEN gen;
236: PetscFunctionBegin;
237: PetscCall(VecSet(X, 0.0));
238: PetscCall(VecSet(localX, 0.0));
240: /* Set initial guess for power subnetwork */
241: PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
242: PetscCall(SetInitialGuess_Power(networkdm, localX, nv, ne, vtx, edges, &appctx_power));
244: /* Set initial guess for water subnetwork */
245: PetscCall(DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges));
246: PetscCall(SetInitialGuess_Water(networkdm, localX, nv, ne, vtx, edges, NULL));
248: /* Set initial guess at the coupling vertex */
249: PetscCall(VecGetArray(localX, &xarr));
250: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
251: for (i = 0; i < nv; i++) {
252: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[i], &ghost));
253: if (ghost) continue;
255: PetscCall(DMNetworkGetNumComponents(networkdm, vtx[i], &ncomp));
256: for (j = 0; j < ncomp; j++) {
257: PetscCall(DMNetworkGetLocalVecOffset(networkdm, vtx[i], j, &offset));
258: PetscCall(DMNetworkGetComponent(networkdm, vtx[i], j, &key, (void **)&component, NULL));
259: if (key == appctx_power.compkey_bus) {
260: bus = (VERTEX_Power)(component);
261: xarr[offset] = bus->va * PETSC_PI / 180.0;
262: xarr[offset + 1] = bus->vm;
263: } else if (key == appctx_power.compkey_gen) {
264: gen = (GEN)(component);
265: if (!gen->status) continue;
266: xarr[offset + 1] = gen->vs;
267: } else if (key == appctx_water.compkey_vtx) {
268: vertex = (VERTEX_Water)(component);
269: if (vertex->type == VERTEX_TYPE_JUNCTION) {
270: xarr[offset] = 100;
271: } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
272: xarr[offset] = vertex->res.head;
273: } else {
274: xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
275: }
276: }
277: }
278: }
279: PetscCall(VecRestoreArray(localX, &xarr));
281: PetscCall(DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X));
282: PetscCall(DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X));
283: PetscFunctionReturn(PETSC_SUCCESS);
284: }
286: /* Set coordinates */
287: static PetscErrorCode CoordinateVecSetUp(DM dm, Vec coords)
288: {
289: DM dmclone;
290: PetscInt i, gidx, offset, v, nv, Nsubnet;
291: const PetscInt *vtx;
292: PetscScalar *carray;
294: PetscFunctionBeginUser;
295: PetscCall(DMGetCoordinateDM(dm, &dmclone));
296: PetscCall(VecGetArrayWrite(coords, &carray));
297: PetscCall(DMNetworkGetNumSubNetworks(dm, NULL, &Nsubnet));
298: for (i = 0; i < Nsubnet; i++) {
299: PetscCall(DMNetworkGetSubnetwork(dm, i, &nv, NULL, &vtx, NULL));
300: for (v = 0; v < nv; v++) {
301: PetscCall(DMNetworkGetGlobalVertexIndex(dm, vtx[v], &gidx));
302: PetscCall(DMNetworkGetLocalVecOffset(dmclone, vtx[v], 0, &offset));
303: switch (gidx) {
304: case 0:
305: carray[offset] = -1.0;
306: carray[offset + 1] = -1.0;
307: break;
308: case 1:
309: carray[offset] = -2.0;
310: carray[offset + 1] = 2.0;
311: break;
312: case 2:
313: carray[offset] = 0.0;
314: carray[offset + 1] = 2.0;
315: break;
316: case 3:
317: carray[offset] = -1.0;
318: carray[offset + 1] = 0.0;
319: break;
320: case 4:
321: carray[offset] = 0.0;
322: carray[offset + 1] = 0.0;
323: break;
324: case 5:
325: carray[offset] = 0.0;
326: carray[offset + 1] = 1.0;
327: break;
328: case 6:
329: carray[offset] = -1.0;
330: carray[offset + 1] = 1.0;
331: break;
332: case 7:
333: carray[offset] = -2.0;
334: carray[offset + 1] = 1.0;
335: break;
336: case 8:
337: carray[offset] = -2.0;
338: carray[offset + 1] = 0.0;
339: break;
340: case 9:
341: carray[offset] = 1.0;
342: carray[offset + 1] = 0.0;
343: break;
344: case 10:
345: carray[offset] = 1.0;
346: carray[offset + 1] = -1.0;
347: break;
348: case 11:
349: carray[offset] = 2.0;
350: carray[offset + 1] = -1.0;
351: break;
352: case 12:
353: carray[offset] = 2.0;
354: carray[offset + 1] = 0.0;
355: break;
356: case 13:
357: carray[offset] = 0.0;
358: carray[offset + 1] = -1.0;
359: break;
360: case 14:
361: carray[offset] = 2.0;
362: carray[offset + 1] = 1.0;
363: break;
364: default:
365: PetscCheck(gidx < 15 && gidx > -1, PETSC_COMM_SELF, PETSC_ERR_ARG_WRONG, "gidx %" PetscInt_FMT "must between 0 and 14", gidx);
366: }
367: }
368: }
369: PetscCall(VecRestoreArrayWrite(coords, &carray));
370: PetscFunctionReturn(PETSC_SUCCESS);
371: }
373: static PetscErrorCode CoordinatePrint(DM dm)
374: {
375: DM dmclone;
376: PetscInt cdim, v, off, vglobal, vStart, vEnd;
377: const PetscScalar *carray;
378: Vec coords;
379: MPI_Comm comm;
380: PetscMPIInt rank;
382: PetscFunctionBegin;
383: PetscCall(PetscObjectGetComm((PetscObject)dm, &comm));
384: PetscCallMPI(MPI_Comm_rank(comm, &rank));
386: PetscCall(DMGetCoordinateDM(dm, &dmclone));
387: PetscCall(DMNetworkGetVertexRange(dm, &vStart, &vEnd));
388: PetscCall(DMGetCoordinatesLocal(dm, &coords));
390: PetscCall(DMGetCoordinateDim(dm, &cdim));
391: PetscCall(VecGetArrayRead(coords, &carray));
393: PetscCall(PetscPrintf(MPI_COMM_WORLD, "\nCoordinatePrint, cdim %" PetscInt_FMT ":\n", cdim));
394: PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "[%i]\n", rank));
395: for (v = vStart; v < vEnd; v++) {
396: PetscCall(DMNetworkGetLocalVecOffset(dmclone, v, 0, &off));
397: PetscCall(DMNetworkGetGlobalVertexIndex(dmclone, v, &vglobal));
398: switch (cdim) {
399: case 2:
400: PetscCall(PetscSynchronizedPrintf(MPI_COMM_WORLD, "Vertex: %" PetscInt_FMT ", x = %f y = %f \n", vglobal, (double)PetscRealPart(carray[off]), (double)PetscRealPart(carray[off + 1])));
401: break;
402: default:
403: PetscCheck(cdim == 2, MPI_COMM_WORLD, PETSC_ERR_SUP, "Only supports Network embedding dimension of 2, not supplied %" PetscInt_FMT, cdim);
404: break;
405: }
406: }
407: PetscCall(PetscSynchronizedFlush(MPI_COMM_WORLD, NULL));
408: PetscCall(VecRestoreArrayRead(coords, &carray));
409: PetscFunctionReturn(PETSC_SUCCESS);
410: }
412: int main(int argc, char **argv)
413: {
414: DM networkdm, dmclone;
415: PetscMPIInt rank, size;
416: PetscInt Nsubnet = 2, numVertices[2], numEdges[2], i, j, nv, ne, it_max = 10;
417: PetscInt vStart, vEnd, compkey;
418: const PetscInt *vtx, *edges;
419: Vec X, F, coords;
420: SNES snes, snes_power, snes_water;
421: Mat Jac;
422: PetscBool ghost, viewJ = PETSC_FALSE, viewX = PETSC_FALSE, test = PETSC_FALSE, distribute = PETSC_TRUE, flg, printCoord = PETSC_FALSE, viewCSV = PETSC_FALSE;
423: UserCtx user;
424: SNESConvergedReason reason;
425: #if defined(PETSC_USE_LOG)
426: PetscLogStage stage[4];
427: #endif
429: /* Power subnetwork */
430: UserCtx_Power *appctx_power = &user.appctx_power;
431: char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
432: PFDATA *pfdata = NULL;
433: PetscInt genj, loadj, *edgelist_power = NULL, power_netnum;
434: PetscScalar Sbase = 0.0;
436: /* Water subnetwork */
437: AppCtx_Water *appctx_water = &user.appctx_water;
438: WATERDATA *waterdata = NULL;
439: char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
440: PetscInt *edgelist_water = NULL, water_netnum;
442: /* Shared vertices between subnetworks */
443: PetscInt power_svtx, water_svtx;
445: PetscFunctionBeginUser;
446: PetscCall(PetscInitialize(&argc, &argv, "ex1options", help));
447: PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &rank));
448: PetscCallMPI(MPI_Comm_size(PETSC_COMM_WORLD, &size));
450: /* (1) Read Data - Only rank 0 reads the data */
451: PetscCall(PetscLogStageRegister("Read Data", &stage[0]));
452: PetscCall(PetscLogStagePush(stage[0]));
454: for (i = 0; i < Nsubnet; i++) {
455: numVertices[i] = 0;
456: numEdges[i] = 0;
457: }
459: /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
460: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
461: PetscCall(PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, PETSC_MAX_PATH_LEN - 1, NULL));
462: PetscCall(PetscNew(&pfdata));
463: PetscCall(PFReadMatPowerData(pfdata, pfdata_file));
464: 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));
465: Sbase = pfdata->sbase;
466: if (rank == 0) { /* proc[0] will create Electric Power Grid */
467: numEdges[0] = pfdata->nbranch;
468: numVertices[0] = pfdata->nbus;
470: PetscCall(PetscMalloc1(2 * numEdges[0], &edgelist_power));
471: PetscCall(GetListofEdges_Power(pfdata, edgelist_power));
472: }
473: /* Broadcast power Sbase to all processors */
474: PetscCallMPI(MPI_Bcast(&Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD));
475: appctx_power->Sbase = Sbase;
476: appctx_power->jac_error = PETSC_FALSE;
477: /* If external option activated. Introduce error in jacobian */
478: PetscCall(PetscOptionsHasName(NULL, NULL, "-jac_error", &appctx_power->jac_error));
480: /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
481: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
482: PetscCall(PetscNew(&waterdata));
483: PetscCall(PetscOptionsGetString(NULL, NULL, "-waterdata", waterdata_file, PETSC_MAX_PATH_LEN - 1, NULL));
484: PetscCall(WaterReadData(waterdata, waterdata_file));
485: if (size == 1 || (size > 1 && rank == 1)) {
486: PetscCall(PetscCalloc1(2 * waterdata->nedge, &edgelist_water));
487: PetscCall(GetListofEdges_Water(waterdata, edgelist_water));
488: numEdges[1] = waterdata->nedge;
489: numVertices[1] = waterdata->nvertex;
490: }
491: PetscCall(PetscLogStagePop());
493: /* (2) Create a network consist of two subnetworks */
494: PetscCall(PetscLogStageRegister("Net Setup", &stage[1]));
495: PetscCall(PetscLogStagePush(stage[1]));
497: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewCSV", &viewCSV, NULL));
498: PetscCall(PetscOptionsGetBool(NULL, NULL, "-printCoord", &printCoord, NULL));
499: PetscCall(PetscOptionsGetBool(NULL, NULL, "-test", &test, NULL));
500: PetscCall(PetscOptionsGetBool(NULL, NULL, "-distribute", &distribute, NULL));
502: /* Create an empty network object */
503: PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));
505: /* Register the components in the network */
506: PetscCall(DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &appctx_power->compkey_branch));
507: PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &appctx_power->compkey_bus));
508: PetscCall(DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &appctx_power->compkey_gen));
509: PetscCall(DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &appctx_power->compkey_load));
511: PetscCall(DMNetworkRegisterComponent(networkdm, "edge_water", sizeof(struct _p_EDGE_Water), &appctx_water->compkey_edge));
512: PetscCall(DMNetworkRegisterComponent(networkdm, "vertex_water", sizeof(struct _p_VERTEX_Water), &appctx_water->compkey_vtx));
514: 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]));
515: PetscCall(PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT));
517: PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, Nsubnet));
518: PetscCall(DMNetworkAddSubnetwork(networkdm, "power", numEdges[0], edgelist_power, &power_netnum));
519: PetscCall(DMNetworkAddSubnetwork(networkdm, "water", numEdges[1], edgelist_water, &water_netnum));
521: /* vertex subnet[0].4 shares with vertex subnet[1].0 */
522: power_svtx = 4;
523: water_svtx = 0;
524: PetscCall(DMNetworkAddSharedVertices(networkdm, power_netnum, water_netnum, 1, &power_svtx, &water_svtx));
526: /* Set up the network layout */
527: PetscCall(DMNetworkLayoutSetUp(networkdm));
529: /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
530: genj = 0;
531: loadj = 0;
532: PetscCall(DMNetworkGetSubnetwork(networkdm, power_netnum, &nv, &ne, &vtx, &edges));
534: for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx_power->compkey_branch, &pfdata->branch[i], 0));
536: for (i = 0; i < nv; i++) {
537: PetscCall(DMNetworkIsSharedVertex(networkdm, vtx[i], &flg));
538: if (flg) continue;
540: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[i], 2));
541: if (pfdata->bus[i].ngen) {
542: for (j = 0; j < pfdata->bus[i].ngen; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_gen, &pfdata->gen[genj++], 0));
543: }
544: if (pfdata->bus[i].nload) {
545: for (j = 0; j < pfdata->bus[i].nload; j++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[loadj++], 0));
546: }
547: }
549: /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
550: PetscCall(DMNetworkGetSubnetwork(networkdm, water_netnum, &nv, &ne, &vtx, &edges));
551: for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx_water->compkey_edge, &waterdata->edge[i], 0));
553: for (i = 0; i < nv; i++) {
554: PetscCall(DMNetworkIsSharedVertex(networkdm, vtx[i], &flg));
555: if (flg) continue;
557: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[i], 1));
558: }
560: /* 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 */
561: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
562: for (i = 0; i < nv; i++) {
563: /* power */
564: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[4], 2));
565: /* bus[4] is a load, add its component */
566: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[0], 0));
568: /* water */
569: PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[0], 1));
570: }
572: /* Set coordinates for visualization */
573: PetscCall(DMSetCoordinateDim(networkdm, 2));
574: PetscCall(DMGetCoordinateDM(networkdm, &dmclone));
575: PetscCall(DMNetworkGetVertexRange(dmclone, &vStart, &vEnd));
576: PetscCall(DMNetworkRegisterComponent(dmclone, "coordinates", 0, &compkey));
577: for (i = vStart; i < vEnd; i++) PetscCall(DMNetworkAddComponent(dmclone, i, compkey, NULL, 2));
578: PetscCall(DMNetworkFinalizeComponents(dmclone));
580: PetscCall(DMCreateLocalVector(dmclone, &coords));
581: PetscCall(DMSetCoordinatesLocal(networkdm, coords)); /* set/get coords to/from networkdm */
582: PetscCall(CoordinateVecSetUp(networkdm, coords));
583: if (printCoord) PetscCall(CoordinatePrint(networkdm));
585: /* Set up DM for use */
586: PetscCall(DMSetUp(networkdm));
588: /* Free user objects */
589: PetscCall(PetscFree(edgelist_power));
590: PetscCall(PetscFree(pfdata->bus));
591: PetscCall(PetscFree(pfdata->gen));
592: PetscCall(PetscFree(pfdata->branch));
593: PetscCall(PetscFree(pfdata->load));
594: PetscCall(PetscFree(pfdata));
596: PetscCall(PetscFree(edgelist_water));
597: PetscCall(PetscFree(waterdata->vertex));
598: PetscCall(PetscFree(waterdata->edge));
599: PetscCall(PetscFree(waterdata));
601: /* Re-distribute networkdm to multiple processes for better job balance */
602: if (distribute) {
603: PetscCall(DMNetworkDistribute(&networkdm, 0));
605: if (printCoord) PetscCall(CoordinatePrint(networkdm));
606: if (viewCSV) { /* CSV View of network with coordinates */
607: PetscCall(PetscViewerPushFormat(PETSC_VIEWER_STDOUT_WORLD, PETSC_VIEWER_ASCII_CSV));
608: PetscCall(DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD));
609: PetscCall(PetscViewerPopFormat(PETSC_VIEWER_STDOUT_WORLD));
610: }
611: }
612: PetscCall(VecDestroy(&coords));
614: /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
615: if (test) {
616: PetscInt v, gidx;
617: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
618: for (i = 0; i < Nsubnet; i++) {
619: PetscCall(DMNetworkGetSubnetwork(networkdm, i, &nv, &ne, &vtx, &edges));
620: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, subnet[%" PetscInt_FMT "] ne %" PetscInt_FMT ", nv %" PetscInt_FMT "\n", rank, i, ne, nv));
621: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
623: for (v = 0; v < nv; v++) {
624: PetscCall(DMNetworkIsGhostVertex(networkdm, vtx[v], &ghost));
625: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx));
626: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] subnet[%" PetscInt_FMT "] v %" PetscInt_FMT " %" PetscInt_FMT "; ghost %d\n", rank, i, vtx[v], gidx, ghost));
627: }
628: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
629: }
630: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
632: PetscCall(DMNetworkGetSharedVertices(networkdm, &nv, &vtx));
633: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, num of shared vertices nsv = %" PetscInt_FMT "\n", rank, nv));
634: for (v = 0; v < nv; v++) {
635: PetscCall(DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx));
636: PetscCall(PetscPrintf(PETSC_COMM_SELF, "[%d] sv %" PetscInt_FMT ", gidx=%" PetscInt_FMT "\n", rank, vtx[v], gidx));
637: }
638: PetscCallMPI(MPI_Barrier(PETSC_COMM_WORLD));
639: }
641: /* Create solution vector X */
642: PetscCall(DMCreateGlobalVector(networkdm, &X));
643: PetscCall(VecDuplicate(X, &F));
644: PetscCall(DMGetLocalVector(networkdm, &user.localXold));
645: PetscCall(PetscLogStagePop());
647: /* (3) Setup Solvers */
648: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewJ", &viewJ, NULL));
649: PetscCall(PetscOptionsGetBool(NULL, NULL, "-viewX", &viewX, NULL));
651: PetscCall(PetscLogStageRegister("SNES Setup", &stage[2]));
652: PetscCall(PetscLogStagePush(stage[2]));
654: PetscCall(SetInitialGuess(networkdm, X, &user));
656: /* Create coupled snes */
657: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_coupled setup ......\n"));
658: user.subsnes_id = Nsubnet;
659: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
660: PetscCall(SNESSetDM(snes, networkdm));
661: PetscCall(SNESSetOptionsPrefix(snes, "coupled_"));
662: PetscCall(SNESSetFunction(snes, F, FormFunction, &user));
663: PetscCall(SNESMonitorSet(snes, UserMonitor, &user, NULL));
664: PetscCall(SNESSetFromOptions(snes));
666: if (viewJ) {
667: /* View Jac structure */
668: PetscCall(SNESGetJacobian(snes, &Jac, NULL, NULL, NULL));
669: PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
670: }
672: if (viewX) {
673: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Solution:\n"));
674: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
675: }
677: if (viewJ) {
678: /* View assembled Jac */
679: PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
680: }
682: /* Create snes_power */
683: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_power setup ......\n"));
685: user.subsnes_id = 0;
686: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_power));
687: PetscCall(SNESSetDM(snes_power, networkdm));
688: PetscCall(SNESSetOptionsPrefix(snes_power, "power_"));
689: PetscCall(SNESSetFunction(snes_power, F, FormFunction, &user));
690: PetscCall(SNESMonitorSet(snes_power, UserMonitor, &user, NULL));
692: /* Use user-provide Jacobian */
693: PetscCall(DMCreateMatrix(networkdm, &Jac));
694: PetscCall(SNESSetJacobian(snes_power, Jac, Jac, FormJacobian_subPower, &user));
695: PetscCall(SNESSetFromOptions(snes_power));
697: if (viewX) {
698: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Power Solution:\n"));
699: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
700: }
701: if (viewJ) {
702: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Power Jac:\n"));
703: PetscCall(SNESGetJacobian(snes_power, &Jac, NULL, NULL, NULL));
704: PetscCall(MatView(Jac, PETSC_VIEWER_DRAW_WORLD));
705: /* PetscCall(MatView(Jac,PETSC_VIEWER_STDOUT_WORLD)); */
706: }
708: /* Create snes_water */
709: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "SNES_water setup......\n"));
711: user.subsnes_id = 1;
712: PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes_water));
713: PetscCall(SNESSetDM(snes_water, networkdm));
714: PetscCall(SNESSetOptionsPrefix(snes_water, "water_"));
715: PetscCall(SNESSetFunction(snes_water, F, FormFunction, &user));
716: PetscCall(SNESMonitorSet(snes_water, UserMonitor, &user, NULL));
717: PetscCall(SNESSetFromOptions(snes_water));
719: if (viewX) {
720: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Water Solution:\n"));
721: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
722: }
723: PetscCall(PetscLogStagePop());
725: /* (4) Solve */
726: PetscCall(PetscLogStageRegister("SNES Solve", &stage[3]));
727: PetscCall(PetscLogStagePush(stage[3]));
728: user.it = 0;
729: reason = SNES_DIVERGED_DTOL;
730: while (user.it < it_max && (PetscInt)reason < 0) {
731: #if 0
732: user.subsnes_id = 0;
733: PetscCall(SNESSolve(snes_power,NULL,X));
735: user.subsnes_id = 1;
736: PetscCall(SNESSolve(snes_water,NULL,X));
737: #endif
738: user.subsnes_id = Nsubnet;
739: PetscCall(SNESSolve(snes, NULL, X));
741: PetscCall(SNESGetConvergedReason(snes, &reason));
742: user.it++;
743: }
744: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Coupled_SNES converged in %" PetscInt_FMT " iterations\n", user.it));
745: if (viewX) {
746: PetscCall(PetscPrintf(PETSC_COMM_WORLD, "Final Solution:\n"));
747: PetscCall(VecView(X, PETSC_VIEWER_STDOUT_WORLD));
748: }
749: PetscCall(PetscLogStagePop());
751: /* Free objects */
752: /* -------------*/
753: PetscCall(VecDestroy(&X));
754: PetscCall(VecDestroy(&F));
755: PetscCall(DMRestoreLocalVector(networkdm, &user.localXold));
757: PetscCall(SNESDestroy(&snes));
758: PetscCall(MatDestroy(&Jac));
759: PetscCall(SNESDestroy(&snes_power));
760: PetscCall(SNESDestroy(&snes_water));
762: PetscCall(DMDestroy(&networkdm));
763: PetscCall(PetscFinalize());
764: return 0;
765: }
767: /*TEST
769: build:
770: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
771: depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c
773: test:
774: args: -coupled_snes_converged_reason -options_left no -dmnetwork_view -fp_trap 0
775: localrunfiles: ex1options power/case9.m water/sample1.inp
776: output_file: output/ex1.out
778: test:
779: suffix: 2
780: nsize: 3
781: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis -fp_trap 0
782: localrunfiles: ex1options power/case9.m water/sample1.inp
783: output_file: output/ex1_2.out
784: requires: parmetis
786: test:
787: suffix: 3
788: nsize: 3
789: args: -coupled_snes_converged_reason -options_left no -distribute false -fp_trap 0
790: localrunfiles: ex1options power/case9.m water/sample1.inp
791: output_file: output/ex1_2.out
793: test:
794: suffix: 4
795: nsize: 4
796: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -dmnetwork_view -dmnetwork_view_distributed -fp_trap 0
797: localrunfiles: ex1options power/case9.m water/sample1.inp
798: output_file: output/ex1_4.out
800: test:
801: suffix: 5
802: args: -coupled_snes_converged_reason -options_left no -viewCSV -fp_trap 0
803: localrunfiles: ex1options power/case9.m water/sample1.inp
804: output_file: output/ex1_5.out
806: test:
807: suffix: 6
808: nsize: 3
809: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis -dmnetwork_view_distributed draw:null -fp_trap 0
810: localrunfiles: ex1options power/case9.m water/sample1.inp
811: output_file: output/ex1_2.out
812: requires: parmetis
814: TEST*/