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: PetscObjectGetComm((PetscObject)snes, &comm);
28: MPI_Comm_rank(comm, &rank);
29: #if 0
30: if (rank == 0) {
31: PetscInt subsnes_id = user->subsnes_id;
32: if (subsnes_id == 2) {
33: PetscPrintf(PETSC_COMM_SELF," it %" PetscInt_FMT ", subsnes_id %" PetscInt_FMT ", fnorm %g\n",user->it,user->subsnes_id,(double)fnorm);
34: } else {
35: PetscPrintf(PETSC_COMM_SELF," subsnes_id %" PetscInt_FMT ", fnorm %g\n",user->subsnes_id,(double)fnorm);
36: }
37: }
38: #endif
39: SNESGetSolution(snes, &X);
40: SNESGetDM(snes, &networkdm);
41: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localXold);
42: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localXold);
43: return 0;
44: }
46: PetscErrorCode FormJacobian_subPower(SNES snes, Vec X, Mat J, Mat Jpre, void *appctx)
47: {
48: DM networkdm;
49: Vec localX;
50: PetscInt nv, ne, i, j, offset, nvar, row;
51: const PetscInt *vtx, *edges;
52: PetscBool ghostvtex;
53: PetscScalar one = 1.0;
54: PetscMPIInt rank;
55: MPI_Comm comm;
57: SNESGetDM(snes, &networkdm);
58: DMGetLocalVector(networkdm, &localX);
60: PetscObjectGetComm((PetscObject)networkdm, &comm);
61: MPI_Comm_rank(comm, &rank);
63: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
64: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
66: MatZeroEntries(J);
68: /* Power subnetwork: copied from power/FormJacobian_Power() */
69: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
70: FormJacobian_Power_private(networkdm, localX, J, nv, ne, vtx, edges, appctx);
72: /* Water subnetwork: Identity */
73: DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges);
74: for (i = 0; i < nv; i++) {
75: DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex);
76: if (ghostvtex) continue;
78: DMNetworkGetGlobalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset);
79: DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar);
80: for (j = 0; j < nvar; j++) {
81: row = offset + j;
82: MatSetValues(J, 1, &row, 1, &row, &one, ADD_VALUES);
83: }
84: }
85: MatAssemblyBegin(J, MAT_FINAL_ASSEMBLY);
86: MatAssemblyEnd(J, MAT_FINAL_ASSEMBLY);
88: DMRestoreLocalVector(networkdm, &localX);
89: return 0;
90: }
92: /* Dummy equation localF(X) = localX - localXold */
93: PetscErrorCode FormFunction_Dummy(DM networkdm, Vec localX, Vec localF, PetscInt nv, PetscInt ne, const PetscInt *vtx, const PetscInt *edges, void *appctx)
94: {
95: const PetscScalar *xarr, *xoldarr;
96: PetscScalar *farr;
97: PetscInt i, j, offset, nvar;
98: PetscBool ghostvtex;
99: UserCtx *user = (UserCtx *)appctx;
100: Vec localXold = user->localXold;
102: VecGetArrayRead(localX, &xarr);
103: VecGetArrayRead(localXold, &xoldarr);
104: VecGetArray(localF, &farr);
106: for (i = 0; i < nv; i++) {
107: DMNetworkIsGhostVertex(networkdm, vtx[i], &ghostvtex);
108: if (ghostvtex) continue;
110: DMNetworkGetLocalVecOffset(networkdm, vtx[i], ALL_COMPONENTS, &offset);
111: DMNetworkGetComponent(networkdm, vtx[i], ALL_COMPONENTS, NULL, NULL, &nvar);
112: for (j = 0; j < nvar; j++) farr[offset + j] = xarr[offset + j] - xoldarr[offset + j];
113: }
115: VecRestoreArrayRead(localX, &xarr);
116: VecRestoreArrayRead(localXold, &xoldarr);
117: VecRestoreArray(localF, &farr);
118: return 0;
119: }
121: PetscErrorCode FormFunction(SNES snes, Vec X, Vec F, void *appctx)
122: {
123: DM networkdm;
124: Vec localX, localF;
125: PetscInt nv, ne, v;
126: const PetscInt *vtx, *edges;
127: PetscMPIInt rank;
128: MPI_Comm comm;
129: UserCtx *user = (UserCtx *)appctx;
130: UserCtx_Power appctx_power = (*user).appctx_power;
131: AppCtx_Water appctx_water = (*user).appctx_water;
133: SNESGetDM(snes, &networkdm);
134: PetscObjectGetComm((PetscObject)networkdm, &comm);
135: MPI_Comm_rank(comm, &rank);
137: DMGetLocalVector(networkdm, &localX);
138: DMGetLocalVector(networkdm, &localF);
139: VecSet(F, 0.0);
140: VecSet(localF, 0.0);
142: DMGlobalToLocalBegin(networkdm, X, INSERT_VALUES, localX);
143: DMGlobalToLocalEnd(networkdm, X, INSERT_VALUES, localX);
145: /* Form Function for power subnetwork */
146: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
147: if (user->subsnes_id == 1) { /* snes_water only */
148: FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user);
149: } else {
150: FormFunction_Power(networkdm, localX, localF, nv, ne, vtx, edges, &appctx_power);
151: }
153: /* Form Function for water subnetwork */
154: DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges);
155: if (user->subsnes_id == 0) { /* snes_power only */
156: FormFunction_Dummy(networkdm, localX, localF, nv, ne, vtx, edges, user);
157: } else {
158: FormFunction_Water(networkdm, localX, localF, nv, ne, vtx, edges, NULL);
159: }
161: /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */
162: DMNetworkGetSharedVertices(networkdm, &nv, &vtx);
163: for (v = 0; v < nv; v++) {
164: PetscInt key, ncomp, nvar, nconnedges, k, e, keye, goffset[3];
165: void *component;
166: const PetscInt *connedges;
168: DMNetworkGetComponent(networkdm, vtx[v], ALL_COMPONENTS, NULL, NULL, &nvar);
169: DMNetworkGetNumComponents(networkdm, vtx[v], &ncomp);
170: /* printf(" [%d] coupling vertex[%" PetscInt_FMT "]: v %" PetscInt_FMT ", ncomp %" PetscInt_FMT "; nvar %" PetscInt_FMT "\n",rank,v,vtx[v], ncomp,nvar); */
172: for (k = 0; k < ncomp; k++) {
173: DMNetworkGetComponent(networkdm, vtx[v], k, &key, &component, &nvar);
174: DMNetworkGetGlobalVecOffset(networkdm, vtx[v], k, &goffset[k]);
176: /* Verify the coupling vertex is a powernet load vertex or a water vertex */
177: switch (k) {
178: case 0:
180: break;
181: case 1:
183: break;
184: case 2:
186: break;
187: default:
188: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %" PetscInt_FMT " is wrong", k);
189: }
190: /* printf(" [%d] coupling vertex[%" PetscInt_FMT "]: key %" PetscInt_FMT "; nvar %" PetscInt_FMT ", goffset %" PetscInt_FMT "\n",rank,v,key,nvar,goffset[k]); */
191: }
193: /* Get its supporting edges */
194: DMNetworkGetSupportingEdges(networkdm, vtx[v], &nconnedges, &connedges);
195: /* printf("\n[%d] coupling vertex: nconnedges %" PetscInt_FMT "\n",rank,nconnedges); */
196: for (k = 0; k < nconnedges; k++) {
197: e = connedges[k];
198: DMNetworkGetNumComponents(networkdm, e, &ncomp);
199: /* printf("\n [%d] connected edge[%" PetscInt_FMT "]=%" PetscInt_FMT " has ncomp %" PetscInt_FMT "\n",rank,k,e,ncomp); */
200: DMNetworkGetComponent(networkdm, e, 0, &keye, &component, NULL);
202: }
203: }
205: DMRestoreLocalVector(networkdm, &localX);
207: DMLocalToGlobalBegin(networkdm, localF, ADD_VALUES, F);
208: DMLocalToGlobalEnd(networkdm, localF, ADD_VALUES, F);
209: DMRestoreLocalVector(networkdm, &localF);
210: #if 0
211: if (rank == 0) printf("F:\n");
212: VecView(F,PETSC_VIEWER_STDOUT_WORLD);
213: #endif
214: return 0;
215: }
217: PetscErrorCode SetInitialGuess(DM networkdm, Vec X, void *appctx)
218: {
219: PetscInt nv, ne, i, j, ncomp, offset, key;
220: const PetscInt *vtx, *edges;
221: UserCtx *user = (UserCtx *)appctx;
222: Vec localX = user->localXold;
223: UserCtx_Power appctx_power = (*user).appctx_power;
224: AppCtx_Water appctx_water = (*user).appctx_water;
225: PetscBool ghost;
226: PetscScalar *xarr;
227: VERTEX_Power bus;
228: VERTEX_Water vertex;
229: void *component;
230: GEN gen;
232: VecSet(X, 0.0);
233: VecSet(localX, 0.0);
235: /* Set initial guess for power subnetwork */
236: DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges);
237: SetInitialGuess_Power(networkdm, localX, nv, ne, vtx, edges, &appctx_power);
239: /* Set initial guess for water subnetwork */
240: DMNetworkGetSubnetwork(networkdm, 1, &nv, &ne, &vtx, &edges);
241: SetInitialGuess_Water(networkdm, localX, nv, ne, vtx, edges, NULL);
243: /* Set initial guess at the coupling vertex */
244: VecGetArray(localX, &xarr);
245: DMNetworkGetSharedVertices(networkdm, &nv, &vtx);
246: for (i = 0; i < nv; i++) {
247: DMNetworkIsGhostVertex(networkdm, vtx[i], &ghost);
248: if (ghost) continue;
250: DMNetworkGetNumComponents(networkdm, vtx[i], &ncomp);
251: for (j = 0; j < ncomp; j++) {
252: DMNetworkGetLocalVecOffset(networkdm, vtx[i], j, &offset);
253: DMNetworkGetComponent(networkdm, vtx[i], j, &key, (void **)&component, NULL);
254: if (key == appctx_power.compkey_bus) {
255: bus = (VERTEX_Power)(component);
256: xarr[offset] = bus->va * PETSC_PI / 180.0;
257: xarr[offset + 1] = bus->vm;
258: } else if (key == appctx_power.compkey_gen) {
259: gen = (GEN)(component);
260: if (!gen->status) continue;
261: xarr[offset + 1] = gen->vs;
262: } else if (key == appctx_water.compkey_vtx) {
263: vertex = (VERTEX_Water)(component);
264: if (vertex->type == VERTEX_TYPE_JUNCTION) {
265: xarr[offset] = 100;
266: } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
267: xarr[offset] = vertex->res.head;
268: } else {
269: xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
270: }
271: }
272: }
273: }
274: VecRestoreArray(localX, &xarr);
276: DMLocalToGlobalBegin(networkdm, localX, ADD_VALUES, X);
277: DMLocalToGlobalEnd(networkdm, localX, ADD_VALUES, X);
278: return 0;
279: }
281: int main(int argc, char **argv)
282: {
283: DM networkdm;
284: PetscLogStage stage[4];
285: PetscMPIInt rank, size;
286: PetscInt Nsubnet = 2, numVertices[2], numEdges[2], i, j, nv, ne, it_max = 10;
287: const PetscInt *vtx, *edges;
288: Vec X, F;
289: SNES snes, snes_power, snes_water;
290: Mat Jac;
291: PetscBool ghost, viewJ = PETSC_FALSE, viewX = PETSC_FALSE, viewDM = PETSC_FALSE, test = PETSC_FALSE, distribute = PETSC_TRUE, flg;
292: UserCtx user;
293: SNESConvergedReason reason;
295: /* Power subnetwork */
296: UserCtx_Power *appctx_power = &user.appctx_power;
297: char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
298: PFDATA *pfdata = NULL;
299: PetscInt genj, loadj, *edgelist_power = NULL, power_netnum;
300: PetscScalar Sbase = 0.0;
302: /* Water subnetwork */
303: AppCtx_Water *appctx_water = &user.appctx_water;
304: WATERDATA *waterdata = NULL;
305: char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
306: PetscInt *edgelist_water = NULL, water_netnum;
308: /* Shared vertices between subnetworks */
309: PetscInt power_svtx, water_svtx;
312: PetscInitialize(&argc, &argv, "ex1options", help);
313: MPI_Comm_rank(PETSC_COMM_WORLD, &rank);
314: MPI_Comm_size(PETSC_COMM_WORLD, &size);
316: /* (1) Read Data - Only rank 0 reads the data */
317: PetscLogStageRegister("Read Data", &stage[0]);
318: PetscLogStagePush(stage[0]);
320: for (i = 0; i < Nsubnet; i++) {
321: numVertices[i] = 0;
322: numEdges[i] = 0;
323: }
325: /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
326: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
327: PetscOptionsGetString(NULL, NULL, "-pfdata", pfdata_file, PETSC_MAX_PATH_LEN - 1, NULL);
328: PetscNew(&pfdata);
329: PFReadMatPowerData(pfdata, pfdata_file);
330: if (rank == 0) 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);
331: Sbase = pfdata->sbase;
332: if (rank == 0) { /* proc[0] will create Electric Power Grid */
333: numEdges[0] = pfdata->nbranch;
334: numVertices[0] = pfdata->nbus;
336: PetscMalloc1(2 * numEdges[0], &edgelist_power);
337: GetListofEdges_Power(pfdata, edgelist_power);
338: }
339: /* Broadcast power Sbase to all processors */
340: MPI_Bcast(&Sbase, 1, MPIU_SCALAR, 0, PETSC_COMM_WORLD);
341: appctx_power->Sbase = Sbase;
342: appctx_power->jac_error = PETSC_FALSE;
343: /* If external option activated. Introduce error in jacobian */
344: PetscOptionsHasName(NULL, NULL, "-jac_error", &appctx_power->jac_error);
346: /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
347: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
348: PetscNew(&waterdata);
349: PetscOptionsGetString(NULL, NULL, "-waterdata", waterdata_file, PETSC_MAX_PATH_LEN - 1, NULL);
350: WaterReadData(waterdata, waterdata_file);
351: if (size == 1 || (size > 1 && rank == 1)) {
352: PetscCalloc1(2 * waterdata->nedge, &edgelist_water);
353: GetListofEdges_Water(waterdata, edgelist_water);
354: numEdges[1] = waterdata->nedge;
355: numVertices[1] = waterdata->nvertex;
356: }
357: PetscLogStagePop();
359: /* (2) Create a network consist of two subnetworks */
360: PetscLogStageRegister("Net Setup", &stage[1]);
361: PetscLogStagePush(stage[1]);
363: PetscOptionsGetBool(NULL, NULL, "-viewDM", &viewDM, NULL);
364: PetscOptionsGetBool(NULL, NULL, "-test", &test, NULL);
365: PetscOptionsGetBool(NULL, NULL, "-distribute", &distribute, NULL);
367: /* Create an empty network object */
368: DMNetworkCreate(PETSC_COMM_WORLD, &networkdm);
370: /* Register the components in the network */
371: DMNetworkRegisterComponent(networkdm, "branchstruct", sizeof(struct _p_EDGE_Power), &appctx_power->compkey_branch);
372: DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Power), &appctx_power->compkey_bus);
373: DMNetworkRegisterComponent(networkdm, "genstruct", sizeof(struct _p_GEN), &appctx_power->compkey_gen);
374: DMNetworkRegisterComponent(networkdm, "loadstruct", sizeof(struct _p_LOAD), &appctx_power->compkey_load);
376: DMNetworkRegisterComponent(networkdm, "edge_water", sizeof(struct _p_EDGE_Water), &appctx_water->compkey_edge);
377: DMNetworkRegisterComponent(networkdm, "vertex_water", sizeof(struct _p_VERTEX_Water), &appctx_water->compkey_vtx);
378: #if 0
379: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_branch %d\n",appctx_power->compkey_branch);
380: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_bus %d\n",appctx_power->compkey_bus);
381: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_gen %d\n",appctx_power->compkey_gen);
382: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_load %d\n",appctx_power->compkey_load);
383: PetscPrintf(PETSC_COMM_WORLD,"water->compkey_edge %d\n",appctx_water->compkey_edge);
384: PetscPrintf(PETSC_COMM_WORLD,"water->compkey_vtx %d\n",appctx_water->compkey_vtx);
385: #endif
386: 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]);
387: PetscSynchronizedFlush(PETSC_COMM_WORLD, PETSC_STDOUT);
389: DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, Nsubnet);
390: DMNetworkAddSubnetwork(networkdm, "power", numEdges[0], edgelist_power, &power_netnum);
391: DMNetworkAddSubnetwork(networkdm, "water", numEdges[1], edgelist_water, &water_netnum);
393: /* vertex subnet[0].4 shares with vertex subnet[1].0 */
394: power_svtx = 4;
395: water_svtx = 0;
396: DMNetworkAddSharedVertices(networkdm, power_netnum, water_netnum, 1, &power_svtx, &water_svtx);
398: /* Set up the network layout */
399: DMNetworkLayoutSetUp(networkdm);
401: /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
402: genj = 0;
403: loadj = 0;
404: DMNetworkGetSubnetwork(networkdm, power_netnum, &nv, &ne, &vtx, &edges);
406: for (i = 0; i < ne; i++) DMNetworkAddComponent(networkdm, edges[i], appctx_power->compkey_branch, &pfdata->branch[i], 0);
408: for (i = 0; i < nv; i++) {
409: DMNetworkIsSharedVertex(networkdm, vtx[i], &flg);
410: if (flg) continue;
412: DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[i], 2);
413: if (pfdata->bus[i].ngen) {
414: for (j = 0; j < pfdata->bus[i].ngen; j++) DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_gen, &pfdata->gen[genj++], 0);
415: }
416: if (pfdata->bus[i].nload) {
417: for (j = 0; j < pfdata->bus[i].nload; j++) DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[loadj++], 0);
418: }
419: }
421: /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
422: DMNetworkGetSubnetwork(networkdm, water_netnum, &nv, &ne, &vtx, &edges);
423: for (i = 0; i < ne; i++) DMNetworkAddComponent(networkdm, edges[i], appctx_water->compkey_edge, &waterdata->edge[i], 0);
425: for (i = 0; i < nv; i++) {
426: DMNetworkIsSharedVertex(networkdm, vtx[i], &flg);
427: if (flg) continue;
429: DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[i], 1);
430: }
432: /* 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 */
433: DMNetworkGetSharedVertices(networkdm, &nv, &vtx);
434: for (i = 0; i < nv; i++) {
435: /* power */
436: DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_bus, &pfdata->bus[4], 2);
437: /* bus[4] is a load, add its component */
438: DMNetworkAddComponent(networkdm, vtx[i], appctx_power->compkey_load, &pfdata->load[0], 0);
440: /* water */
441: DMNetworkAddComponent(networkdm, vtx[i], appctx_water->compkey_vtx, &waterdata->vertex[0], 1);
442: }
444: /* Set up DM for use */
445: DMSetUp(networkdm);
446: if (viewDM) {
447: PetscPrintf(PETSC_COMM_WORLD, "\nAfter DMSetUp, DMView:\n");
448: DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD);
449: }
451: /* Free user objects */
452: PetscFree(edgelist_power);
453: PetscFree(pfdata->bus);
454: PetscFree(pfdata->gen);
455: PetscFree(pfdata->branch);
456: PetscFree(pfdata->load);
457: PetscFree(pfdata);
459: PetscFree(edgelist_water);
460: PetscFree(waterdata->vertex);
461: PetscFree(waterdata->edge);
462: PetscFree(waterdata);
464: /* Re-distribute networkdm to multiple processes for better job balance */
465: if (size > 1 && distribute) {
466: DMNetworkDistribute(&networkdm, 0);
467: if (viewDM) {
468: PetscPrintf(PETSC_COMM_WORLD, "\nAfter DMNetworkDistribute, DMView:\n");
469: DMView(networkdm, PETSC_VIEWER_STDOUT_WORLD);
470: }
471: }
473: /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
474: if (test) {
475: PetscInt v, gidx;
476: MPI_Barrier(PETSC_COMM_WORLD);
477: for (i = 0; i < Nsubnet; i++) {
478: DMNetworkGetSubnetwork(networkdm, i, &nv, &ne, &vtx, &edges);
479: PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, subnet[%" PetscInt_FMT "] ne %" PetscInt_FMT ", nv %" PetscInt_FMT "\n", rank, i, ne, nv);
480: MPI_Barrier(PETSC_COMM_WORLD);
482: for (v = 0; v < nv; v++) {
483: DMNetworkIsGhostVertex(networkdm, vtx[v], &ghost);
484: DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx);
485: PetscPrintf(PETSC_COMM_SELF, "[%d] subnet[%" PetscInt_FMT "] v %" PetscInt_FMT " %" PetscInt_FMT "; ghost %d\n", rank, i, vtx[v], gidx, ghost);
486: }
487: MPI_Barrier(PETSC_COMM_WORLD);
488: }
489: MPI_Barrier(PETSC_COMM_WORLD);
491: DMNetworkGetSharedVertices(networkdm, &nv, &vtx);
492: PetscPrintf(PETSC_COMM_SELF, "[%d] After distribute, num of shared vertices nsv = %" PetscInt_FMT "\n", rank, nv);
493: for (v = 0; v < nv; v++) {
494: DMNetworkGetGlobalVertexIndex(networkdm, vtx[v], &gidx);
495: PetscPrintf(PETSC_COMM_SELF, "[%d] sv %" PetscInt_FMT ", gidx=%" PetscInt_FMT "\n", rank, vtx[v], gidx);
496: }
497: MPI_Barrier(PETSC_COMM_WORLD);
498: }
500: /* Create solution vector X */
501: DMCreateGlobalVector(networkdm, &X);
502: VecDuplicate(X, &F);
503: DMGetLocalVector(networkdm, &user.localXold);
504: PetscLogStagePop();
506: /* (3) Setup Solvers */
507: PetscOptionsGetBool(NULL, NULL, "-viewJ", &viewJ, NULL);
508: PetscOptionsGetBool(NULL, NULL, "-viewX", &viewX, NULL);
510: PetscLogStageRegister("SNES Setup", &stage[2]);
511: PetscLogStagePush(stage[2]);
513: SetInitialGuess(networkdm, X, &user);
515: /* Create coupled snes */
516: PetscPrintf(PETSC_COMM_WORLD, "SNES_coupled setup ......\n");
517: user.subsnes_id = Nsubnet;
518: SNESCreate(PETSC_COMM_WORLD, &snes);
519: SNESSetDM(snes, networkdm);
520: SNESSetOptionsPrefix(snes, "coupled_");
521: SNESSetFunction(snes, F, FormFunction, &user);
522: SNESMonitorSet(snes, UserMonitor, &user, NULL);
523: SNESSetFromOptions(snes);
525: if (viewJ) {
526: /* View Jac structure */
527: SNESGetJacobian(snes, &Jac, NULL, NULL, NULL);
528: MatView(Jac, PETSC_VIEWER_DRAW_WORLD);
529: }
531: if (viewX) {
532: PetscPrintf(PETSC_COMM_WORLD, "Solution:\n");
533: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
534: }
536: if (viewJ) {
537: /* View assembled Jac */
538: MatView(Jac, PETSC_VIEWER_DRAW_WORLD);
539: }
541: /* Create snes_power */
542: PetscPrintf(PETSC_COMM_WORLD, "SNES_power setup ......\n");
544: user.subsnes_id = 0;
545: SNESCreate(PETSC_COMM_WORLD, &snes_power);
546: SNESSetDM(snes_power, networkdm);
547: SNESSetOptionsPrefix(snes_power, "power_");
548: SNESSetFunction(snes_power, F, FormFunction, &user);
549: SNESMonitorSet(snes_power, UserMonitor, &user, NULL);
551: /* Use user-provide Jacobian */
552: DMCreateMatrix(networkdm, &Jac);
553: SNESSetJacobian(snes_power, Jac, Jac, FormJacobian_subPower, &user);
554: SNESSetFromOptions(snes_power);
556: if (viewX) {
557: PetscPrintf(PETSC_COMM_WORLD, "Power Solution:\n");
558: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
559: }
560: if (viewJ) {
561: PetscPrintf(PETSC_COMM_WORLD, "Power Jac:\n");
562: SNESGetJacobian(snes_power, &Jac, NULL, NULL, NULL);
563: MatView(Jac, PETSC_VIEWER_DRAW_WORLD);
564: /* MatView(Jac,PETSC_VIEWER_STDOUT_WORLD); */
565: }
567: /* Create snes_water */
568: PetscPrintf(PETSC_COMM_WORLD, "SNES_water setup......\n");
570: user.subsnes_id = 1;
571: SNESCreate(PETSC_COMM_WORLD, &snes_water);
572: SNESSetDM(snes_water, networkdm);
573: SNESSetOptionsPrefix(snes_water, "water_");
574: SNESSetFunction(snes_water, F, FormFunction, &user);
575: SNESMonitorSet(snes_water, UserMonitor, &user, NULL);
576: SNESSetFromOptions(snes_water);
578: if (viewX) {
579: PetscPrintf(PETSC_COMM_WORLD, "Water Solution:\n");
580: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
581: }
582: PetscLogStagePop();
584: /* (4) Solve */
585: PetscLogStageRegister("SNES Solve", &stage[3]);
586: PetscLogStagePush(stage[3]);
587: user.it = 0;
588: reason = SNES_DIVERGED_DTOL;
589: while (user.it < it_max && (PetscInt)reason < 0) {
590: #if 0
591: user.subsnes_id = 0;
592: SNESSolve(snes_power,NULL,X);
594: user.subsnes_id = 1;
595: SNESSolve(snes_water,NULL,X);
596: #endif
597: user.subsnes_id = Nsubnet;
598: SNESSolve(snes, NULL, X);
600: SNESGetConvergedReason(snes, &reason);
601: user.it++;
602: }
603: PetscPrintf(PETSC_COMM_WORLD, "Coupled_SNES converged in %" PetscInt_FMT " iterations\n", user.it);
604: if (viewX) {
605: PetscPrintf(PETSC_COMM_WORLD, "Final Solution:\n");
606: VecView(X, PETSC_VIEWER_STDOUT_WORLD);
607: }
608: PetscLogStagePop();
610: /* Free objects */
611: /* -------------*/
612: VecDestroy(&X);
613: VecDestroy(&F);
614: DMRestoreLocalVector(networkdm, &user.localXold);
616: SNESDestroy(&snes);
617: MatDestroy(&Jac);
618: SNESDestroy(&snes_power);
619: SNESDestroy(&snes_water);
621: DMDestroy(&networkdm);
622: PetscFinalize();
623: return 0;
624: }
626: /*TEST
628: build:
629: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
630: depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c
632: test:
633: args: -coupled_snes_converged_reason -options_left no -viewDM
634: localrunfiles: ex1options power/case9.m water/sample1.inp
635: output_file: output/ex1.out
637: test:
638: suffix: 2
639: nsize: 3
640: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis
641: localrunfiles: ex1options power/case9.m water/sample1.inp
642: output_file: output/ex1_2.out
643: requires: parmetis
645: # test:
646: # suffix: 3
647: # nsize: 3
648: # args: -coupled_snes_converged_reason -options_left no -distribute false
649: # localrunfiles: ex1options power/case9.m water/sample1.inp
650: # output_file: output/ex1_2.out
652: test:
653: suffix: 4
654: nsize: 4
655: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -viewDM
656: localrunfiles: ex1options power/case9.m water/sample1.inp
657: output_file: output/ex1_4.out
659: TEST*/