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 %D, subsnes_id %D, fnorm %g\n",user->it,user->subsnes_id,(double)fnorm);
34: } else {
35: PetscPrintf(PETSC_COMM_SELF," subsnes_id %D, 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++) {
113: farr[offset+j] = xarr[offset+j] - xoldarr[offset+j];
114: }
115: }
117: VecRestoreArrayRead(localX,&xarr);
118: VecRestoreArrayRead(localXold,&xoldarr);
119: VecRestoreArray(localF,&farr);
120: return 0;
121: }
123: PetscErrorCode FormFunction(SNES snes,Vec X,Vec F,void *appctx)
124: {
125: DM networkdm;
126: Vec localX,localF;
127: PetscInt nv,ne,v;
128: const PetscInt *vtx,*edges;
129: PetscMPIInt rank;
130: MPI_Comm comm;
131: UserCtx *user = (UserCtx*)appctx;
132: UserCtx_Power appctx_power = (*user).appctx_power;
133: AppCtx_Water appctx_water = (*user).appctx_water;
135: SNESGetDM(snes,&networkdm);
136: PetscObjectGetComm((PetscObject)networkdm,&comm);
137: MPI_Comm_rank(comm,&rank);
139: DMGetLocalVector(networkdm,&localX);
140: DMGetLocalVector(networkdm,&localF);
141: VecSet(F,0.0);
142: VecSet(localF,0.0);
144: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
145: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
147: /* Form Function for power subnetwork */
148: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
149: if (user->subsnes_id == 1) { /* snes_water only */
150: FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
151: } else {
152: FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,&appctx_power);
153: }
155: /* Form Function for water subnetwork */
156: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
157: if (user->subsnes_id == 0) { /* snes_power only */
158: FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
159: } else {
160: FormFunction_Water(networkdm,localX,localF,nv,ne,vtx,edges,NULL);
161: }
163: /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */
164: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
165: for (v=0; v<nv; v++) {
166: PetscInt key,ncomp,nvar,nconnedges,k,e,keye,goffset[3];
167: void* component;
168: const PetscInt *connedges;
170: DMNetworkGetComponent(networkdm,vtx[v],ALL_COMPONENTS,NULL,NULL,&nvar);
171: DMNetworkGetNumComponents(networkdm,vtx[v],&ncomp);
172: /* printf(" [%d] coupling vertex[%D]: v %D, ncomp %D; nvar %D\n",rank,v,vtx[v], ncomp,nvar); */
174: for (k=0; k<ncomp; k++) {
175: DMNetworkGetComponent(networkdm,vtx[v],k,&key,&component,&nvar);
176: DMNetworkGetGlobalVecOffset(networkdm,vtx[v],k,&goffset[k]);
178: /* Verify the coupling vertex is a powernet load vertex or a water vertex */
179: switch (k) {
180: case 0:
182: break;
183: case 1:
185: break;
186: case 2:
188: break;
189: default: SETERRQ(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %D is wrong",k);
190: }
191: /* printf(" [%d] coupling vertex[%D]: key %D; nvar %D, goffset %D\n",rank,v,key,nvar,goffset[k]); */
192: }
194: /* Get its supporting edges */
195: DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
196: /* printf("\n[%d] coupling vertex: nconnedges %D\n",rank,nconnedges); */
197: for (k=0; k<nconnedges; k++) {
198: e = connedges[k];
199: DMNetworkGetNumComponents(networkdm,e,&ncomp);
200: /* printf("\n [%d] connected edge[%D]=%D has ncomp %D\n",rank,k,e,ncomp); */
201: DMNetworkGetComponent(networkdm,e,0,&keye,&component,NULL);
202: if (keye == appctx_water.compkey_edge) { /* water_compkey_edge */
203: EDGE_Water edge=(EDGE_Water)component;
204: if (edge->type == EDGE_TYPE_PUMP) {
205: /* printf(" connected edge[%D]=%D has keye=%D, is appctx_water.compkey_edge with EDGE_TYPE_PUMP\n",k,e,keye); */
206: }
207: } else { /* ower->compkey_branch */
209: }
210: }
211: }
213: DMRestoreLocalVector(networkdm,&localX);
215: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
216: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
217: DMRestoreLocalVector(networkdm,&localF);
218: #if 0
219: if (rank == 0) printf("F:\n");
220: VecView(F,PETSC_VIEWER_STDOUT_WORLD);
221: #endif
222: return 0;
223: }
225: PetscErrorCode SetInitialGuess(DM networkdm,Vec X,void* appctx)
226: {
227: PetscInt nv,ne,i,j,ncomp,offset,key;
228: const PetscInt *vtx,*edges;
229: UserCtx *user = (UserCtx*)appctx;
230: Vec localX = user->localXold;
231: UserCtx_Power appctx_power = (*user).appctx_power;
232: AppCtx_Water appctx_water = (*user).appctx_water;
233: PetscBool ghost;
234: PetscScalar *xarr;
235: VERTEX_Power bus;
236: VERTEX_Water vertex;
237: void* component;
238: GEN gen;
240: VecSet(X,0.0);
241: VecSet(localX,0.0);
243: /* Set initial guess for power subnetwork */
244: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
245: SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,&appctx_power);
247: /* Set initial guess for water subnetwork */
248: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
249: SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);
251: /* Set initial guess at the coupling vertex */
252: VecGetArray(localX,&xarr);
253: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
254: for (i=0; i<nv; i++) {
255: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);
256: if (ghost) continue;
258: DMNetworkGetNumComponents(networkdm,vtx[i],&ncomp);
259: for (j=0; j<ncomp; j++) {
260: DMNetworkGetLocalVecOffset(networkdm,vtx[i],j,&offset);
261: DMNetworkGetComponent(networkdm,vtx[i],j,&key,(void**)&component,NULL);
262: if (key == appctx_power.compkey_bus) {
263: bus = (VERTEX_Power)(component);
264: xarr[offset] = bus->va*PETSC_PI/180.0;
265: xarr[offset+1] = bus->vm;
266: } else if (key == appctx_power.compkey_gen) {
267: gen = (GEN)(component);
268: if (!gen->status) continue;
269: xarr[offset+1] = gen->vs;
270: } else if (key == appctx_water.compkey_vtx) {
271: vertex = (VERTEX_Water)(component);
272: if (vertex->type == VERTEX_TYPE_JUNCTION) {
273: xarr[offset] = 100;
274: } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
275: xarr[offset] = vertex->res.head;
276: } else {
277: xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
278: }
279: }
280: }
281: }
282: VecRestoreArray(localX,&xarr);
284: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
285: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
286: return 0;
287: }
289: int main(int argc,char **argv)
290: {
291: DM networkdm;
292: PetscLogStage stage[4];
293: PetscMPIInt rank,size;
294: PetscInt Nsubnet=2,numVertices[2],numEdges[2],i,j,nv,ne,it_max=10;
295: const PetscInt *vtx,*edges;
296: Vec X,F;
297: SNES snes,snes_power,snes_water;
298: Mat Jac;
299: PetscBool ghost,viewJ=PETSC_FALSE,viewX=PETSC_FALSE,viewDM=PETSC_FALSE,test=PETSC_FALSE,distribute=PETSC_TRUE,flg;
300: UserCtx user;
301: SNESConvergedReason reason;
303: /* Power subnetwork */
304: UserCtx_Power *appctx_power = &user.appctx_power;
305: char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
306: PFDATA *pfdata = NULL;
307: PetscInt genj,loadj,*edgelist_power = NULL,power_netnum;
308: PetscScalar Sbase = 0.0;
310: /* Water subnetwork */
311: AppCtx_Water *appctx_water = &user.appctx_water;
312: WATERDATA *waterdata = NULL;
313: char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
314: PetscInt *edgelist_water = NULL,water_netnum;
316: /* Shared vertices between subnetworks */
317: PetscInt power_svtx,water_svtx;
319: PetscInitialize(&argc,&argv,"ex1options",help);
320: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
321: MPI_Comm_size(PETSC_COMM_WORLD,&size);
323: /* (1) Read Data - Only rank 0 reads the data */
324: /*--------------------------------------------*/
325: PetscLogStageRegister("Read Data",&stage[0]);
326: PetscLogStagePush(stage[0]);
328: for (i=0; i<Nsubnet; i++) {
329: numVertices[i] = 0;
330: numEdges[i] = 0;
331: }
333: /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
334: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
335: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
336: PetscNew(&pfdata);
337: PFReadMatPowerData(pfdata,pfdata_file);
338: if (rank == 0) {
339: PetscPrintf(PETSC_COMM_SELF,"Power network: nb = %D, ngen = %D, nload = %D, nbranch = %D\n",pfdata->nbus,pfdata->ngen,pfdata->nload,pfdata->nbranch);
340: }
341: Sbase = pfdata->sbase;
342: if (rank == 0) { /* proc[0] will create Electric Power Grid */
343: numEdges[0] = pfdata->nbranch;
344: numVertices[0] = pfdata->nbus;
346: PetscMalloc1(2*numEdges[0],&edgelist_power);
347: GetListofEdges_Power(pfdata,edgelist_power);
348: }
349: /* Broadcast power Sbase to all processors */
350: MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
351: appctx_power->Sbase = Sbase;
352: appctx_power->jac_error = PETSC_FALSE;
353: /* If external option activated. Introduce error in jacobian */
354: PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);
356: /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
357: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
358: PetscNew(&waterdata);
359: PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,PETSC_MAX_PATH_LEN-1,NULL);
360: WaterReadData(waterdata,waterdata_file);
361: if (size == 1 || (size > 1 && rank == 1)) {
362: PetscCalloc1(2*waterdata->nedge,&edgelist_water);
363: GetListofEdges_Water(waterdata,edgelist_water);
364: numEdges[1] = waterdata->nedge;
365: numVertices[1] = waterdata->nvertex;
366: }
367: PetscLogStagePop();
369: /* (2) Create a network consist of two subnetworks */
370: /*-------------------------------------------------*/
371: PetscLogStageRegister("Net Setup",&stage[1]);
372: PetscLogStagePush(stage[1]);
374: PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);
375: PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);
376: PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);
378: /* Create an empty network object */
379: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
381: /* Register the components in the network */
382: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);
383: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);
384: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);
385: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);
387: DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);
388: DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);
389: #if 0
390: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_branch %d\n",appctx_power->compkey_branch);
391: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_bus %d\n",appctx_power->compkey_bus);
392: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_gen %d\n",appctx_power->compkey_gen);
393: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_load %d\n",appctx_power->compkey_load);
394: PetscPrintf(PETSC_COMM_WORLD,"water->compkey_edge %d\n",appctx_water->compkey_edge);
395: PetscPrintf(PETSC_COMM_WORLD,"water->compkey_vtx %d\n",appctx_water->compkey_vtx);
396: #endif
397: PetscSynchronizedPrintf(PETSC_COMM_WORLD,"[%d] Total local nvertices %D + %D = %D, nedges %D + %D = %D\n",rank,numVertices[0],numVertices[1],numVertices[0]+numVertices[1],numEdges[0],numEdges[1],numEdges[0]+numEdges[1]);
398: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
400: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,Nsubnet);
401: DMNetworkAddSubnetwork(networkdm,"power",numEdges[0],edgelist_power,&power_netnum);
402: DMNetworkAddSubnetwork(networkdm,"water",numEdges[1],edgelist_water,&water_netnum);
404: /* vertex subnet[0].4 shares with vertex subnet[1].0 */
405: power_svtx = 4; water_svtx = 0;
406: DMNetworkAddSharedVertices(networkdm,power_netnum,water_netnum,1,&power_svtx,&water_svtx);
408: /* Set up the network layout */
409: DMNetworkLayoutSetUp(networkdm);
411: /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
412: /*-------------------------------------------------------*/
413: genj = 0; loadj = 0;
414: DMNetworkGetSubnetwork(networkdm,power_netnum,&nv,&ne,&vtx,&edges);
416: for (i = 0; i < ne; i++) {
417: DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i],0);
418: }
420: for (i = 0; i < nv; i++) {
421: DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
422: if (flg) continue;
424: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i],2);
425: if (pfdata->bus[i].ngen) {
426: for (j = 0; j < pfdata->bus[i].ngen; j++) {
427: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++],0);
428: }
429: }
430: if (pfdata->bus[i].nload) {
431: for (j=0; j < pfdata->bus[i].nload; j++) {
432: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++],0);
433: }
434: }
435: }
437: /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
438: /*-------------------------------------------------------*/
439: DMNetworkGetSubnetwork(networkdm,water_netnum,&nv,&ne,&vtx,&edges);
440: for (i = 0; i < ne; i++) {
441: DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i],0);
442: }
444: for (i = 0; i < nv; i++) {
445: DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
446: if (flg) continue;
448: DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[i],1);
449: }
451: /* 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 */
452: /*----------------------------------------------------------------------------------------------------------------------------*/
453: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
454: for (i = 0; i < nv; i++) {
455: /* power */
456: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[4],2);
457: /* bus[4] is a load, add its component */
458: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[0],0);
460: /* water */
461: DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[0],1);
462: }
464: /* Set up DM for use */
465: DMSetUp(networkdm);
466: if (viewDM) {
467: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
468: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
469: }
471: /* Free user objects */
472: PetscFree(edgelist_power);
473: PetscFree(pfdata->bus);
474: PetscFree(pfdata->gen);
475: PetscFree(pfdata->branch);
476: PetscFree(pfdata->load);
477: PetscFree(pfdata);
479: PetscFree(edgelist_water);
480: PetscFree(waterdata->vertex);
481: PetscFree(waterdata->edge);
482: PetscFree(waterdata);
484: /* Re-distribute networkdm to multiple processes for better job balance */
485: if (size >1 && distribute) {
486: DMNetworkDistribute(&networkdm,0);
487: if (viewDM) {
488: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
489: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
490: }
491: }
493: /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
494: if (test) {
495: PetscInt v,gidx;
496: MPI_Barrier(PETSC_COMM_WORLD);
497: for (i=0; i<Nsubnet; i++) {
498: DMNetworkGetSubnetwork(networkdm,i,&nv,&ne,&vtx,&edges);
499: PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, subnet[%d] ne %d, nv %d\n",rank,i,ne,nv);
500: MPI_Barrier(PETSC_COMM_WORLD);
502: for (v=0; v<nv; v++) {
503: DMNetworkIsGhostVertex(networkdm,vtx[v],&ghost);
504: DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
505: PetscPrintf(PETSC_COMM_SELF,"[%d] subnet[%d] v %d %d; ghost %d\n",rank,i,vtx[v],gidx,ghost);
506: }
507: MPI_Barrier(PETSC_COMM_WORLD);
508: }
509: MPI_Barrier(PETSC_COMM_WORLD);
511: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
512: PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, num of shared vertices nsv = %d\n",rank,nv);
513: for (v=0; v<nv; v++) {
514: DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
515: PetscPrintf(PETSC_COMM_SELF,"[%d] sv %d, gidx=%d\n",rank,vtx[v],gidx);
516: }
517: MPI_Barrier(PETSC_COMM_WORLD);
518: }
520: /* Create solution vector X */
521: DMCreateGlobalVector(networkdm,&X);
522: VecDuplicate(X,&F);
523: DMGetLocalVector(networkdm,&user.localXold);
524: PetscLogStagePop();
526: /* (3) Setup Solvers */
527: /*-------------------*/
528: PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);
529: PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
531: PetscLogStageRegister("SNES Setup",&stage[2]);
532: PetscLogStagePush(stage[2]);
534: SetInitialGuess(networkdm,X,&user);
536: /* Create coupled snes */
537: /*-------------------- */
538: PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");
539: user.subsnes_id = Nsubnet;
540: SNESCreate(PETSC_COMM_WORLD,&snes);
541: SNESSetDM(snes,networkdm);
542: SNESSetOptionsPrefix(snes,"coupled_");
543: SNESSetFunction(snes,F,FormFunction,&user);
544: SNESMonitorSet(snes,UserMonitor,&user,NULL);
545: SNESSetFromOptions(snes);
547: if (viewJ) {
548: /* View Jac structure */
549: SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
550: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
551: }
553: if (viewX) {
554: PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");
555: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
556: }
558: if (viewJ) {
559: /* View assembled Jac */
560: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
561: }
563: /* Create snes_power */
564: /*-------------------*/
565: PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");
567: user.subsnes_id = 0;
568: SNESCreate(PETSC_COMM_WORLD,&snes_power);
569: SNESSetDM(snes_power,networkdm);
570: SNESSetOptionsPrefix(snes_power,"power_");
571: SNESSetFunction(snes_power,F,FormFunction,&user);
572: SNESMonitorSet(snes_power,UserMonitor,&user,NULL);
574: /* Use user-provide Jacobian */
575: DMCreateMatrix(networkdm,&Jac);
576: SNESSetJacobian(snes_power,Jac,Jac,FormJacobian_subPower,&user);
577: SNESSetFromOptions(snes_power);
579: if (viewX) {
580: PetscPrintf(PETSC_COMM_WORLD,"Power Solution:\n");
581: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
582: }
583: if (viewJ) {
584: PetscPrintf(PETSC_COMM_WORLD,"Power Jac:\n");
585: SNESGetJacobian(snes_power,&Jac,NULL,NULL,NULL);
586: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
587: /* MatView(Jac,PETSC_VIEWER_STDOUT_WORLD); */
588: }
590: /* Create snes_water */
591: /*-------------------*/
592: PetscPrintf(PETSC_COMM_WORLD,"SNES_water setup......\n");
594: user.subsnes_id = 1;
595: SNESCreate(PETSC_COMM_WORLD,&snes_water);
596: SNESSetDM(snes_water,networkdm);
597: SNESSetOptionsPrefix(snes_water,"water_");
598: SNESSetFunction(snes_water,F,FormFunction,&user);
599: SNESMonitorSet(snes_water,UserMonitor,&user,NULL);
600: SNESSetFromOptions(snes_water);
602: if (viewX) {
603: PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");
604: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
605: }
606: PetscLogStagePop();
608: /* (4) Solve */
609: /*-----------*/
610: PetscLogStageRegister("SNES Solve",&stage[3]);
611: PetscLogStagePush(stage[3]);
612: user.it = 0;
613: reason = SNES_DIVERGED_DTOL;
614: while (user.it < it_max && (PetscInt)reason<0) {
615: #if 0
616: user.subsnes_id = 0;
617: SNESSolve(snes_power,NULL,X);
619: user.subsnes_id = 1;
620: SNESSolve(snes_water,NULL,X);
621: #endif
622: user.subsnes_id = Nsubnet;
623: SNESSolve(snes,NULL,X);
625: SNESGetConvergedReason(snes,&reason);
626: user.it++;
627: }
628: PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);
629: if (viewX) {
630: PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");
631: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
632: }
633: PetscLogStagePop();
635: /* Free objects */
636: /* -------------*/
637: VecDestroy(&X);
638: VecDestroy(&F);
639: DMRestoreLocalVector(networkdm,&user.localXold);
641: SNESDestroy(&snes);
642: MatDestroy(&Jac);
643: SNESDestroy(&snes_power);
644: SNESDestroy(&snes_water);
646: DMDestroy(&networkdm);
647: PetscFinalize();
648: return 0;
649: }
651: /*TEST
653: build:
654: requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
655: depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c
657: test:
658: args: -coupled_snes_converged_reason -options_left no -viewDM
659: localrunfiles: ex1options power/case9.m water/sample1.inp
660: output_file: output/ex1.out
662: test:
663: suffix: 2
664: nsize: 3
665: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis
666: localrunfiles: ex1options power/case9.m water/sample1.inp
667: output_file: output/ex1_2.out
668: requires: parmetis
670: # test:
671: # suffix: 3
672: # nsize: 3
673: # args: -coupled_snes_converged_reason -options_left no -distribute false
674: # localrunfiles: ex1options power/case9.m water/sample1.inp
675: # output_file: output/ex1_2.out
677: test:
678: suffix: 4
679: nsize: 4
680: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -viewDM
681: localrunfiles: ex1options power/case9.m water/sample1.inp
682: output_file: output/ex1_4.out
684: TEST*/