1 #include "water.h" 2 #include <string.h> 3 #include <ctype.h> 4 5 PetscErrorCode PumpHeadCurveResidual(SNES snes,Vec X, Vec F,void *ctx) 6 { 7 const PetscScalar *x; 8 PetscScalar *f; 9 Pump *pump=(Pump*)ctx; 10 PetscScalar *head=pump->headcurve.head,*flow=pump->headcurve.flow; 11 PetscInt i; 12 13 PetscFunctionBegin; 14 PetscCall(VecGetArrayRead(X,&x)); 15 PetscCall(VecGetArray(F,&f)); 16 17 f[0] = f[1] = f[2] = 0; 18 for (i=0; i < pump->headcurve.npt;i++) { 19 f[0] += x[0] - x[1]*PetscPowScalar(flow[i],x[2]) - head[i]; /* Partial w.r.t x[0] */ 20 f[1] += (x[0] - x[1]*PetscPowScalar(flow[i],x[2]) - head[i])*-1*PetscPowScalar(flow[i],x[2]); /*Partial w.r.t x[1] */ 21 f[2] += (x[0] - x[1]*PetscPowScalar(flow[i],x[2]) - head[i])*-1*x[1]*x[2]*PetscPowScalar(flow[i],x[2]-1); /*Partial w.r.t x[2] */ 22 } 23 24 PetscCall(VecRestoreArrayRead(X,&x)); 25 PetscCall(VecRestoreArray(F,&f)); 26 27 PetscFunctionReturn(0); 28 } 29 30 PetscErrorCode SetPumpHeadCurveParams(Pump *pump) 31 { 32 SNES snes; 33 Vec X,F; 34 PetscScalar *head,*flow,*x; 35 SNESConvergedReason reason; 36 37 PetscFunctionBegin; 38 head = pump->headcurve.head; 39 flow = pump->headcurve.flow; 40 if (pump->headcurve.npt == 1) { 41 /* Single point head curve, set the other two data points */ 42 flow[1] = 0; 43 head[1] = 1.33*head[0]; /* 133% of design head -- From EPANET manual */ 44 flow[2] = 2*flow[0]; /* 200% of design flow -- From EPANET manual */ 45 head[2] = 0; 46 pump->headcurve.npt += 2; 47 } 48 49 PetscCall(SNESCreate(PETSC_COMM_SELF,&snes)); 50 51 PetscCall(VecCreate(PETSC_COMM_SELF,&X)); 52 PetscCall(VecSetSizes(X,3,3)); 53 PetscCall(VecSetFromOptions(X)); 54 PetscCall(VecDuplicate(X,&F)); 55 56 PetscCall(SNESSetFunction(snes,F,PumpHeadCurveResidual,(void*)pump)); 57 PetscCall(SNESSetJacobian(snes,NULL,NULL,SNESComputeJacobianDefault,NULL)); 58 PetscCall(SNESSetFromOptions(snes)); 59 60 PetscCall(VecGetArray(X,&x)); 61 x[0] = head[1]; x[1] = 10; x[2] = 3; 62 PetscCall(VecRestoreArray(X,&x)); 63 64 PetscCall(SNESSolve(snes,NULL,X)); 65 66 PetscCall(SNESGetConvergedReason(snes,&reason)); 67 PetscCheck(reason >= 0,PETSC_COMM_SELF,PETSC_ERR_CONV_FAILED,"Pump head curve did not converge"); 68 69 PetscCall(VecGetArray(X,&x)); 70 pump->h0 = x[0]; 71 pump->r = x[1]; 72 pump->n = x[2]; 73 PetscCall(VecRestoreArray(X,&x)); 74 75 PetscCall(VecDestroy(&X)); 76 PetscCall(VecDestroy(&F)); 77 PetscCall(SNESDestroy(&snes)); 78 PetscFunctionReturn(0); 79 } 80 81 int LineStartsWith(const char *a, const char *b) 82 { 83 if (strncmp(a, b, strlen(b)) == 0) return 1; 84 return 0; 85 } 86 87 int CheckDataSegmentEnd(const char *line) 88 { 89 if (LineStartsWith(line,"[JUNCTIONS]") || \ 90 LineStartsWith(line,"[RESERVOIRS]") || \ 91 LineStartsWith(line,"[TANKS]") || \ 92 LineStartsWith(line,"[PIPES]") || \ 93 LineStartsWith(line,"[PUMPS]") || \ 94 LineStartsWith(line,"[CURVES]") || \ 95 LineStartsWith(line,"[VALVES]") || \ 96 LineStartsWith(line,"[PATTERNS]") || \ 97 LineStartsWith(line,"[VALVES]") || \ 98 LineStartsWith(line,"[QUALITY]") || \ 99 LineStartsWith(line,"\n") || LineStartsWith(line,"\r\n")) { 100 return 1; 101 } 102 return 0; 103 } 104 105 /* Gets the file pointer positiion for the start of the data segment and the 106 number of data segments (lines) read 107 */ 108 PetscErrorCode GetDataSegment(FILE *fp,char *line,fpos_t *data_segment_start_pos,PetscInt *ndatalines) 109 { 110 PetscInt data_segment_end; 111 PetscInt nlines=0; 112 113 PetscFunctionBegin; 114 data_segment_end = 0; 115 fgetpos(fp,data_segment_start_pos); 116 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read data segment from file"); 117 while (LineStartsWith(line,";")) { 118 fgetpos(fp,data_segment_start_pos); 119 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read data segment from file"); 120 } 121 while (!data_segment_end) { 122 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read data segment from file"); 123 nlines++; 124 data_segment_end = CheckDataSegmentEnd(line); 125 } 126 *ndatalines = nlines; 127 PetscFunctionReturn(0); 128 } 129 130 PetscErrorCode WaterReadData(WATERDATA *water,char *filename) 131 { 132 FILE *fp=NULL; 133 VERTEX_Water vert; 134 EDGE_Water edge; 135 fpos_t junc_start_pos,res_start_pos,tank_start_pos,pipe_start_pos,pump_start_pos; 136 fpos_t curve_start_pos,title_start_pos; 137 char line[MAXLINE]; 138 PetscInt i,j,nv=0,ne=0,ncurve=0,ntitle=0,nlines,ndata,curve_id; 139 Junction *junction=NULL; 140 Reservoir *reservoir=NULL; 141 Tank *tank=NULL; 142 Pipe *pipe=NULL; 143 Pump *pump=NULL; 144 PetscScalar curve_x,curve_y; 145 double v1,v2,v3,v4,v5,v6; 146 147 PetscFunctionBegin; 148 water->nvertex = water->nedge = 0; 149 fp = fopen(filename,"rb"); 150 /* Check for valid file */ 151 PetscCheck(fp,PETSC_COMM_SELF,PETSC_ERR_FILE_OPEN,"Can't open EPANET data file %s",filename); 152 153 /* Read file and get line numbers for different data segments */ 154 while (fgets(line,MAXLINE,fp)) { 155 156 if (strstr(line,"[TITLE]")) { 157 GetDataSegment(fp,line,&title_start_pos,&ntitle); 158 } 159 160 if (strstr(line,"[JUNCTIONS]")) { 161 GetDataSegment(fp,line,&junc_start_pos,&nlines); 162 water->nvertex += nlines; 163 water->njunction = nlines; 164 } 165 166 if (strstr(line,"[RESERVOIRS]")) { 167 GetDataSegment(fp,line,&res_start_pos,&nlines); 168 water->nvertex += nlines; 169 water->nreservoir = nlines; 170 } 171 172 if (strstr(line,"[TANKS]")) { 173 GetDataSegment(fp,line,&tank_start_pos,&nlines); 174 water->nvertex += nlines; 175 water->ntank = nlines; 176 } 177 178 if (strstr(line,"[PIPES]")) { 179 GetDataSegment(fp,line,&pipe_start_pos,&nlines); 180 water->nedge += nlines; 181 water->npipe = nlines; 182 } 183 184 if (strstr(line,"[PUMPS]")) { 185 GetDataSegment(fp,line,&pump_start_pos,&nlines); 186 water->nedge += nlines; 187 water->npump = nlines; 188 } 189 190 if (strstr(line,"[CURVES]")) { 191 GetDataSegment(fp,line,&curve_start_pos,&ncurve); 192 } 193 } 194 195 /* Allocate vertex and edge data structs */ 196 PetscCall(PetscCalloc1(water->nvertex,&water->vertex)); 197 PetscCall(PetscCalloc1(water->nedge,&water->edge)); 198 vert = water->vertex; 199 edge = water->edge; 200 201 /* Junctions */ 202 fsetpos(fp,&junc_start_pos); 203 for (i=0; i < water->njunction; i++) { 204 int id=0,pattern=0; 205 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read junction from file"); 206 vert[nv].type = VERTEX_TYPE_JUNCTION; 207 junction = &vert[nv].junc; 208 ndata = sscanf(line,"%d %lf %lf %d",&id,&v1,&v2,&pattern);PetscCheck(ndata >= 3,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read junction data"); 209 vert[nv].id = id; 210 junction->dempattern = pattern; 211 junction->elev = (PetscScalar)v1; 212 junction->demand = (PetscScalar)v2; 213 junction->demand *= GPM_CFS; 214 junction->id = vert[nv].id; 215 nv++; 216 } 217 218 /* Reservoirs */ 219 fsetpos(fp,&res_start_pos); 220 for (i=0; i < water->nreservoir; i++) { 221 int id=0,pattern=0; 222 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read reservoir from file"); 223 vert[nv].type = VERTEX_TYPE_RESERVOIR; 224 reservoir = &vert[nv].res; 225 ndata = sscanf(line,"%d %lf %d",&id,&v1,&pattern);PetscCheck(ndata >= 2,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read reservoir data"); 226 vert[nv].id = id; 227 reservoir->headpattern = pattern; 228 reservoir->head = (PetscScalar)v1; 229 reservoir->id = vert[nv].id; 230 nv++; 231 } 232 233 /* Tanks */ 234 fsetpos(fp,&tank_start_pos); 235 for (i=0; i < water->ntank; i++) { 236 int id=0,curve=0; 237 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read data tank from file"); 238 vert[nv].type = VERTEX_TYPE_TANK; 239 tank = &vert[nv].tank; 240 ndata = sscanf(line,"%d %lf %lf %lf %lf %lf %lf %d",&id,&v1,&v2,&v3,&v4,&v5,&v6,&curve);PetscCheck(ndata >= 7,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read tank data"); 241 vert[nv].id = id; 242 tank->volumecurve = curve; 243 tank->elev = (PetscScalar)v1; 244 tank->initlvl = (PetscScalar)v2; 245 tank->minlvl = (PetscScalar)v3; 246 tank->maxlvl = (PetscScalar)v4; 247 tank->diam = (PetscScalar)v5; 248 tank->minvolume = (PetscScalar)v6; 249 tank->id = vert[nv].id; 250 nv++; 251 } 252 253 /* Pipes */ 254 fsetpos(fp,&pipe_start_pos); 255 for (i=0; i < water->npipe; i++) { 256 int id=0,node1=0,node2=0; 257 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read data pipe from file"); 258 edge[ne].type = EDGE_TYPE_PIPE; 259 pipe = &edge[ne].pipe; 260 ndata = sscanf(line,"%d %d %d %lf %lf %lf %lf %s",&id,&node1,&node2,&v1,&v2,&v3,&v4,pipe->stat); 261 pipe->id = id; 262 pipe->node1 = node1; 263 pipe->node2 = node2; 264 pipe->length = (PetscScalar)v1; 265 pipe->diam = (PetscScalar)v2; 266 pipe->roughness = (PetscScalar)v3; 267 pipe->minorloss = (PetscScalar)v4; 268 edge[ne].id = pipe->id; 269 if (strcmp(pipe->stat,"OPEN") == 0) pipe->status = PIPE_STATUS_OPEN; 270 if (ndata < 8) { 271 strcpy(pipe->stat,"OPEN"); /* default OPEN */ 272 pipe->status = PIPE_STATUS_OPEN; 273 } 274 if (ndata < 7) pipe->minorloss = 0.; 275 pipe->n = 1.85; 276 pipe->k = 4.72*pipe->length/(PetscPowScalar(pipe->roughness,pipe->n)*PetscPowScalar(0.0833333*pipe->diam,4.87)); 277 ne++; 278 } 279 280 /* Pumps */ 281 fsetpos(fp,&pump_start_pos); 282 for (i=0; i < water->npump; i++) { 283 int id=0,node1=0,node2=0,paramid=0; 284 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read data pump from file"); 285 edge[ne].type = EDGE_TYPE_PUMP; 286 pump = &edge[ne].pump; 287 ndata = sscanf(line,"%d %d %d %s %d",&id,&node1,&node2,pump->param,¶mid);PetscCheck(ndata == 5,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read pump data"); 288 pump->id = id; 289 pump->node1 = node1; 290 pump->node2 = node2; 291 pump->paramid = paramid; 292 edge[ne].id = pump->id; 293 ne++; 294 } 295 296 /* Curves */ 297 fsetpos(fp,&curve_start_pos); 298 for (i=0; i < ncurve; i++) { 299 int icurve_id=0; 300 PetscCheck(fgets(line,MAXLINE,fp),PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Cannot read data curve from file"); 301 ndata = sscanf(line,"%d %lf %lf",&icurve_id,&v1,&v2);PetscCheck(ndata == 3,PETSC_COMM_SELF,PETSC_ERR_FILE_READ,"Unable to read curve data"); 302 curve_id = icurve_id; 303 curve_x = (PetscScalar)v1; 304 curve_y = (PetscScalar)v2; 305 /* Check for pump with the curve_id */ 306 for (j=water->npipe;j < water->npipe+water->npump;j++) { 307 if (water->edge[j].pump.paramid == curve_id) { 308 PetscCheck(pump->headcurve.npt != 3,PETSC_COMM_SELF,PETSC_ERR_SUP,"Pump %" PetscInt_FMT " [%" PetscInt_FMT " --> %" PetscInt_FMT "]: No support for more than 3-pt head-flow curve",pump->id,pump->node1,pump->node2); 309 pump = &water->edge[j].pump; 310 pump->headcurve.flow[pump->headcurve.npt] = curve_x*GPM_CFS; 311 pump->headcurve.head[pump->headcurve.npt] = curve_y; 312 pump->headcurve.npt++; 313 break; 314 } 315 } 316 } 317 318 fclose(fp); 319 320 /* Get pump curve parameters */ 321 for (j=water->npipe;j < water->npipe+water->npump;j++) { 322 pump = &water->edge[j].pump; 323 if (strcmp(pump->param,"HEAD") == 0) { 324 /* Head-flow curve */ 325 PetscCall(SetPumpHeadCurveParams(pump)); 326 } 327 } 328 PetscFunctionReturn(0); 329 } 330