diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index f6b2f68..1f922fc 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -3,14 +3,9 @@ This document describes the changes and updates that have been made in version 2.3 of EPANET. - - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files can be opened by the toolkit. + - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files to be opened by the toolkit. - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). - A `EN_setvertex` function was added to allow API clients to change the coordinates of a link's vertex. - The index of a General Purpose Valve's (GPV's) head loss curve was added to the list of editable Link Properties using the symbolic constant name `EN_GPV_CURVE`. - The `EN_getlinkvalue` and `EN_setlinkvalue` functions were updated to get and set the value of `EN_GPV_CURVE`. - - For `EN_CUSTOM` type pump curves the maximum head value is now extrapolated to the y-axis intercept instead of being based on the first curve data point. Similarly, the maximum flow value is extrapolated to the x-axis intercept. - - Status checking for a pump not able to deliver enough head has been replaced by adding a penalty term to the pump's operating curve that prevents it from having negative flow (i.e., from crossing the y-axis). - - Status checking for Flow Control Valves has been eliminated by using a continuous head v. flow function. If the current flow is below the valve setting then the normal open head loss relation is used; otherwise a linear penalty function is applied to any flow in excess of the setting. Warnings are no longer issued when the valve operates fully opened at flows below the setting. - - diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index b461946..9f4ac3d 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/07/2020 + Last Updated: 10/04/2019 ****************************************************************************** */ @@ -59,9 +59,9 @@ static void valvecoeff(Project *pr, int k); static void gpvcoeff(Project *pr, int k); static void pbvcoeff(Project *pr, int k); static void tcvcoeff(Project *pr, int k); -static void fcvcoeff(Project *pr, int k); static void prvcoeff(Project *pr, int k, int n1, int n2); static void psvcoeff(Project *pr, int k, int n1, int n2); +static void fcvcoeff(Project *pr, int k, int n1, int n2); void resistcoeff(Project *pr, int k) @@ -152,8 +152,6 @@ void headlosscoeffs(Project *pr) gpvcoeff(pr, k); break; case FCV: - fcvcoeff(pr, k); - break; case PRV: case PSV: if (hyd->LinkSetting[k] == MISSING) valvecoeff(pr, k); @@ -287,8 +285,8 @@ void valvecoeffs(Project *pr) ** Input: none ** Output: none ** Purpose: computes coeffs. of the linearized hydraulic eqns. -** contributed by PRVs & PSVs whose status is not -** fixed to OPEN/CLOSED +** contributed by PRVs, PSVs & FCVs whose status is +** not fixed to OPEN/CLOSED **-------------------------------------------------------------- */ { @@ -315,8 +313,19 @@ void valvecoeffs(Project *pr) n2 = link->N2; // Call valve-specific function - if (link->Type == PRV) prvcoeff(pr, k, n1, n2); - if (link->Type == PSV) psvcoeff(pr, k, n1, n2); + switch (link->Type) + { + case PRV: + prvcoeff(pr, k, n1, n2); + break; + case PSV: + psvcoeff(pr, k, n1, n2); + break; + case FCV: + fcvcoeff(pr, k, n1, n2); + break; + default: continue; + } } } @@ -692,31 +701,20 @@ void pumpcoeff(Project *pr, int k) } // Obtain reference to pump object + q = ABS(hyd->LinkFlow[k]); p = findpump(&pr->network, k); pump = &pr->network.Pump[p]; - - // Prevent negative flow - q = hyd->LinkFlow[k]; - if (q < 0.0) - { - hloss = -(SQR(setting) * pump->Hmax) + CBIG * q; - hgrad = CBIG; - hyd->P[k] = 1.0 / hgrad; - hyd->Y[k] = hloss / hgrad; - return; - } // If no pump curve treat pump as an open valve if (pump->Ptype == NOCURVE) { hyd->P[k] = 1.0 / CSMALL; - hyd->Y[k] = q; + hyd->Y[k] = hyd->LinkFlow[k]; return; } // Get pump curve coefficients for custom pump curve // (Other pump types have pre-determined coeffs.) - q = ABS(q); if (pump->Ptype == CUSTOM) { // Find intercept (h0) & slope (r) of pump curve @@ -1046,10 +1044,12 @@ void psvcoeff(Project *pr, int k, int n1, int n2) } -void fcvcoeff(Project *pr, int k) +void fcvcoeff(Project *pr, int k, int n1, int n2) /* **-------------------------------------------------------------- ** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve ** Output: none ** Purpose: computes solution matrix coeffs. for flow control ** valve @@ -1059,28 +1059,40 @@ void fcvcoeff(Project *pr, int k) Hydraul *hyd = &pr->hydraul; Smatrix *sm = &hyd->smatrix; - double qset; // Valve flow setting - double flow; // Current valve flow - double hloss, hgrad; // Head loss & gradient + int i, j; // Rows in solution matrix + double q; // Valve flow setting - // Treat as a regular valve if status fixed or flow below setting - qset = hyd->LinkSetting[k]; - flow = hyd->LinkFlow[k]; - if (qset == MISSING || hyd->LinkStatus[k] <= CLOSED || flow < qset) + q = hyd->LinkSetting[k]; + i = sm->Row[n1]; + j = sm->Row[n2]; + + // If valve active, break network at valve and treat + // flow setting as external demand at upstream node + // and external supply at downstream node. + + if (hyd->LinkStatus[k] == ACTIVE) { - valvecoeff(pr, k); + hyd->Xflow[n1] -= q; + hyd->Xflow[n2] += q; + hyd->Y[k] = hyd->LinkFlow[k] - q; + sm->F[i] -= q; + sm->F[j] += q; + hyd->P[k] = 1.0 / CBIG; + sm->Aij[sm->Ndx[k]] -= hyd->P[k]; + sm->Aii[i] += hyd->P[k]; + sm->Aii[j] += hyd->P[k]; } - - // Otherwise prevent flow from exceeding the setting + + // Otherwise treat valve as an open pipe + else { - hyd->LinkFlow[k] = qset; valvecoeff(pr, k); - hloss = hyd->Y[k] / hyd->P[k] + CBIG * (flow - qset); - hgrad = CBIG; - hyd->P[k] = 1.0 / hgrad; - hyd->Y[k] = hloss / hgrad; - hyd->LinkFlow[k] = flow; + sm->Aij[sm->Ndx[k]] -= hyd->P[k]; + sm->Aii[i] += hyd->P[k]; + sm->Aii[j] += hyd->P[k]; + sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]); + sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]); } } diff --git a/src/hydsolver.c b/src/hydsolver.c index 53057c2..fcb8258 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -8,7 +8,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/07/2020 + Last Updated: 07/15/2019 ****************************************************************************** */ @@ -111,7 +111,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) maxtrials = hyd->MaxIter; if (hyd->ExtraIter > 0) maxtrials += hyd->ExtraIter; *iter = 1; - headlosscoeffs(pr); while (*iter <= maxtrials) { // Compute coefficient matrices A & F and solve A*H = F @@ -119,6 +118,7 @@ int hydsolve(Project *pr, int *iter, double *relerr) // head loss gradients, & F = flow correction terms. // Solution for H is returned in F from call to linsolve(). + headlosscoeffs(pr); matrixcoeffs(pr); errcode = linsolve(sm, net->Njuncs); @@ -138,9 +138,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) } newerr = newflows(pr, &hydbal); // Update flows *relerr = newerr; - - // Compute head loss coeffs. for new flows - headlosscoeffs(pr); // Write convergence error to status report if called for if (rpt->Statflag == FULL) @@ -246,7 +243,7 @@ int badvalve(Project *pr, int n) if (n == n1 || n == n2) { t = link->Type; - if (t == PRV || t == PSV) + if (t == PRV || t == PSV || t == FCV) { if (hyd->LinkStatus[k] == ACTIVE) { @@ -256,7 +253,8 @@ int badvalve(Project *pr, int n) clocktime(rpt->Atime, time->Htime), link->ID); writeline(pr, pr->Msg); } - hyd->LinkStatus[k] = XPRESSURE; + if (link->Type == FCV) hyd->LinkStatus[k] = XFCV; + else hyd->LinkStatus[k] = XPRESSURE; return 1; } } diff --git a/src/hydstatus.c b/src/hydstatus.c index 76ca823..87d45c6 100644 --- a/src/hydstatus.c +++ b/src/hydstatus.c @@ -7,7 +7,7 @@ Description: updates hydraulic status of network elements Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 02/07/2020 +Last Updated: 05/15/2019 ****************************************************************************** */ @@ -141,6 +141,18 @@ int linkstatus(Project *pr) hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh, hyd->LinkFlow[k]); } + if (link->Type == PUMP && hyd->LinkStatus[k] >= OPEN && + hyd->LinkSetting[k] > 0.0) + { + hyd->LinkStatus[k] = pumpstatus(pr, k, -dh); + } + + // Check for status changes in non-fixed FCVs + if (link->Type == FCV && hyd->LinkSetting[k] != MISSING) + { + hyd->LinkStatus[k] = fcvstatus(pr, k, status, hyd->NodeHead[n1], + hyd->NodeHead[n2]); + } // Check for flow into (out of) full (empty) tanks if (n1 > net->Njuncs || n2 > net->Njuncs) diff --git a/src/input2.c b/src/input2.c index 136d297..a4d7fc8 100644 --- a/src/input2.c +++ b/src/input2.c @@ -7,7 +7,7 @@ Description: reads and interprets network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 02/07/2020 +Last Updated: 02/03/2020 ****************************************************************************** */ @@ -429,14 +429,9 @@ int updatepumpparams(Project *pr, int pumpindex) { if (curve->Y[m] >= curve->Y[m - 1]) return 227; } - pump->Q0 = (curve->X[0] + curve->X[npts-1]) / 2.0; - - // Extend curve to find Hmax (at 0 flow) and Qmax (at 0 head) - b = (curve->Y[1] - curve->Y[0]) / (curve->X[1] - curve->X[0]); - pump->Hmax = curve->Y[0] + b * curve->X[0]; - b = (curve->Y[npts-1] - curve->Y[npts-2]) / - (curve->X[npts-1] - curve->X[npts-2]); - pump->Qmax = curve->X[npts-1] - curve->Y[npts-1] / b; + pump->Qmax = curve->X[npts - 1]; + pump->Q0 = (curve->X[0] + pump->Qmax) / 2.0; + pump->Hmax = curve->Y[0]; } // Compute shape factors & limits of power function curves