Actual source code: water.c

      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";

      6: #include "water.h"
      7: #include <petscdmnetwork.h>

      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:   PetscLogStage       stage1, stage2;
     15:   PetscMPIInt         crank;
     16:   DM                  networkdm;
     17:   PetscInt           *edgelist = NULL;
     18:   PetscInt            nv, ne, i;
     19:   const PetscInt     *vtx, *edges;
     20:   Vec                 X, F;
     21:   SNES                snes;
     22:   SNESConvergedReason reason;

     24:   PetscFunctionBeginUser;
     25:   PetscCall(PetscInitialize(&argc, &argv, "wateroptions", help));
     26:   PetscCallMPI(MPI_Comm_rank(PETSC_COMM_WORLD, &crank));

     28:   /* Create an empty network object */
     29:   PetscCall(DMNetworkCreate(PETSC_COMM_WORLD, &networkdm));

     31:   /* Register the components in the network */
     32:   PetscCall(DMNetworkRegisterComponent(networkdm, "edgestruct", sizeof(struct _p_EDGE_Water), &appctx.compkey_edge));
     33:   PetscCall(DMNetworkRegisterComponent(networkdm, "busstruct", sizeof(struct _p_VERTEX_Water), &appctx.compkey_vtx));

     35:   PetscCall(PetscLogStageRegister("Read Data", &stage1));
     36:   PetscCall(PetscLogStagePush(stage1));
     37:   PetscCall(PetscNew(&waterdata));

     39:   /* READ THE DATA */
     40:   if (!crank) {
     41:     /* READ DATA. Only rank 0 reads the data */
     42:     PetscCall(PetscOptionsGetString(NULL, NULL, "-waterdata", waterdata_file, sizeof(waterdata_file), NULL));
     43:     PetscCall(WaterReadData(waterdata, waterdata_file));

     45:     PetscCall(PetscCalloc1(2 * waterdata->nedge, &edgelist));
     46:     PetscCall(GetListofEdges_Water(waterdata, edgelist));
     47:   }
     48:   PetscCall(PetscLogStagePop());

     50:   PetscCall(PetscLogStageRegister("Create network", &stage2));
     51:   PetscCall(PetscLogStagePush(stage2));

     53:   /* Set numbers of nodes and edges */
     54:   PetscCall(DMNetworkSetNumSubNetworks(networkdm, PETSC_DECIDE, 1));
     55:   PetscCall(DMNetworkAddSubnetwork(networkdm, "", waterdata->nedge, edgelist, NULL));
     56:   if (!crank) PetscCall(PetscPrintf(PETSC_COMM_SELF, "water nvertices %" PetscInt_FMT ", nedges %" PetscInt_FMT "\n", waterdata->nvertex, waterdata->nedge));

     58:   /* Set up the network layout */
     59:   PetscCall(DMNetworkLayoutSetUp(networkdm));

     61:   if (!crank) PetscCall(PetscFree(edgelist));

     63:   /* ADD VARIABLES AND COMPONENTS FOR THE NETWORK */
     64:   PetscCall(DMNetworkGetSubnetwork(networkdm, 0, &nv, &ne, &vtx, &edges));

     66:   for (i = 0; i < ne; i++) PetscCall(DMNetworkAddComponent(networkdm, edges[i], appctx.compkey_edge, &waterdata->edge[i], 0));

     68:   for (i = 0; i < nv; i++) PetscCall(DMNetworkAddComponent(networkdm, vtx[i], appctx.compkey_vtx, &waterdata->vertex[i], 1));

     70:   /* Set up DM for use */
     71:   PetscCall(DMSetUp(networkdm));

     73:   if (!crank) {
     74:     PetscCall(PetscFree(waterdata->vertex));
     75:     PetscCall(PetscFree(waterdata->edge));
     76:   }
     77:   PetscCall(PetscFree(waterdata));

     79:   /* Distribute networkdm to multiple processes */
     80:   PetscCall(DMNetworkDistribute(&networkdm, 0));

     82:   PetscCall(PetscLogStagePop());

     84:   PetscCall(DMCreateGlobalVector(networkdm, &X));
     85:   PetscCall(VecDuplicate(X, &F));

     87:   /* HOOK UP SOLVER */
     88:   PetscCall(SNESCreate(PETSC_COMM_WORLD, &snes));
     89:   PetscCall(SNESSetDM(snes, networkdm));
     90:   PetscCall(SNESSetOptionsPrefix(snes, "water_"));
     91:   PetscCall(SNESSetFunction(snes, F, WaterFormFunction, NULL));
     92:   PetscCall(SNESSetFromOptions(snes));

     94:   PetscCall(WaterSetInitialGuess(networkdm, X));
     95:   /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */

     97:   PetscCall(SNESSolve(snes, NULL, X));
     98:   PetscCall(SNESGetConvergedReason(snes, &reason));

    100:   PetscCheck(reason >= 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "No solution found for the water network");
    101:   /* PetscCall(VecView(X,PETSC_VIEWER_STDOUT_WORLD)); */

    103:   PetscCall(VecDestroy(&X));
    104:   PetscCall(VecDestroy(&F));
    105:   PetscCall(SNESDestroy(&snes));
    106:   PetscCall(DMDestroy(&networkdm));
    107:   PetscCall(PetscFinalize());
    108:   return 0;
    109: }

    111: /*TEST

    113:    build:
    114:       depends: waterreaddata.c waterfunctions.c
    115:       requires: !complex double defined(PETSC_HAVE_ATTRIBUTEALIGNED)

    117:    test:
    118:       args: -water_snes_converged_reason -options_left no
    119:       localrunfiles: wateroptions sample1.inp
    120:       output_file: output/water.out
    121:       requires: double !complex defined(PETSC_HAVE_ATTRIBUTEALIGNED)

    123:    test:
    124:       suffix: 2
    125:       nsize: 3
    126:       args: -water_snes_converged_reason -options_left no
    127:       localrunfiles: wateroptions sample1.inp
    128:       output_file: output/water.out
    129:       requires: double !complex defined(PETSC_HAVE_ATTRIBUTEALIGNED)

    131: TEST*/