xref: /petsc/src/snes/tutorials/network/water/water.c (revision bef158480efac06de457f7a665168877ab3c2fd7)
1 static char help[] = "This example demonstrates the use of DMNetwork interface for solving a steady-state water network model.\n\
2                       The water network equations follow those used for the package EPANET\n\
3                       The data file format used is from the EPANET package (https://www.epa.gov/water-research/epanet).\n\
4                       Run this program: mpiexec -n <n> ./water\n\\n";
5 
6 /* T
7    Concepts: DMNetwork
8    Concepts: PETSc SNES solver
9 */
10 
11 #include "water.h"
12 #include <petscdmnetwork.h>
13 
14 int main(int argc,char ** argv)
15 {
16   PetscErrorCode   ierr;
17   char             waterdata_file[PETSC_MAX_PATH_LEN]="sample1.inp";
18   WATERDATA        *waterdata;
19   AppCtx_Water     appctx;
20   PetscLogStage    stage1,stage2;
21   PetscMPIInt      crank;
22   DM               networkdm;
23   PetscInt         *edgelist = NULL;
24   PetscInt         nv,ne,i;
25   const PetscInt   *vtx,*edges;
26   Vec              X,F;
27   SNES             snes;
28   SNESConvergedReason reason;
29 
30   ierr = PetscInitialize(&argc,&argv,"wateroptions",help);if (ierr) return ierr;
31   ierr = MPI_Comm_rank(PETSC_COMM_WORLD,&crank);CHKERRQ(ierr);
32 
33   /* Create an empty network object */
34   ierr = DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);CHKERRQ(ierr);
35 
36   /* Register the components in the network */
37   ierr = DMNetworkRegisterComponent(networkdm,"edgestruct",sizeof(struct _p_EDGE_Water),&appctx.compkey_edge);CHKERRQ(ierr);
38   ierr = DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Water),&appctx.compkey_vtx);CHKERRQ(ierr);
39 
40   ierr = PetscLogStageRegister("Read Data",&stage1);CHKERRQ(ierr);
41   PetscLogStagePush(stage1);
42   ierr = PetscNew(&waterdata);CHKERRQ(ierr);
43 
44   /* READ THE DATA */
45   if (!crank) {
46     /* READ DATA. Only rank 0 reads the data */
47     ierr = PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,sizeof(waterdata_file),NULL);CHKERRQ(ierr);
48     ierr = WaterReadData(waterdata,waterdata_file);CHKERRQ(ierr);
49 
50     ierr = PetscCalloc1(2*waterdata->nedge,&edgelist);CHKERRQ(ierr);
51     ierr = GetListofEdges_Water(waterdata,edgelist);CHKERRQ(ierr);
52   }
53   PetscLogStagePop();
54 
55   ierr = PetscLogStageRegister("Create network",&stage2);CHKERRQ(ierr);
56   PetscLogStagePush(stage2);
57 
58   /* Set numbers of nodes and edges */
59   ierr = DMNetworkSetSizes(networkdm,1,&waterdata->nvertex,&waterdata->nedge,0,NULL);CHKERRQ(ierr);
60   if (!crank) {
61     ierr = PetscPrintf(PETSC_COMM_SELF,"water nvertices %D, nedges %D\n",waterdata->nvertex,waterdata->nedge);CHKERRQ(ierr);
62   }
63 
64   /* Add edge connectivity */
65   ierr = DMNetworkSetEdgeList(networkdm,&edgelist,NULL);CHKERRQ(ierr);
66 
67   /* Set up the network layout */
68   ierr = DMNetworkLayoutSetUp(networkdm);CHKERRQ(ierr);
69 
70   if (!crank) {
71     ierr = PetscFree(edgelist);CHKERRQ(ierr);
72   }
73 
74   /* ADD VARIABLES AND COMPONENTS FOR THE NETWORK */
75   ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr);
76 
77   for (i = 0; i < ne; i++) {
78     ierr = DMNetworkAddComponent(networkdm,edges[i],appctx.compkey_edge,&waterdata->edge[i]);CHKERRQ(ierr);
79   }
80 
81   for (i = 0; i < nv; i++) {
82     ierr = DMNetworkAddComponent(networkdm,vtx[i],appctx.compkey_vtx,&waterdata->vertex[i]);CHKERRQ(ierr);
83     /* Add number of variables */
84     ierr = DMNetworkAddNumVariables(networkdm,vtx[i],1);CHKERRQ(ierr);
85   }
86 
87   /* Set up DM for use */
88   ierr = DMSetUp(networkdm);CHKERRQ(ierr);
89 
90   if (!crank) {
91     ierr = PetscFree(waterdata->vertex);CHKERRQ(ierr);
92     ierr = PetscFree(waterdata->edge);CHKERRQ(ierr);
93   }
94   ierr = PetscFree(waterdata);CHKERRQ(ierr);
95 
96   /* Distribute networkdm to multiple processes */
97   ierr = DMNetworkDistribute(&networkdm,0);CHKERRQ(ierr);
98 
99   PetscLogStagePop();
100 
101   ierr = DMCreateGlobalVector(networkdm,&X);CHKERRQ(ierr);
102   ierr = VecDuplicate(X,&F);CHKERRQ(ierr);
103 
104   /* HOOK UP SOLVER */
105   ierr = SNESCreate(PETSC_COMM_WORLD,&snes);CHKERRQ(ierr);
106   ierr = SNESSetDM(snes,networkdm);CHKERRQ(ierr);
107   ierr = SNESSetOptionsPrefix(snes,"water_");CHKERRQ(ierr);
108   ierr = SNESSetFunction(snes,F,WaterFormFunction,NULL);CHKERRQ(ierr);
109   ierr = SNESSetFromOptions(snes);CHKERRQ(ierr);
110 
111   ierr = WaterSetInitialGuess(networkdm,X);CHKERRQ(ierr);
112   /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
113 
114   ierr = SNESSolve(snes,NULL,X);CHKERRQ(ierr);
115   ierr = SNESGetConvergedReason(snes,&reason);CHKERRQ(ierr);
116   if (reason < 0) {
117     SETERRQ(PETSC_COMM_SELF,0,"No solution found for the water network");
118   }
119   /* ierr = VecView(X,PETSC_VIEWER_STDOUT_WORLD);CHKERRQ(ierr); */
120 
121   ierr = VecDestroy(&X);CHKERRQ(ierr);
122   ierr = VecDestroy(&F);CHKERRQ(ierr);
123   ierr = SNESDestroy(&snes);CHKERRQ(ierr);
124   ierr = DMDestroy(&networkdm);CHKERRQ(ierr);
125   ierr = PetscFinalize();
126   return ierr;
127 }
128 
129 /*TEST
130 
131    build:
132       depends: waterreaddata.c waterfunctions.c
133       requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)
134 
135    test:
136       args: -water_snes_converged_reason -options_left no
137       localrunfiles: wateroptions sample1.inp
138       output_file: output/water.out
139       requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
140 
141    test:
142       suffix: 2
143       nsize: 3
144       args: -water_snes_converged_reason -options_left no
145       localrunfiles: wateroptions sample1.inp
146       output_file: output/water.out
147       requires: double !complex define(PETSC_HAVE_ATTRIBUTEALIGNED)
148 
149 TEST*/
150