Actual source code: ex1.c
petsc-3.15.0 2021-03-30
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: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
79: FormJacobian_Power_private(networkdm,localX,J,nv,ne,vtx,edges,appctx);
81: /* Water subnetwork: Identity */
82: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
83: for (i=0; i<nv; i++) {
84: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);
85: if (ghostvtex) continue;
87: DMNetworkGetGlobalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
88: DMNetworkGetComponent(networkdm,vtx[i],ALL_COMPONENTS,NULL,NULL,&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: DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);
122: DMNetworkGetComponent(networkdm,vtx[i],ALL_COMPONENTS,NULL,NULL,&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,v;
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;
145: AppCtx_Water appctx_water = (*user).appctx_water;
148: SNESGetDM(snes,&networkdm);
149: PetscObjectGetComm((PetscObject)networkdm,&comm);
150: MPI_Comm_rank(comm,&rank);
152: DMGetLocalVector(networkdm,&localX);
153: DMGetLocalVector(networkdm,&localF);
154: VecSet(F,0.0);
155: VecSet(localF,0.0);
157: DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);
158: DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);
160: /* Form Function for power subnetwork */
161: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
162: if (user->subsnes_id == 1) { /* snes_water only */
163: FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
164: } else {
165: FormFunction_Power(networkdm,localX,localF,nv,ne,vtx,edges,&appctx_power);
166: }
168: /* Form Function for water subnetwork */
169: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
170: if (user->subsnes_id == 0) { /* snes_power only */
171: FormFunction_Dummy(networkdm,localX,localF,nv,ne,vtx,edges,user);
172: } else {
173: FormFunction_Water(networkdm,localX,localF,nv,ne,vtx,edges,NULL);
174: }
176: /* Illustrate how to access the coupling vertex of the subnetworks without doing anything to F yet */
177: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
178: for (v=0; v<nv; v++) {
179: PetscInt key,ncomp,nvar,nconnedges,k,e,keye,goffset[3];
180: void* component;
181: const PetscInt *connedges;
183: DMNetworkGetComponent(networkdm,vtx[v],ALL_COMPONENTS,NULL,NULL,&nvar);
184: DMNetworkGetNumComponents(networkdm,vtx[v],&ncomp);
185: /* printf(" [%d] coupling vertex[%D]: v %D, ncomp %D; nvar %D\n",rank,v,vtx[v], ncomp,nvar); */
187: for (k=0; k<ncomp; k++) {
188: DMNetworkGetComponent(networkdm,vtx[v],k,&key,&component,&nvar);
189: DMNetworkGetGlobalVecOffset(networkdm,vtx[v],k,&goffset[k]);
191: /* Verify the coupling vertex is a powernet load vertex or a water vertex */
192: switch (k) {
193: case 0:
194: if (key != appctx_power.compkey_bus || nvar != 2) SETERRQ2(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"key %D not a power bus vertex or nvar %D != 2",key,nvar);
195: break;
196: case 1:
197: if (key != appctx_power.compkey_load || nvar != 0 || goffset[1] != goffset[0]+2) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power load vertex");
198: break;
199: case 2:
200: if (key != appctx_water.compkey_vtx || nvar != 1 || goffset[2] != goffset[1]) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a water vertex");
201: break;
202: default: SETERRQ1(PETSC_COMM_SELF, PETSC_ERR_ARG_OUTOFRANGE, "k %D is wrong",k);
203: }
204: /* printf(" [%d] coupling vertex[%D]: key %D; nvar %D, goffset %D\n",rank,v,key,nvar,goffset[k]); */
205: }
207: /* Get its supporting edges */
208: DMNetworkGetSupportingEdges(networkdm,vtx[v],&nconnedges,&connedges);
209: /* printf("\n[%d] coupling vertex: nconnedges %D\n",rank,nconnedges); */
210: for (k=0; k<nconnedges; k++) {
211: e = connedges[k];
212: DMNetworkGetNumComponents(networkdm,e,&ncomp);
213: /* printf("\n [%d] connected edge[%D]=%D has ncomp %D\n",rank,k,e,ncomp); */
214: DMNetworkGetComponent(networkdm,e,0,&keye,&component,NULL);
215: if (keye == appctx_water.compkey_edge) { /* water_compkey_edge */
216: EDGE_Water edge=(EDGE_Water)component;
217: if (edge->type == EDGE_TYPE_PUMP) {
218: /* printf(" connected edge[%D]=%D has keye=%D, is appctx_water.compkey_edge with EDGE_TYPE_PUMP\n",k,e,keye); */
219: }
220: } else { /* ower->compkey_branch */
221: if (keye != appctx_power.compkey_branch) SETERRQ(PETSC_COMM_SELF,PETSC_ERR_ARG_WRONG,"Not a power branch");
222: }
223: }
224: }
226: DMRestoreLocalVector(networkdm,&localX);
228: DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);
229: DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);
230: DMRestoreLocalVector(networkdm,&localF);
231: #if 0
232: if (!rank) printf("F:\n");
233: VecView(F,PETSC_VIEWER_STDOUT_WORLD);
234: #endif
235: return(0);
236: }
238: PetscErrorCode SetInitialGuess(DM networkdm,Vec X,void* appctx)
239: {
241: PetscInt nv,ne,i,j,ncomp,offset,key;
242: const PetscInt *vtx,*edges;
243: UserCtx *user = (UserCtx*)appctx;
244: Vec localX = user->localXold;
245: UserCtx_Power appctx_power = (*user).appctx_power;
246: AppCtx_Water appctx_water = (*user).appctx_water;
247: PetscBool ghost;
248: PetscScalar *xarr;
249: VERTEX_Power bus;
250: VERTEX_Water vertex;
251: void* component;
252: GEN gen;
255: VecSet(X,0.0);
256: VecSet(localX,0.0);
258: /* Set initial guess for power subnetwork */
259: DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);
260: SetInitialGuess_Power(networkdm,localX,nv,ne,vtx,edges,&appctx_power);
262: /* Set initial guess for water subnetwork */
263: DMNetworkGetSubnetwork(networkdm,1,&nv,&ne,&vtx,&edges);
264: SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);
266: /* Set initial guess at the coupling vertex */
267: VecGetArray(localX,&xarr);
268: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
269: for (i=0; i<nv; i++) {
270: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);
271: if (ghost) continue;
273: DMNetworkGetNumComponents(networkdm,vtx[i],&ncomp);
274: for (j=0; j<ncomp; j++) {
275: DMNetworkGetLocalVecOffset(networkdm,vtx[i],j,&offset);
276: DMNetworkGetComponent(networkdm,vtx[i],j,&key,(void**)&component,NULL);
277: if (key == appctx_power.compkey_bus) {
278: bus = (VERTEX_Power)(component);
279: xarr[offset] = bus->va*PETSC_PI/180.0;
280: xarr[offset+1] = bus->vm;
281: } else if (key == appctx_power.compkey_gen) {
282: gen = (GEN)(component);
283: if (!gen->status) continue;
284: xarr[offset+1] = gen->vs;
285: } else if (key == appctx_water.compkey_vtx) {
286: vertex = (VERTEX_Water)(component);
287: if (vertex->type == VERTEX_TYPE_JUNCTION) {
288: xarr[offset] = 100;
289: } else if (vertex->type == VERTEX_TYPE_RESERVOIR) {
290: xarr[offset] = vertex->res.head;
291: } else {
292: xarr[offset] = vertex->tank.initlvl + vertex->tank.elev;
293: }
294: }
295: }
296: }
297: VecRestoreArray(localX,&xarr);
299: DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);
300: DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);
301: return(0);
302: }
304: int main(int argc,char **argv)
305: {
306: PetscErrorCode ierr;
307: DM networkdm;
308: PetscLogStage stage[4];
309: PetscMPIInt rank,size;
310: PetscInt Nsubnet=2,numVertices[2],numEdges[2],i,j,nv,ne,it_max=10;
311: const PetscInt *vtx,*edges;
312: Vec X,F;
313: SNES snes,snes_power,snes_water;
314: Mat Jac;
315: PetscBool ghost,viewJ=PETSC_FALSE,viewX=PETSC_FALSE,viewDM=PETSC_FALSE,test=PETSC_FALSE,distribute=PETSC_TRUE,flg;
316: UserCtx user;
317: SNESConvergedReason reason;
319: /* Power subnetwork */
320: UserCtx_Power *appctx_power = &user.appctx_power;
321: char pfdata_file[PETSC_MAX_PATH_LEN] = "power/case9.m";
322: PFDATA *pfdata = NULL;
323: PetscInt genj,loadj,*edgelist_power = NULL,power_netnum;
324: PetscScalar Sbase = 0.0;
326: /* Water subnetwork */
327: AppCtx_Water *appctx_water = &user.appctx_water;
328: WATERDATA *waterdata = NULL;
329: char waterdata_file[PETSC_MAX_PATH_LEN] = "water/sample1.inp";
330: PetscInt *edgelist_water = NULL,water_netnum;
332: /* Shared vertices between subnetworks */
333: PetscInt power_svtx,water_svtx;
335: PetscInitialize(&argc,&argv,"ex1options",help);if (ierr) return ierr;
336: MPI_Comm_rank(PETSC_COMM_WORLD,&rank);
337: MPI_Comm_size(PETSC_COMM_WORLD,&size);
339: /* (1) Read Data - Only rank 0 reads the data */
340: /*--------------------------------------------*/
341: PetscLogStageRegister("Read Data",&stage[0]);
342: PetscLogStagePush(stage[0]);
344: for (i=0; i<Nsubnet; i++) {
345: numVertices[i] = 0;
346: numEdges[i] = 0;
347: }
349: /* All processes READ THE DATA FOR THE FIRST SUBNETWORK: Electric Power Grid */
350: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
351: PetscOptionsGetString(NULL,NULL,"-pfdata",pfdata_file,PETSC_MAX_PATH_LEN-1,NULL);
352: PetscNew(&pfdata);
353: PFReadMatPowerData(pfdata,pfdata_file);
354: if (!rank) {
355: PetscPrintf(PETSC_COMM_SELF,"Power network: nb = %D, ngen = %D, nload = %D, nbranch = %D\n",pfdata->nbus,pfdata->ngen,pfdata->nload,pfdata->nbranch);
356: }
357: Sbase = pfdata->sbase;
358: if (rank == 0) { /* proc[0] will create Electric Power Grid */
359: numEdges[0] = pfdata->nbranch;
360: numVertices[0] = pfdata->nbus;
362: PetscMalloc1(2*numEdges[0],&edgelist_power);
363: GetListofEdges_Power(pfdata,edgelist_power);
364: }
365: /* Broadcast power Sbase to all processors */
366: MPI_Bcast(&Sbase,1,MPIU_SCALAR,0,PETSC_COMM_WORLD);
367: appctx_power->Sbase = Sbase;
368: appctx_power->jac_error = PETSC_FALSE;
369: /* If external option activated. Introduce error in jacobian */
370: PetscOptionsHasName(NULL,NULL, "-jac_error", &appctx_power->jac_error);
372: /* All processes READ THE DATA FOR THE SECOND SUBNETWORK: Water */
373: /* Used for shared vertex, because currently the coupling info must be available in all processes!!! */
374: PetscNew(&waterdata);
375: PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,PETSC_MAX_PATH_LEN-1,NULL);
376: WaterReadData(waterdata,waterdata_file);
377: if (size == 1 || (size > 1 && rank == 1)) {
378: PetscCalloc1(2*waterdata->nedge,&edgelist_water);
379: GetListofEdges_Water(waterdata,edgelist_water);
380: numEdges[1] = waterdata->nedge;
381: numVertices[1] = waterdata->nvertex;
382: }
383: PetscLogStagePop();
385: /* (2) Create a network consist of two subnetworks */
386: /*-------------------------------------------------*/
387: PetscLogStageRegister("Net Setup",&stage[1]);
388: PetscLogStagePush(stage[1]);
390: PetscOptionsGetBool(NULL,NULL,"-viewDM",&viewDM,NULL);
391: PetscOptionsGetBool(NULL,NULL,"-test",&test,NULL);
392: PetscOptionsGetBool(NULL,NULL,"-distribute",&distribute,NULL);
394: /* Create an empty network object */
395: DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);
397: /* Register the components in the network */
398: DMNetworkRegisterComponent(networkdm,"branchstruct",sizeof(struct _p_EDGE_Power),&appctx_power->compkey_branch);
399: DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Power),&appctx_power->compkey_bus);
400: DMNetworkRegisterComponent(networkdm,"genstruct",sizeof(struct _p_GEN),&appctx_power->compkey_gen);
401: DMNetworkRegisterComponent(networkdm,"loadstruct",sizeof(struct _p_LOAD),&appctx_power->compkey_load);
403: DMNetworkRegisterComponent(networkdm,"edge_water",sizeof(struct _p_EDGE_Water),&appctx_water->compkey_edge);
404: DMNetworkRegisterComponent(networkdm,"vertex_water",sizeof(struct _p_VERTEX_Water),&appctx_water->compkey_vtx);
405: #if 0
406: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_branch %d\n",appctx_power->compkey_branch);
407: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_bus %d\n",appctx_power->compkey_bus);
408: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_gen %d\n",appctx_power->compkey_gen);
409: PetscPrintf(PETSC_COMM_WORLD,"power->compkey_load %d\n",appctx_power->compkey_load);
410: PetscPrintf(PETSC_COMM_WORLD,"water->compkey_edge %d\n",appctx_water->compkey_edge);
411: PetscPrintf(PETSC_COMM_WORLD,"water->compkey_vtx %d\n",appctx_water->compkey_vtx);
412: #endif
413: 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]);
414: PetscSynchronizedFlush(PETSC_COMM_WORLD,PETSC_STDOUT);
416: DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,Nsubnet);
417: DMNetworkAddSubnetwork(networkdm,"power",numVertices[0],numEdges[0],edgelist_power,&power_netnum);
418: DMNetworkAddSubnetwork(networkdm,"water",numVertices[1],numEdges[1],edgelist_water,&water_netnum);
420: /* vertex subnet[0].4 shares with vertex subnet[1].0 */
421: power_svtx = 4; water_svtx = 0;
422: DMNetworkAddSharedVertices(networkdm,power_netnum,water_netnum,1,&power_svtx,&water_svtx);
424: /* Set up the network layout */
425: DMNetworkLayoutSetUp(networkdm);
427: /* ADD VARIABLES AND COMPONENTS FOR THE POWER SUBNETWORK */
428: /*-------------------------------------------------------*/
429: genj = 0; loadj = 0;
430: DMNetworkGetSubnetwork(networkdm,power_netnum,&nv,&ne,&vtx,&edges);
432: for (i = 0; i < ne; i++) {
433: DMNetworkAddComponent(networkdm,edges[i],appctx_power->compkey_branch,&pfdata->branch[i],0);
434: }
436: for (i = 0; i < nv; i++) {
437: DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
438: if (flg) continue;
440: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[i],2);
441: if (pfdata->bus[i].ngen) {
442: for (j = 0; j < pfdata->bus[i].ngen; j++) {
443: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_gen,&pfdata->gen[genj++],0);
444: }
445: }
446: if (pfdata->bus[i].nload) {
447: for (j=0; j < pfdata->bus[i].nload; j++) {
448: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[loadj++],0);
449: }
450: }
451: }
453: /* ADD VARIABLES AND COMPONENTS FOR THE WATER SUBNETWORK */
454: /*-------------------------------------------------------*/
455: DMNetworkGetSubnetwork(networkdm,water_netnum,&nv,&ne,&vtx,&edges);
456: for (i = 0; i < ne; i++) {
457: DMNetworkAddComponent(networkdm,edges[i],appctx_water->compkey_edge,&waterdata->edge[i],0);
458: }
460: for (i = 0; i < nv; i++) {
461: DMNetworkIsSharedVertex(networkdm,vtx[i],&flg);
462: if (flg) continue;
464: DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[i],1);
465: }
467: /* ADD VARIABLES AND COMPONENTS AT THE SHARED VERTEX: net[0].4 coupls with net[1].0 -- only the owner of the vertex does this */
468: /*----------------------------------------------------------------------------------------------------------------------------*/
469: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
470: for (i = 0; i < nv; i++) {
471: DMNetworkIsGhostVertex(networkdm,vtx[i],&ghost);
472: /* printf("[%d] coupling info: nv %d; sv[0] %d; ghost %d\n",rank,nv,vtx[0],ghost); */
473: if (ghost) continue;
475: /* power */
476: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_bus,&pfdata->bus[4],2);
477: /* bus[4] is a load, add its component */
478: DMNetworkAddComponent(networkdm,vtx[i],appctx_power->compkey_load,&pfdata->load[0],0);
480: /* water */
481: DMNetworkAddComponent(networkdm,vtx[i],appctx_water->compkey_vtx,&waterdata->vertex[0],1);
482: }
484: /* Set up DM for use */
485: DMSetUp(networkdm);
486: if (viewDM) {
487: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMSetUp, DMView:\n");
488: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
489: }
491: /* Free user objects */
492: PetscFree(edgelist_power);
493: PetscFree(pfdata->bus);
494: PetscFree(pfdata->gen);
495: PetscFree(pfdata->branch);
496: PetscFree(pfdata->load);
497: PetscFree(pfdata);
499: PetscFree(edgelist_water);
500: PetscFree(waterdata->vertex);
501: PetscFree(waterdata->edge);
502: PetscFree(waterdata);
504: /* Re-distribute networkdm to multiple processes for better job balance */
505: if (size >1 && distribute) {
506: DMNetworkDistribute(&networkdm,0);
507: if (viewDM) {
508: PetscPrintf(PETSC_COMM_WORLD,"\nAfter DMNetworkDistribute, DMView:\n");
509: DMView(networkdm,PETSC_VIEWER_STDOUT_WORLD);
510: }
511: }
513: /* Test DMNetworkGetSubnetwork() and DMNetworkGetSubnetworkSharedVertices() */
514: if (test) {
515: PetscInt v,gidx;
516: MPI_Barrier(PETSC_COMM_WORLD);
517: for (i=0; i<Nsubnet; i++) {
518: DMNetworkGetSubnetwork(networkdm,i,&nv,&ne,&vtx,&edges);
519: PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, subnet[%d] ne %d, nv %d\n",rank,i,ne,nv);
520: MPI_Barrier(PETSC_COMM_WORLD);
522: for (v=0; v<nv; v++) {
523: DMNetworkIsGhostVertex(networkdm,vtx[v],&ghost);
524: DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
525: PetscPrintf(PETSC_COMM_SELF,"[%d] subnet[%d] v %d %d; ghost %d\n",rank,i,vtx[v],gidx,ghost);
526: }
527: MPI_Barrier(PETSC_COMM_WORLD);
528: }
529: MPI_Barrier(PETSC_COMM_WORLD);
531: DMNetworkGetSharedVertices(networkdm,&nv,&vtx);
532: PetscPrintf(PETSC_COMM_SELF,"[%d] After distribute, num of shared vertices nsv = %d\n",rank,nv);
533: for (v=0; v<nv; v++) {
534: DMNetworkGetGlobalVertexIndex(networkdm,vtx[v],&gidx);
535: PetscPrintf(PETSC_COMM_SELF,"[%d] sv %d, gidx=%d\n",rank,vtx[v],gidx);
536: }
537: MPI_Barrier(PETSC_COMM_WORLD);
538: }
540: /* Create solution vector X */
541: DMCreateGlobalVector(networkdm,&X);
542: VecDuplicate(X,&F);
543: DMGetLocalVector(networkdm,&user.localXold);
544: PetscLogStagePop();
546: /* (3) Setup Solvers */
547: /*-------------------*/
548: PetscOptionsGetBool(NULL,NULL,"-viewJ",&viewJ,NULL);
549: PetscOptionsGetBool(NULL,NULL,"-viewX",&viewX,NULL);
551: PetscLogStageRegister("SNES Setup",&stage[2]);
552: PetscLogStagePush(stage[2]);
554: SetInitialGuess(networkdm,X,&user);
556: /* Create coupled snes */
557: /*-------------------- */
558: PetscPrintf(PETSC_COMM_WORLD,"SNES_coupled setup ......\n");
559: user.subsnes_id = Nsubnet;
560: SNESCreate(PETSC_COMM_WORLD,&snes);
561: SNESSetDM(snes,networkdm);
562: SNESSetOptionsPrefix(snes,"coupled_");
563: SNESSetFunction(snes,F,FormFunction,&user);
564: SNESMonitorSet(snes,UserMonitor,&user,NULL);
565: SNESSetFromOptions(snes);
567: if (viewJ) {
568: /* View Jac structure */
569: SNESGetJacobian(snes,&Jac,NULL,NULL,NULL);
570: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
571: }
573: if (viewX) {
574: PetscPrintf(PETSC_COMM_WORLD,"Solution:\n");
575: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
576: }
578: if (viewJ) {
579: /* View assembled Jac */
580: MatView(Jac,PETSC_VIEWER_DRAW_WORLD);
581: }
583: /* Create snes_power */
584: /*-------------------*/
585: PetscPrintf(PETSC_COMM_WORLD,"SNES_power setup ......\n");
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");
614: user.subsnes_id = 1;
615: SNESCreate(PETSC_COMM_WORLD,&snes_water);
616: SNESSetDM(snes_water,networkdm);
617: SNESSetOptionsPrefix(snes_water,"water_");
618: SNESSetFunction(snes_water,F,FormFunction,&user);
619: SNESMonitorSet(snes_water,UserMonitor,&user,NULL);
620: SNESSetFromOptions(snes_water);
622: if (viewX) {
623: PetscPrintf(PETSC_COMM_WORLD,"Water Solution:\n");
624: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
625: }
626: PetscLogStagePop();
628: /* (4) Solve */
629: /*-----------*/
630: PetscLogStageRegister("SNES Solve",&stage[3]);
631: PetscLogStagePush(stage[3]);
632: user.it = 0;
633: reason = SNES_DIVERGED_DTOL;
634: while (user.it < it_max && (PetscInt)reason<0) {
635: #if 0
636: user.subsnes_id = 0;
637: SNESSolve(snes_power,NULL,X);
639: user.subsnes_id = 1;
640: SNESSolve(snes_water,NULL,X);
641: #endif
642: user.subsnes_id = Nsubnet;
643: SNESSolve(snes,NULL,X);
645: SNESGetConvergedReason(snes,&reason);
646: user.it++;
647: }
648: PetscPrintf(PETSC_COMM_WORLD,"Coupled_SNES converged in %D iterations\n",user.it);
649: if (viewX) {
650: PetscPrintf(PETSC_COMM_WORLD,"Final Solution:\n");
651: VecView(X,PETSC_VIEWER_STDOUT_WORLD);
652: }
653: PetscLogStagePop();
655: /* Free objects */
656: /* -------------*/
657: VecDestroy(&X);
658: VecDestroy(&F);
659: DMRestoreLocalVector(networkdm,&user.localXold);
661: SNESDestroy(&snes);
662: MatDestroy(&Jac);
663: SNESDestroy(&snes_power);
664: SNESDestroy(&snes_water);
666: DMDestroy(&networkdm);
667: PetscFinalize();
668: return ierr;
669: }
671: /*TEST
673: build:
674: requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
675: depends: power/PFReadData.c power/pffunctions.c water/waterreaddata.c water/waterfunctions.c
677: test:
678: args: -coupled_snes_converged_reason -options_left no -viewDM
679: localrunfiles: ex1options power/case9.m water/sample1.inp
680: output_file: output/ex1.out
682: test:
683: suffix: 2
684: nsize: 3
685: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type parmetis
686: localrunfiles: ex1options power/case9.m water/sample1.inp
687: output_file: output/ex1_2.out
688: requires: parmetis
690: # test:
691: # suffix: 3
692: # nsize: 3
693: # args: -coupled_snes_converged_reason -options_left no -distribute false
694: # localrunfiles: ex1options power/case9.m water/sample1.inp
695: # output_file: output/ex1_2.out
697: test:
698: suffix: 4
699: nsize: 4
700: args: -coupled_snes_converged_reason -options_left no -petscpartitioner_type simple -viewDM
701: localrunfiles: ex1options power/case9.m water/sample1.inp
702: output_file: output/ex1_4.out
704: TEST*/