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