xref: /petsc/src/snes/tutorials/network/water/water.c (revision fbf9dbe564678ed6eff1806adbc4c4f01b9743f4)
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) PetscCall(PetscPrintf(PETSC_COMM_SELF, "water nvertices %" PetscInt_FMT ", nedges %" PetscInt_FMT "\n", waterdata->nvertex, waterdata->nedge));
59 
60   /* Set up the network layout */
61   PetscCall(DMNetworkLayoutSetUp(networkdm));
62 
63   if (!crank) PetscCall(PetscFree(edgelist));
64 
65   /* ADD VARIABLES AND COMPONENTS FOR THE NETWORK */
66   PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));
67 
68   for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx.compkey_edge, &waterdata->edge[i], 0));
69 
70   for (i = 0; i < nv; i++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx.compkey_vtx, &waterdata->vertex[i], 1));
71 
72   /* Set up DM for use */
73   PetscCall(DMSetUp(networkdm));
74 
75   if (!crank) {
76     PetscCall(PetscFree(waterdata->vertex));
77     PetscCall(PetscFree(waterdata->edge));
78   }
79   PetscCall(PetscFree(waterdata));
80 
81   /* Distribute networkdm to multiple processes */
82   PetscCall(DMNetworkDistribute(&networkdm, 0));
83 
84   PetscCall(PetscLogStagePop());
85 
86   PetscCall(DMCreateGlobalVector(networkdm, &X));
87   PetscCall(VecDuplicate(X, &F));
88 
89   /* HOOK UP SOLVER */
90   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
91   PetscCall(SNESSetDM(snes, networkdm));
92   PetscCall(SNESSetOptionsPrefix(snes, "water_"));
93   PetscCall(SNESSetFunction(snes, F, WaterFormFunction, NULL));
94   PetscCall(SNESSetFromOptions(snes));
95 
96   PetscCall(WaterSetInitialGuess(networkdm, X));
97   /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
98 
99   PetscCall(SNESSolve(snes, NULL, X));
100   PetscCall(SNESGetConvergedReason(snes, &reason));
101 
102   PetscCheck(reason >= 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "No solution found for the water network");
103   /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */
104 
105   PetscCall(VecDestroy(&X));
106   PetscCall(VecDestroy(&F));
107   PetscCall(SNESDestroy(&snes));
108   PetscCall(DMDestroy(&networkdm));
109   PetscCall(PetscFinalize());
110   return 0;
111 }
112 
113 /*TEST
114 
115    build:
116       depends: waterreaddata.c waterfunctions.c
117       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)
118 
119    test:
120       args: -water_snes_converged_reason -options_left no -fp_trap 0
121       localrunfiles: wateroptions sample1.inp
122       output_file: output/water.out
123       requires: double !complex defined(PETSC_HAVE_ATTRIBUTEALIGNED)
124 
125    test:
126       suffix: 2
127       nsize: 3
128       args: -water_snes_converged_reason -options_left no -fp_trap 0
129       localrunfiles: wateroptions sample1.inp
130       output_file: output/water.out
131       requires: double !complex defined(PETSC_HAVE_ATTRIBUTEALIGNED)
132 
133 TEST*/
134