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