Actual source code: ex1.c
petsc-3.14.4 2021-02-03
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: /* T
9: Concepts: DMNetwork
10: Concepts: PETSc SNES solver
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: } UserCtx;
24: PetscErrorCode UserMonitor(SNES snes,PetscInt its,PetscReal fnorm ,void *appctx)
25: {
27: UserCtx *user = (UserCtx*)appctx;
28: Vec X,localXold = user->localXold;
29: DM networkdm;
30: PetscMPIInt rank;
31: MPI_Comm comm;
34: PetscObjectGetComm((PetscObject)snes,&comm);
35: MPI_Comm_rank(comm,&rank);
36: #if 0
37: if (!rank) {
38: PetscInt subsnes_id = user->subsnes_id;
39: if (subsnes_id == 2) {
40: PetscPrintf(PETSC_COMM_SELF," it %D, subsnes_id %D, fnorm %g\n",user->it,user->subsnes_id,(double)fnorm);
41: } else {
42: PetscPrintf(PETSC_COMM_SELF," subsnes_id %D, fnorm %g\n",user->subsnes_id,(double)fnorm);
43: }
44: }
45: #endif
46: SNESGetSolution(snes,&X);
47: SNESGetDM(snes,&networkdm);
48: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localXold);
49: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localXold);
50: return(0);
51: }
53: PetscErrorCode FormJacobian_subPower(SNES snes,Vec X, Mat J,Mat Jpre,void *appctx)
54: {
56: DM networkdm;
57: Vec localX;
58: PetscInt nv,ne,i,j,offset,nvar,row;
59: const PetscInt *vtx,*edges;
60: PetscBool ghostvtex;
61: PetscScalar one = 1.0;
62: PetscMPIInt rank;
63: MPI_Comm comm;
66: SNESGetDM(snes,&networkdm);
67: DMGetLocalVector(networkdm,&localX);
69: PetscObjectGetComm((PetscObject)networkdm,&comm);
70: MPI_Comm_rank(comm,&rank);
72: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
73: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
75: MatZeroEntries(J);
77: /* Power subnetwork: copied from power/FormJacobian_Power() */
78: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
79: FormJacobian_Power_private(networkdm,localX,J,nv,ne,vtx,edges,appctx);
81: /* Water subnetwork: Identity */
82: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
83: for (i=0; i<nv; i++) {
84: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
85: if (ghostvtex) continue;
87: DMNetworkGetVariableGlobalOffset(networkdm,vtx[i],&offset);
88: DMNetworkGetNumVariables(networkdm,vtx[i],&nvar);
89: for (j=0; j<nvar; j++) {
90: row = offset + j;
91: MatSetValues(J,1,&row,1,&row,&one,ADD_VALUES);
92: }
93: }
94: MatAssemblyBegin(J,MAT_FINAL_ASSEMBLY);
95: MatAssemblyEnd(J,MAT_FINAL_ASSEMBLY);
97: DMRestoreLocalVector(networkdm,&localX);
98: return(0);
99: }
101: /* Dummy equation localF(X) = localX - localXold */
102: PetscErrorCode FormFunction_Dummy(DM networkdm,Vec localX, Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx)
103: {
104: PetscErrorCode ierr;
105: const PetscScalar *xarr,*xoldarr;
106: PetscScalar *farr;
107: PetscInt i,j,offset,nvar;
108: PetscBool ghostvtex;
109: UserCtx *user = (UserCtx*)appctx;
110: Vec localXold = user->localXold;
113: VecGetArrayRead(localX,&xarr);
114: VecGetArrayRead(localXold,&xoldarr);
115: VecGetArray(localF,&farr);
117: for (i=0; i<nv; i++) {
118: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
119: if (ghostvtex) continue;
121: DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);
122: DMNetworkGetNumVariables(networkdm,vtx[i],&nvar);
123: for (j=0; j<nvar; j++) {
124: farr[offset+j] = xarr[offset+j] - xoldarr[offset+j];
125: }
126: }
128: VecRestoreArrayRead(localX,&xarr);
129: VecRestoreArrayRead(localXold,&xoldarr);
130: VecRestoreArray(localF,&farr);
131: return(0);
132: }
134: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *appctx)
135: {
137: DM networkdm;
138: Vec localX,localF;
139: PetscInt nv,ne;
140: const PetscInt *vtx,*edges;
141: PetscMPIInt rank;
142: MPI_Comm comm;
143: UserCtx *user = (UserCtx*)appctx;
144: UserCtx_Power appctx_power = (*user).appctx_power;
147: SNESGetDM(snes,&networkdm);
148: PetscObjectGetComm((PetscObject)networkdm,&comm);
149: MPI_Comm_rank(comm,&rank);
151: DMGetLocalVector(networkdm,&localX);
152: DMGetLocalVector(networkdm,&localF);
153: VecSet(F,0.0);
154: VecSet(localF,0.0);
156: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
157: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
159: /* Form Function for power subnetwork */
160: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
161: if (user->subsnes_id == 1) { /* snes_water only */
162: FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
163: } else {
164: FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,&appctx_power);
165: }
167: /* Form Function for water subnetwork */
168: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
169: if (user->subsnes_id == 0) { /* snes_power only */
170: FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
171: } else {
172: FormFunction_Water(networkdm,localX,localF,nv,ne,vtx,edges,NULL);
173: }
175: #if 0
176: /* Form Function for the coupling subnetwork -- experimental, not done yet */
177: DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);
178: if (ne) {
179: const PetscInt *cone,*connedges,*econe;
180: PetscInt key,vid[2],nconnedges,k,e,keye;
181: void* component;
182: AppCtx_Water appctx_water = (*user).appctx_water;
184: DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);
186: DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);
187: DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);
188: PetscPrintf(PETSC_COMM_SELF,"[%d] Formfunction, coupling subnetwork: ne %d; connected vertices %d %d\n",rank,ne,vid[0],vid[1]);
190: /* Get coupling powernet load vertex */
191: DMNetworkGetComponent(networkdm,cone[0],1,&key,&component);
192: if (key != appctx_power.compkey_load) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power load vertex");
194: /* Get coupling water vertex and pump edge */
195: DMNetworkGetComponent(networkdm,cone[1],0,&key,&component);
196: if (key != appctx_water.compkey_vtx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex");
198: /* get its supporting edges */
199: DMNetworkGetSupportingEdges(networkdm,cone[1],&nconnedges,&connedges);
201: for (k=0; k<nconnedges; k++) {
202: e = connedges[k];
203: DMNetworkGetComponent(networkdm,e,0,&keye,&component);
205: if (keye == appctx_water.compkey_edge) { /* water_compkey_edge */
206: EDGE_Water edge=(EDGE_Water)component;
207: if (edge->type == EDGE_TYPE_PUMP) {
208: /* compute flow_pump */
209: PetscInt offsetnode1,offsetnode2,key_0,key_1;
210: const PetscScalar *xarr;
211: PetscScalar *farr;
212: VERTEX_Water vertexnode1,vertexnode2;
214: DMNetworkGetConnectedVertices(networkdm,e,&econe);
215: DMNetworkGetGlobalVertexIndex(networkdm,econe[0],&vid[0]);
216: DMNetworkGetGlobalVertexIndex(networkdm,econe[1],&vid[1]);
218: VecGetArray(localF,&farr);
219: DMNetworkGetVariableOffset(networkdm,econe[0],&offsetnode1);
220: DMNetworkGetVariableOffset(networkdm,econe[1],&offsetnode2);
222: VecGetArrayRead(localX,&xarr);
223: #if 0
224: PetscScalar hf,ht;
225: Pump *pump;
226: pump = &edge->pump;
227: hf = xarr[offsetnode1];
228: ht = xarr[offsetnode2];
230: PetscScalar flow = Flow_Pump(pump,hf,ht);
231: PetscScalar Hp = 0.1; /* load->pl */
232: PetscScalar flow_couple = 8.81*Hp*1.e6/(ht-hf); /* pump->h0; */
233: /* PetscPrintf(PETSC_COMM_SELF,"pump %d: connected vtx %d %d; flow_pump %g flow_couple %g; offset %d %d\n",e,vid[0],vid[1],flow,flow_couple,offsetnode1,offsetnode2); */
234: #endif
235: /* Get the components at the two vertices */
236: DMNetworkGetComponent(networkdm,econe[0],0,&key_0,(void**)&vertexnode1);
237: if (key_0 != appctx_water.compkey_vtx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex");
238: DMNetworkGetComponent(networkdm,econe[1],0,&key_1,(void**)&vertexnode2);
239: if (key_1 != appctx_water.compkey_vtx) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex");
240: #if 0
241: /* subtract flow_pump computed in FormFunction_Water() and add flow_couple to connected nodes */
242: if (vertexnode1->type == VERTEX_TYPE_JUNCTION) {
243: farr[offsetnode1] += flow;
244: farr[offsetnode1] -= flow_couple;
245: }
246: if (vertexnode2->type == VERTEX_TYPE_JUNCTION) {
247: farr[offsetnode2] -= flow;
248: farr[offsetnode2] += flow_couple;
249: }
250: #endif
251: VecRestoreArrayRead(localX,&xarr);
252: VecRestoreArray(localF,&farr);
253: }
254: }
255: }
256: }
257: #endif
259: DMRestoreLocalVector(networkdm,&localX);
261: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
262: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
263: DMRestoreLocalVector(networkdm,&localF);
264: /* VecView(F,PETSC_VIEWER_STDOUT_WORLD); */
265: return(0);
266: }
268: PetscErrorCode SetInitialGuess(DM networkdm,Vec X,void* appctx)
269: {
271: PetscInt nv,ne;
272: const PetscInt *vtx,*edges;
273: UserCtx *user = (UserCtx*)appctx;
274: Vec localX = user->localXold;
275: UserCtx_Power appctx_power = (*user).appctx_power;
278: VecSet(X,0.0);
279: VecSet(localX,0.0);
281: /* Set initial guess for power subnetwork */
282: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
283: SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,&appctx_power);
285: /* Set initial guess for water subnetwork */
286: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
287: SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);
289: #if 0
290: /* Set initial guess for the coupling subnet */
291: DMNetworkGetSubnetworkInfo(networkdm,2,&nv,&ne,&vtx,&edges);
292: if (ne) {
293: const PetscInt *cone;
294: DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);
295: }
296: #endif
298: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
299: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
300: return(0);
301: }
303: int main(int argc,char **argv)
304: {
305: PetscErrorCode ierr;
306: DM networkdm;
307: #if defined(PETSC_USE_LOG)
308: PetscLogStage stage[4];
309: #endif
310: PetscMPIInt rank,size;
311: PetscInt nsubnet = 2,nsubnetCouple = 0,numVertices[2],numEdges[2],numEdgesCouple[1];
312: PetscInt i,j,nv,ne;
313: PetscInt *edgelist[2];
314: const PetscInt *vtx,*edges;
315: Vec X,F;
316: SNES snes,snes_power,snes_water;
317: Mat Jac;
318: PetscBool viewJ = PETSC_FALSE,viewX = PETSC_FALSE,viewDM = PETSC_FALSE,test = PETSC_FALSE,distribute = PETSC_TRUE;
319: UserCtx user;
320: PetscInt it_max = 10;
321: SNESConvergedReason reason;
323: /* Power subnetwork */
324: UserCtx_Power *appctx_power = &user.appctx_power;
325: char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
326: PFDATA *pfdata = NULL;
327: PetscInt genj,loadj;
328: PetscInt *edgelist_power = NULL;
329: PetscScalar Sbase = 0.0;
331: /* Water subnetwork */
332: AppCtx_Water *appctx_water = &user.appctx_water;
333: WATERDATA *waterdata = NULL;
334: char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
335: PetscInt *edgelist_water = NULL;
337: /* Coupling subnetwork */
338: PetscInt *edgelist_couple = NULL;
340: PetscInitialize(&argc,&argv,"ex1options",help);if (ierr) return ierr;
341: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
342: MPI_Comm_size(PETSC_COMM_WORLD,&size);
344: /* (1) Read Data - Only rank 0 reads the data */
345: /*--------------------------------------------*/
346: PetscLogStageRegister("Read Data",&stage[0]);
347: PetscLogStagePush(stage[0]);
349: for (i=0; i<nsubnet; i++) {
350: numVertices[i] = 0;
351: numEdges[i] = 0;
352: }
353: numEdgesCouple[0] = 0;
355: /* proc[0] READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
356: if (rank == 0) {
357: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,sizeof(pfdata_file),NULL);
358: PetscNew(&pfdata);
359: PFReadMatPowerData(pfdata,pfdata_file);
360: Sbase = pfdata->sbase;
362: numEdges[0] = pfdata->nbranch;
363: numVertices[0] = pfdata->nbus;
365: PetscMalloc1(2*numEdges[0],&edgelist_power);
366: GetListofEdges_Power(pfdata,edgelist_power);
367: }
368: /* Broadcast power Sbase to all processors */
369: MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
370: appctx_power->Sbase = Sbase;
371: appctx_power->jac_error = PETSC_FALSE;
372: /* If external option activated. Introduce error in jacobian */
373: PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);
375: /* proc[1] GET DATA FOR THE SECOND SUBNETWORK: Water */
376: if (size == 1 || (size > 1 && rank == 1)) {
377: PetscNew(&waterdata);
378: PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,sizeof(waterdata_file),NULL);
379: WaterReadData(waterdata,waterdata_file);
381: PetscCalloc1(2*waterdata->nedge,&edgelist_water);
382: GetListofEdges_Water(waterdata,edgelist_water);
383: numEdges[1] = waterdata->nedge;
384: numVertices[1] = waterdata->nvertex;
385: }
387: /* Get data for the coupling subnetwork */
388: if (size == 1) { /* TODO: for size > 1, parallel processing coupling is buggy */
389: nsubnetCouple = 1;
390: numEdgesCouple[0] = 1;
392: PetscMalloc1(4*numEdgesCouple[0],&edgelist_couple);
393: edgelist_couple[0] = 0; edgelist_couple[1] = 4; /* from node: net[0] vertex[4] */
394: edgelist_couple[2] = 1; edgelist_couple[3] = 0; /* to node: net[1] vertex[0] */
395: }
396: PetscLogStagePop();
398: /* (2) Create network */
399: /*--------------------*/
400: MPI_Barrier(PETSC_COMM_WORLD);
401: PetscLogStageRegister("Net Setup",&stage[1]);
402: PetscLogStagePush(stage[1]);
404: PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);
405: PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);
406: PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);
408: /* Create an empty network object */
409: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
411: /* Register the components in the network */
412: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);
413: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);
414: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);
415: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);
417: DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);
418: DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);
420: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Total local nvertices %D + %D = %D, nedges %D + %D + %D = %D\n",rank,numVertices[0],numVertices[1],numVertices[0]+numVertices[1],numEdges[0],numEdges[1],numEdgesCouple[0],numEdges[0]+numEdges[1]+numEdgesCouple[0]);
421: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
423: DMNetworkSetSizes(networkdm,nsubnet,numVertices,numEdges,nsubnetCouple,numEdgesCouple);
425: /* Add edge connectivity */
426: edgelist[0] = edgelist_power;
427: edgelist[1] = edgelist_water;
428: DMNetworkSetEdgeList(networkdm,edgelist,&edgelist_couple);
430: /* Set up the network layout */
431: DMNetworkLayoutSetUp(networkdm);
433: /* Add network components - only process[0] has any data to add */
434: /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
435: genj = 0; loadj = 0;
436: if (rank == 0) {
437: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
438: /* PetscPrintf(PETSC_COMM_SELF,"Power network: nv %D, ne %D\n",nv,ne); */
439: for (i = 0; i < ne; i++) {
440: DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i]);
441: }
443: for (i = 0; i < nv; i++) {
444: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i]);
445: if (pfdata->bus[i].ngen) {
446: for (j = 0; j < pfdata->bus[i].ngen; j++) {
447: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++]);
448: }
449: }
450: if (pfdata->bus[i].nload) {
451: for (j=0; j < pfdata->bus[i].nload; j++) {
452: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++]);
453: }
454: }
455: /* Add number of variables */
456: DMNetworkAddNumVariables(networkdm,vtx[i],2);
457: }
458: }
460: /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
461: if (size == 1 || (size > 1 && rank == 1)) {
462: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
463: /* PetscPrintf(PETSC_COMM_SELF,"Water network: nv %D, ne %D\n",nv,ne); */
464: for (i = 0; i < ne; i++) {
465: DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i]);
466: }
468: for (i = 0; i < nv; i++) {
469: DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[i]);
470: /* Add number of variables */
471: DMNetworkAddNumVariables(networkdm,vtx[i],1);
472: }
473: }
475: /* Set up DM for use */
476: DMSetUp(networkdm);
477: if (viewDM) {
478: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
479: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
480: }
482: DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);
483: if (ne && viewDM) {
484: const PetscInt *cone;
485: PetscInt vid[2];
486: DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);
488: DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);
489: DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);
490: PetscPrintf(PETSC_COMM_SELF,"[%d] Coupling subnetwork: ne %d; connected vertices %d %d\n",rank,ne,vid[0],vid[1]);
491: }
493: /* Free user objects */
494: if (!rank) {
495: PetscFree(edgelist_power);
496: PetscFree(pfdata->bus);
497: PetscFree(pfdata->gen);
498: PetscFree(pfdata->branch);
499: PetscFree(pfdata->load);
500: PetscFree(pfdata);
501: }
502: if (size == 1 || (size > 1 && rank == 1)) {
503: PetscFree(edgelist_water);
504: PetscFree(waterdata->vertex);
505: PetscFree(waterdata->edge);
506: PetscFree(waterdata);
507: }
508: if (size == 1) {
509: PetscFree(edgelist_couple);
510: }
512: /* Re-distribute networkdm to multiple processes for better job balance */
513: if (distribute) {
514: DMNetworkDistribute(&networkdm,0);
515: if (viewDM) {
516: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
517: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
518: }
519: }
521: /* Test DMNetworkGetSubnetworkCoupleInfo() */
522: MPI_Barrier(PETSC_COMM_WORLD);
523: if (test) {
524: for (i=0; i<nsubnetCouple; i++) {
525: DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);
526: if (ne && viewDM) {
527: const PetscInt *cone;
528: PetscInt vid[2];
529: DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);
531: DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);
532: DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);
533: PetscPrintf(PETSC_COMM_SELF,"[%d] After DMNetworkDistribute(), subnetworkCouple[%d]: ne %d; connected vertices %d %d\n",rank,i,ne,vid[0],vid[1]);
534: }
535: }
536: }
538: DMCreateGlobalVector(networkdm,&X);
539: VecDuplicate(X,&F);
540: DMGetLocalVector(networkdm,&user.localXold);
542: PetscLogStagePop();
544: /* (3) Setup Solvers */
545: /*-------------------*/
546: PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);
547: PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
549: PetscLogStageRegister("SNES Setup",&stage[2]);
550: PetscLogStagePush(stage[2]);
552: SetInitialGuess(networkdm,X,&user);
553: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
555: /* Create coupled snes */
556: /*-------------------- */
557: PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");
558: user.subsnes_id = nsubnet;
559: SNESCreate(PETSC_COMM_WORLD,&snes);
560: SNESSetDM(snes,networkdm);
561: SNESSetOptionsPrefix(snes,"coupled_");
562: SNESSetFunction(snes,F,FormFunction,&user);
563: SNESMonitorSet(snes,UserMonitor,&user,NULL);
564: SNESSetFromOptions(snes);
566: if (viewJ) {
567: /* View Jac structure */
568: SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
569: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
570: }
572: if (viewX) {
573: PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");
574: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
575: }
577: if (viewJ) {
578: /* View assembled Jac */
579: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
580: }
582: /* Create snes_power */
583: /*-------------------*/
584: PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");
585: SetInitialGuess(networkdm,X,&user);
587: user.subsnes_id = 0;
588: SNESCreate(PETSC_COMM_WORLD,&snes_power);
589: SNESSetDM(snes_power,networkdm);
590: SNESSetOptionsPrefix(snes_power,"power_");
591: SNESSetFunction(snes_power,F,FormFunction,&user);
592: SNESMonitorSet(snes_power,UserMonitor,&user,NULL);
594: /* Use user-provide Jacobian */
595: DMCreateMatrix(networkdm,&Jac);
596: SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);
597: SNESSetFromOptions(snes_power);
599: if (viewX) {
600: PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");
601: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
602: }
603: if (viewJ) {
604: PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");
605: SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);
606: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
607: /* MatView(Jac,PETSC_VIEWER_STDOUT_WORLD); */
608: }
610: /* Create snes_water */
611: /*-------------------*/
612: PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");
613: SetInitialGuess(networkdm,X,&user);
615: user.subsnes_id = 1;
616: SNESCreate(PETSC_COMM_WORLD,&snes_water);
617: SNESSetDM(snes_water,networkdm);
618: SNESSetOptionsPrefix(snes_water,"water_");
619: SNESSetFunction(snes_water,F,FormFunction,&user);
620: SNESMonitorSet(snes_water,UserMonitor,&user,NULL);
621: SNESSetFromOptions(snes_water);
623: if (viewX) {
624: PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");
625: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
626: }
627: PetscLogStagePop();
629: /* (4) Solve */
630: /*-----------*/
631: PetscLogStageRegister("SNES Solve",&stage[3]);
632: PetscLogStagePush(stage[3]);
633: user.it = 0;
634: reason = SNES_DIVERGED_DTOL;
635: while (user.it < it_max && (PetscInt)reason<0) {
636: #if 0
637: user.subsnes_id = 0;
638: SNESSolve(snes_power,NULL,X);
640: user.subsnes_id = 1;
641: SNESSolve(snes_water,NULL,X);
642: #endif
643: user.subsnes_id = nsubnet;
644: SNESSolve(snes,NULL,X);
646: SNESGetConvergedReason(snes,&reason);
647: user.it++;
648: }
649: PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);
650: if (viewX) {
651: PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");
652: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
653: }
654: PetscLogStagePop();
656: /* Free objects */
657: /* -------------*/
658: VecDestroy(&X);
659: VecDestroy(&F);
660: DMRestoreLocalVector(networkdm,&user.localXold);
662: SNESDestroy(&snes);
663: MatDestroy(&Jac);
664: SNESDestroy(&snes_power);
665: SNESDestroy(&snes_water);
667: DMDestroy(&networkdm);
668: PetscFinalize();
669: return ierr;
670: }
672: /*TEST
674: build:
675: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
676: depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c
678: test:
679: args: -coupled_snes_converged_reason -options_left no
680: localrunfiles: ex1options power/case9.m water/sample1.inp
681: output_file: output/ex1.out
682: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
684: test:
685: suffix: 2
686: nsize: 3
687: args: -coupled_snes_converged_reason -options_left no
688: localrunfiles: ex1options power/case9.m water/sample1.inp
689: output_file: output/ex1_2.out
690: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
692: test:
693: suffix: 3
694: nsize: 3
695: args: -coupled_snes_converged_reason -options_left no -distribute false
696: localrunfiles: ex1options power/case9.m water/sample1.inp
697: output_file: output/ex1_2.out
698: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
700: TEST*/