xref: /petsc/src/snes/tutorials/network/water/waterreaddata.c (revision f2c6b1a247e0aba1e6cff92019aae48a2a13617a) 
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 
27   PetscFunctionReturn(0);
28 }
29 
30 PetscErrorCode SetPumpHeadCurveParams(Pump *pump)
31 {
32   SNES                snes;
33   Vec                 X, F;
34   PetscScalar        *head, *flow, *x;
35   SNESConvergedReason reason;
36 
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   }
48 
49   PetscCall(SNESCreate(PETSC_COMM_SELF, &snes));
50 
51   PetscCall(VecCreate(PETSC_COMM_SELF, &X));
52   PetscCall(VecSetSizes(X, 3, 3));
53   PetscCall(VecSetFromOptions(X));
54   PetscCall(VecDuplicate(X, &F));
55 
56   PetscCall(SNESSetFunction(snes, F, PumpHeadCurveResidual, (void *)pump));
57   PetscCall(SNESSetJacobian(snes, NULL, NULL, SNESComputeJacobianDefault, NULL));
58   PetscCall(SNESSetFromOptions(snes));
59 
60   PetscCall(VecGetArray(X, &x));
61   x[0] = head[1];
62   x[1] = 10;
63   x[2] = 3;
64   PetscCall(VecRestoreArray(X, &x));
65 
66   PetscCall(SNESSolve(snes, NULL, X));
67 
68   PetscCall(SNESGetConvergedReason(snes, &reason));
69   PetscCheck(reason >= 0, PETSC_COMM_SELF, PETSC_ERR_CONV_FAILED, "Pump head curve did not converge");
70 
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));
76 
77   PetscCall(VecDestroy(&X));
78   PetscCall(VecDestroy(&F));
79   PetscCall(SNESDestroy(&snes));
80   PetscFunctionReturn(0);
81 }
82 
83 int LineStartsWith(const char *a, const char *b)
84 {
85   if (strncmp(a, b, strlen(b)) == 0) return 1;
86   return 0;
87 }
88 
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 }
96 
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;
104 
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(0);
120 }
121 
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;
138 
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);
144 
145   /* Read file and get line numbers for different data segments */
146   while (fgets(line, MAXLINE, fp)) {
147     if (strstr(line, "[TITLE]")) GetDataSegment(fp, line, &title_start_pos, &ntitle);
148 
149     if (strstr(line, "[JUNCTIONS]")) {
150       GetDataSegment(fp, line, &junc_start_pos, &nlines);
151       water->nvertex += nlines;
152       water->njunction = nlines;
153     }
154 
155     if (strstr(line, "[RESERVOIRS]")) {
156       GetDataSegment(fp, line, &res_start_pos, &nlines);
157       water->nvertex += nlines;
158       water->nreservoir = nlines;
159     }
160 
161     if (strstr(line, "[TANKS]")) {
162       GetDataSegment(fp, line, &tank_start_pos, &nlines);
163       water->nvertex += nlines;
164       water->ntank = nlines;
165     }
166 
167     if (strstr(line, "[PIPES]")) {
168       GetDataSegment(fp, line, &pipe_start_pos, &nlines);
169       water->nedge += nlines;
170       water->npipe = nlines;
171     }
172 
173     if (strstr(line, "[PUMPS]")) {
174       GetDataSegment(fp, line, &pump_start_pos, &nlines);
175       water->nedge += nlines;
176       water->npump = nlines;
177     }
178 
179     if (strstr(line, "[CURVES]")) GetDataSegment(fp, line, &curve_start_pos, &ncurve);
180   }
181 
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;
187 
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   }
205 
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   }
221 
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   }
242 
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   }
269 
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, &paramid);
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   }
286 
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   }
309 
310   fclose(fp);
311 
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(0);
321 }
322