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