1 #include "water.h"
2 #include <string.h>
3 #include <ctype.h>
4
PumpHeadCurveResidual(SNES snes,Vec X,Vec F,PetscCtx ctx)5 PetscErrorCode PumpHeadCurveResidual(SNES snes, Vec X, Vec F, PetscCtx 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
SetPumpHeadCurveParams(Pump * pump)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
LineStartsWith(const char * a,const char * b)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
CheckDataSegmentEnd(const char * line)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 */
GetDataSegment(FILE * fp,char * line,fpos_t * data_segment_start_pos,PetscInt * ndatalines)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
WaterReadData(WATERDATA * water,char * filename)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