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