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: /* T
      7:    Concepts: DMNetwork
      8:    Concepts: PETSc SNES solver
      9: */

     11: #include "water.h"
     12: #include <petscdmnetwork.h>

     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;

     32:   PetscInitialize(&argc,&argv,"wateroptions",help);if (ierr) return ierr;
     33:   MPI_Comm_rank(PETSC_COMM_WORLD,&crank);

     35:   /* Create an empty network object */
     36:   DMNetworkCreate(PETSC_COMM_WORLD,&networkdm);

     38:   /* Register the components in the network */
     39:   DMNetworkRegisterComponent(networkdm,"edgestruct",sizeof(struct _p_EDGE_Water),&appctx.compkey_edge);
     40:   DMNetworkRegisterComponent(networkdm,"busstruct",sizeof(struct _p_VERTEX_Water),&appctx.compkey_vtx);

     42:   PetscLogStageRegister("Read Data",&stage1);
     43:   PetscLogStagePush(stage1);
     44:   PetscNew(&waterdata);

     46:   /* READ THE DATA */
     47:   if (!crank) {
     48:     /* READ DATA. Only rank 0 reads the data */
     49:     PetscOptionsGetString(NULL,NULL,"-waterdata",waterdata_file,sizeof(waterdata_file),NULL);
     50:     WaterReadData(waterdata,waterdata_file);

     52:     PetscCalloc1(2*waterdata->nedge,&edgelist);
     53:     GetListofEdges_Water(waterdata,edgelist);
     54:   }
     55:   PetscLogStagePop();

     57:   PetscLogStageRegister("Create network",&stage2);
     58:   PetscLogStagePush(stage2);

     60:   /* Set numbers of nodes and edges */
     61:   DMNetworkSetNumSubNetworks(networkdm,PETSC_DECIDE,1);
     62:   DMNetworkAddSubnetwork(networkdm,"",waterdata->nvertex,waterdata->nedge,edgelist,NULL);
     63:   if (!crank) {
     64:     PetscPrintf(PETSC_COMM_SELF,"water nvertices %D, nedges %D\n",waterdata->nvertex,waterdata->nedge);
     65:   }

     67:   /* Set up the network layout */
     68:   DMNetworkLayoutSetUp(networkdm);

     70:   if (!crank) {
     71:     PetscFree(edgelist);
     72:   }

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

     77:   for (i = 0; i < ne; i++) {
     78:     DMNetworkAddComponent(networkdm,edges[i],appctx.compkey_edge,&waterdata->edge[i],0);
     79:   }

     81:   for (i = 0; i < nv; i++) {
     82:     DMNetworkAddComponent(networkdm,vtx[i],appctx.compkey_vtx,&waterdata->vertex[i],1);
     83:   }

     85:   /* Set up DM for use */
     86:   DMSetUp(networkdm);

     88:   if (!crank) {
     89:     PetscFree(waterdata->vertex);
     90:     PetscFree(waterdata->edge);
     91:   }
     92:   PetscFree(waterdata);

     94:   /* Distribute networkdm to multiple processes */
     95:   DMNetworkDistribute(&networkdm,0);

     97:   PetscLogStagePop();

     99:   DMCreateGlobalVector(networkdm,&X);
    100:   VecDuplicate(X,&F);

    102:   /* HOOK UP SOLVER */
    103:   SNESCreate(PETSC_COMM_WORLD,&snes);
    104:   SNESSetDM(snes,networkdm);
    105:   SNESSetOptionsPrefix(snes,"water_");
    106:   SNESSetFunction(snes,F,WaterFormFunction,NULL);
    107:   SNESSetFromOptions(snes);

    109:   WaterSetInitialGuess(networkdm,X);
    110:   /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

    112:   SNESSolve(snes,NULL,X);
    113:   SNESGetConvergedReason(snes,&reason);

    115:   if (reason < 0) SETERRQ(PETSC_COMM_SELF,0,"No solution found for the water network");
    116:   /* VecView(X,PETSC_VIEWER_STDOUT_WORLD); */

    118:   VecDestroy(&X);
    119:   VecDestroy(&F);
    120:   SNESDestroy(&snes);
    121:   DMDestroy(&networkdm);
    122:   PetscFinalize();
    123:   return ierr;
    124: }

    126: /*TEST

    128:    build:
    129:       depends: waterreaddata.c waterfunctions.c
    130:       requires: !complex double define(PETSC_HAVE_ATTRIBUTEALIGNED)

    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 define(PETSC_HAVE_ATTRIBUTEALIGNED)

    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 define(PETSC_HAVE_ATTRIBUTEALIGNED)

    146: TEST*/