Merge pull request #282 from LRossman/lrossman-dev

Fixes zero flow (#164) and D-W eqn. (#199) issues

README.md

@@ -18,4 +18,3 @@ A step-by-step tutorial on how to contribute to EPANET using GitHub is also [ava
 __Note:__ This repository is not affiliated with, or endorsed by, the USEPA. For the last "official" release of EPANET (2.00.12 UI and Toolkit) please go to the [EPA's GitHub repo](https://github.com/USEPA/Water-Distribution-Network-Model) or [the USEPA website](http://www2.epa.gov/water-research/epanet). It is also not the graphical user interface version. This is the hydraulic and water quality solver engine.
 
 However, if you are interested in extending EPANET for academic, personal, or commercial use, then you've come to the right place. For community discussion, FAQ, and roadmapping of the project, go to the [Community Forum](http://community.wateranalytics.org/category/epanet). 
-

src/funcs.h

@@ -178,11 +178,12 @@ int     hydsolve(EN_Project *pr, int *,double *);   /* Solves network equations
 void    resistcoeff(EN_Project *pr, int k);         /* Finds pipe flow resistance */
 void    headlosscoeffs(EN_Project *pr);             // Finds link head loss coeffs.
 void    matrixcoeffs(EN_Project *pr);               /* Finds hyd. matrix coeffs.  */
-double  emitflowchange(EN_Project *pr, int i);      /* Change in emitter outflow  */
+void    emitheadloss(EN_Project *pr, int,           // Finds emitter head loss
-double  demandflowchange(EN_Project *pr, int i,     // Change in demand outflow
+                     double *, double *);           
-                         double dp, double n);
+double  demandflowchange(EN_Project *pr, int,       // Change in demand outflow
-void    demandparams(EN_Project *pr, double *dp,    // PDA function parameters
+                         double, double);
-                     double *n); 
+void    demandparams(EN_Project *pr, double *,      // PDA function parameters
+                     double *); 
 
 /* ----------- SMATRIX.C ---------------*/
 int     createsparse(EN_Project *pr);               /* Creates sparse matrix      */

src/hydcoeffs.c

@@ -18,21 +18,21 @@ HYDCOEFFS.C --  hydraulic coefficients for the EPANET Program
 #include "funcs.h"
 
 // Constants used for computing Darcy-Weisbach friction factor
-const double A1 = 0.314159265359e04;   // 1000*PI
+const double A1 = 3.14159265358979323850e+03;   // 1000*PI
-const double A2 = 0.157079632679e04;   // 500*PI
+const double A2 = 1.57079632679489661930e+03;   // 500*PI
-const double A3 = 0.502654824574e02;   // 16*PI
+const double A3 = 5.02654824574366918160e+01;   // 16*PI
-const double A4 = 6.283185307;         // 2*PI
+const double A4 = 6.28318530717958647700e+00;   // 2*PI
-const double A8 = 4.61841319859;       // 5.74*(PI/4)^.9
+const double A8 = 4.61841319859066668690e+00;   // 5.74*(PI/4)^.9
-const double A9 = -8.685889638e-01;    // -2/ln(10)
+const double A9 = -8.68588963806503655300e-01;  // -2/ln(10)
-const double AA = -1.5634601348;       // -2*.9*2/ln(10)
+const double AA = -1.5634601348517065795e+00;   // -2*.9*2/ln(10)
-const double AB = 3.28895476345e-03;   // 5.74/(4000^.9)
+const double AB = 3.28895476345399058690e-03;   // 5.74/(4000^.9)
-const double AC = -5.14214965799e-03;  // AA*AB
+const double AC = -5.14214965799093883760e-03;  // AA*AB
 
 // External functions
 //void   resistcoeff(EN_Project *pr, int k);
 //void   headlosscoeffs(EN_Project *pr);
 //void   matrixcoeffs(EN_Project *pr);
-//double emitflowchange(EN_Project *pr, int i);
+//void   emitheadloss(EN_Project *pr, int i, double *hloss, double *dhdq);
 //double demandflowchange(EN_Project *pr, int i, double dp, double n);
 //void   demandparams(EN_Project *pr, double *dp, double *n);
 
@@ -47,7 +47,7 @@ static void    demandheadloss(double d, double dfull, double dp,
 
 static void    pipecoeff(EN_Project *pr, int k);
 static void    DWpipecoeff(EN_Project *pr, int k);
-static double  frictionFactor(EN_Project *pr, int k, double *dfdq);
+static double  frictionFactor(double q, double e, double s, double *dfdq);
 
 static void    pumpcoeff(EN_Project *pr, int k);
 static void    curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r);
@@ -71,10 +71,12 @@ void  resistcoeff(EN_Project *pr, int k)
 */
 {
     double e, d, L;
+
     EN_Network   *net = &pr->network;
     hydraulics_t *hyd = &pr->hydraulics;
     Slink *link = &net->Link[k];
     
+    link->Qa = 0.0;
     switch (link->Type) {
 
     // ... Link is a pipe. Compute resistance based on headloss formula.
@@ -88,14 +90,20 @@ void  resistcoeff(EN_Project *pr, int k)
         switch (hyd->Formflag)
         {
         case HW:
+            // ... resistance coeff.
             link->R = 4.727 * L / pow(e, hyd->Hexp) / pow(d, 4.871);
+            // ... flow below which linear head loss applies
+            link->Qa = pow(hyd->RQtol / hyd->Hexp / link->R, 1.17371);
             break;
         case DW:
             link->R = L / 2.0 / 32.2 / d / SQR(PI * SQR(d) / 4.0);
             break;
         case CM:
+            // ... resistance coeff.
             link->R = SQR(4.0 * e / (1.49 * PI * SQR(d))) *
                       pow((d / 4.0), -1.333) * L;
+            // ... flow below which linear head loss applies
+            link->Qa = hyd->RQtol / 2.0 / link->R;
         }
         break;
 
@@ -195,7 +203,8 @@ void  linkcoeffs(EN_Project *pr)
 **--------------------------------------------------------------
 **   Input:   none
 **   Output:  none
-**   Purpose: computes matrix coefficients for links
+**   Purpose: computes coefficients contributed by links to the
+**            linearized system of hydraulic equations.
 **--------------------------------------------------------------
 */
 {
@@ -214,32 +223,34 @@ void  linkcoeffs(EN_Project *pr)
         n1 = link->N1;           // Start node of link
         n2 = link->N2;           // End node of link
 
-        // Update net nodal inflows (X_tmp), solution matrix (A) and RHS array (F)
+        // Update nodal flow balance (X_tmp)
-        // (Use covention that flow out of node is (-), flow into node is (+))
+        // (Flow out of node is (-), flow into node is (+))
         hyd->X_tmp[n1] -= hyd->LinkFlows[k];
         hyd->X_tmp[n2] += hyd->LinkFlows[k];
 
-        // Off-diagonal coeff.
+        // Add to off-diagonal coeff. of linear system matrix
         sol->Aij[sol->Ndx[k]] -= sol->P[k];
 
-        // Node n1 is junction
+        // Update linear system coeffs. associated with start node n1
+        // ... node n1 is junction
         if (n1 <= net->Njuncs)                     
         {
             sol->Aii[sol->Row[n1]] += sol->P[k];   // Diagonal coeff.
             sol->F[sol->Row[n1]] += sol->Y[k];     // RHS coeff.
         }
 
-        // Node n1 is a tank/reservoir
+        // ... node n1 is a tank/reservoir
         else sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); 
 
-        // Node n2 is junction
+        // Update linear system coeffs. associated with end node n2
+        // ... node n2 is junction
         if (n2 <= net->Njuncs)
         {
             sol->Aii[sol->Row[n2]] += sol->P[k];   // Diagonal coeff.
             sol->F[sol->Row[n2]] -= sol->Y[k];     // RHS coeff.
         }
 
-        // Node n2 is a tank/reservoir
+        // ... node n2 is a tank/reservoir
         else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); 
     }
 }
@@ -250,8 +261,8 @@ void  nodecoeffs(EN_Project *pr)
 **----------------------------------------------------------------
 **  Input:   none
 **  Output:  none
-**  Purpose: completes calculation of nodal flow imbalance (X_tmp)
+**  Purpose: completes calculation of nodal flow balance array
-**           & flow correction (F) arrays
+**           (X_tmp) & r.h.s. (F) of linearized hydraulic eqns.
 **----------------------------------------------------------------
 */
 {
@@ -261,7 +272,7 @@ void  nodecoeffs(EN_Project *pr)
     EN_Network   *net = &pr->network;
 
     // For junction nodes, subtract demand flow from net
-    // flow imbalance & add imbalance to RHS array F.
+    // flow balance & add flow balance to RHS array F
     for (i = 1; i <= net->Njuncs; i++)
     {
         hyd->X_tmp[i] -= hyd->DemandFlows[i];
@@ -275,8 +286,9 @@ void  valvecoeffs(EN_Project *pr)
 **--------------------------------------------------------------
 **   Input:   none
 **   Output:  none
-**   Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs
+**   Purpose: computes coeffs. of the linearized hydraulic eqns.
-**            whose status is not fixed to OPEN/CLOSED
+**            contributed by PRVs, PSVs & FCVs whose status is
+**            not fixed to OPEN/CLOSED
 **--------------------------------------------------------------
 */
 {
@@ -325,7 +337,8 @@ void  emittercoeffs(EN_Project *pr)
 **--------------------------------------------------------------
 **   Input:   none
 **   Output:  none
-**   Purpose: computes matrix coeffs. for emitters
+**   Purpose: computes coeffs. of the linearized hydraulic eqns.
+**            contributed by emitters.
 **
 **   Note: Emitters consist of a fictitious pipe connected to
 **         a fictitious reservoir whose elevation equals that
@@ -334,12 +347,8 @@ void  emittercoeffs(EN_Project *pr)
 **--------------------------------------------------------------
 */
 {
-    int   i;
+    int     i, row;
-    double  ke;
+    double  hloss, hgrad;
-    double  p;
-    double  q;
-    double  y;
-    double  z;
 
     hydraulics_t *hyd = &pr->hydraulics;
     solver_t     *sol = &hyd->solver;
@@ -348,46 +357,61 @@ void  emittercoeffs(EN_Project *pr)
 
     for (i = 1; i <= net->Njuncs; i++)
     {
+        // Skip junctions without emitters
         node = &net->Node[i];
         if (node->Ke == 0.0) continue;
-        ke = MAX(CSMALL, node->Ke);              // emitter coeff.
+
-        q = hyd->EmitterFlows[i];                // emitter flow
+        // Find emitter head loss and gradient
-        z = ke * pow(ABS(q), hyd->Qexp);         // emitter head loss
+        emitheadloss(pr, i, &hloss, &hgrad);
-        p = hyd->Qexp * z / ABS(q);              // head loss gradient
+
-        if (p < hyd->RQtol)
+        // Row of solution matrix
-        {
+        row = sol->Row[i];
-            p = 1.0 / hyd->RQtol;
+
-        }
+        // Addition to matrix diagonal & r.h.s
-        else
+        sol->Aii[row] += 1.0 / hgrad;
-        {
+        sol->F[row] += (hloss + node->El) / hgrad;
-            p = 1.0 / p;                         // inverse head loss gradient
+
-        }
+        // Update to node flow balance
-        y = SGN(q)*z*p;                          // head loss / gradient
+        hyd->X_tmp[i] -= hyd->EmitterFlows[i];
-        sol->Aii[sol->Row[i]] += p;              // addition to main diagonal
-        sol->F[sol->Row[i]] += y + p * node->El; // addition to r.h.s.
-        hyd->X_tmp[i] -= q;                      // addition to net node inflow
     }
 }
 
 
-double emitflowchange(EN_Project *pr, int i)
+void emitheadloss(EN_Project *pr, int i, double *hloss, double *hgrad)
 /*
-**--------------------------------------------------------------
+**-------------------------------------------------------------
 **   Input:   i = node index
-**   Output:  returns change in flow at an emitter node
+**   Output:  hloss = head loss across node's emitter
-**   Purpose: computes flow change at an emitter node
+**            hgrad = head loss gradient
-**--------------------------------------------------------------
+**   Purpose: computes an emitters's head loss and gradient.
+**-------------------------------------------------------------
 */
 {
-    double ke, p;
+    double  ke;
+    double  q;
+    double  qa;
     hydraulics_t *hyd = &pr->hydraulics;
-    Snode *node = &pr->network.Node[i];
 
-    ke = MAX(CSMALL, node->Ke);
+    // Set adjusted emitter coeff.
-    p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0));
+    ke = MAX(CSMALL, pr->network.Node[i].Ke);
-    if (p < hyd->RQtol) p = 1 / hyd->RQtol;
+
-    else                p = 1.0 / p;
+    // Find flow below which head loss is linear
-    return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El));
+    qa = pow(hyd->RQtol / ke / hyd->Qexp, 1.0 / (hyd->Qexp - 1.0));
+
+    // Use linear head loss relation for small flow
+    q = hyd->EmitterFlows[i];
+    if (fabs(q) <= qa)
+    {
+        *hgrad = hyd->RQtol;
+        *hloss = (*hgrad) * q;
+    }
+
+    // Otherwise use normal emitter function
+    else
+    {
+        *hgrad = hyd->Qexp * ke * pow(fabs(q), hyd->Qexp - 1.0);
+        *hloss = (*hgrad) * q / hyd->Qexp;
+    }
 }
 
 
@@ -397,8 +421,8 @@ void demandparams(EN_Project *pr, double *dp, double *n)
 **   Input:   none
 **   Output:  dp = pressure range over which demands can vary
 **            n = exponent in head loss v. demand function
-**   Purpose: retrieves parameters that define a pressure dependent
+**   Purpose: retrieves parameters that define a pressure
-**            demand function.
+**            dependent demand function.
 **--------------------------------------------------------------
 */
 {
@@ -427,12 +451,13 @@ void  demandcoeffs(EN_Project *pr)
 **--------------------------------------------------------------
 **   Input:   none
 **   Output:  none
-**   Purpose: computes matrix coeffs. for pressure dependent demands
+**   Purpose: computes coeffs. of the linearized hydraulic eqns.
+**            contributed by pressure dependent demands.
 **
 **   Note: Pressure dependent demands are modelled like emitters
-**         with Hloss = Pserv * (D / Dfull)^(1/Pexp)
+**         with Hloss = Preq * (D / Dfull)^(1/Pexp)
 **         where D (actual demand) is zero for negative pressure
-**         and is Dfull above pressure Pserv.
+**         and is Dfull above pressure Preq.
 **--------------------------------------------------------------
 */
 {
@@ -524,7 +549,7 @@ void demandheadloss(double d, double dfull, double dp, double n,
     // Use linear head loss function for near zero demand
     else if (r < EPS)
     {
-        *hgrad = dp * pow(EPS, n - 1.0) / dfull;
+        *hgrad = dp * pow(EPS, n) / dfull / EPS;
         *hloss = (*hgrad) * d;
     }
 
@@ -542,25 +567,23 @@ void  pipecoeff(EN_Project *pr, int k)
 **--------------------------------------------------------------
 **   Input:   k = link index
 **   Output:  none
-**  Purpose:  computes P & Y coefficients for pipe k
+**   Purpose:  computes P & Y coefficients for pipe k.
 **
-**    P = inverse head loss gradient = 1/(dh/dQ)
+**    P = inverse head loss gradient = 1/hgrad
-**    Y = flow correction term = h*P
+**    Y = flow correction term = hloss / hgrad
 **--------------------------------------------------------------
 */
 {
-    double  hpipe,     // Normal head loss
+    double  hloss,     // Head loss
-            hml,       // Minor head loss
+            hgrad,     // Head loss gradient
             ml,        // Minor loss coeff.
-            p,         // q*(dh/dq)
             q,         // Abs. value of flow
             r;         // Resistance coeff.
 
     hydraulics_t *hyd = &pr->hydraulics;
     solver_t     *sol = &hyd->solver;
-    Slink *link = &pr->network.Link[k];
 
-    // For closed pipe use headloss formula: h = CBIG*q
+    // For closed pipe use headloss formula: hloss = CBIG*q
     if (hyd->LinkStatus[k] <= CLOSED)
     {
         sol->P[k] = 1.0 / CBIG;
@@ -575,31 +598,37 @@ void  pipecoeff(EN_Project *pr, int k)
         return;
     }
 
-    // Evaluate headloss coefficients
+    q = ABS(hyd->LinkFlows[k]);
-    q = ABS(hyd->LinkFlows[k]);                  // Absolute flow
+    ml = pr->network.Link[k].Km;
-    ml = link->Km;                               // Minor loss coeff.
+    r = pr->network.Link[k].R;
-    r = link->R;                                 // Resistance coeff.
 
-    // Use large P coefficient for small flow resistance product
+    // Friction head loss
-    if ( (r+ml)*q < hyd->RQtol)
+    // ... use linear relation for small flows
+    if (q <= pr->network.Link[k].Qa)
     {
-        sol->P[k] = 1.0 / hyd->RQtol;
+        hgrad = hyd->RQtol;
-        sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp;
+        hloss = hgrad * q;
-        return;
+    }
+    // ... use original formula for other flows
+    else
+    {
+        hgrad = hyd->Hexp * r * pow(q, hyd->Hexp - 1.0);
+        hloss = hgrad * q / hyd->Hexp;
     }
 
-    // Compute P and Y coefficients
+    // Contribution of minor head loss
-    hpipe = r*pow(q, hyd->Hexp);             // Friction head loss
-    p = hyd->Hexp*hpipe;                     // Q*dh(friction)/dQ 
     if (ml > 0.0)
     {
-        hml = ml*q*q;                        // Minor head loss
+        hloss += ml * q * q;
-        p += 2.0*hml;                        // Q*dh(Total)/dQ
+        hgrad += 2.0 * ml * q;
     }
-    else  hml = 0.0;
+
-    p = hyd->LinkFlows[k] / p;               // 1 / (dh/dQ)
+    // Adjust head loss sign for flow direction
-    sol->P[k] = ABS(p);
+    hloss *= SGN(hyd->LinkFlows[k]);
-    sol->Y[k] = p*(hpipe + hml);
+
+    // P and Y coeffs.
+    sol->P[k] = 1.0 / hgrad;
+    sol->Y[k] = hloss / hgrad;
 }
 
 
@@ -618,102 +647,85 @@ void DWpipecoeff(EN_Project *pr, int k)
     Slink *link = &pr->network.Link[k];
 
     double q = ABS(hyd->LinkFlows[k]);
-    double dfdq = 0.0;
+    double r = link->R;                         // Resistance coeff. 
-    double r, r1, f, ml, p, hloss;
+    double ml = link->Km;                       // Minor loss coeff. 
+    double e = link->Kc / link->Diam;           // Relative roughness
+    double s = hyd->Viscos * link->Diam;        // Viscosity / diameter
     
-    ml = link->Km;                       // Minor loss coeff. 
+    double hloss, hgrad, f, dfdq, r1;
-    r = link->R;                         // Resistance coeff. 
-    f =  frictionFactor(pr,k,&dfdq);     // D-W friction factor
-    r1 = f*r+ml;
     
-    // Use large P coefficient for small flow resistance product
+    // Compute head loss and its derivative
-    if (r1*q < hyd->RQtol)
+    // ... use Hagen-Poiseuille formula for laminar flow (Re <= 2000)
+    if (q <= A2 * s)
     {
-        sol->P[k] = 1.0/hyd->RQtol;
+        r = 16.0 * PI * s * r;
-        sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp;
+        hloss = hyd->LinkFlows[k] * (r + ml * q);
-        return;
+        hgrad  = r + 2.0 * ml * q;
+    }
+    
+    // ... otherwise use Darcy-Weisbach formula with friction factor
+    else
+    {
+        dfdq = 0.0;
+        f = frictionFactor(q, e, s, &dfdq);
+        r1 = f * r + ml;
+        hloss = r1 * q * hyd->LinkFlows[k];
+        hgrad = (2.0 * r1 * q) + (dfdq * r * q * q);
     }
     
     // Compute P and Y coefficients
-    hloss = r1*SQR(q);                   // Total head loss
+    sol->P[k] = 1.0 / hgrad;
-    p = 2.0*r1*q;                        // |dHloss/dQ|
+    sol->Y[k] = hloss / hgrad;
-    // + dfdq*r*q*q;                     // Ignore df/dQ term 
-    p = 1.0 / p;
-    sol->P[k] = p;
-    sol->Y[k] = SGN(hyd->LinkFlows[k]) * hloss * p;
 }
 
 
-double frictionFactor(EN_Project *pr, int k, double *dfdq)
+double frictionFactor(double q, double e, double s, double *dfdq)
 /*
 **--------------------------------------------------------------
-**   Input:   k = link index
+**   Input:   q = |pipe flow|
-**   Output:  returns friction factor and
+**            e = pipe roughness  / diameter
-**            replaces dfdq (derivative of f w.r.t. flow)
+**            s = viscosity * pipe diameter
+**   Output:  dfdq = derivative of friction factor w.r.t. flow
+**   Returns: pipe's friction factor
 **   Purpose: computes Darcy-Weisbach friction factor and its
 **            derivative as a function of Reynolds Number (Re).
-**
-**   Note:    Current formulas for dfdq need to be corrected
-**            so dfdq returned as 0.
 **--------------------------------------------------------------
 */
 {
-   double q,             // Abs. value of flow
+    double f;                // friction factor
-          f;             // Friction factor 
+    double x1, x2, x3, x4,
-   double x1,x2,x3,x4,
+           y1, y2, y3,
-          y1,y2,y3,
+           fa, fb, r;
-          fa,fb,r;
+    double w = q / s;        // Re*Pi/4
-   double s,w;
 
-  hydraulics_t *hyd = &pr->hydraulics;
+    //   For Re >= 4000 use Swamee & Jain approximation
-  Slink *link = &pr->network.Link[k];
+    //   of the Colebrook-White Formula
-
+    if ( w >= A1 )
-   *dfdq = 0.0;
-   if (link->Type > EN_PIPE)
-     return(1.0);                 // Only apply to pipes
-   q = ABS(hyd->LinkFlows[k]);
-   s = hyd->Viscos * link->Diam;
-   w = q/s;                       // w = Re(Pi/4)
-   
-   // For Re >= 4000 use Colebrook Formula
-   if (w >= A1)                   
     {
-      y1 = A8/pow(w,0.9);
+        y1 = A8 / pow(w, 0.9);
-      y2 = link->Kc/(3.7*link->Diam) + y1;
+        y2 = e / 3.7 + y1;
-      y3 = A9*log(y2);
+        y3 = A9 * log(y2);
-      f = 1.0/SQR(y3);
+        f = 1.0 / (y3*y3);
-      /*  *dfdq = (2.0+AA*y1/(y2*y3))*f; */   /* df/dq */
+        *dfdq = 1.8 * f * y1 * A9 / y2 / y3 / q;
     }
 
-   // For Re > 2000 use Interpolation Formula
+    //   Use interpolating polynomials developed by
-   else if (w > A2)              
+    //   E. Dunlop for transition flow from 2000 < Re < 4000.
-   {
-      y2 = link->Kc/(3.7*link->Diam) + AB;
-      y3 = A9*log(y2);
-      fa = 1.0/SQR(y3);
-      fb = (2.0+AC/(y2*y3))*fa;
-      r = w/A2;
-      x1 = 7.0*fa - fb;
-      x2 = 0.128 - 17.0*fa + 2.5*fb;
-      x3 = -0.128 + 13.0*fa - (fb+fb);
-      x4 = r*(0.032 - 3.0*fa + 0.5*fb);
-      f = x1 + r*(x2 + r*(x3+x4));
-      /*  *dfdq = (x1 + x1 + r*(3.0*x2 + r*(4.0*x3 + 5.0*x4)));  */
-   }
-   
-   // For Re > 8 (Laminar Flow) use Hagen-Poiseuille Formula
-   else if (w > A4)
-   {
-      f = A3*s/q;  // 16 * PI * Viscos * Diam / Flow
-      /*  *dfdq = A3*s; */
-   }
     else
     {
-      f = 8.0;
+        y2 = e / 3.7 + AB;
-      *dfdq = 0.0;
+        y3 = A9 * log(y2);
+        fa = 1.0 / (y3*y3);
+        fb = (2.0 + AC / (y2*y3)) * fa;
+        r = w / A2;
+        x1 = 7.0 * fa - fb;
+        x2 = 0.128 - 17.0 * fa + 2.5 * fb;
+        x3 = -0.128 + 13.0 * fa - (fb + fb);
+        x4 = 0.032 - 3.0 * fa + 0.5 *fb;
+        f = x1 + r * (x2 + r * (x3 + r * x4));
+        *dfdq = (x2 + r * (2.0 * x3 + r * 3.0 * x4)) / s / A2;
     }
-   return(f);
+    return f;
-
 }
 
 
@@ -730,13 +742,18 @@ void  pumpcoeff(EN_Project *pr, int k)
     double h0,               // Shutoff head
            q,                // Abs. value of flow 
            r,                // Flow resistance coeff.
-           n;                // Flow exponent coeff.
+           n,                // Flow exponent coeff.
+           setting,          // Pump speed setting
+           qa,               // Flow limit for linear head loss
+           hloss,            // Head loss across pump
+           hgrad;            // Head loss gradient
+
     hydraulics_t *hyd = &pr->hydraulics;
     solver_t     *sol = &hyd->solver;
-    double setting = hyd->LinkSetting[k];
     Spump        *pump;
 
     // Use high resistance pipe if pump closed or cannot deliver head
+    setting = hyd->LinkSetting[k];
     if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0)
     {
         sol->P[k] = 1.0 / CBIG;
@@ -744,35 +761,55 @@ void  pumpcoeff(EN_Project *pr, int k)
         return;
     }
 
+    // Obtain reference to pump object & its speed setting
     q = ABS(hyd->LinkFlows[k]);
-    q = MAX(q, TINY);
-
-    // Obtain reference to pump object
     p = findpump(&pr->network, k);
     pump = &pr->network.Pump[p];
 
-    // Get pump curve coefficients for custom pump curve.
+    // Get pump curve coefficients for custom pump curve
+    // (Other pump types have pre-determined coeffs.)
     if (pump->Ptype == CUSTOM)
     {
         // Find intercept (h0) & slope (r) of pump curve
         // line segment which contains speed-adjusted flow.
         curvecoeff(pr, pump->Hcurve, q / setting, &h0, &r);
 
-        // Determine head loss coefficients.
+        // Determine head loss coefficients (negative sign
+        // converts from pump curve's head gain to head loss)
         pump->H0 = -h0;
         pump->R = -r;
         pump->N = 1.0;
-    }
 
-    // Adjust head loss coefficients for pump speed.
+        // Compute head loss and its gradient
+        hgrad = pump->R * setting ;
+        hloss = pump->H0 * SQR(setting) + hgrad * hyd->LinkFlows[k];
+    }
+    else
+    {
+        // Adjust head loss coefficients for pump speed
         h0 = SQR(setting) * pump->H0;
         n = pump->N;
         r = pump->R * pow(setting, 2.0 - n);
-    if (n != 1.0) r = n * r * pow(q, n - 1.0);
 
-    // Compute inverse headloss gradient (P) and flow correction factor (Y)
+        // Compute head loss and its gradient
-    sol->P[k] = 1.0 / MAX(r, hyd->RQtol);
+        // ... use linear approx. to pump curve for small flows
-    sol->Y[k] = hyd->LinkFlows[k] / n + sol->P[k] * h0;
+        qa = pow(hyd->RQtol / n / r, 1.0 / (n - 1.0));
+        if (q <= qa)
+        {
+            hgrad = hyd->RQtol;
+            hloss = h0 + hgrad * hyd->LinkFlows[k];
+        }
+        // ... use original pump curve for normal flows
+        else
+        {
+            hgrad = n * r * pow(q, n - 1.0);
+            hloss = h0 + hgrad * hyd->LinkFlows[k] / n;
+        }
+    }
+
+    // P and Y coeffs.
+    sol->P[k] = 1.0 / hgrad;
+    sol->Y[k] = hloss / hgrad;
 }
 
 
@@ -825,9 +862,10 @@ void  gpvcoeff(EN_Project *pr, int k)
 **--------------------------------------------------------------
 */
 {
-    double h0,        // Headloss curve intercept
+    int    i;
-           q,         // Abs. value of flow
+    double h0,        // Intercept of head loss curve segment
-           r;         // Flow resistance coeff.
+           r,         // Slope of head loss curve segment
+           q;         // Abs. value of flow
 
     hydraulics_t *hyd = &pr->hydraulics;
     solver_t     *sol = &hyd->solver;
@@ -835,19 +873,25 @@ void  gpvcoeff(EN_Project *pr, int k)
     // Treat as a pipe if valve closed
     if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k);                          
 
-    // Otherwise utilize headloss curve
+    // Otherwise utilize segment of head loss curve
-    // whose index is stored in K
+    // bracketing current flow (curve index is stored
+    // in valve's setting)
     else
     {
-        // Find slope & intercept of headloss curve.
+        // Index of valve's head loss curve
+        i = (int)ROUND(hyd->LinkSetting[k]);
+
+        // Adjusted flow rate
         q = ABS(hyd->LinkFlows[k]);
         q = MAX(q, TINY);
-        curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r);
 
-        // Compute inverse headloss gradient (P)
+        // Intercept and slope of curve segment containing q
-        // and flow correction factor (Y).
+        curvecoeff(pr, i, q, &h0, &r);
-        sol->P[k] = 1.0 / MAX(r, hyd->RQtol);
+        r = MAX(r, TINY);
-        sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]);
+
+        // Resulting P and Y coeffs.
+        sol->P[k] = 1.0 / r;
+        sol->Y[k] = (h0 / r + q) * SGN(hyd->LinkFlows[k]);
     }
 }
 
@@ -1084,14 +1128,14 @@ void valvecoeff(EN_Project *pr, int k)
 **--------------------------------------------------------------
 */
 {
-    double p;
+    double flow, q, y, qa, hgrad;
     
     EN_Network   *net = &pr->network;
     hydraulics_t *hyd = &pr->hydraulics;
     solver_t     *sol = &hyd->solver;
     Slink *link = &net->Link[k];
 
-    double flow = hyd->LinkFlows[k];
+    flow = hyd->LinkFlows[k];
 
     // Valve is closed. Use a very small matrix coeff.
     if (hyd->LinkStatus[k] <= CLOSED)
@@ -1104,14 +1148,30 @@ void valvecoeff(EN_Project *pr, int k)
     // Account for any minor headloss through the valve
     if (link->Km > 0.0)
     {
-        p = 2.0 * link->Km * fabs(flow);
+        // Adjust for very small flow
-        if (p < hyd->RQtol) p = hyd->RQtol;
+        q = fabs(flow);
-        sol->P[k] = 1.0 / p;
+        qa = hyd->RQtol / 2.0 / link->Km;
-        sol->Y[k] = flow / 2.0;
+        if (q <= qa)
+        {
+            hgrad = hyd->RQtol;
+            y = flow;
         }
         else
         {
-        sol->P[k] = 1.0 / hyd->RQtol;
+            hgrad = 2.0 * link->Km * q;
+            y = flow / 2.0;
+        }
+
+        // P and Y coeffs.
+        sol->P[k] = 1.0 / hgrad;
+        sol->Y[k] = y;
+    }
+
+    // If no minor loss coeff. specified use a
+    // low resistance linear head loss relation
+    else
+    {
+        sol->P[k] = 1.0 / CSMALL;
         sol->Y[k] = flow;
     }
 }

src/hydsolver.c

@@ -494,30 +494,34 @@ void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum,
 **----------------------------------------------------------------
 */
 {
-    double  dq;
+    int     i;
-    int     k;
+    double  hloss, hgrad, dh, dq;
     EN_Network   *net = &pr->network;
     hydraulics_t *hyd = &pr->hydraulics;
 
     // Examine each network junction
-    for (k = 1; k <= net->Njuncs; k++)
+    for (i = 1; i <= net->Njuncs; i++)
     {
         // Skip junction if it does not have an emitter
-        if (net->Node[k].Ke == 0.0) continue;
+        if (net->Node[i].Ke == 0.0) continue;
 
-        // Find emitter flow change (see hydcoeffs.c)
+        // Find emitter head loss and gradient 
-        dq = emitflowchange(pr, k);
+        emitheadloss(pr, i, &hloss, &hgrad);
-        hyd->EmitterFlows[k] -= dq;
+
+        // Find emitter flow change
+        dh = hyd->NodeHead[i] - net->Node[i].El;
+        dq = (hloss - dh) / hgrad;
+        hyd->EmitterFlows[i] -= dq;
 
         // Update system flow summation
-        *qsum += ABS(hyd->EmitterFlows[k]);
+        *qsum += ABS(hyd->EmitterFlows[i]);
         *dqsum += ABS(dq);
 
         // Update identity of element with max. flow change
         if (ABS(dq) > hbal->maxflowchange)
         {
             hbal->maxflowchange = ABS(dq);
-            hbal->maxflownode = k;
+            hbal->maxflownode = i;
             hbal->maxflowlink = -1;
         }
     }

src/types.h

@@ -410,7 +410,8 @@ typedef struct            /* LINK OBJECT */
    double  Kb;             /* Bulk react. coeff */
    double  Kw;             /* Wall react. coeff */
    double  R;              /* Flow resistance   */
-   double  Rc;             /* Reaction cal      */
+   double  Rc;             /* Reaction coeff.   */
+   double  Qa;             // Low flow limit
    EN_LinkType Type;       /* Link type         */
    StatType Stat;          /* Initial status    */
    char Rpt;            /* Reporting flag    */

tools/before-test.cmd

@@ -25,8 +25,8 @@ set SCRIPT_HOME=%~dp0
 set TEST_HOME=%~1
 
 
-set EXAMPLES_VER=1.0.2-dev.2
+set EXAMPLES_VER=1.0.2-dev.3
-set BENCHMARK_VER=220dev2
+set BENCHMARK_VER=220dev3
 
 
 set TESTFILES_URL=https://github.com/OpenWaterAnalytics/epanet-example-networks/archive/v%EXAMPLES_VER%.zip

tools/nrtest-epanet/report_diff.py

@@ -51,16 +51,22 @@ def _log_relative_error(q, c):
     '''
     diff = np.subtract(q, c)
     tmp_c = np.copy(c)
+
     # If ref value is small compute absolute error
     tmp_c[np.fabs(tmp_c) < 1.0e-6] = 1.0    
-    
     re = np.fabs(diff)/np.fabs(tmp_c)
+
     # If re is tiny set lre to number of digits
     re[re < 1.0e-7] = 1.0e-7
     # If re is very large set lre to zero
     re[re > 2.0] = 1.0
 
-    return np.negative(np.log10(re))
+    lre = np.negative(np.log10(re))
+
+    # If lre is negative set to zero
+    lre[lre < 1.0] = 0.0
+
+    return lre
     
     
 def _print_diff(idx, lre, test, ref):
@@ -71,7 +77,7 @@ def _print_diff(idx, lre, test, ref):
     diff_val = (test_val - ref_val)
     lre_val = (lre[idx[0]])
     
-    print("Idx: %s\nSut: %f  Ref: %f  Diff: %f  LRE: %f\n" 
+    print("Idx: %s\nSut: %e  Ref: %e  Diff: %e  LRE: %.2f\n" 
           % (idx_val, test_val, ref_val, diff_val, lre_val))
 
 

tools/run-nrtest.cmd

@@ -18,7 +18,7 @@ setlocal
 set NRTEST_SCRIPT_PATH=%~1
 set TEST_SUITE_PATH=%~2
 
-set BENCHMARK_VER=220dev2
+set BENCHMARK_VER=220dev3
 
 
 set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute