Replaced status checking for pumps & FCVs

See ReleaseNotes2_3.md.

ReleaseNotes2_3.md

@@ -3,9 +3,14 @@
 
 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 to 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 can 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.
+ 
+
 

src/hydcoeffs.c

@@ -7,7 +7,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 10/04/2019
+ Last Updated: 02/07/2020
  ******************************************************************************
 */
 
@@ -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,6 +152,8 @@ 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);
@@ -285,8 +287,8 @@ void  valvecoeffs(Project *pr)
 **   Input:   none
 **   Output:  none
 **   Purpose: computes coeffs. of the linearized hydraulic eqns.
-**            contributed by PRVs, PSVs & FCVs whose status is
+**            contributed by PRVs & PSVs whose status is not
-**            not fixed to OPEN/CLOSED
+**            fixed to OPEN/CLOSED
 **--------------------------------------------------------------
 */
 {
@@ -313,19 +315,8 @@ void  valvecoeffs(Project *pr)
         n2 = link->N2;
 
         // Call valve-specific function
-        switch (link->Type)
+        if (link->Type == PRV) prvcoeff(pr, k, n1, n2);
-        {
+        if (link->Type == PSV) psvcoeff(pr, k, n1, n2);
-        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;
-        }
     }
 }
 
@@ -701,20 +692,31 @@ 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] = hyd->LinkFlow[k];
+        hyd->Y[k] = q;
         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
@@ -1044,12 +1046,10 @@ void  psvcoeff(Project *pr, int k, int n1, int n2)
 }
 
 
-void  fcvcoeff(Project *pr, int k, int n1, int n2)
+void  fcvcoeff(Project *pr, int k)
 /*
 **--------------------------------------------------------------
 **   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,40 +1059,28 @@ void  fcvcoeff(Project *pr, int k, int n1, int n2)
     Hydraul *hyd = &pr->hydraul;
     Smatrix *sm = &hyd->smatrix;
 
-    int   i, j;                   // Rows in solution matrix
+    double qset;                  // Valve flow setting
-    double q;                     // Valve flow setting
+    double flow;                  // Current valve flow
+    double hloss, hgrad;          // Head loss & gradient
 
-    q = hyd->LinkSetting[k];
+    // Treat as a regular valve if status fixed or flow below setting
-    i = sm->Row[n1];
+    qset = hyd->LinkSetting[k];
-    j = sm->Row[n2];
+    flow = hyd->LinkFlow[k];
-
+    if (qset == MISSING  || hyd->LinkStatus[k] <= CLOSED || flow < qset)
-    // 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)
-    {
-        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 treat valve as an open pipe
-
-    else
     {
         valvecoeff(pr, k);
-        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
-        sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]);
+    else
-        sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]);
+    {
+        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;
     }
 }
 

src/hydsolver.c

@@ -8,7 +8,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 07/15/2019
+ Last Updated: 02/07/2020
  ******************************************************************************
 */
 
@@ -111,6 +111,7 @@ 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
@@ -118,7 +119,6 @@ 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);
 
@@ -139,6 +139,9 @@ 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)
         {
@@ -243,7 +246,7 @@ int  badvalve(Project *pr, int n)
         if (n == n1 || n == n2)
         {
             t = link->Type;
-            if (t == PRV || t == PSV || t == FCV)
+            if (t == PRV || t == PSV)
             {
                 if (hyd->LinkStatus[k] == ACTIVE)
                 {
@@ -253,8 +256,7 @@ int  badvalve(Project *pr, int n)
                                 clocktime(rpt->Atime, time->Htime), link->ID);
                         writeline(pr, pr->Msg);
                     }
-                    if (link->Type == FCV) hyd->LinkStatus[k] = XFCV;
+                    hyd->LinkStatus[k] = XPRESSURE;
-                    else                   hyd->LinkStatus[k] = XPRESSURE;
                     return 1;
                 }
             }

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: 05/15/2019
+Last Updated: 02/07/2020
 ******************************************************************************
 */
 
@@ -141,18 +141,6 @@ 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)

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/03/2020
+Last Updated: 02/07/2020
 ******************************************************************************
 */
 
@@ -429,9 +429,14 @@ int updatepumpparams(Project *pr, int pumpindex)
             {
                 if (curve->Y[m] >= curve->Y[m - 1]) return 227;
             }
-            pump->Qmax = curve->X[npts - 1];
+            pump->Q0 = (curve->X[0] + curve->X[npts-1]) / 2.0;
-            pump->Q0 = (curve->X[0] + pump->Qmax) / 2.0;
+            
-            pump->Hmax = curve->Y[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;
         }
 
         // Compute shape factors & limits of power function curves