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