Actual source code: ex1.c
petsc-3.13.4 2020-08-01
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,fnorm);
41: } else {
42: PetscPrintf(PETSC_COMM_SELF," subsnes_id %d, fnorm %g\n",user->subsnes_id,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: PetscLogStage stage[4];
308: PetscMPIInt rank;
309: PetscInt nsubnet=2,nsubnetCouple=1,numVertices[2],numEdges[2],numEdgesCouple[1];
310: PetscInt i,j,nv,ne;
311: PetscInt *edgelist[2];
312: const PetscInt *vtx,*edges;
313: Vec X,F;
314: SNES snes,snes_power,snes_water;
315: Mat Jac;
316: PetscBool viewJ=PETSC_FALSE,viewX=PETSC_FALSE,viewDM=PETSC_FALSE,test=PETSC_FALSE;
317: UserCtx user;
318: PetscInt it_max=10;
319: SNESConvergedReason reason;
321: /* Power subnetwork */
322: UserCtx_Power *appctx_power = &user.appctx_power;
323: char pfdata_file[PETSC_MAX_PATH_LEN]="power/case9.m";
324: PFDATA *pfdata=NULL;
325: PetscInt genj,loadj;
326: PetscInt *edgelist_power=NULL;
327: PetscScalar Sbase=0.0;
329: /* Water subnetwork */
330: AppCtx_Water *appctx_water = &user.appctx_water;
331: WATERDATA *waterdata=NULL;
332: char waterdata_file[PETSC_MAX_PATH_LEN]="water/sample1.inp";
333: PetscInt *edgelist_water=NULL;
335: /* Coupling subnetwork */
336: PetscInt *edgelist_couple=NULL;
338: PetscInitialize(&argc,&argv,"ex1options",help);if (ierr) return ierr;
339: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
341: /* (1) Read Data - Only rank 0 reads the data */
342: /*--------------------------------------------*/
343: PetscLogStageRegister("Read Data",&stage[0]);
344: PetscLogStagePush(stage[0]);
346: for (i=0; i<nsubnet; i++) {
347: numVertices[i] = 0;
348: numEdges[i] = 0;
349: }
350: for (i=0; i<nsubnetCouple; i++) {
351: numEdgesCouple[0] = 0;
352: }
354: /* READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
355: if (!rank) {
356: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
357: PetscNew(&pfdata);
358: PFReadMatPowerData(pfdata,pfdata_file);
359: Sbase = pfdata->sbase;
361: numEdges[0] = pfdata->nbranch;
362: numVertices[0] = pfdata->nbus;
364: PetscMalloc1(2*numEdges[0],&edgelist_power);
365: GetListofEdges_Power(pfdata,edgelist_power);
366: #if 0
367: PetscPrintf(PETSC_COMM_WORLD,"edgelist_power:\n");
368: for (i=0; i<numEdges[0]; i++) {
369: PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edgelist_power[2*i],edgelist_power[2*i+1]);
370: }
371: PetscPrintf(PETSC_COMM_WORLD,"\n");
372: #endif
373: }
374: /* Broadcast power Sbase to all processors */
375: MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
376: appctx_power->Sbase = Sbase;
377: appctx_power->jac_error = PETSC_FALSE;
378: /* If external option activated. Introduce error in jacobian */
379: PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);
381: /* GET DATA FOR THE SECOND SUBNETWORK: Water */
382: if (!rank) {
383: PetscNew(&waterdata);
384: PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,PETSC_MAX_PATH_LEN-1,NULL);
385: WaterReadData(waterdata,waterdata_file);
387: PetscCalloc1(2*waterdata->nedge,&edgelist_water);
388: GetListofEdges_Water(waterdata,edgelist_water);
389: numEdges[1] = waterdata->nedge;
390: numVertices[1] = waterdata->nvertex;
391: #if 0
392: PetscPrintf(PETSC_COMM_WORLD,"edgelist_water:\n");
393: for (i=0; i<numEdges[1]; i++) {
394: PetscPrintf(PETSC_COMM_SELF,"[%D %D]",edgelist_water[2*i],edgelist_water[2*i+1]);
395: }
396: PetscPrintf(PETSC_COMM_WORLD,("\n");
397: #endif
398: }
400: /* Get data for the coupling subnetwork */
401: if (!rank) {
402: numEdgesCouple[0] = 1;
404: PetscMalloc1(4*numEdgesCouple[0],&edgelist_couple);
405: edgelist_couple[0] = 0; edgelist_couple[1] = 4; /* from node: net[0] vertex[4] */
406: edgelist_couple[2] = 1; edgelist_couple[3] = 0; /* to node: net[1] vertex[0] */
407: }
408: PetscLogStagePop();
410: /* (2) Create network */
411: /*--------------------*/
412: MPI_Barrier(PETSC_COMM_WORLD);
413: PetscLogStageRegister("Net Setup",&stage[1]);
414: PetscLogStagePush(stage[1]);
416: PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);
417: PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);
419: /* Create an empty network object */
420: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
422: /* Register the components in the network */
423: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);
424: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);
425: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);
426: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);
428: DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);
429: DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);
431: if (!rank) {
432: PetscPrintf(PETSC_COMM_SELF,"[%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]);
433: }
435: DMNetworkSetSizes(networkdm,nsubnet,numVertices,numEdges,nsubnetCouple,numEdgesCouple);
437: /* Add edge connectivity */
438: edgelist[0] = edgelist_power;
439: edgelist[1] = edgelist_water;
440: DMNetworkSetEdgeList(networkdm,edgelist,&edgelist_couple);
442: /* Set up the network layout */
443: DMNetworkLayoutSetUp(networkdm);
445: /* Add network components - only process[0] has any data to add */
446: /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
447: genj = 0; loadj = 0;
448: if (!rank) {
449: DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);
450: PetscPrintf(PETSC_COMM_SELF,"Power network: nv %D, ne %D\n",nv,ne);
451: for (i = 0; i < ne; i++) {
452: DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i]);
453: }
455: for (i = 0; i < nv; i++) {
456: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i]);
457: if (pfdata->bus[i].ngen) {
458: for (j = 0; j < pfdata->bus[i].ngen; j++) {
459: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++]);
460: }
461: }
462: if (pfdata->bus[i].nload) {
463: for (j=0; j < pfdata->bus[i].nload; j++) {
464: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++]);
465: }
466: }
467: /* Add number of variables */
468: DMNetworkAddNumVariables(networkdm,vtx[i],2);
469: }
470: }
472: /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
473: if (!rank) {
474: DMNetworkGetSubnetworkInfo(networkdm,1,&nv,&ne,&vtx,&edges);
475: PetscPrintf(PETSC_COMM_SELF,"Water network: nv %D, ne %D\n",nv,ne);
476: for (i = 0; i < ne; i++) {
477: DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i]);
478: }
480: for (i = 0; i < nv; i++) {
481: DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[i]);
482: /* Add number of variables */
483: DMNetworkAddNumVariables(networkdm,vtx[i],1);
484: }
485: }
487: /* Set up DM for use */
488: DMSetUp(networkdm);
489: if (viewDM) {
490: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
491: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
492: }
494: DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);
495: if (ne && viewDM) {
496: const PetscInt *cone;
497: PetscInt vid[2];
498: DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);
500: DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);
501: DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);
502: PetscPrintf(PETSC_COMM_SELF,"[%d] Coupling subnetwork: ne %d; connected vertices %d %d\n",rank,ne,vid[0],vid[1]);
503: }
505: /* Free user objects */
506: if (!rank) {
507: PetscFree(edgelist_power);
508: PetscFree(pfdata->bus);
509: PetscFree(pfdata->gen);
510: PetscFree(pfdata->branch);
511: PetscFree(pfdata->load);
512: PetscFree(pfdata);
514: PetscFree(edgelist_water);
515: PetscFree(waterdata->vertex);
516: PetscFree(waterdata->edge);
518: PetscFree(edgelist_couple);
519: PetscFree(waterdata);
520: }
522: /* Distribute networkdm to multiple processes */
523: DMNetworkDistribute(&networkdm,0);
524: if (viewDM) {
525: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
526: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
527: }
529: /* Test DMNetworkGetSubnetworkCoupleInfo() */
530: MPI_Barrier(PETSC_COMM_WORLD);
531: if (test) {
532: for (i=0; i<nsubnetCouple; i++) {
533: DMNetworkGetSubnetworkCoupleInfo(networkdm,0,&ne,&edges);
534: if (ne && viewDM) {
535: const PetscInt *cone;
536: PetscInt vid[2];
537: DMNetworkGetConnectedVertices(networkdm,edges[0],&cone);
539: DMNetworkGetGlobalVertexIndex(networkdm,cone[0],&vid[0]);
540: DMNetworkGetGlobalVertexIndex(networkdm,cone[1],&vid[1]);
541: PetscPrintf(PETSC_COMM_SELF,"[%d] After DMNetworkDistribute(), subnetworkCouple[%d]: ne %d; connected vertices %d %d\n",rank,i,ne,vid[0],vid[1]);
542: }
543: }
544: }
546: DMCreateGlobalVector(networkdm,&X);
547: VecDuplicate(X,&F);
548: DMGetLocalVector(networkdm,&user.localXold);
550: PetscLogStagePop();
552: /* (3) Setup Solvers */
553: /*-------------------*/
554: PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);
555: PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
557: PetscLogStageRegister("SNES Setup",&stage[2]);
558: PetscLogStagePush(stage[2]);
560: SetInitialGuess(networkdm,X,&user);
561: /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */
563: /* Create coupled snes */
564: /*-------------------- */
565: PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");
566: user.subsnes_id = nsubnet;
567: SNESCreate(PETSC_COMM_WORLD,&snes);
568: SNESSetDM(snes,networkdm);
569: SNESSetOptionsPrefix(snes,"coupled_");
570: SNESSetFunction(snes,F,FormFunction,&user);
571: SNESMonitorSet(snes,UserMonitor,&user,NULL);
572: SNESSetFromOptions(snes);
574: if (viewJ) {
575: /* View Jac structure */
576: SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
577: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
578: }
580: if (viewX) {
581: PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");
582: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
583: }
585: if (viewJ) {
586: /* View assembled Jac */
587: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
588: }
590: /* Create snes_power */
591: /*-------------------*/
592: PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");
593: SetInitialGuess(networkdm,X,&user);
595: user.subsnes_id = 0;
596: SNESCreate(PETSC_COMM_WORLD,&snes_power);
597: SNESSetDM(snes_power,networkdm);
598: SNESSetOptionsPrefix(snes_power,"power_");
599: SNESSetFunction(snes_power,F,FormFunction,&user);
600: SNESMonitorSet(snes_power,UserMonitor,&user,NULL);
602: /* Use user-provide Jacobian */
603: DMCreateMatrix(networkdm,&Jac);
604: SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);
605: SNESSetFromOptions(snes_power);
607: if (viewX) {
608: PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");
609: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
610: }
611: if (viewJ) {
612: PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");
613: SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);
614: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
615: /* MatView(Jac,PETSC_VIEWER_STDOUT_WORLD); */
616: }
618: /* Create snes_water */
619: /*-------------------*/
620: PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");
621: SetInitialGuess(networkdm,X,&user);
623: user.subsnes_id = 1;
624: SNESCreate(PETSC_COMM_WORLD,&snes_water);
625: SNESSetDM(snes_water,networkdm);
626: SNESSetOptionsPrefix(snes_water,"water_");
627: SNESSetFunction(snes_water,F,FormFunction,&user);
628: SNESMonitorSet(snes_water,UserMonitor,&user,NULL);
629: SNESSetFromOptions(snes_water);
631: if (viewX) {
632: PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");
633: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
634: }
635: PetscLogStagePop();
637: /* (4) Solve */
638: /*-----------*/
639: PetscLogStageRegister("SNES Solve",&stage[3]);
640: PetscLogStagePush(stage[3]);
641: user.it = 0;
642: reason = SNES_DIVERGED_DTOL;
643: while (user.it < it_max && (PetscInt)reason<0) {
644: #if 0
645: user.subsnes_id = 0;
646: SNESSolve(snes_power,NULL,X);
648: user.subsnes_id = 1;
649: SNESSolve(snes_water,NULL,X);
650: #endif
651: user.subsnes_id = nsubnet;
652: SNESSolve(snes,NULL,X);
654: SNESGetConvergedReason(snes,&reason);
655: user.it++;
656: }
657: PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);
658: if (viewX) {
659: PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");
660: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
661: }
662: PetscLogStagePop();
664: /* Free objects */
665: /* -------------*/
666: VecDestroy(&X);
667: VecDestroy(&F);
668: DMRestoreLocalVector(networkdm,&user.localXold);
670: SNESDestroy(&snes);
671: MatDestroy(&Jac);
672: SNESDestroy(&snes_power);
673: SNESDestroy(&snes_water);
675: DMDestroy(&networkdm);
676: PetscFinalize();
677: return ierr;
678: }
680: /*TEST
682: build:
683: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
684: depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c
686: test:
687: args: -coupled_snes_converged_reason -options_left no
688: localrunfiles: ex1options power/case9.m water/sample1.inp
689: output_file: output/ex1.out
690: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
692: test:
693: suffix: 2
694: nsize: 3
695: args: -coupled_snes_converged_reason -options_left no
696: localrunfiles: ex1options power/case9.m water/sample1.inp
697: output_file: output/ex1.out
698: requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
700: TEST*/