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));
26: PetscFunctionReturn(PETSC_SUCCESS);
27: }
29: PetscErrorCode SetPumpHeadCurveParams(Pump *pump)
30: {
31: SNES snes;
32: Vec X, F;
33: PetscScalar *head, *flow, *x;
34: SNESConvergedReason reason;
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: }
48: PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
50: PetscCall(VecCreate(PETSC_COMM_SELF, &X));
51: PetscCall(VecSetSizes(X, 3, 3));
52: PetscCall(VecSetFromOptions(X));
53: PetscCall(VecDuplicate(X, &F));
55: PetscCall(SNESSetFunction(snes, F, PumpHeadCurveResidual, (void *)pump));
56: PetscCall(SNESSetJacobian(snes, NULL, NULL, SNESComputeJacobianDefault, NULL));
57: PetscCall(SNESSetFromOptions(snes));
59: PetscCall(VecGetArray(X, &x));
60: x[0] = head[1];
61: x[1] = 10;
62: x[2] = 3;
63: PetscCall(VecRestoreArray(X, &x));
65: PetscCall(SNESSolve(snes, NULL, X));
67: PetscCall(SNESGetConvergedReason(snes, &reason));
68: PetscCheck(reason >= 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Pump head curve did not converge");
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));
76: PetscCall(VecDestroy(&X));
77: PetscCall(VecDestroy(&F));
78: PetscCall(SNESDestroy(&snes));
79: PetscFunctionReturn(PETSC_SUCCESS);
80: }
82: int LineStartsWith(const char *a, const char *b)
83: {
84: if (strncmp(a, b, strlen(b)) == 0) return 1;
85: return 0;
86: }
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: }
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;
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: }
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;
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);
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));
148: if (strstr(line, "[JUNCTIONS]")) {
149: PetscCall(GetDataSegment(fp, line, &junc_start_pos, &nlines));
150: water->nvertex += nlines;
151: water->njunction = nlines;
152: }
154: if (strstr(line, "[RESERVOIRS]")) {
155: PetscCall(GetDataSegment(fp, line, &res_start_pos, &nlines));
156: water->nvertex += nlines;
157: water->nreservoir = nlines;
158: }
160: if (strstr(line, "[TANKS]")) {
161: PetscCall(GetDataSegment(fp, line, &tank_start_pos, &nlines));
162: water->nvertex += nlines;
163: water->ntank = nlines;
164: }
166: if (strstr(line, "[PIPES]")) {
167: PetscCall(GetDataSegment(fp, line, &pipe_start_pos, &nlines));
168: water->nedge += nlines;
169: water->npipe = nlines;
170: }
172: if (strstr(line, "[PUMPS]")) {
173: PetscCall(GetDataSegment(fp, line, &pump_start_pos, &nlines));
174: water->nedge += nlines;
175: water->npump = nlines;
176: }
178: if (strstr(line, "[CURVES]")) PetscCall(GetDataSegment(fp, line, &curve_start_pos, &ncurve));
179: }
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;
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: }
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: }
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: }
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: }
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: }
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: }
309: fclose(fp);
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: }