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,ncomp; 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,NULL);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, their variable offsets */ 51 ierr = DMNetworkGetNumComponents(networkdm,vnode1,&ncomp);CHKERRQ(ierr); 52 ierr = DMNetworkGetComponent(networkdm,vnode1,ncomp-1,&key,(void**)&vertexnode1,NULL);CHKERRQ(ierr); 53 ierr = DMNetworkGetLocalVecOffset(networkdm,vnode1,ncomp-1,&offsetnode1);CHKERRQ(ierr); 54 55 ierr = DMNetworkGetNumComponents(networkdm,vnode2,&ncomp);CHKERRQ(ierr); 56 ierr = DMNetworkGetComponent(networkdm,vnode2,ncomp-1,&key,(void**)&vertexnode2,NULL);CHKERRQ(ierr); 57 ierr = DMNetworkGetLocalVecOffset(networkdm,vnode2,ncomp-1,&offsetnode2);CHKERRQ(ierr); 58 59 /* Variables at node1 and node 2 */ 60 hf = xarr[offsetnode1]; 61 ht = xarr[offsetnode2]; 62 63 flow = 0.0; 64 if (edge->type == EDGE_TYPE_PIPE) { 65 pipe = &edge->pipe; 66 flow = Flow_Pipe(pipe,hf,ht); 67 } else if (edge->type == EDGE_TYPE_PUMP) { 68 pump = &edge->pump; 69 flow = Flow_Pump(pump,hf,ht); 70 } 71 /* Convention: Node 1 has outgoing flow and Node 2 has incoming flow */ 72 if (vertexnode1->type == VERTEX_TYPE_JUNCTION) farr[offsetnode1] -= flow; 73 if (vertexnode2->type == VERTEX_TYPE_JUNCTION) farr[offsetnode2] += flow; 74 } 75 76 /* Subtract Demand flows at the vertices */ 77 for (i=0; i<nv; i++) { 78 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr); 79 if (ghostvtex) continue; 80 81 ierr = DMNetworkGetLocalVecOffset(networkdm,vtx[i],ALL_COMPONENTS,&offset);CHKERRQ(ierr); 82 ierr = DMNetworkGetNumComponents(networkdm,vtx[i],&ncomp);CHKERRQ(ierr); 83 ierr = DMNetworkGetComponent(networkdm,vtx[i],ncomp-1,&key,(void**)&vertex,NULL);CHKERRQ(ierr); 84 85 if (vertex->type == VERTEX_TYPE_JUNCTION) { 86 farr[offset] -= vertex->junc.demand; 87 } else if (vertex->type == VERTEX_TYPE_RESERVOIR) { 88 reservoir = &vertex->res; 89 farr[offset] = xarr[offset] - reservoir->head; 90 } else if (vertex->type == VERTEX_TYPE_TANK) { 91 tank = &vertex->tank; 92 farr[offset] = xarr[offset] - (tank->elev + tank->initlvl); 93 } 94 } 95 96 ierr = VecRestoreArrayRead(localX,&xarr);CHKERRQ(ierr); 97 ierr = VecRestoreArray(localF,&farr);CHKERRQ(ierr); 98 PetscFunctionReturn(0); 99 } 100 101 PetscErrorCode WaterFormFunction(SNES snes,Vec X, Vec F, void *user) 102 { 103 PetscErrorCode ierr; 104 DM networkdm; 105 Vec localX,localF; 106 const PetscInt *v,*e; 107 PetscInt nv,ne; 108 109 PetscFunctionBegin; 110 /* Get the DM attached with the SNES */ 111 ierr = SNESGetDM(snes,&networkdm);CHKERRQ(ierr); 112 113 /* Get local vertices and edges */ 114 ierr = DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&v,&e);CHKERRQ(ierr); 115 116 /* Get local vectors */ 117 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 118 ierr = DMGetLocalVector(networkdm,&localF);CHKERRQ(ierr); 119 120 /* Scatter values from global to local vector */ 121 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 122 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 123 124 /* Initialize residual */ 125 ierr = VecSet(F,0.0);CHKERRQ(ierr); 126 ierr = VecSet(localF,0.0);CHKERRQ(ierr); 127 128 ierr = FormFunction_Water(networkdm,localX,localF,nv,ne,v,e,NULL);CHKERRQ(ierr); 129 130 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 131 ierr = DMLocalToGlobalBegin(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 132 ierr = DMLocalToGlobalEnd(networkdm,localF,ADD_VALUES,F);CHKERRQ(ierr); 133 134 ierr = DMRestoreLocalVector(networkdm,&localF);CHKERRQ(ierr); 135 PetscFunctionReturn(0); 136 } 137 138 PetscErrorCode WaterSetInitialGuess(DM networkdm,Vec X) 139 { 140 PetscErrorCode ierr; 141 PetscInt nv,ne; 142 const PetscInt *vtx,*edges; 143 Vec localX; 144 145 PetscFunctionBegin; 146 ierr = DMGetLocalVector(networkdm,&localX);CHKERRQ(ierr); 147 148 ierr = VecSet(X,0.0);CHKERRQ(ierr); 149 ierr = DMGlobalToLocalBegin(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 150 ierr = DMGlobalToLocalEnd(networkdm,X,INSERT_VALUES,localX);CHKERRQ(ierr); 151 152 /* Get subnetwork */ 153 ierr = DMNetworkGetSubnetwork(networkdm,0,&nv,&ne,&vtx,&edges);CHKERRQ(ierr); 154 ierr = SetInitialGuess_Water(networkdm,localX,nv,ne,vtx,edges,NULL);CHKERRQ(ierr); 155 156 ierr = DMLocalToGlobalBegin(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 157 ierr = DMLocalToGlobalEnd(networkdm,localX,ADD_VALUES,X);CHKERRQ(ierr); 158 ierr = DMRestoreLocalVector(networkdm,&localX);CHKERRQ(ierr); 159 PetscFunctionReturn(0); 160 } 161 162 PetscErrorCode GetListofEdges_Water(WATERDATA *water,PetscInt *edgelist) 163 { 164 PetscErrorCode ierr; 165 PetscInt i,j,node1,node2; 166 Pipe *pipe; 167 Pump *pump; 168 PetscBool netview=PETSC_FALSE; 169 170 PetscFunctionBegin; 171 ierr = PetscOptionsHasName(NULL,NULL, "-water_view",&netview);CHKERRQ(ierr); 172 for (i=0; i < water->nedge; i++) { 173 if (water->edge[i].type == EDGE_TYPE_PIPE) { 174 pipe = &water->edge[i].pipe; 175 node1 = pipe->node1; 176 node2 = pipe->node2; 177 if (netview) { 178 ierr = PetscPrintf(PETSC_COMM_SELF,"edge %d, pipe v[%d] -> v[%d]\n",i,node1,node2);CHKERRQ(ierr); 179 } 180 } else { 181 pump = &water->edge[i].pump; 182 node1 = pump->node1; 183 node2 = pump->node2; 184 if (netview) { 185 ierr = PetscPrintf(PETSC_COMM_SELF,"edge %d, pump v[%d] -> v[%d]\n",i,node1,node2);CHKERRQ(ierr); 186 } 187 } 188 189 for (j=0; j < water->nvertex; j++) { 190 if (water->vertex[j].id == node1) { 191 edgelist[2*i] = j; 192 break; 193 } 194 } 195 196 for (j=0; j < water->nvertex; j++) { 197 if (water->vertex[j].id == node2) { 198 edgelist[2*i+1] = j; 199 break; 200 } 201 } 202 } 203 PetscFunctionReturn(0); 204 } 205 206 PetscErrorCode SetInitialGuess_Water(DM networkdm,Vec localX,PetscInt nv,PetscInt ne, const PetscInt *vtx, const PetscInt *edges,void* appctx) 207 { 208 PetscErrorCode ierr; 209 PetscInt i,offset,key; 210 PetscBool ghostvtex,sharedv; 211 VERTEX_Water vertex; 212 PetscScalar *xarr; 213 214 PetscFunctionBegin; 215 ierr = VecGetArray(localX,&xarr);CHKERRQ(ierr); 216 for (i=0; i < nv; i++) { 217 ierr = DMNetworkIsGhostVertex(networkdm,vtx[i],&ghostvtex);CHKERRQ(ierr); 218 ierr = DMNetworkIsSharedVertex(networkdm,vtx[i],&sharedv);CHKERRQ(ierr); 219 if (ghostvtex || sharedv) continue; 220 221 ierr = DMNetworkGetComponent(networkdm,vtx[i],0,&key,(void**)&vertex,NULL);CHKERRQ(ierr); 222 ierr = DMNetworkGetLocalVecOffset(networkdm,vtx[i],0,&offset);CHKERRQ(ierr); 223 if (vertex->type == VERTEX_TYPE_JUNCTION) { 224 xarr[offset] = 100; 225 } else if (vertex->type == VERTEX_TYPE_RESERVOIR) { 226 xarr[offset] = vertex->res.head; 227 } else { 228 xarr[offset] = vertex->tank.initlvl + vertex->tank.elev; 229 } 230 } 231 ierr = VecRestoreArray(localX,&xarr);CHKERRQ(ierr); 232 PetscFunctionReturn(0); 233 } 234