Merge pull request #195 from LRossman/contributor-lr

Feature Update - more rigorous hydraulic convergence criteria (#161)

include/epanet2.bas

@@ -128,6 +128,8 @@ Public Const EN_ACCURACY = 1
 Public Const EN_TOLERANCE = 2
 Public Const EN_EMITEXPON = 3
 Public Const EN_DEMANDMULT = 4
+Public Const EN_HEADERROR = 5
+Public Const EN_FLOWCHANGE = 6
 
 Public Const EN_LOWLEVEL = 0     ' Control types
 Public Const EN_HILEVEL = 1

include/epanet2.h

@@ -212,7 +212,9 @@ typedef enum {
   EN_ACCURACY     = 1,
   EN_TOLERANCE    = 2,
   EN_EMITEXPON    = 3,
-  EN_DEMANDMULT   = 4
+  EN_DEMANDMULT   = 4,
+  EN_HEADERROR    = 5,
+  EN_FLOWCHANGE   = 6
 } EN_Option;
 
 typedef enum {

include/epanet2.vb

@@ -121,6 +121,8 @@ Public Const EN_ACCURACY = 1
 Public Const EN_TOLERANCE = 2
 Public Const EN_EMITEXPON = 3
 Public Const EN_DEMANDMULT = 4
+Public Const EN_HEADERROR = 5
+Public Const EN_FLOWCHANGE = 6
 
 Public Const EN_LOWLEVEL = 0     ' Control types
 Public Const EN_HILEVEL = 1

src/epanet.c

@@ -71,11 +71,12 @@ This module calls the following functions that reside in other modules:
      runhyd()
      nexthyd()
      closehyd()
-     resistance()
      tankvolume()
      getenergy()
      setlinkstatus()
      setlinksetting()
+   HYDCOEFFS
+     resistcoeff()
    QUALITY.C
      openqual()
      initqual()
@@ -129,6 +130,7 @@ execute function x and set the error code equal to its return value.
 #include "text.h"
 #include "types.h"
 #define EXTERN
+////////////////////////////////////////////#include "epanet2.h"
 #include "vars.h"
 
 /****************************************************************
@@ -1341,6 +1343,14 @@ int DLLEXPORT EN_getoption(EN_Project *pr, EN_Option code,
   case EN_DEMANDMULT:
     v = hyd->Dmult;
     break;
+
+  case EN_HEADERROR:
+    v = hyd->HeadErrorLimit * Ucf[HEAD];
+    break;
+  case EN_FLOWCHANGE:
+    v = hyd->FlowChangeLimit * Ucf[FLOW];
+    break;
+
   default:
     return (251);
   }
@@ -2648,7 +2658,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code,
       r = Link[index].Diam / value;      /* Ratio of old to new diam */
       Link[index].Km *= SQR(r) * SQR(r); /* Adjust minor loss factor */
       Link[index].Diam = value;          /* Update diameter */
-      resistance(p, index);                 /* Update resistance factor */
+      resistcoeff(p, index);             /* Update resistance coeff. */
     }
     break;
 
@@ -2657,7 +2667,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code,
       if (value <= 0.0)
         return (202);
       Link[index].Len = value / Ucf[ELEV];
-      resistance(p, index);
+      resistcoeff(p, index);
     }
     break;
 
@@ -2668,7 +2678,7 @@ int DLLEXPORT EN_setlinkvalue(EN_Project *p, int index, int code,
       Link[index].Kc = value;
       if (hyd->Formflag == DW)
         Link[index].Kc /= (1000.0 * Ucf[ELEV]);
-      resistance(p, index);
+      resistcoeff(p, index);
     }
     break;
 
@@ -3171,6 +3181,18 @@ int DLLEXPORT EN_setoption(EN_Project *p, int code, EN_API_FLOAT_TYPE v)
       return (202);
     hyd->Dmult = value;
     break;
+
+  case EN_HEADERROR:
+    if (value < 0.0)
+        return (202);
+    hyd->HeadErrorLimit = value / Ucf[HEAD];
+    break;
+  case EN_FLOWCHANGE:
+      if (value < 0.0)
+          return (202);
+      hyd->FlowChangeLimit = value / Ucf[FLOW];
+      break;
+
   default:
     return (251);
   }
@@ -3419,7 +3441,7 @@ int openhydfile(EN_Project *p)
   out->HydFile = NULL;
   switch (out->Hydflag) {
   case SCRATCH:
-    getTmpName(out->HydFname); 
+    getTmpName(p, out->HydFname); 
     out->HydFile = fopen(out->HydFname, "w+b"); 
     break;
   case SAVE:
@@ -3512,7 +3534,7 @@ int openoutfile(EN_Project *p)
   // else if ( (OutFile = tmpfile()) == NULL) 
   else 
   {
-    getTmpName(out->OutFname);                           
+    getTmpName(p, out->OutFname);                           
     if ((out->OutFile = fopen(out->OutFname, "w+b")) == NULL) 
     {
       writecon(FMT08);
@@ -3530,7 +3552,7 @@ int openoutfile(EN_Project *p)
   if (!errcode) {
     if (rep->Tstatflag != SERIES) {
       // if ( (TmpOutFile = tmpfile()) == NULL) errcode = 304; 
-      getTmpName(out->TmpFname);                
+      getTmpName(p, out->TmpFname);                
       out->TmpOutFile = fopen(out->TmpFname, "w+b"); 
       if (out->TmpOutFile == NULL)
         errcode = 304; 
@@ -3837,7 +3859,7 @@ void freedata(EN_Project *p)
 ----------------------------------------------------------------
 */
 
-char *getTmpName(char *fname)
+char *getTmpName(EN_Project *p, char *fname)
 //
 //  Input:   fname = file name string
 //  Output:  returns pointer to file name
@@ -3865,7 +3887,6 @@ char *getTmpName(char *fname)
 ///////////////////  DEPRECATED   /////////////////////////////////////
 // --- for Windows systems:
 #ifdef WINDOWS
-  out_file_t *out = &p->out_files;
   // --- use system function tmpnam() to create a temporary file name
   char name[MAXFNAME + 1];
   int n;

src/funcs.h

@@ -31,17 +31,18 @@ AUTHOR:     L. Rossman
 
 #include "types.h"
 
-void    initpointers(EN_Project *p);                /* Initializes pointers       */
+void    initpointers(EN_Project *pr);              /* Initializes pointers       */
-int     allocdata(EN_Project *p);                   /* Allocates memory           */
+int     allocdata(EN_Project *pr);                 /* Allocates memory           */
 void    freeTmplist(STmplist *);                   /* Frees items in linked list */
 void    freeFloatlist(SFloatlist *);               /* Frees list of floats       */
-void    freedata(EN_Project *p);                    /* Frees allocated memory     */
+void    freedata(EN_Project *pr);                  /* Frees allocated memory     */
-int     openfiles(EN_Project *p, char *,char *,char *);  /* Opens input & report files */
+int     openfiles(EN_Project *pr, char *,char *,char *);  /* Opens input & report files */
-int     openhydfile(EN_Project *p);                 /* Opens hydraulics file      */
+int     openhydfile(EN_Project *pr);               /* Opens hydraulics file      */
-int     openoutfile(EN_Project *p);                 /* Opens binary output file   */
+int     openoutfile(EN_Project *pr);               /* Opens binary output file   */
 int     strcomp(char *, char *);                   /* Compares two strings       */
-char*   getTmpName(char* fname);                    /* Gets temporary file name   */     
+char*   getTmpName(EN_Project *p, char* fname);    /* Gets temporary file name   */     
-double  interp(int, double *,double *, double);     /* Interpolates a data curve  */
+double  interp(int n, double x[], double y[],
+        double xx);                                /* Interpolates a data curve  */
                
 int     findnode(EN_Network *n, char *);           /* Finds node's index from ID */
 int     findlink(EN_Network *n, char *);           /* Finds link's index from ID */
@@ -77,7 +78,8 @@ int     getpatterns(EN_Project *pr);                /* Gets pattern data from li
 int     getcurves(EN_Project *pr);                  /* Gets curve data from list  */
 int     findmatch(char *, char *[]);                /* Finds keyword in line      */
 int     match(const char *, const char *);          /* Checks for word match      */
-int     gettokens(char *s, char** Tok, int maxToks, char *comment); /* Tokenizes input line       */
+int     gettokens(char *s, char** Tok, int maxToks,
+                  char *comment);                   /* Tokenizes input line       */
 int     getfloat(char *, double *);                 /* Converts string to double  */
 double  hour(char *, char *);                       /* Converts time to hours     */
 int     setreport(EN_Project *pr, char *);          /* Processes reporting command*/
@@ -111,7 +113,8 @@ int     powercurve(double, double, double,          /* Coeffs. of power pump cur
                    double, double, double *,
                    double *, double *);
 int     valvecheck(EN_Project *pr, int, int, int);  /* Checks valve placement     */
-void    changestatus(EN_Network *net, int, StatType, double);  /* Changes status of a link   */
+void    changestatus(EN_Network *net, int, StatType,
+                     double);                       /* Changes status of a link   */
 
 /* -------------- RULES.C --------------*/
 void    initrules(rules_t *rules);                  /* Initializes rule base      */
@@ -142,7 +145,8 @@ int     disconnected(EN_Project *pr);                     /* Checks for disconne
 void    marknodes(EN_Project *pr, int, int *, char *);    /* Identifies connected nodes */
 void    getclosedlink(EN_Project *pr, int, char *);       /* Finds a disconnecting link */
 void    writelimits(EN_Project *pr, int,int);             /* Writes reporting limits    */
-int     checklimits(report_options_t *rep, double *,int,int);     /* Checks variable limits     */
+int     checklimits(report_options_t *rep, double *,
+                    int,int);                             /* Checks variable limits     */
 void    writetime(EN_Project *pr, char *);                /* Writes current clock time  */
 char    *clocktime(char *, long);                         /* Converts time to hrs:min   */
 char    *fillstr(char *, char, int);                      /* Fills string with character*/
@@ -150,10 +154,7 @@ int     getnodetype(EN_Network *net, int);                /* Determines node typ
 
 /* --------- HYDRAUL.C -----------------*/
 int     openhyd(EN_Project *pr);                    /* Opens hydraulics solver    */
-
-/*** Updated 3/1/01 ***/
 void    inithyd(EN_Project *pr, int initFlags);     /* Re-sets initial conditions */
-
 int     runhyd(EN_Project *pr, long *);             /* Solves 1-period hydraulics */
 int     nexthyd(EN_Project *pr, long *);            /* Moves to next time period  */
 void    closehyd(EN_Project *pr);                   /* Closes hydraulics solver   */
@@ -164,64 +165,31 @@ void    initlinkflow(EN_Project *pr, int, char,
 void    setlinkflow(EN_Project *pr, int, double);   /* Sets link flow via headloss*/
 void    setlinkstatus(EN_Project *pr, int, char,
                       StatType *, double *);        /* Sets link status           */
-                      
 void    setlinksetting(EN_Project *pr, int, double,
                        StatType *, double *);       /* Sets pump/valve setting    */
                        
-void    resistance(EN_Project *pr, int);            /* Computes resistance coeff. */
 void    demands(EN_Project *pr);                    /* Computes current demands   */
 int     controls(EN_Project *pr);                   /* Controls link settings     */
 long    timestep(EN_Project *pr);                   /* Computes new time step     */
 int     tanktimestep(EN_Project *pr, long *);       /* Time till tanks fill/drain */
 void    controltimestep(EN_Project *pr, long *);    /* Time till control action   */
 void    ruletimestep(EN_Project *pr, long *);       /* Time till rule action      */
+
 void    addenergy(EN_Project *pr, long);            /* Accumulates energy usage   */
 void    getenergy(EN_Project *pr, int, double *,
                   double *);                        /* Computes link energy use   */
+
 void    tanklevels(EN_Project *pr, long);           /* Computes new tank levels   */
 double  tankvolume(EN_Project *pr, int,double);     /* Finds tank vol. from grade */
 double  tankgrade(EN_Project *pr, int,double);      /* Finds tank grade from vol. */
-int     netsolve(EN_Project *pr, int *,double *);   /* Solves network equations   */
-int     badvalve(EN_Project *pr, int);              /* Checks for bad valve       */
-int     valvestatus(EN_Project *pr);                /* Updates valve status       */
-int     linkstatus(EN_Project *pr);                 /* Updates link status        */
-StatType    cvstatus(EN_Project *pr, StatType,
-                     double,double);                /* Updates CV status          */
-StatType    pumpstatus(EN_Project *pr, int,double); /* Updates pump status        */
-StatType    prvstatus(EN_Project *pr, int,StatType,
-                      double,double,double);        /* Updates PRV status         */
 
-StatType    psvstatus(EN_Project *pr, int,StatType,
+/* ----------- HYDSOLVER.C -  ----------*/
-                      double,double,double);        /* Updates PSV status         */
+int     hydsolve(EN_Project *pr, int *,double *);   /* Solves network equations   */
 
-StatType    fcvstatus(EN_Project *pr, int,StatType,
+/* ----------- HYDCOEFFS.C --------------*/
-                      double,double);               /* Updates FCV status         */
+void    resistcoeff(EN_Project *pr, int k);         /* Finds pipe flow resistance */
-                  
+void    hlosscoeff(EN_Project *pr, int k);          /* Finds link head loss coeff */
-void    tankstatus(EN_Project *pr, int,int,int);    /* Checks if tank full/empty  */
+void    matrixcoeffs(EN_Project *pr);               /* Finds hyd. matrix coeffs.  */
-int     pswitch(EN_Project *pr);                    /* Pressure switch controls   */
-double  newflows(EN_Project *pr);                   /* Updates link flows         */
-void    newcoeffs(EN_Project *pr);                  /* Computes matrix coeffs.    */
-void    linkcoeffs(EN_Project *pr);                 /* Computes link coeffs.      */
-void    nodecoeffs(EN_Project *pr);                 /* Computes node coeffs.      */
-void    valvecoeffs(EN_Project *pr);                /* Computes valve coeffs.     */
-void    pipecoeff(EN_Project *pr, int);             /* Computes pipe coeff.       */
-double  DWcoeff(EN_Project *pr, int, double *);     /* Computes D-W coeff.        */
-void    pumpcoeff(EN_Project *pr, int);             /* Computes pump coeff.       */
-
-/*** Updated 10/25/00 ***/
-/*** Updated 12/29/00 ***/
-void    curvecoeff(EN_Project *pr, int,double,
-                   double *,double *);              /* Computes curve coeffs.     */
-                             
-
-void    gpvcoeff(EN_Project *pr, int iLink);        /* Computes GPV coeff.        */
-void    pbvcoeff(EN_Project *pr, int iLink);        /* Computes PBV coeff.        */
-void    tcvcoeff(EN_Project *pr, int iLink);        /* Computes TCV coeff.        */
-void    prvcoeff(EN_Project *pr, int iLink, int n1, int n2);  /* Computes PRV coeff.        */
-void    psvcoeff(EN_Project *pr, int iLink, int n1, int n2);  /* Computes PSV coeff.        */
-void    fcvcoeff(EN_Project *pr, int iLink, int n1, int n2);  /* Computes FCV coeff.        */
-void    emittercoeffs(EN_Project *pr);              /* Computes emitter coeffs.   */
-double  emitflowchange(EN_Project *pr, int);        /* Computes new emitter flow  */
 
 /* ----------- SMATRIX.C ---------------*/
 int     createsparse(EN_Project *pr);               /* Creates sparse matrix      */
@@ -242,7 +210,7 @@ int     storesparse(EN_Project *pr, int);           /* Stores sparse matrix
 int     ordersparse(hydraulics_t *h, int);          /* Orders matrix storage      */
 void    transpose(int,int *,int *,                  /* Transposes sparse matrix   */
         int *,int *,int *,int *,int *);
-int     linsolve(solver_t *s, int);                 /* Solves set of linear eqns. */
+int     linsolve(solver_t *s, int);                 /* Solves set of inear eqns.  */
 
 /* ----------- QUALITY.C ---------------*/
 int     openqual(EN_Project *pr);                   /* Opens WQ solver system     */
@@ -282,7 +250,6 @@ double  bulkrate(EN_Project *pr, double,double,
 double  wallrate(EN_Project *pr, double,double,
                  double,double);                    /* Finds wall reaction rate   */
 
-
 /* ------------ OUTPUT.C ---------------*/
 int     savenetdata(EN_Project *pr);                /* Saves basic data to file   */
 int     savehyd(EN_Project *pr, long *);            /* Saves hydraulic solution   */

src/hydcoeffs.c

@@ -0,0 +1,933 @@
+/*
+*********************************************************************
+
+HYDCOEFFS.C --  hydraulic coefficients for the EPANET Program
+
+*******************************************************************
+*/
+
+#include <stdio.h>
+#include <string.h>
+#ifndef __APPLE__
+#include <malloc.h>
+#else
+#include <stdlib.h>
+#endif
+#include <math.h>
+#include "types.h"
+#include "funcs.h"
+
+// Constants used for computing Darcy-Weisbach friction factor
+const double A1 = 0.314159265359e04;   // 1000*PI
+const double A2 = 0.157079632679e04;   // 500*PI
+const double A3 = 0.502654824574e02;   // 16*PI
+const double A4 = 6.283185307;         // 2*PI
+const double A8 = 4.61841319859;       // 5.74*(PI/4)^.9
+const double A9 = -8.685889638e-01;    // -2/ln(10)
+const double AA = -1.5634601348;       // -2*.9*2/ln(10)
+const double AB = 3.28895476345e-03;   // 5.74/(4000^.9)
+const double AC = -5.14214965799e-03;  // AA*AB
+
+// External functions
+//void  resistcoeff(EN_Project *pr, int k);
+//void  hlosscoeff(EN_Project *pr, int k);
+//void  matrixcoeffs(EN_Project *pr);
+
+// Local functions
+static void    linkcoeffs(EN_Project *pr);
+static void    nodecoeffs(EN_Project *pr);
+static void    valvecoeffs(EN_Project *pr);
+static void    emittercoeffs(EN_Project *pr);
+
+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 void    pumpcoeff(EN_Project *pr, int k);
+static void    curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r);
+
+static void    valvecoeff(EN_Project *pr, int k);
+static void    gpvcoeff(EN_Project *pr, int k);
+static void    pbvcoeff(EN_Project *pr, int k);
+static void    tcvcoeff(EN_Project *pr, int k);
+static void    prvcoeff(EN_Project *pr, int k, int n1, int n2);
+static void    psvcoeff(EN_Project *pr, int k, int n1, int n2);
+static void    fcvcoeff(EN_Project *pr, int k, int n1, int n2);
+
+
+
+void  resistcoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------------
+**  Input:   k = link index
+**  Output:  none
+**  Purpose: computes link flow resistance coefficient
+**--------------------------------------------------------------------
+*/
+{
+    double e, d, L;
+    EN_Network *net = &pr->network;
+    hydraulics_t *hyd = &pr->hydraulics;
+    Slink *link = &net->Link[k];
+    
+    switch (link->Type) {
+
+    // ... Link is a pipe. Compute resistance based on headloss formula.
+    //     Friction factor for D-W formula gets included during head loss 
+    //     calculation.
+    case EN_CVPIPE:
+    case EN_PIPE:
+        e = link->Kc;                 // Roughness coeff.
+        d = link->Diam;               // Diameter
+        L = link->Len;                // Length
+        switch (hyd->Formflag)
+        {
+        case HW:
+            link->R = 4.727*L / pow(e, hyd->Hexp) / pow(d, 4.871);
+            break;
+        case DW:
+            link->R = L / 2.0 / 32.2 / d / SQR(PI*SQR(d) / 4.0);
+            break;
+        case CM:
+            link->R = SQR(4.0*e / (1.49*PI*d*d)) * pow((d / 4.0), -1.333)*L;
+        }
+        break;
+
+    // ... Link is a pump. Use huge resistance.
+    case EN_PUMP:
+        link->R = CBIG;
+        break;
+
+    // ... For all other links (e.g. valves) use a small resistance
+    default:
+        link->R = CSMALL;
+        break;
+    }
+}
+
+
+void  hlosscoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  none
+**   Purpose: computes P and Y coefficients for a link
+**--------------------------------------------------------------
+*/
+{
+    EN_Network   *net = &pr->network;
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    Slink        *link = &net->Link[k];
+
+    switch (link->Type)
+    {
+    case EN_CVPIPE:
+    case EN_PIPE:
+        pipecoeff(pr, k);
+        break;
+    case EN_PUMP:
+        pumpcoeff(pr, k);
+        break;
+    case EN_PBV:
+        pbvcoeff(pr, k);
+        break;
+    case EN_TCV:
+        tcvcoeff(pr, k);
+        break;
+    case EN_GPV:
+        gpvcoeff(pr, k);
+        break;
+    case EN_FCV:
+    case EN_PRV:
+    case EN_PSV:
+        if (hyd->LinkSetting[k] == MISSING) {
+            valvecoeff(pr, k);
+        }
+        else sol->P[k] = 0.0;
+    }
+}
+
+
+void   matrixcoeffs(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+**  Input:   none
+**  Output:  none
+**  Purpose: computes coefficients of linearized network eqns.
+**--------------------------------------------------------------
+*/
+{
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    EN_Network   *net = &pr->network;
+
+    memset(sol->Aii, 0, (net->Nnodes + 1) * sizeof(double));
+    memset(sol->Aij, 0, (hyd->Ncoeffs + 1) * sizeof(double));
+    memset(sol->F, 0, (net->Nnodes + 1) * sizeof(double));
+    memset(hyd->X_tmp, 0, (net->Nnodes + 1) * sizeof(double));
+    linkcoeffs(pr);
+    emittercoeffs(pr);
+    nodecoeffs(pr);
+    valvecoeffs(pr);
+}
+
+
+void  linkcoeffs(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+**   Input:   none
+**   Output:  none
+**   Purpose: computes matrix coefficients for links
+**--------------------------------------------------------------
+*/
+{
+    int   k, n1, n2;
+
+    EN_Network   *net = &pr->network;
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    Slink *link;
+
+    // Examine each link of network */
+    for (k = 1; k <= net->Nlinks; k++)
+    {
+        if (sol->P[k] == 0.0) continue;
+        link = &net->Link[k];
+        n1 = link->N1;           // Start node of link
+        n2 = link->N2;           // End node of link
+
+        // Update net nodal inflows (X), solution matrix (A) and RHS array (F)
+        // (Use covention that 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.
+        sol->Aij[sol->Ndx[k]] -= sol->P[k];
+
+        // 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
+        else {
+            sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); 
+        }
+
+        // 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
+        else {
+            sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); 
+        }
+    }
+}
+
+
+void  nodecoeffs(EN_Project *pr)
+/*
+**----------------------------------------------------------------
+**  Input:   none
+**  Output:  none
+**  Purpose: completes calculation of nodal flow imbalance (X)
+**           & flow correction (F) arrays
+**----------------------------------------------------------------
+*/
+{
+    int   i;
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    EN_Network   *net = &pr->network;
+
+    // For junction nodes, subtract demand flow from net
+    // flow imbalance & add imbalance to RHS array F.
+    for (i = 1; i <= net->Njuncs; i++)
+    {
+        hyd->X_tmp[i] -= hyd->NodeDemand[i];
+        sol->F[sol->Row[i]] += hyd->X_tmp[i];
+    }
+}
+
+
+void  valvecoeffs(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+**   Input:   none
+**   Output:  none
+**   Purpose: computes matrix coeffs. for PRVs, PSVs & FCVs
+**            whose status is not fixed to OPEN/CLOSED
+**--------------------------------------------------------------
+*/
+{
+    int i, k, n1, n2;
+
+    hydraulics_t *hyd = &pr->hydraulics;
+    EN_Network   *net = &pr->network;
+    Slink *link;
+    Svalve *valve;
+
+    // Examine each valve
+    for (i = 1; i <= net->Nvalves; i++)
+    {
+        valve = &net->Valve[i];
+        k = valve->Link;                         // Link index of valve
+        link = &net->Link[k];
+        if (hyd->LinkSetting[k] == MISSING) {
+            continue;                            // Valve status fixed
+        }
+        n1 = link->N1;                           // Start & end nodes
+        n2 = link->N2;
+        switch (link->Type)                      // Call valve-specific function 
+        {
+        case EN_PRV:
+            prvcoeff(pr, k, n1, n2);
+            break;
+        case EN_PSV:
+            psvcoeff(pr, k, n1, n2);
+            break;
+        case EN_FCV:
+            fcvcoeff(pr, k, n1, n2);
+            break;
+        default:   continue;
+        }
+    }
+}
+
+
+void  emittercoeffs(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+**   Input:   none
+**   Output:  none
+**   Purpose: computes matrix coeffs. for emitters
+**
+**   Note: Emitters consist of a fictitious pipe connected to
+**         a fictitious reservoir whose elevation equals that
+**         of the junction. The headloss through this pipe is
+**         Ke*(Flow)^hyd->Qexp, where Ke = emitter headloss coeff.
+**--------------------------------------------------------------
+*/
+{
+    int   i;
+    double  ke;
+    double  p;
+    double  q;
+    double  y;
+    double  z;
+
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    EN_Network   *net = &pr->network;
+    Snode *node;
+
+    for (i = 1; i <= net->Njuncs; i++)
+    {
+        node = &net->Node[i];
+        if (node->Ke == 0.0) {
+            continue;
+        }
+        ke = MAX(CSMALL, node->Ke);              // emitter coeff.
+        q = hyd->EmitterFlows[i];                // emitter flow
+        z = ke * pow(ABS(q), hyd->Qexp);         // emitter head loss
+        p = hyd->Qexp * z / ABS(q);              // head loss gradient
+        if (p < hyd->RQtol) {
+            p = 1.0 / hyd->RQtol;
+        }
+        else {
+            p = 1.0 / p;                         // inverse head loss gradient
+        }
+        y = SGN(q)*z*p;                          // head loss / gradient
+        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
+    }
+}
+
+
+void  pipecoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  none
+**  Purpose:  computes P & Y coefficients for pipe k
+**
+**    P = inverse head loss gradient = 1/(dh/dQ)
+**    Y = flow correction term = h*P
+**--------------------------------------------------------------
+*/
+{
+    double  hpipe,     // Normal head loss
+            hml,       // Minor head loss
+            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
+    if (hyd->LinkStatus[k] <= CLOSED)
+    {
+        sol->P[k] = 1.0 / CBIG;
+        sol->Y[k] = hyd->LinkFlows[k];
+        return;
+    }
+
+    // ... head loss formula is Darcy-Weisbach
+    if (hyd->Formflag == DW) {
+        DWpipecoeff(pr, k);
+        return;
+    }
+
+    // Evaluate headloss coefficients
+    q = ABS(hyd->LinkFlows[k]);                  // Absolute flow
+    ml = link->Km;                               // Minor loss coeff.
+    r = link->R;                                 // Resistance coeff.
+
+    // Use large P coefficient for small flow resistance product
+    if ( (r+ml)*q < hyd->RQtol)
+    {
+        sol->P[k] = 1.0 / hyd->RQtol;
+        sol->Y[k] = hyd->LinkFlows[k] / hyd->Hexp;
+        return;
+    }
+
+    // Compute P and Y coefficients
+    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
+        p += 2.0*hml;                        // Q*dh(Total)/dQ
+    }
+    else  hml = 0.0;
+    p = hyd->LinkFlows[k] / p;               // 1 / (dh/dQ)
+    sol->P[k] = ABS(p);
+    sol->Y[k] = p*(hpipe + hml);
+}
+
+
+void DWpipecoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  none
+**   Purpose: computes pipe head loss coeffs. for Darcy-Weisbach
+**            formula.
+**--------------------------------------------------------------
+*/
+{
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    Slink *link = &pr->network.Link[k];
+
+    double q = ABS(hyd->LinkFlows[k]);
+    double dfdq = 0.0;
+    double r, r1, f, ml, p, hloss;
+
+    ml = link->Km;                       // Minor loss coeff. 
+    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
+    if (r1*q < hyd->RQtol)
+    {
+        sol->P[k] = 1.0/hyd->RQtol;
+        sol->Y[k] = hyd->LinkFlows[k]/hyd->Hexp;
+        return;
+    }
+
+    // Compute P and Y coefficients
+    hloss = r1*SQR(q);                   // Total head loss
+    p = 2.0*r1*q;                        // |dHloss/dQ|
+    // + 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)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  returns friction factor and
+**            replaces dfdq (derivative of f w.r.t. flow)
+**   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
+          f;             // Friction factor 
+   double x1,x2,x3,x4,
+          y1,y2,y3,
+          fa,fb,r;
+   double s,w;
+  
+  hydraulics_t *hyd = &pr->hydraulics;
+  Slink *link = &pr->network.Link[k];
+
+   *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);
+      y2 = link->Kc/(3.7*link->Diam) + y1;
+      y3 = A9*log(y2);
+      f = 1.0/SQR(y3);
+      /*  *dfdq = (2.0+AA*y1/(y2*y3))*f; */   /* df/dq */
+   }
+   
+   // For Re > 2000 use Interpolation Formula
+   else if (w > A2)              
+   {
+      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;
+      *dfdq = 0.0;
+   }
+   return(f);
+
+}
+
+
+void  pumpcoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  none
+**   Purpose: computes P & Y coeffs. for pump in link k
+**--------------------------------------------------------------
+*/
+{
+    int    p;                // Pump index
+    double h0,               // Shutoff head
+           q,                // Abs. value of flow 
+           r,                // Flow resistance coeff.
+           n;                // Flow exponent coeff.
+    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
+    if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) {
+        sol->P[k] = 1.0 / CBIG;
+        sol->Y[k] = hyd->LinkFlows[k];
+        return;
+    }
+
+    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.
+    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.
+        pump->H0 = -h0;
+        pump->R = -r;
+        pump->N = 1.0;
+    }
+
+    // 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)
+    sol->P[k] = 1.0 / MAX(r, hyd->RQtol);
+    sol->Y[k] = hyd->LinkFlows[k] / n + sol->P[k] * h0;
+}
+
+
+void  curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r)
+/*
+**-------------------------------------------------------------------
+**   Input:   i   = curve index
+**            q   = flow rate
+**   Output:  *h0  = head at zero flow (y-intercept)
+**            *r  = dHead/dFlow (slope)
+**   Purpose: computes intercept and slope of head v. flow curve
+**            at current flow.
+**-------------------------------------------------------------------
+*/
+{
+    int   k1, k2, npts;
+    double *x, *y;
+    Scurve *curve;
+
+    // Remember that curve is stored in untransformed units
+    q *= pr->Ucf[FLOW];
+    curve = &pr->network.Curve[i];
+    x = curve->X;                      // x = flow
+    y = curve->Y;                      // y = head
+    npts = curve->Npts;
+
+    // Find linear segment of curve that brackets flow q
+    k2 = 0;
+    while (k2 < npts && x[k2] < q)
+        k2++;
+
+    if (k2 == 0)
+        k2++;
+
+    else if (k2 == npts)
+        k2--;
+
+    k1 = k2 - 1;
+
+    // Compute slope and intercept of this segment
+    *r = (y[k2] - y[k1]) / (x[k2] - x[k1]);
+    *h0 = y[k1] - (*r)*x[k1];
+
+    // Convert units
+    *h0 = (*h0) / pr->Ucf[HEAD];
+    *r = (*r) * pr->Ucf[FLOW] / pr->Ucf[HEAD];
+}
+
+
+void  gpvcoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  none
+**   Purpose: computes P & Y coeffs. for general purpose valve
+**--------------------------------------------------------------
+*/
+{
+    double h0,        // Headloss curve intercept
+           q,         // Abs. value of flow
+           r;         // Flow resistance coeff.
+
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+
+    // Treat as a pipe if valve closed
+    if (hyd->LinkStatus[k] == CLOSED) {
+        valvecoeff(pr, k);                          
+    }
+
+    // Otherwise utilize headloss curve
+    // whose index is stored in K
+    else {
+        // Find slope & intercept of headloss curve.
+        q = ABS(hyd->LinkFlows[k]);
+        q = MAX(q, TINY);
+        curvecoeff(pr, (int)ROUND(hyd->LinkSetting[k]), q, &h0, &r);
+
+        // Compute inverse headloss gradient (P)
+        // and flow correction factor (Y).
+        sol->P[k] = 1.0 / MAX(r, hyd->RQtol);
+        sol->Y[k] = sol->P[k] * (h0 + r*q) * SGN(hyd->LinkFlows[k]);
+    }
+}
+
+
+void  pbvcoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  none
+**   Purpose: computes P & Y coeffs. for pressure breaker valve
+**--------------------------------------------------------------
+*/
+{
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    Slink *link = &pr->network.Link[k];
+
+    // If valve fixed OPEN or CLOSED then treat as a pipe
+    if (hyd->LinkSetting[k] == MISSING || hyd->LinkSetting[k] == 0.0) {
+        valvecoeff(pr, k);         
+    }
+
+    // If valve is active
+    else {
+        // Treat as a pipe if minor loss > valve setting
+        if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k]) {
+            valvecoeff(pr, k);        
+        }
+        // Otherwise force headloss across valve to be equal to setting
+        else {
+            sol->P[k] = CBIG;
+            sol->Y[k] = hyd->LinkSetting[k] * CBIG;
+        }
+    }
+}
+
+
+void  tcvcoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k = link index
+**   Output:  none
+**   Purpose: computes P & Y coeffs. for throttle control valve
+**--------------------------------------------------------------
+*/
+{
+    double km;
+    hydraulics_t *hyd = &pr->hydraulics;
+    Slink *link = &pr->network.Link[k];
+
+    // Save original loss coeff. for open valve
+    km = link->Km;
+
+    // If valve not fixed OPEN or CLOSED, compute its loss coeff.
+    if (hyd->LinkSetting[k] != MISSING) {
+        link->Km = 0.02517 * hyd->LinkSetting[k] / (SQR(link->Diam)*SQR(link->Diam));
+    }
+
+    // Then apply usual valve formula
+    valvecoeff(pr, k);                                             
+
+    // Restore original loss coeff.
+    link->Km = km;
+}
+
+
+void  prvcoeff(EN_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 pressure
+**            reducing valves
+**--------------------------------------------------------------
+*/
+{
+    int   i, j;                        // Rows of solution matrix
+    double hset;                       // Valve head setting
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+
+    i = sol->Row[n1];                  // Matrix rows of nodes
+    j = sol->Row[n2];
+    hset = pr->network.Node[n2].El +
+           hyd->LinkSetting[k];        // Valve setting
+
+    if (hyd->LinkStatus[k] == ACTIVE)
+    {
+
+        // Set coeffs. to force head at downstream
+        // node equal to valve setting & force flow 
+        // to equal to flow imbalance at downstream node.
+
+        sol->P[k] = 0.0;
+        sol->Y[k] = hyd->LinkFlows[k] + hyd->X_tmp[n2];    // Force flow balance
+        sol->F[j] += (hset * CBIG);                        // Force head = hset
+        sol->Aii[j] += CBIG;                               // at downstream node
+        if (hyd->X_tmp[n2] < 0.0) {
+            sol->F[i] += hyd->X_tmp[n2];
+        }
+        return;
+    }
+
+    // For OPEN, CLOSED, or XPRESSURE valve
+    // compute matrix coeffs. using the valvecoeff() function.
+
+    valvecoeff(pr, k);
+    sol->Aij[sol->Ndx[k]] -= sol->P[k];
+    sol->Aii[i] += sol->P[k];
+    sol->Aii[j] += sol->P[k];
+    sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]);
+    sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]);
+}
+
+
+void  psvcoeff(EN_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 pressure
+**            sustaining valve
+**--------------------------------------------------------------
+*/
+{
+    int   i, j;                        // Rows of solution matrix
+    double hset;                       // Valve head setting
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+
+    i = sol->Row[n1];                  // Matrix rows of nodes
+    j = sol->Row[n2];
+    hset = pr->network.Node[n1].El +
+           hyd->LinkSetting[k];        // Valve setting
+
+    if (hyd->LinkStatus[k] == ACTIVE)
+    {
+        // Set coeffs. to force head at upstream
+        // node equal to valve setting & force flow 
+        // equal to flow imbalance at upstream node.
+
+        sol->P[k] = 0.0;
+        sol->Y[k] = hyd->LinkFlows[k] - hyd->X_tmp[n1];    // Force flow balance
+        sol->F[i] += (hset * CBIG);                        // Force head = hset
+        sol->Aii[i] += CBIG;                               // at upstream node
+        if (hyd->X_tmp[n1] > 0.0) {
+            sol->F[j] += hyd->X_tmp[n1];
+        }
+        return;
+    }
+
+    // For OPEN, CLOSED, or XPRESSURE valve
+    // compute matrix coeffs. using the valvecoeff() function.
+
+    valvecoeff(pr, k);
+    sol->Aij[sol->Ndx[k]] -= sol->P[k];
+    sol->Aii[i] += sol->P[k];
+    sol->Aii[j] += sol->P[k];
+    sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]);
+    sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]);
+}
+
+
+void  fcvcoeff(EN_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
+**--------------------------------------------------------------
+*/
+{
+    int   i, j;                   // Rows in solution matrix
+    double q;                     // Valve flow setting
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+
+    q = hyd->LinkSetting[k];
+    i = hyd->solver.Row[n1];
+    j = hyd->solver.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)
+    {
+        hyd->X_tmp[n1] -= q;
+        sol->F[i] -= q;
+        hyd->X_tmp[n2] += q;
+        sol->F[j] += q;
+        sol->P[k] = 1.0 / CBIG;
+        sol->Aij[sol->Ndx[k]] -= sol->P[k];
+        sol->Aii[i] += sol->P[k];
+        sol->Aii[j] += sol->P[k];
+        sol->Y[k] = hyd->LinkFlows[k] - q;
+    }
+
+    // Otherwise treat valve as an open pipe
+
+    else
+    {
+        valvecoeff(pr, k);                                          
+        sol->Aij[sol->Ndx[k]] -= sol->P[k];
+        sol->Aii[i] += sol->P[k];
+        sol->Aii[j] += sol->P[k];
+        sol->F[i] += (sol->Y[k] - hyd->LinkFlows[k]);
+        sol->F[j] -= (sol->Y[k] - hyd->LinkFlows[k]);
+    }
+}
+
+
+void valvecoeff(EN_Project *pr, int k)
+/*
+**--------------------------------------------------------------
+**   Input:   k    = link index
+**   Output:  none
+**   Purpose: computes solution matrix coeffs. for a completely
+**            open, closed, or throttled control valve.
+**--------------------------------------------------------------
+*/
+{
+    double p;
+
+    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];
+
+    // Valve is closed. Use a very small matrix coeff.
+    if (hyd->LinkStatus[k] <= CLOSED)
+    {
+        sol->P[k] = 1.0 / CBIG;
+        sol->Y[k] = flow;
+        return;
+    }
+
+    // Account for any minor headloss through the valve
+    if (link->Km > 0.0)
+    {
+        p = 2.0 * link->Km * fabs(flow);
+        if (p < hyd->RQtol) {
+            p = hyd->RQtol;
+        }
+        sol->P[k] = 1.0 / p;
+        sol->Y[k] = flow / 2.0;
+    }
+    else
+    {
+        sol->P[k] = 1.0 / hyd->RQtol;
+        sol->Y[k] = flow;
+    }
+}

src/inpfile.c

@@ -618,6 +618,12 @@ int saveinpfile(EN_Project *pr, char *fname)
   fprintf(f, "\n CHECKFREQ           %-d", hyd->CheckFreq);
   fprintf(f, "\n MAXCHECK            %-d", hyd->MaxCheck);
   fprintf(f, "\n DAMPLIMIT           %-.8f", hyd->DampLimit);
+  if (hyd->HeadErrorLimit > 0.0) {
+      fprintf(f, "\n HEADERROR           %-.8f", hyd->HeadErrorLimit * pr->Ucf[HEAD]);
+  }
+  if (hyd->FlowChangeLimit > 0.0) {
+      fprintf(f, "\n FLOWCHANGE          %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]);
+  }
 
   /* Write [REPORT] section */
 

src/input1.c

@@ -127,6 +127,10 @@ void setdefaults(EN_Project *pr)
   hyd->Htol = HTOL;        /* Default head tolerance         */
   hyd->Qtol = QTOL;        /* Default flow tolerance         */
   hyd->Hacc = HACC;        /* Default hydraulic accuracy     */
+
+  hyd->FlowChangeLimit = 0.0; /* Default flow change limit   */
+  hyd->HeadErrorLimit = 0.0;  /* Default head error limit    */
+
   qu->Ctol = MISSING;      /* No pre-set quality tolerance   */
   hyd->MaxIter = MAXITER;  /* Default max. hydraulic trials  */
   hyd->ExtraIter = -1;     /* Stop if network unbalanced     */
@@ -604,6 +608,10 @@ void convertunits(EN_Project *pr)
     tank->V1max *= tank->Vmax;
   }
 
+  /* Convert hydraulic convergence criteria */
+  hyd->FlowChangeLimit /= pr->Ucf[FLOW];
+  hyd->HeadErrorLimit  /= pr->Ucf[HEAD];
+
   /* Convert WQ option concentration units */
   qu->Climit /= pr->Ucf[QUALITY];
   qu->Ctol /= pr->Ucf[QUALITY];
@@ -676,8 +684,9 @@ void convertunits(EN_Project *pr)
       
     }
 
+//////  Moved to inithyd() in hydraul.c  ///////
     /* Compute flow resistances */
-    resistance(pr, k);
+    //resistance(pr, k);
   }
 
   /* Convert units on control settings */

src/input3.c

@@ -1913,6 +1913,10 @@ int optionvalue(EN_Project *pr, int n)
 **    SPECIFIC GRAVITY    value
 **    TRIALS              value
 **    ACCURACY            value
+
+**    HEADLIMIT           value
+**    FLOWLIMIT           value
+
 **    TOLERANCE           value
 **    SEGMENTS            value  (not used)
 **  ------ Undocumented Options -----
@@ -1928,18 +1932,20 @@ int optionvalue(EN_Project *pr, int n)
   hydraulics_t *hyd = &pr->hydraulics;
   quality_t *qu = &pr->quality;
   parser_data_t *par = &pr->parser;
+  char* tok0 = par->Tok[0];
 
   
   int nvalue = 1; /* Index of token with numerical value */
   double y;
 
   /* Check for obsolete SEGMENTS keyword */
-  if (match(par->Tok[0], w_SEGMENTS))
+  //if (match(par->Tok[0], w_SEGMENTS))
+  if (match(tok0, w_SEGMENTS))
           return (0);
 
   /* Check for missing value (which is permissible) */
-  if (match(par->Tok[0], w_SPECGRAV) || match(par->Tok[0], w_EMITTER) ||
+  if (match(tok0, w_SPECGRAV) || match(tok0, w_EMITTER) ||
-      match(par->Tok[0], w_DEMAND))
+      match(tok0, w_DEMAND))
     nvalue = 2;
   if (n < nvalue)
     return (0);
@@ -1949,7 +1955,7 @@ int optionvalue(EN_Project *pr, int n)
     return (213);
 
   /* Check for WQ tolerance option (which can be 0) */
-  if (match(par->Tok[0], w_TOLERANCE)) {
+  if (match(tok0, w_TOLERANCE)) {
     if (y < 0.0)
       return (213);
     qu->Ctol = y; /* Quality tolerance*/
@@ -1957,7 +1963,7 @@ int optionvalue(EN_Project *pr, int n)
   }
 
   /* Check for Diffusivity option */
-  if (match(par->Tok[0], w_DIFFUSIVITY)) {
+  if (match(tok0, w_DIFFUSIVITY)) {
     if (y < 0.0)
       return (213);
     qu->Diffus = y;
@@ -1965,42 +1971,59 @@ int optionvalue(EN_Project *pr, int n)
   }
 
   /* Check for Damping Limit option */ 
-  if (match(par->Tok[0], w_DAMPLIMIT)) {
+  if (match(tok0, w_DAMPLIMIT)) {
     hyd->DampLimit = y;
     return (0);
   }
   
+  /* Check for flow change limit*/
+  else if (match(tok0, w_FLOWCHANGE))
+  {
+      if (y < 0.0) return 213;
+      hyd->FlowChangeLimit = y;
+      return 0;
+  }
+
+  /* Check for head error limit*/
+  else if (match(tok0, w_HEADERROR))
+  {
+      if (y < 0.0) return 213;
+      hyd->HeadErrorLimit = y;
+      return 0;
+  }
+
   /* All other options must be > 0 */
   if (y <= 0.0)
     return (213);
 
   /* Assign value to specified option */
-  if (match(par->Tok[0], w_VISCOSITY))
+  if (match(tok0, w_VISCOSITY))
     hyd->Viscos = y; /* Viscosity */
-  else if (match(par->Tok[0], w_SPECGRAV))
+  else if (match(tok0, w_SPECGRAV))
     hyd->SpGrav = y; /* Spec. gravity */
-  else if (match(par->Tok[0], w_TRIALS))
+  else if (match(tok0, w_TRIALS))
     hyd->MaxIter = (int)y;                   /* Max. trials */
-  else if (match(par->Tok[0], w_ACCURACY)) /* Accuracy */
+  else if (match(tok0, w_ACCURACY)) /* Accuracy */
   {
     y = MAX(y, 1.e-5);
     y = MIN(y, 1.e-1);
     hyd->Hacc = y;
-  } else if (match(par->Tok[0], w_HTOL))
+  }
+  else if (match(tok0, w_HTOL))
     hyd->Htol = y;
-  else if (match(par->Tok[0], w_QTOL))
+  else if (match(tok0, w_QTOL))
     hyd->Qtol = y;
-  else if (match(par->Tok[0], w_RQTOL)) {
+  else if (match(tok0, w_RQTOL)) {
     if (y >= 1.0)
       return (213);
     hyd->RQtol = y;
-  } else if (match(par->Tok[0], w_CHECKFREQ))
+  } else if (match(tok0, w_CHECKFREQ))
     hyd->CheckFreq = (int)y;
-  else if (match(par->Tok[0], w_MAXCHECK))
+  else if (match(tok0, w_MAXCHECK))
     hyd->MaxCheck = (int)y;
-  else if (match(par->Tok[0], w_EMITTER))
+  else if (match(tok0, w_EMITTER))
     hyd->Qexp = 1.0 / y;
-  else if (match(par->Tok[0], w_DEMAND))
+  else if (match(tok0, w_DEMAND))
     hyd->Dmult = y;
   else
     return (201);

src/output.c

@@ -258,54 +258,69 @@ int saveenergy(EN_Project *pr)
   parser_data_t  *par = &pr->parser;
   time_options_t *time = &pr->time_options;
   FILE           *outFile = out->OutFile;
+  Spump *pump;
   
   int i, j;
   INT4 index;
-  REAL4 x[6];        /* work array */
+  REAL4 x[MAX_ENERGY_STATS]; // work array
-  double hdur, /* total time->Duration in hours */
+  double hdur,               // total simulation duration in hours
-      t;             /* pumping time->Duration */
+         t;                  // total pumping time duration
 
   hdur = time->Dur / 3600.0;
   for (i = 1; i <= net->Npumps; i++) {
-    Spump *pump = &net->Pump[i];
+    pump = &net->Pump[i];
     if (hdur == 0.0) {
-      for (j = 0; j < 5; j++)
+        pump->Energy[TOTAL_COST] *= 24.0;
-        x[j] = (REAL4)pump->Energy[j];
+    }
-      x[5] = (REAL4)(pump->Energy[5] * 24.0);
+    else {
-    } else {
+        // ... convert total hrs. online to fraction of total time online
-      t = pump->Energy[0];
+        t = pump->Energy[PCNT_ONLINE];  //currently holds total hrs. online
-      x[0] = (REAL4)(t / hdur);
+        pump->Energy[PCNT_ONLINE] = t / hdur;
-      x[1] = 0.0f;
+
-      x[2] = 0.0f;
+        // ... convert cumulative values to time-averaged ones
-      x[3] = 0.0f;
-      x[4] = 0.0f;
         if (t > 0.0) {
-        x[1] = (REAL4)(pump->Energy[1] / t);
+            pump->Energy[PCNT_EFFIC]   /= t;
-        x[2] = (REAL4)(pump->Energy[2] / t);
+            pump->Energy[KWH_PER_FLOW] /= t;
-        x[3] = (REAL4)(pump->Energy[3] / t);
+            pump->Energy[TOTAL_KWH]    /= t;
         }
-      x[4] = (REAL4)pump->Energy[4];
+
-      x[5] = (REAL4)(pump->Energy[5] * 24.0 / hdur);
+        // ... convert total cost to cost per day
+        pump->Energy[TOTAL_COST] *= 24.0 / hdur;
     }
-    x[0] *= 100.0f;
+
-    x[1] *= 100.0f;
+    // ... express time online and avg. efficiency as percentages
-    /* Compute Kw-hr per MilGal (or per cubic meter) */
+    pump->Energy[PCNT_ONLINE] *= 100.0;
-    if (par->Unitsflag == SI)
+    pump->Energy[PCNT_EFFIC] *= 100.0;
-      x[2] *= (REAL4)(1000.0 / LPSperCFS / 3600.0);
+
-    else
+    // ... compute KWH per Million Gallons or per Cubic Meter
-      x[2] *= (REAL4)(1.0e6 / GPMperCFS / 60.0);
+    if (par->Unitsflag == SI) {
-    for (j = 0; j < 6; j++)
+        pump->Energy[KWH_PER_FLOW] *= (1000. / LPSperCFS / 3600.);
-      pump->Energy[j] = x[j];
+    }
+    else {
+        pump->Energy[KWH_PER_FLOW] *= (1.0e6 / GPMperCFS / 60.);
+    }
+
+    // ... save energy stats to REAL4 work array
+    for (j = 0; j < MAX_ENERGY_STATS; j++) {
+        x[j] = (REAL4)pump->Energy[j];
+    }
+
+    // ... save energy results to output file
     index = pump->Link;
-    if (fwrite(&index, sizeof(INT4), 1, outFile) < 1)
+    if (fwrite(&index, sizeof(INT4), 1, outFile) < 1) {
-      return (308);
-    if (fwrite(x, sizeof(REAL4), 6, outFile) < 6)
       return (308);
     }
+    if (fwrite(x, sizeof(REAL4), MAX_ENERGY_STATS, outFile) < MAX_ENERGY_STATS) {
+      return (308);
+    }
+  }
+
+  // ... compute and save demand charge
   hyd->Emax = hyd->Emax * hyd->Dcost;
   x[0] = (REAL4)hyd->Emax;
-  if (fwrite(&x[0], sizeof(REAL4), 1, outFile) < 1)
+  if (fwrite(&x[0], sizeof(REAL4), 1, outFile) < 1) {
     return (308);
+  }
   return (0);
 }
 

src/report.c

@@ -229,6 +229,15 @@ void writesummary(EN_Project *pr)
   sprintf(s, FMT27, hyd->Hacc);
   writeline(pr, s);
 
+  if (hyd->HeadErrorLimit > 0.0) {
+      sprintf(s, FMT27d, hyd->HeadErrorLimit*pr->Ucf[HEAD], rep->Field[HEAD].Units);
+      writeline(pr, s);
+  }
+  if (hyd->FlowChangeLimit > 0.0) {
+      sprintf(s, FMT27e, hyd->FlowChangeLimit*pr->Ucf[FLOW], rep->Field[FLOW].Units);
+      writeline(pr, s);
+  }
+
   sprintf(s, FMT27a, hyd->CheckFreq); 
   writeline(pr, s);                   
   sprintf(s, FMT27b, hyd->MaxCheck);       
@@ -378,6 +387,7 @@ void writeenergy(EN_Project *pr)
   EN_Network *net = &pr->network;
   hydraulics_t *hyd = &pr->hydraulics;
   report_options_t *rep = &pr->report;
+  Spump *pump;
 
   int j;
   double csum;
@@ -388,14 +398,14 @@ void writeenergy(EN_Project *pr)
   writeheader(pr,ENERHDR, 0);
   csum = 0.0;
   for (j = 1; j <= net->Npumps; j++) {
-    Spump *pump = &net->Pump[j];
+    pump = &net->Pump[j];
-    csum += pump->Energy[5];
+    csum += pump->Energy[TOTAL_COST];
     if (rep->LineNum == (long)rep->PageSize)
       writeheader(pr, ENERHDR, 1);
     sprintf(s, "%-8s  %6.2f %6.2f %9.2f %9.2f %9.2f %9.2f",
-            net->Link[pump->Link].ID, pump->Energy[0], pump->Energy[1],
+            net->Link[pump->Link].ID, pump->Energy[PCNT_ONLINE], pump->Energy[PCNT_EFFIC],
-            pump->Energy[2], pump->Energy[3], pump->Energy[4],
+            pump->Energy[KWH_PER_FLOW], pump->Energy[TOTAL_KWH], pump->Energy[MAX_KW],
-            pump->Energy[5]);
+            pump->Energy[TOTAL_COST]);
     writeline(pr, s);
   }
   fillstr(s, '-', 63);

src/text.h

@@ -142,6 +142,9 @@ AUTHOR:     L. Rossman
 #define   w_MAXCHECK    "MAXCHECK"
 #define   w_DAMPLIMIT   "DAMPLIMIT"
 
+#define   w_FLOWCHANGE  "FLOWCHANGE"
+#define   w_HEADERROR   "HEADERROR"
+
 #define   w_SECONDS     "SEC"
 #define   w_MINUTES     "MIN"
 #define   w_HOURS       "HOU"
@@ -366,6 +369,9 @@ AUTHOR:     L. Rossman
 #define FMT27b "    Maximum Trials Checked ............ %-d"                   
 #define FMT27c "    Damping Limit Threshold ........... %-.6f"
 
+#define FMT27d "    Headloss Error Limit .............. %-.6f %s"
+#define FMT27e "    Flow Change Limit ................. %-.6f %s"
+
 #define FMT28  "    Maximum Trials .................... %-d"
 #define FMT29  "    Quality Analysis .................. None"
 #define FMT30  "    Quality Analysis .................. %s"
@@ -409,10 +415,13 @@ AUTHOR:     L. Rossman
 #define FMT61  "%10s: Valve %s caused ill-conditioning"
 #define FMT62  "%10s: System ill-conditioned at node %s"
 #define FMT63  "%10s: %s %s changed by rule %s"
-#define FMT64  "%10s: Balancing the network:"
+#define FMT64  "%10s: Balancing the network:\n"
 #define FMT65  "            Trial %2d: relative flow change = %-.6f"
 /*** End of update ***/
 
+#define FMT66  "                      maximum  flow change = %.4f for Link %s"
+#define FMT67  "                      maximum  head error  = %.4f for Link %s\n"
+
 /* -------------------- Energy Report Table ------------------- */
 #define FMT71  "Energy Usage:"
 #define FMT72  \

src/types.h

@@ -59,6 +59,9 @@ typedef  int          INT4;
 #define   TINY      1.E-6
 #define   MISSING   -1.E10
 
+#define   CBIG      1.e8     /* Big coefficient         */
+#define   CSMALL    1.e-6    /* Small coefficient       */
+
 #ifdef M_PI
   #define   PI        M_PI
 #else
@@ -169,8 +172,8 @@ typedef enum {
   LOWLEVEL,      /*    act when grade below set level   */
   HILEVEL,       /*    act when grade above set level   */
   TIMER,         /*    act when set time reached        */
-  TIMEOFDAY
+  TIMEOFDAY      /*    act when time of day occurs      */
-} ControlType;    /*    act when time of day occurs      */
+} ControlType;    
 
 typedef enum {
   XHEAD,        /*   pump cannot deliver head (closed) */
@@ -188,13 +191,13 @@ typedef enum {
 typedef enum {
   HW,           /*   Hazen-Williams                    */
   DW,           /*   Darcy-Weisbach                    */
-  CM
+  CM            /*   Chezy-Manning                     */
-} FormType;          /*   Chezy-Manning                     */
+} FormType;    
 
 typedef enum {
   US,           /*   US                                */
-  SI
+  SI            /*   SI (metric)                       */
-} UnitsType;          /*   SI (metric)                       */
+} UnitsType;          
 
 typedef enum {
   CFS,          /*   cubic feet per second             */
@@ -206,35 +209,36 @@ typedef enum {
   LPM,          /*   liters per minute                 */
   MLD,          /*   megaliters per day                */
   CMH,          /*   cubic meters per hour             */
-  CMD
+  CMD           /*   cubic meters per day              */
-} FlowUnitsType;         /*   cubic meters per day              */
+} FlowUnitsType;         
 
 typedef enum {
   PSI,          /*   pounds per square inch            */
   KPA,          /*   kiloPascals                       */
-  METERS
+  METERS        /*   meters                            */
-} PressUnitsType;      /*   meters                            */
 
+} PressUnitsType;      
 typedef enum {
   LOW,          /*   lower limit                       */
   HI,           /*   upper limit                       */
-  PREC
+  PREC          /*   precision                         */
-} RangeType;        /*   precision                         */
 
+} RangeType;        
 typedef enum {
   MIX1,         /*   1-compartment model               */
   MIX2,         /*   2-compartment model               */
   FIFO,         /*   First in, first out model         */
-  LIFO
+  LIFO          /*   Last in, first out model          */
-} MixType;        /*   Last in, first out model          */
+} MixType;       
 
 typedef enum {
   SERIES,       /*   none                              */
   AVG,          /*   time-averages                     */
   MIN,          /*   minimum values                    */
   MAX,          /*   maximum values                    */
-  RANGE
+  RANGE         /*   max - min values                  */
-} TstatType;       /*   max - min values                  */
+} TstatType;
+
 
 #define MAXVAR   21 /* Max. # types of network variables */        
                     /* (equals # items enumed below)   */
@@ -273,17 +277,27 @@ typedef enum {
 } SectType;
 
 typedef enum {
-  STATHDR,      /*  Hydraulic Status       */
+  STATHDR,      /*  Hydraulic Status header  */
-  ENERHDR,      /*  Energy Usage           */
+  ENERHDR,      /*  Energy Usage header      */
-  NODEHDR,      /*  Node Results           */
+  NODEHDR,      /*  Node Results header      */
-  LINKHDR
+  LINKHDR       /*  Link Results header      */
-} HdrType;     /*  Link Results           */
+} HdrType;    
 
 typedef enum {
   POSITIVE,
   NEGATIVE
 } FlowDirection;
 
+typedef enum {
+    PCNT_ONLINE,
+    PCNT_EFFIC,
+    KWH_PER_FLOW,
+    TOTAL_KWH,
+    MAX_KW,
+    TOTAL_COST,
+    MAX_ENERGY_STATS
+} EnergyStats;
+
 /*
 ------------------------------------------------------
    Global Data Structures                             
@@ -344,6 +358,16 @@ struct Sdemand            /* DEMAND CATEGORY OBJECT */
 };
 typedef struct Sdemand *Pdemand; /* Pointer to demand object */
 
+typedef struct
+{
+    double hrsOnLine;        // hours pump is online
+    double efficiency;       // total time wtd. efficiency
+    double kwHrsPerCFS;      // total kw-hrs per cfs of flow
+    double kwHrs;            // total kw-hrs consumed
+    double maxKwatts;        // max. kw consumed
+    double totalCost;        // total pumping cost
+} Senergy;
+
 struct Ssource     /* WQ SOURCE OBJECT */
 {
  /*int   Node;*/     /* Node index of source     */
@@ -422,13 +446,7 @@ typedef struct     /* PUMP OBJECT */
    int    Upat;     /* Utilization pattern index   */
    int    Epat;     /* Energy cost pattern index   */
    double Ecost;    /* Unit energy cost            */
-   double Energy[6];  /* Energy usage statistics:  */
+   double Energy[MAX_ENERGY_STATS];  /* Energy usage statistics     */
-                     /* 0 = pump utilization      */
-                     /* 1 = avg. efficiency       */
-                     /* 2 = avg. kW/flow          */
-                     /* 3 = avg. kwatts           */
-                     /* 4 = peak kwatts           */
-                     /* 5 = cost/day              */
 }  Spump;
 
 typedef struct     /* VALVE OBJECT */
@@ -699,7 +717,6 @@ typedef struct {
   Premise *Plast;     /* Previous premise clause */
   Action *Tlast;      /* Previous true action */
   Action *Flast;      /* Previous false action */
-  
 } rules_t;
 
 /*
@@ -754,6 +771,10 @@ typedef struct {
   Qexp,                  /* Exponent in orifice formula  */
   Dmult,                 /* Demand multiplier            */
   Hacc,                  /* Hydraulics solution accuracy */
+
+  FlowChangeLimit,       /* Hydraulics flow change limit */
+  HeadErrorLimit,        /* Hydraulics head error limit  */
+
   DampLimit,             /* Solution damping threshold   */
   Viscos,                /* Kin. viscosity (sq ft/sec)   */
   SpGrav,                /* Specific gravity             */
@@ -819,7 +840,8 @@ typedef struct {
   Scurve   *Curve;       /* Curve data                   */
   Scoord   *Coord;       /* Coordinate data              */
   Scontrol *Control;     /* Control data                 */
-  ENHashTable  *NodeHashTable, *LinkHashTable;            /* Hash tables for ID labels    */
+  ENHashTable *NodeHashTable,
+              *LinkHashTable; /* Hash tables for ID labels    */
   Padjlist *Adjlist;          /* Node adjacency lists         */
   
 } EN_Network;
@@ -852,7 +874,4 @@ struct EN_Project {
   
 };
 
-
-
-
 #endif

tools/before-test.cmd

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

tools/run-nrtest.cmd

@@ -18,7 +18,7 @@ setlocal
 set NRTEST_SCRIPT_PATH=%~1
 set TEST_SUITE_PATH=%~2
 
-set BENCHMARK_VER=2012vs10
+set BENCHMARK_VER=220dev-vs17
 
 
 set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute

win_build/WinSDK/Makefile.bat

@@ -21,9 +21,9 @@ Find /i "x86" < checkOS.tmp > StringCheck.tmp
 If %ERRORLEVEL% == 1 (
 	CALL "%SDK_PATH%bin\"SetEnv.cmd /x64 /release
 	rem : create EPANET2.DLL
-	cl -o epanet2.dll epanet.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL 
+	cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL 
 	rem : create EPANET2.EXE
-	cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
+	cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
 	md "%Build_PATH%"\64bit
 	move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\64bit
 	move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\64bit
@@ -35,9 +35,9 @@ rem : 32 bit with DEF
 CALL "%SDK_PATH%bin\"SetEnv.cmd /x86 /release
 echo "32 bit with epanet2.def mapping"
 rem : create EPANET2.DLL
-cl -o epanet2.dll epanet.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
+cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /link /DLL /def:..\win_build\WinSDK\epanet2.def /MAP
 rem : create EPANET2.EXE
-cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
+cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c /I ..\include /I ..\run /I ..\src /link
 md "%Build_PATH%"\32bit
 move /y "%SRC_PATH%"\*.dll "%Build_PATH%"\32bit
 move /y "%SRC_PATH%"\*.exe "%Build_PATH%"\32bit