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 const PetscScalar *xarr; 24 const PetscInt *cone; 25 PetscScalar *farr,hf,ht,flow; 26 PetscInt i,key,vnode1,vnode2,offsetnode1,offsetnode2,offset,ncomp; 27 PetscBool ghostvtex; 28 VERTEX_Water vertex,vertexnode1,vertexnode2; 29 EDGE_Water edge; 30 Pipe *pipe; 31 Pump *pump; 32 Reservoir *reservoir; 33 Tank *tank; 34 35 PetscFunctionBegin; 36 /* Get arrays for the vectors */ 37 PetscCall(VecGetArrayRead(localX,&xarr)); 38 PetscCall(VecGetArray(localF,&farr)); 39 40 for (i=0; i<ne; i++) { /* for each edge */ 41 /* Get the offset and the key for the component for edge number e[i] */ 42 PetscCall(DMNetworkGetComponent(networkdm,edges[i],0,&key,(void**)&edge,NULL)); 43 44 /* Get the numbers for the vertices covering this edge */ 45 PetscCall(DMNetworkGetConnectedVertices(networkdm,edges[i],&cone)); 46 vnode1 = cone[0]; 47 vnode2 = cone[1]; 48 49 /* Get the components at the two vertices, their variable offsets */ 50 PetscCall(DMNetworkGetNumComponents(networkdm,vnode1,&ncomp)); 51 PetscCall(DMNetworkGetComponent(networkdm,vnode1,ncomp-1,&key,(void**)&vertexnode1,NULL)); 52 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vnode1,ncomp-1,&offsetnode1)); 53 54 PetscCall(DMNetworkGetNumComponents(networkdm,vnode2,&ncomp)); 55 PetscCall(DMNetworkGetComponent(networkdm,vnode2,ncomp-1,&key,(void**)&vertexnode2,NULL)); 56 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vnode2,ncomp-1,&offsetnode2)); 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 PetscCall(DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex)); 78 if (ghostvtex) continue; 79 80 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset)); 81 PetscCall(DMNetworkGetNumComponents(networkdm,vtx[i],&ncomp)); 82 PetscCall(DMNetworkGetComponent(networkdm,vtx[i],ncomp-1,&key,(void**)&vertex,NULL)); 83 84 if (vertex->type == VERTEX_TYPE_JUNCTION) { 85 farr[offset] -= vertex->junc.demand; 86 } else if (vertex->type == VERTEX_TYPE_RESERVOIR) { 87 reservoir = &vertex->res; 88 farr[offset] = xarr[offset] - reservoir->head; 89 } else if (vertex->type == VERTEX_TYPE_TANK) { 90 tank = &vertex->tank; 91 farr[offset] = xarr[offset] - (tank->elev + tank->initlvl); 92 } 93 } 94 95 PetscCall(VecRestoreArrayRead(localX,&xarr)); 96 PetscCall(VecRestoreArray(localF,&farr)); 97 PetscFunctionReturn(0); 98 } 99 100 PetscErrorCode WaterFormFunction(SNES snes,Vec X, Vec F, void *user) 101 { 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 PetscCall(SNESGetDM(snes,&networkdm)); 110 111 /* Get local vertices and edges */ 112 PetscCall(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&v,&e)); 113 114 /* Get local vectors */ 115 PetscCall(DMGetLocalVector(networkdm,&localX)); 116 PetscCall(DMGetLocalVector(networkdm,&localF)); 117 118 /* Scatter values from global to local vector */ 119 PetscCall(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 120 PetscCall(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 121 122 /* Initialize residual */ 123 PetscCall(VecSet(F,0.0)); 124 PetscCall(VecSet(localF,0.0)); 125 126 PetscCall(FormFunction_Water(networkdm,localX,localF,nv,ne,v,e,NULL)); 127 128 PetscCall(DMRestoreLocalVector(networkdm,&localX)); 129 PetscCall(DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F)); 130 PetscCall(DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F)); 131 132 PetscCall(DMRestoreLocalVector(networkdm,&localF)); 133 PetscFunctionReturn(0); 134 } 135 136 PetscErrorCode WaterSetInitialGuess(DM networkdm,Vec X) 137 { 138 PetscInt nv,ne; 139 const PetscInt *vtx,*edges; 140 Vec localX; 141 142 PetscFunctionBegin; 143 PetscCall(DMGetLocalVector(networkdm,&localX)); 144 145 PetscCall(VecSet(X,0.0)); 146 PetscCall(DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX)); 147 PetscCall(DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX)); 148 149 /* Get subnetwork */ 150 PetscCall(DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges)); 151 PetscCall(SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL)); 152 153 PetscCall(DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X)); 154 PetscCall(DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X)); 155 PetscCall(DMRestoreLocalVector(networkdm,&localX)); 156 PetscFunctionReturn(0); 157 } 158 159 PetscErrorCode GetListofEdges_Water(WATERDATA *water,PetscInt *edgelist) 160 { 161 PetscInt i,j,node1,node2; 162 Pipe *pipe; 163 Pump *pump; 164 PetscBool netview=PETSC_FALSE; 165 166 PetscFunctionBegin; 167 PetscCall(PetscOptionsHasName(NULL,NULL, "-water_view",&netview)); 168 for (i=0; i < water->nedge; i++) { 169 if (water->edge[i].type == EDGE_TYPE_PIPE) { 170 pipe = &water->edge[i].pipe; 171 node1 = pipe->node1; 172 node2 = pipe->node2; 173 if (netview) { 174 PetscCall(PetscPrintf(PETSC_COMM_SELF,"edge %" PetscInt_FMT ", pipe v[%" PetscInt_FMT "] -> v[%" PetscInt_FMT "]\n",i,node1,node2)); 175 } 176 } else { 177 pump = &water->edge[i].pump; 178 node1 = pump->node1; 179 node2 = pump->node2; 180 if (netview) { 181 PetscCall(PetscPrintf(PETSC_COMM_SELF,"edge %" PetscInt_FMT ", pump v[%" PetscInt_FMT "] -> v[%" PetscInt_FMT "]\n",i,node1,node2)); 182 } 183 } 184 185 for (j=0; j < water->nvertex; j++) { 186 if (water->vertex[j].id == node1) { 187 edgelist[2*i] = j; 188 break; 189 } 190 } 191 192 for (j=0; j < water->nvertex; j++) { 193 if (water->vertex[j].id == node2) { 194 edgelist[2*i+1] = j; 195 break; 196 } 197 } 198 } 199 PetscFunctionReturn(0); 200 } 201 202 PetscErrorCode SetInitialGuess_Water(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx) 203 { 204 PetscInt i,offset,key; 205 PetscBool ghostvtex,sharedv; 206 VERTEX_Water vertex; 207 PetscScalar *xarr; 208 209 PetscFunctionBegin; 210 PetscCall(VecGetArray(localX,&xarr)); 211 for (i=0; i < nv; i++) { 212 PetscCall(DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex)); 213 PetscCall(DMNetworkIsSharedVertex(networkdm,vtx[i],&sharedv)); 214 if (ghostvtex || sharedv) continue; 215 216 PetscCall(DMNetworkGetComponent(networkdm,vtx[i],0,&key,(void**)&vertex,NULL)); 217 PetscCall(DMNetworkGetLocalVecOffset(networkdm,vtx[i],0,&offset)); 218 if (vertex->type == VERTEX_TYPE_JUNCTION) { 219 xarr[offset] = 100; 220 } else if (vertex->type == VERTEX_TYPE_RESERVOIR) { 221 xarr[offset] = vertex->res.head; 222 } else { 223 xarr[offset] = vertex->tank.initlvl + vertex->tank.elev; 224 } 225 } 226 PetscCall(VecRestoreArray(localX,&xarr)); 227 PetscFunctionReturn(0); 228 } 229