xref: /petsc/src/snes/tutorials/network/water/water.c (revision 40badf4fbc550ac1f60bd080eaff6de6d55b946d)
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 #if defined(PETSC_USE_LOG)
21   PetscLogStage    stage1,stage2;
22 #endif
23   PetscMPIInt      crank;
24   DM               networkdm;
25   PetscInt         *edgelist = NULL;
26   PetscInt         nv,ne,i;
27   const PetscInt   *vtx,*edges;
28   Vec              X,F;
29   SNES             snes;
30   SNESConvergedReason reason;
31 
32   ierr = PetscInitialize(&argc,&argv,"wateroptions",help);if (ierr) return ierr;
33   CHKERRMPI(MPI_Comm_rank(PETSC_COMM_WORLD,&crank));
34 
35   /* Create an empty network object */
36   CHKERRQ(DMNetworkCreate(PETSC_COMM_WORLD,&networkdm));
37 
38   /* Register the components in the network */
39   CHKERRQ(DMNetworkRegisterComponent(networkdm,"edgestruct",sizeof(struct _p_EDGE_Water),&appctx.compkey_edge));
40   CHKERRQ(DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Water),&appctx.compkey_vtx));
41 
42   CHKERRQ(PetscLogStageRegister("Read Data",&stage1));
43   CHKERRQ(PetscLogStagePush(stage1));
44   CHKERRQ(PetscNew(&waterdata));
45 
46   /* READ THE DATA */
47   if (!crank) {
48     /* READ DATA. Only rank 0 reads the data */
49     CHKERRQ(PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,sizeof(waterdata_file),NULL));
50     CHKERRQ(WaterReadData(waterdata,waterdata_file));
51 
52     CHKERRQ(PetscCalloc1(2*waterdata->nedge,&edgelist));
53     CHKERRQ(GetListofEdges_Water(waterdata,edgelist));
54   }
55   CHKERRQ(PetscLogStagePop());
56 
57   CHKERRQ(PetscLogStageRegister("Create network",&stage2));
58   CHKERRQ(PetscLogStagePush(stage2));
59 
60   /* Set numbers of nodes and edges */
61   CHKERRQ(DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1));
62   CHKERRQ(DMNetworkAddSubnetwork(networkdm,"",waterdata->nedge,edgelist,NULL));
63   if (!crank) {
64     CHKERRQ(PetscPrintf(PETSC_COMM_SELF,"water nvertices %D, nedges %D\n",waterdata->nvertex,waterdata->nedge));
65   }
66 
67   /* Set up the network layout */
68   CHKERRQ(DMNetworkLayoutSetUp(networkdm));
69 
70   if (!crank) {
71     CHKERRQ(PetscFree(edgelist));
72   }
73 
74   /* ADD VARIABLES AND COMPONENTS FOR THE NETWORK */
75   CHKERRQ(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges));
76 
77   for (i = 0; i < ne; i++) {
78     CHKERRQ(DMNetworkAddComponent(networkdm,edges[i],appctx.compkey_edge,&waterdata->edge[i],0));
79   }
80 
81   for (i = 0; i < nv; i++) {
82     CHKERRQ(DMNetworkAddComponent(networkdm,vtx[i],appctx.compkey_vtx,&waterdata->vertex[i],1));
83   }
84 
85   /* Set up DM for use */
86   CHKERRQ(DMSetUp(networkdm));
87 
88   if (!crank) {
89     CHKERRQ(PetscFree(waterdata->vertex));
90     CHKERRQ(PetscFree(waterdata->edge));
91   }
92   CHKERRQ(PetscFree(waterdata));
93 
94   /* Distribute networkdm to multiple processes */
95   CHKERRQ(DMNetworkDistribute(&networkdm,0));
96 
97   CHKERRQ(PetscLogStagePop());
98 
99   CHKERRQ(DMCreateGlobalVector(networkdm,&X));
100   CHKERRQ(VecDuplicate(X,&F));
101 
102   /* HOOK UP SOLVER */
103   CHKERRQ(SNESCreate(PETSC_COMM_WORLD,&snes));
104   CHKERRQ(SNESSetDM(snes,networkdm));
105   CHKERRQ(SNESSetOptionsPrefix(snes,"water_"));
106   CHKERRQ(SNESSetFunction(snes,F,WaterFormFunction,NULL));
107   CHKERRQ(SNESSetFromOptions(snes));
108 
109   CHKERRQ(WaterSetInitialGuess(networkdm,X));
110   /* CHKERRQ(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
111 
112   CHKERRQ(SNESSolve(snes,NULL,X));
113   CHKERRQ(SNESGetConvergedReason(snes,&reason));
114 
115   PetscCheckFalse(reason < 0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"No solution found for the water network");
116   /* CHKERRQ(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
117 
118   CHKERRQ(VecDestroy(&X));
119   CHKERRQ(VecDestroy(&F));
120   CHKERRQ(SNESDestroy(&snes));
121   CHKERRQ(DMDestroy(&networkdm));
122   ierr = PetscFinalize();
123   return ierr;
124 }
125 
126 /*TEST
127 
128    build:
129       depends: waterreaddata.c waterfunctions.c
130       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
131 
132    test:
133       args: -water_snes_converged_reason -options_left no
134       localrunfiles: wateroptions sample1.inp
135       output_file: output/water.out
136       requires: double !complex defined(PETSC_HAVE_ATTRIBUTEALIGNED)
137 
138    test:
139       suffix: 2
140       nsize: 3
141       args: -water_snes_converged_reason -options_left no
142       localrunfiles: wateroptions sample1.inp
143       output_file: output/water.out
144       requires: double !complex defined(PETSC_HAVE_ATTRIBUTEALIGNED)
145 
146 TEST*/
147