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