1 /* function subroutines used by water.c */ 2 3 #include "water.h" 4 #include <petscdmnetwork.h> 5 6 PetscScalar Flow_Pipe(Pipe *pipe,PetscScalar hf,PetscScalar ht) 7 { 8 PetscScalar flow_pipe; 9 10 flow_pipe = PetscSign(hf-ht)*PetscPowScalar(PetscAbsScalar(hf - ht)/pipe->k,(1/pipe->n)); 11 return flow_pipe; 12 } 13 14 PetscScalar Flow_Pump(Pump *pump,PetscScalar hf, PetscScalar ht) 15 { 16 PetscScalar flow_pump; 17 flow_pump = PetscPowScalar((hf - ht + pump->h0)/pump->r,(1/pump->n)); 18 return flow_pump; 19 } 20 21 PetscErrorCode FormFunction_Water(DM networkdm,Vec localX,Vec localF,PetscInt nv,PetscInt ne,const PetscInt* vtx,const PetscInt* edges,void* appctx) 22 { 23 PetscErrorCode ierr; 24 const PetscScalar *xarr; 25 const PetscInt *cone; 26 PetscScalar *farr,hf,ht,flow; 27 PetscInt i,key,vnode1,vnode2,offsetnode1,offsetnode2,offset; 28 PetscBool ghostvtex; 29 VERTEX_Water vertex,vertexnode1,vertexnode2; 30 EDGE_Water edge; 31 Pipe *pipe; 32 Pump *pump; 33 Reservoir *reservoir; 34 Tank *tank; 35 36 PetscFunctionBegin; 37 /* Get arrays for the vectors */ 38 ierr = VecGetArrayRead(localX,&xarr);CHKERRQ(ierr); 39 ierr = VecGetArray(localF,&farr);CHKERRQ(ierr); 40 41 for (i=0; i<ne; i++) { /* for each edge */ 42 /* Get the offset and the key for the component for edge number e[i] */ 43 ierr = DMNetworkGetComponent(networkdm,edges[i],0,&key,(void**)&edge);CHKERRQ(ierr); 44 45 /* Get the numbers for the vertices covering this edge */ 46 ierr = DMNetworkGetConnectedVertices(networkdm,edges[i],&cone);CHKERRQ(ierr); 47 vnode1 = cone[0]; 48 vnode2 = cone[1]; 49 50 /* Get the components at the two vertices */ 51 ierr = DMNetworkGetComponent(networkdm,vnode1,0,&key,(void**)&vertexnode1);CHKERRQ(ierr); 52 ierr = DMNetworkGetComponent(networkdm,vnode2,0,&key,(void**)&vertexnode2);CHKERRQ(ierr); 53 54 /* Get the variable offset (the starting location for the variables in the farr array) for node1 and node2 */ 55 ierr = DMNetworkGetVariableOffset(networkdm,vnode1,&offsetnode1);CHKERRQ(ierr); 56 ierr = DMNetworkGetVariableOffset(networkdm,vnode2,&offsetnode2);CHKERRQ(ierr); 57 58 /* Variables at node1 and node 2 */ 59 hf = xarr[offsetnode1]; 60 ht = xarr[offsetnode2]; 61 62 flow = 0.0; 63 if (edge->type == EDGE_TYPE_PIPE) { 64 pipe = &edge->pipe; 65 flow = Flow_Pipe(pipe,hf,ht); 66 } else if (edge->type == EDGE_TYPE_PUMP) { 67 pump = &edge->pump; 68 flow = Flow_Pump(pump,hf,ht); 69 } 70 /* Convention: Node 1 has outgoing flow and Node 2 has incoming flow */ 71 if (vertexnode1->type == VERTEX_TYPE_JUNCTION) farr[offsetnode1] -= flow; 72 if (vertexnode2->type == VERTEX_TYPE_JUNCTION) farr[offsetnode2] += flow; 73 } 74 75 /* Subtract Demand flows at the vertices */ 76 for (i=0; i<nv; i++) { 77 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr); 78 if(ghostvtex) continue; 79 80 ierr = DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);CHKERRQ(ierr); 81 ierr = DMNetworkGetComponent(networkdm,vtx[i],0,&key,(void**)&vertex);CHKERRQ(ierr); 82 83 if (vertex->type == VERTEX_TYPE_JUNCTION) { 84 farr[offset] -= vertex->junc.demand; 85 } else if (vertex->type == VERTEX_TYPE_RESERVOIR) { 86 reservoir = &vertex->res; 87 farr[offset] = xarr[offset] - reservoir->head; 88 } else if(vertex->type == VERTEX_TYPE_TANK) { 89 tank = &vertex->tank; 90 farr[offset] = xarr[offset] - (tank->elev + tank->initlvl); 91 } 92 } 93 94 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 95 ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); 96 PetscFunctionReturn(0); 97 } 98 99 PetscErrorCode WaterFormFunction(SNES snes,Vec X, Vec F, void *user) 100 { 101 PetscErrorCode ierr; 102 DM networkdm; 103 Vec localX,localF; 104 const PetscInt *v,*e; 105 PetscInt nv,ne; 106 107 PetscFunctionBegin; 108 /* Get the DM attached with the SNES */ 109 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 110 111 /* Get local vertices and edges */ 112 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&v,&e);CHKERRQ(ierr); 113 114 /* Get local vectors */ 115 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 116 ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr); 117 118 /* Scatter values from global to local vector */ 119 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 120 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 121 122 /* Initialize residual */ 123 ierr = VecSet(F,0.0);CHKERRQ(ierr); 124 ierr = VecSet(localF,0.0);CHKERRQ(ierr); 125 126 ierr = FormFunction_Water(networkdm,localX,localF,nv,ne,v,e,NULL);CHKERRQ(ierr); 127 128 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 129 ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 130 ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 131 132 ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr); 133 PetscFunctionReturn(0); 134 } 135 136 PetscErrorCode WaterSetInitialGuess(DM networkdm,Vec X) 137 { 138 PetscErrorCode ierr; 139 PetscInt nv,ne; 140 const PetscInt *vtx,*edges; 141 Vec localX; 142 143 PetscFunctionBegin; 144 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 145 146 ierr = VecSet(X,0.0);CHKERRQ(ierr); 147 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 148 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 149 150 /* Get subnetwork info */ 151 ierr = DMNetworkGetSubnetworkInfo(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 152 ierr = SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);CHKERRQ(ierr); 153 154 ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 155 ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 156 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 157 PetscFunctionReturn(0); 158 } 159 160 PetscErrorCode GetListofEdges_Water(WATERDATA *water,PetscInt *edgelist) 161 { 162 PetscErrorCode ierr; 163 PetscInt i,j,node1,node2; 164 Pipe *pipe; 165 Pump *pump; 166 PetscBool netview=PETSC_FALSE; 167 168 PetscFunctionBegin; 169 ierr = PetscOptionsHasName(NULL,NULL, "-water_view",&netview);CHKERRQ(ierr); 170 for (i=0; i < water->nedge; i++) { 171 if (water->edge[i].type == EDGE_TYPE_PIPE) { 172 pipe = &water->edge[i].pipe; 173 node1 = pipe->node1; 174 node2 = pipe->node2; 175 if (netview) { 176 ierr = PetscPrintf(PETSC_COMM_SELF,"edge %d, pipe v[%d] -> v[%d]\n",i,node1,node2);CHKERRQ(ierr); 177 } 178 } else { 179 pump = &water->edge[i].pump; 180 node1 = pump->node1; 181 node2 = pump->node2; 182 if (netview) { 183 ierr = PetscPrintf(PETSC_COMM_SELF,"edge %d, pump v[%d] -> v[%d]\n",i,node1,node2);CHKERRQ(ierr); 184 } 185 } 186 187 for (j=0; j < water->nvertex; j++) { 188 if (water->vertex[j].id == node1) { 189 edgelist[2*i] = j; 190 break; 191 } 192 } 193 194 for (j=0; j < water->nvertex; j++) { 195 if (water->vertex[j].id == node2) { 196 edgelist[2*i+1] = j; 197 break; 198 } 199 } 200 } 201 PetscFunctionReturn(0); 202 } 203 204 PetscErrorCode SetInitialGuess_Water(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx) 205 { 206 PetscErrorCode ierr; 207 PetscInt i,offset,key; 208 PetscBool ghostvtex; 209 VERTEX_Water vertex; 210 PetscScalar *xarr; 211 212 PetscFunctionBegin; 213 ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr); 214 for (i=0; i < nv; i++) { 215 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr); 216 if (ghostvtex) continue; 217 ierr = DMNetworkGetVariableOffset(networkdm,vtx[i],&offset);CHKERRQ(ierr); 218 ierr = DMNetworkGetComponent(networkdm,vtx[i],0,&key,(void**)&vertex);CHKERRQ(ierr); 219 220 if (vertex->type == VERTEX_TYPE_JUNCTION) { 221 xarr[offset] = 100; 222 } else if (vertex->type == VERTEX_TYPE_RESERVOIR) { 223 xarr[offset] = vertex->res.head; 224 } else { 225 xarr[offset] = vertex->tank.initlvl + vertex->tank.elev; 226 } 227 } 228 ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr); 229 PetscFunctionReturn(0); 230 } 231