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