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, &paramid);
    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: }