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