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