Pressure Dependent Demands added to address issue 163

ReleaseNotes2_2.md

@@ -0,0 +1,96 @@
+
+
+
+Release Notes for EPANET 2.2 (Draft)
+============================
+
+This document describes the changes and updates that have been made to version 2.2 of EPANET.
+
+## Thread-Safe API Functions
+
+A duplicate set of the version 2.1 API functions has been provided that allow multiple EPANET projects to be analyzed concurrently in a thread-safe manner. These functions maintain the same name as the original but use a `EN_` prefix instead of `EN`. In addition, the first argument to each of these functions is a pointer to an `EN_Project` structure that encapsulates the network data for the particular project being analyzed. For example, instead of writing:
+
+`ENgetnodevalue(nodeIndex, EN_ELEVATION, &elev)`
+
+one would use:
+
+`EN_getnodevalue(pr, nodeIndex, EN_ELEVATION, &elev)`
+
+where `pr` is the pointer to an `EN_Project`.
+
+Two new functions have been added to the API to manage the creation and deletion of project pointers. `EN_createproject` creates a new project along with a pointer to it, while `EN_deleteproject` deletes a project. An example of using the thread-safe version of the API is shown below:
+```
+#include "epanet2.h"
+int runEpanet(char *finp, char *frpt)
+{
+    EN_Project *pr = NULL;
+    int err;
+    err = EN_createproject(&pr);
+    if (err) return err;
+    err = EN_open(pr, finp, frpt, "");
+    if (!err) err = EN_solveH(pr);
+    if (!err) err = EN_report(pr);
+    EN_close(pr);
+    EN_deleteproject(pr);
+    return err;
+}
+```
+
+## Additional Convergence Parameters
+
+Two new analysis options have been added to provide more rigorous convergence criteria for EPANET's hydraulic solver. In the API they are named `EN_HEADERROR` and `EN_FLOWCHANGE` while in the `[OPTIONS]` section of an EPANET input file they are named `HEADERROR` and `FLOWCHANGE`, respectively.
+
+`EN_HEADERROR` is the maximum head loss error that any network link can have for hydraulic convergence to occur. A link's head loss error is the difference between the head loss found as a function of computed flow in the link (such as by the Hazen-Williams equation for a pipe) and the difference in computed heads for the link's end nodes. The units of this parameter are feet (or meters for SI units). The default value of 0 indicates that no head error limit applies.
+
+`EN_FLOWCHANGE` is the largest change in flow that any network element (link, emitter, or pressure-dependent demand) can have for hydraulic convergence to occur. It is specified in whatever flow units the project is using. The default value of 0 indicates that no flow change limit applies.
+
+These new parameters augment the current `EN_ACCURACY` option which always remains in effect. In addition, both `EN_HEADERROR` and `EN_FLOWCHANGE` can be used as parameters in the `ENgetstatistic` (or `EN_getstatistic`) function to retrieve their computed values (even when their option values are 0) after a hydraulic solution has been completed.  
+
+## Improved Linear Solver Routine
+EPANET's hydraulic solver requires solving a system of linear equations over a series of iterations until a set of convergence criteria are met. The coefficient matrix of this linear system is square and symmetric. It has a row for each network node and a non-zero off-diagonal coefficient for each link. The numerical effort needed to solve the linear system can be reduced if the nodes are re-ordered so that the non-zero coefficients cluster more tightly around the diagonal.
+
+EPANET's original node re-ordering scheme has been replaced by the more powerful **Multiple Minimum Degree (MMD)** algorithm. On a series of eight networks ranging in size from 7,700 to 100,000 nodes **MMD** reduced the solution time for a single period (steady state) hydraulic analysis by an average of more than 50%.
+
+## Pressure Dependent Demands
+EPANET has always employed a Demand Driven Analysis (**DDA**) when modeling network hydraulics. Under this approach nodal demands at a given point in time are fixed values that must be delivered no matter what nodal heads and link flows are produced by a hydraulic solution. This can result in situations where required demands are satisfied at nodes that have negative pressures - a physical impossibility. 
+
+To address this issue EPANET has been extended to use a Pressure Driven Analysis (**PDA**) if so desired. Under **PDA**, the demand *D* delivered at a node depends on the node's available pressure *P* according to:
+$$D =D_f\left(\frac{P-P_{min}}{P_{req}-P_{min}}\right)^{P_{exp}} for P_{0}<=P<=P_{req}$$where *D<sub>f</sub>* is the full demand required, *P<sub>min</sub>* is the pressure below which demand is zero, *P<sub>req</sub>* is the pressure required to deliver the full required demand and *P<sub>exp</sub>* is an exponent. When *P < P<sub>min</sub>* demand is 0 and when *P > P<sub>req</sub>* demand equals *D<sub>f</sub>*.
+
+To implement pressure dependent analysis four new parameters have been added to the [OPTIONS] section of the EPANET input file:
+
+
+| Parameter | Description  | Default |
+|--|--|--|
+| DEMAND MODEL | either DDA or PDA | DDA |
+| MINIMUM PRESSURE | value for *P<sub>min</sub>* | 0
+| REQUIRED PRESSURE | value for *P<sub>req</sub>* | 0
+| PRESSURE EXPONENT | value for *P<sub>exp</sub>* | 0.5 |
+
+These parameters can also be set and retrieved in code using the following API functions
+```
+int ENsetdemandmodel(int modelType, double pMin, double pReq, double pExp);
+int ENgetdemandmodel(int *modelType, double *pMin, double *pReq, double *pExp);
+```
+for the legacy API and
+```
+int EN_setdemandmodel(EN_Project *pr, int modelType, double pMin, double pReq, double pExp);
+int EN_getdemandmodel(EN_Project *pr, int *modelType, double *pMin, double *pReq, double *pExp);
+```
+for the thread-safe API. Some additional points regarding the new **PDA** option are:
+
+ - If no DEMAND  MODEL and its parameters are specified then the analysis defaults to being demand driven (**DDA**).
+ - This implementation of **PDA** assumes that the same parameters apply to all nodes in the network. Extending the framework to allow different parameters for specific nodes is straightforward to do but is left as a future feature to implement.
+ - *P<sub>0</sub>* is allowed to equal to *P<sub>req</sub>*. This condition can be used to find a solution that results in the smallest amount of demand reductions needed to insure that no node delivers positive demand at a pressure below *P<sub>min</min>*.
+
+## Code Changes
+
+ - The header file `vars.h` containing global variables has been eliminated. Instead a number of new structures incorporating these variables has been added to `types.h`. These structures have been incorporated into the new `EN_Project` structure, also defined in `types.h`, which gets passed into each of the thread-safe API functions as a pointer.
+ - Each of the legacy API functions now simply calls its thread-safe counterpart passing in a pointer to a default global`EN_Project` variable that is declared in `types.h`.
+ -  Throughout all code modules, global variables that were previously accessed through `vars.h` are now accessed using the `EN_Project` pointer that is passed into the functions where the variables appear.
+ - The exceedingly long `hydraul.c` file has been split into four separate files:
+     - `hydraul.c` now contains just the code needed to initialize a hydraulic analysis, set demands and control actions at each time step, and determine the length of the next time step to take.
+     - `hydsolver.c` implements EPANET's hydraulic solver at a single point in time.
+     - `hydcoeffs.c` computes values of the matrix coefficients (derived from link head losses and their gradients) used by the hydraulic solver.
+     - `hydstatus.c` checks for status changes in valves and pumps as requested by the hydraulic solver.
+ - The Multiple Minimum Degree re-ordering algorithm appears in a new file named `genmmd.c`. This is 1990's legacy code that is readily available on the web and can be found in several linear equation solver libraries.

include/epanet2.h

@@ -138,10 +138,11 @@ typedef enum {
   EN_NEXTEVENTIDX = 15
 } EN_TimeProperty;
 
-
 typedef enum {
   EN_ITERATIONS    = 0,
-  EN_RELATIVEERROR = 1
+  EN_RELATIVEERROR = 1,
+  EN_MAXHEADERROR  = 2,
+  EN_MAXFLOWCHANGE = 3
 } EN_AnalysisStatistic;
 
 typedef enum {
@@ -150,7 +151,7 @@ typedef enum {
   EN_LINKCOUNT    = 2,   /**< Number of Links (Pipes + Pumps + Valves) */
   EN_PATCOUNT     = 3,   /**< Number of Time Patterns */
   EN_CURVECOUNT   = 4,   /**< Number of Curves */
-  EN_CONTROLCOUNT = 5,    /**< Number of Control Statements */
+  EN_CONTROLCOUNT = 5,   /**< Number of Control Statements */
   EN_RULECOUNT	  = 6    /**< Number of Rule-based Control Statements */
 } EN_CountType;
 
@@ -205,6 +206,10 @@ typedef enum {
   EN_CMD         = 9
 } EN_FlowUnits;
 
+typedef enum {           /* Demand model types. */
+  EN_DDA         = 0,   /**< Demand driven analysis */
+  EN_PDA         = 1    /**< Pressure driven analysis */
+} EN_DemandModel;
 
 /// Simulation Option codes
 typedef enum {
@@ -224,8 +229,6 @@ typedef enum {
   EN_TIMEOFDAY   = 3
 } EN_ControlType;
 
-
-
 typedef enum {
   EN_AVERAGE     = 1,   /* Time statistic types.    */
   EN_MINIMUM     = 2,   /* See TstatType in TYPES.H */
@@ -233,8 +236,6 @@ typedef enum {
   EN_RANGE       = 4
 } EN_StatisticType;
 
-
-
 typedef enum {
   EN_MIX1        = 0,   /* Tank mixing models */
   EN_MIX2        = 1,
@@ -242,8 +243,6 @@ typedef enum {
   EN_LIFO        = 3
 } EN_MixingModel;
 
-
-
 typedef enum {
   EN_NOSAVE        = 0,
   EN_SAVE          = 1,
@@ -251,8 +250,6 @@ typedef enum {
   EN_SAVE_AND_INIT = 11
 } EN_SaveOption;
 
-
-
 typedef enum {
   EN_CONST_HP    = 0,   /* constant horsepower       */
   EN_POWER_FUNC  = 1,   /* power function            */
@@ -260,7 +257,6 @@ typedef enum {
 } EN_CurveType;
 
 
-
 // --- Declare the EPANET toolkit functions
 #if defined(__cplusplus)
 extern "C" {
@@ -270,8 +266,6 @@ extern "C" {
    @brief The EPANET Project wrapper object
    */
   typedef struct EN_Project EN_Project;
-  typedef struct EN_Pattern EN_Pattern;
-  typedef struct EN_Curve EN_Curve;
   
   /**
    @brief runs a complete EPANET simulation
@@ -522,7 +516,29 @@ extern "C" {
    @return Error code
    */
   int  DLLEXPORT ENsetflowunits(int code);
-  
+
+  /**
+   @brief Retrieves the type of demand model in use and its parameters
+   @param[out] type  Type of demand model (EN_DDA or EN_PDA)
+   @param[out] pmin  Pressure below which there is no demand
+   @param[out] preq  Pressure required to deliver full demand
+   @param[out] pexp  Pressure exponent in demand function
+   @return Error code
+  */
+  int DLLEXPORT ENgetdemandmodel(int *type, EN_API_FLOAT_TYPE *pmin,
+      EN_API_FLOAT_TYPE *preq, EN_API_FLOAT_TYPE *pexp);
+
+  /**
+  @brief Sets the type of demand model to use and its parameters
+  @param type  Type of demand model (EN_DDA or EN_PDA)
+  @param pmin  Pressure below which there is no demand
+  @param preq  Pressure required to deliver full demand
+  @param pexp  Pressure exponent in demand function
+  @return Error code
+  */
+  int DLLEXPORT ENsetdemandmodel(int type, EN_API_FLOAT_TYPE pmin,
+      EN_API_FLOAT_TYPE preq, EN_API_FLOAT_TYPE pexp);
+
   /**
    @brief Retrieves the index of the time pattern with specified ID
    @param id String ID of the time pattern
@@ -795,7 +811,7 @@ extern "C" {
   int  DLLEXPORT ENsetnodevalue(int index, int code, EN_API_FLOAT_TYPE v);
   
   /**
-   @brief Set a proprty value for a link.
+   @brief Set a property value for a link.
    @param index The index of a link. First link is index 1.
    @param code The code for the property to set.
    @param v The value to set for this link and property.
@@ -1120,15 +1136,14 @@ extern "C" {
   int DLLEXPORT ENdeletelink(int linkIndex);
   
   
-  
-  
   /***************************************************
    
    Threadsafe versions of all epanet functions
    
    ***************************************************/
-  int DLLEXPORT EN_alloc(EN_Project **p);
+  int  DLLEXPORT EN_createproject(EN_Project **p);
-  int DLLEXPORT EN_free(EN_Project *p);
+  int  DLLEXPORT EN_deleteproject(EN_Project *p);
+
   int DLLEXPORT EN_epanet(char *inpFile, char *rptFile, char *binOutFile, void (*callback) (char *));
   int DLLEXPORT EN_init(EN_Project *p, char *rptFile, char *binOutFile, EN_FlowUnits UnitsType, EN_FormType HeadlossFormula);
   int DLLEXPORT EN_open(EN_Project *p, char *inpFile, char *rptFile, char *binOutFile);
@@ -1160,6 +1175,10 @@ extern "C" {
   int DLLEXPORT EN_gettimeparam(EN_Project *p, int code, long *value);
   int DLLEXPORT EN_getflowunits(EN_Project *p, int *code);
   int DLLEXPORT EN_setflowunits(EN_Project *p, int code);
+  int DLLEXPORT EN_getdemandmodel(EN_Project *p, int *type, EN_API_FLOAT_TYPE *pmin, EN_API_FLOAT_TYPE *preq,
+                                  EN_API_FLOAT_TYPE *pexp);
+  int DLLEXPORT EN_setdemandmodel(EN_Project *p, int type, EN_API_FLOAT_TYPE pmin, EN_API_FLOAT_TYPE preq,
+                                  EN_API_FLOAT_TYPE pexp);
   int DLLEXPORT EN_getpatternindex(EN_Project *p, char *id, int *index);
   int DLLEXPORT EN_getpatternid(EN_Project *p, int index, char *id);
   int DLLEXPORT EN_getpatternlen(EN_Project *p, int index, int *len);
@@ -1224,11 +1243,6 @@ extern "C" {
   int DLLEXPORT EN_deletenode(EN_Project *p, int nodeIndex);
   int DLLEXPORT EN_deletelink(EN_Project *p, int linkIndex);
   
-  
-  
-  
-  
-  
 #if defined(__cplusplus)
 }
 #endif

src/enumstxt.h

@@ -74,6 +74,10 @@ char *PressUnitsTxt[]   = {w_PSI,
                            w_KPA,
                            w_METERS};
 
+char *DemandModelTxt[] = { w_DDA,
+                           w_PDA,
+                           NULL };
+
 char *QualTxt[]         = {w_NONE,
                            w_CHEM,
                            w_AGE,

src/epanet.c

@@ -129,9 +129,10 @@ execute function x and set the error code equal to its return value.
 #include "funcs.h"
 #include "text.h"
 #include "types.h"
-#define EXTERN
+
-////////////////////////////////////////////#include "epanet2.h"
+// This single global variable is used only when the library is called
-#include "vars.h"
+// in "legacy mode" with the 2.1-style API. 
+EN_Project *_defaultModel;
 
 /****************************************************************
 
@@ -159,7 +160,7 @@ execute function x and set the error code equal to its return value.
 int DLLEXPORT ENepanet(char *f1, char *f2, char *f3,
                        void (*pviewprog)(char *)) {
   int errcode = 0;
-  ERRCODE(EN_alloc(&_defaultModel));
+  ERRCODE(EN_createproject(&_defaultModel));
   ERRCODE(EN_open(_defaultModel, f1, f2, f3));
   _defaultModel->viewprog = pviewprog;
   if (_defaultModel->out_files.Hydflag != USE) {
@@ -168,13 +169,13 @@ int DLLEXPORT ENepanet(char *f1, char *f2, char *f3,
   ERRCODE(EN_solveQ(_defaultModel));
   ERRCODE(EN_report(_defaultModel));
   EN_close(_defaultModel);
-  EN_free(_defaultModel);
+  EN_deleteproject(_defaultModel);
   return (errcode);
 }
 
 int DLLEXPORT ENopen(char *f1, char *f2, char *f3) {
   int errcode = 0;
-  ERRCODE(EN_alloc(&_defaultModel));
+  ERRCODE(EN_createproject(&_defaultModel));
   EN_open(_defaultModel, f1, f2, f3);
   return (errcode);
 }
@@ -183,7 +184,10 @@ int DLLEXPORT ENsaveinpfile(char *filename) {
   return EN_saveinpfile(_defaultModel, filename);
 }
 
-int DLLEXPORT ENclose() { return EN_close(_defaultModel); }
+int DLLEXPORT ENclose() {
+    EN_close(_defaultModel);
+    return EN_deleteproject(_defaultModel);
+}
 
 int DLLEXPORT ENsolveH() { return EN_solveH(_defaultModel); }
 
@@ -262,6 +266,16 @@ int DLLEXPORT ENsetflowunits(int code) {
   return EN_setflowunits(_defaultModel, code);
 }
 
+int DLLEXPORT ENgetdemandmodel(int *type, EN_API_FLOAT_TYPE *pmin, EN_API_FLOAT_TYPE *preq,
+                               EN_API_FLOAT_TYPE *pexp) {
+    return EN_getdemandmodel(_defaultModel, type, pmin, preq, pexp);
+}
+
+int DLLEXPORT ENsetdemandmodel(int type, EN_API_FLOAT_TYPE pmin, EN_API_FLOAT_TYPE preq,
+                               EN_API_FLOAT_TYPE pexp) {
+    return EN_setdemandmodel(_defaultModel, type, pmin, preq, pexp);
+}
+
 int DLLEXPORT ENgetpatternindex(char *id, int *index) {
   return EN_getpatternindex(_defaultModel, id, index);
 }
@@ -532,23 +546,32 @@ int DLLEXPORT ENdeletenode(int index) {
   return EN_deletenode(_defaultModel, index);
 }
 
+
+int DLLEXPORT EN_epanet(char *inpFile, char *rptFile, char *binOutFile, void(*callback) (char *))
+{
+    return ENepanet(inpFile, rptFile, binOutFile, callback);
+}
+
 /*
 ----------------------------------------------------------------
    Functions for opening & closing the EPANET system
 ----------------------------------------------------------------
 */
 
-/// allocate a project pointer
+/// Create an EPANET project
-int DLLEXPORT EN_alloc(EN_Project **p)
+int DLLEXPORT EN_createproject(EN_Project **p)
 {
-  EN_Project *project = calloc(1, sizeof(EN_Project));
+    EN_Project *project = calloc(1, sizeof(EN_Project));
-  *p = project;
+    if (project == NULL) return 101;
-  return 0;
+    *p = project;
+    return 0;
 }
 
-int DLLEXPORT EN_free(EN_Project *p) 
+/// Delete an EPANET project
+int DLLEXPORT EN_deleteproject(EN_Project *p) 
 {
-  free(p);
+  if (p) free(p);
+  p = NULL;
   return 0;
 }
 
@@ -558,9 +581,8 @@ int DLLEXPORT EN_init(EN_Project *pr, char *f2, char *f3,
  **  Input:
  **           f2               = pointer to name of report file
  **           f3               = pointer to name of binary output file
- UnitsType        = flow units flag
+ **           UnitsType        = flow units flag
- HeadlossFormula  = headloss formula flag
+ **           HeadlossFormula  = headloss formula flag
-
  **  Output:  none
  **  Returns: error code
  **  Purpose: opens EPANET
@@ -1504,6 +1526,30 @@ int DLLEXPORT EN_setflowunits(EN_Project *p, int code) {
   return(0);
 }
 
+
+int DLLEXPORT EN_getdemandmodel(EN_Project *p, int *type, EN_API_FLOAT_TYPE *pmin,
+              EN_API_FLOAT_TYPE *preq, EN_API_FLOAT_TYPE *pexp)
+{
+    *type = p->hydraulics.DemandModel;
+    *pmin = (EN_API_FLOAT_TYPE)(p->hydraulics.Pmin * p->Ucf[PRESSURE]);
+    *preq = (EN_API_FLOAT_TYPE)(p->hydraulics.Preq * p->Ucf[PRESSURE]);
+    *pexp = (EN_API_FLOAT_TYPE)(p->hydraulics.Pexp);
+    return 0;
+}
+
+int DLLEXPORT EN_setdemandmodel(EN_Project *p, int type, EN_API_FLOAT_TYPE pmin,
+              EN_API_FLOAT_TYPE preq, EN_API_FLOAT_TYPE pexp)
+{
+    if (type < 0 || type > EN_PDA) return 251;
+    if (pmin > preq || pexp <= 0.0) return 202;
+    p->hydraulics.DemandModel = type;
+    p->hydraulics.Pmin = pmin / p->Ucf[PRESSURE];
+    p->hydraulics.Preq = preq / p->Ucf[PRESSURE];
+    p->hydraulics.Pexp = pexp;
+    return 0;
+}
+
+
 int DLLEXPORT EN_getpatternindex(EN_Project *p, char *id, int *index) {
   int i;
   *index = 0;
@@ -1657,11 +1703,17 @@ int DLLEXPORT EN_geterror(int errcode, char *errmsg, int n) {
 int DLLEXPORT EN_getstatistic(EN_Project *p, int code, EN_API_FLOAT_TYPE *value) {
   switch (code) {
   case EN_ITERATIONS:
-    *value = (EN_API_FLOAT_TYPE)p->hydraulics.iterations;
+    *value = (EN_API_FLOAT_TYPE)p->hydraulics.Iterations;
     break;
   case EN_RELATIVEERROR:
-    *value = (EN_API_FLOAT_TYPE)p->hydraulics.relativeError;
+    *value = (EN_API_FLOAT_TYPE)p->hydraulics.RelativeError;
     break;
+  case EN_MAXHEADERROR:
+      *value = (EN_API_FLOAT_TYPE)(p->hydraulics.MaxHeadError * p->Ucf[HEAD]);
+      break;
+  case EN_MAXFLOWCHANGE:
+      *value = (EN_API_FLOAT_TYPE)(p->hydraulics.MaxFlowChange * p->Ucf[FLOW]);
+      break;
   default:
     break;
   }
@@ -1725,8 +1777,8 @@ int DLLEXPORT EN_getcoord(EN_Project *p, int index, EN_API_FLOAT_TYPE *x,
   if (p->network.Coord[index].HaveCoords == FALSE)
     return (254);
 
-  *x = p->network.Coord[index].X;
+  *x = (EN_API_FLOAT_TYPE)p->network.Coord[index].X;
-  *y = p->network.Coord[index].Y;
+  *y = (EN_API_FLOAT_TYPE)p->network.Coord[index].Y;
   return 0;
 }
 
@@ -4215,7 +4267,7 @@ int DLLEXPORT EN_getaveragepatternvalue(EN_Project *p, int index, EN_API_FLOAT_T
     return (205);
   // if (period < 1 || period > Pattern[index].Length) return(251);
   for (i = 0; i < Pattern[index].Length; i++) {
-    *value += Pattern[index].F[i];
+    *value += (EN_API_FLOAT_TYPE)Pattern[index].F[i];
   }
   *value /= (EN_API_FLOAT_TYPE)Pattern[index].Length;
   return (0);
@@ -4686,7 +4738,7 @@ int DLLEXPORT EN_getpremise(EN_Project *pr, int indexRule, int idxPremise, int *
   *variable = p->variable;
   *relop = p->relop;
   *status = p->status;
-  *value = p[0].value;
+  *value = (EN_API_FLOAT_TYPE)p[0].value;
   return (0);
 }
 
@@ -4820,7 +4872,7 @@ int DLLEXPORT EN_gettrueaction(EN_Project *pr, int indexRule, int indexAction, i
   }
   *indexLink = a->link;
   *status = a->status;
-  *setting = a->setting;
+  *setting = (EN_API_FLOAT_TYPE)a->setting;
   return (0);
 }
 
@@ -4868,7 +4920,7 @@ int DLLEXPORT EN_getfalseaction(EN_Project *pr, int indexRule, int indexAction,
   }
   *indexLink = a->link;
   *status = a->status;
-  *setting = a->setting;
+  *setting = (EN_API_FLOAT_TYPE)a->setting;
   return (0);
 }
 

src/funcs.h

@@ -158,28 +158,13 @@ void    inithyd(EN_Project *pr, int initFlags);     /* Re-sets initial condition
 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   */
-int     allocmatrix(EN_Project *pr);                /* Allocates matrix coeffs.   */
-void    freematrix(EN_Project *pr);                 /* Frees matrix coeffs.       */
-void    initlinkflow(EN_Project *pr, int, char,
-                     double);                       /* Initializes link flow      */
-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    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. */
 
@@ -188,8 +173,13 @@ int     hydsolve(EN_Project *pr, int *,double *);   /* Solves network equations
 
 /* ----------- HYDCOEFFS.C --------------*/
 void    resistcoeff(EN_Project *pr, int k);         /* Finds pipe flow resistance */
-void    hlosscoeff(EN_Project *pr, int k);          /* Finds link head loss coeff */
+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  */
+double  demandflowchange(EN_Project *pr, int i,     // Change in demand outflow
+                         double dp, double n);
+void    demandparams(EN_Project *pr, double *dp,    // PDA function parameters
+                     double *n); 
 
 /* ----------- SMATRIX.C ---------------*/
 int     createsparse(EN_Project *pr);               /* Creates sparse matrix      */

src/hydcoeffs.c

@@ -29,15 +29,21 @@ 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   resistcoeff(EN_Project *pr, int k);
-//void  hlosscoeff(EN_Project *pr, int k);
+//void   headlosscoeffs(EN_Project *pr);
-//void  matrixcoeffs(EN_Project *pr);
+//void   matrixcoeffs(EN_Project *pr);
+//double emitflowchange(EN_Project *pr, int i);
+//double demandflowchange(EN_Project *pr, int i, double dp, double n);
+//void   demandparams(EN_Project *pr, double *dp, double *n);
 
 // 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    demandcoeffs(EN_Project *pr);
+static void    demandheadloss(double d, double dfull, double dp,
+               double n, double *hloss, double *hgrad);
 
 static void    pipecoeff(EN_Project *pr, int k);
 static void    DWpipecoeff(EN_Project *pr, int k);
@@ -55,7 +61,6 @@ 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)
 /*
 **--------------------------------------------------------------------
@@ -83,13 +88,14 @@ void  resistcoeff(EN_Project *pr, int k)
         switch (hyd->Formflag)
         {
         case HW:
-            link->R = 4.727*L / pow(e, hyd->Hexp) / pow(d, 4.871);
+            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);
+            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;
+            link->R = SQR(4.0 * e / (1.49 * PI * SQR(d))) *
+                      pow((d / 4.0), -1.333) * L;
         }
         break;
 
@@ -106,45 +112,46 @@ void  resistcoeff(EN_Project *pr, int k)
 }
 
 
-void  hlosscoeff(EN_Project *pr, int k)
+void headlosscoeffs(EN_Project *pr)
 /*
 **--------------------------------------------------------------
-**   Input:   k = link index
+**   Input:   none
 **   Output:  none
-**   Purpose: computes P and Y coefficients for a link
+**   Purpose: computes coefficients P (1 / head loss gradient)
+**            and Y (head loss / gradient) for all links.
 **--------------------------------------------------------------
 */
 {
+    int k;
     EN_Network   *net = &pr->network;
     hydraulics_t *hyd = &pr->hydraulics;
-    solver_t     *sol = &hyd->solver;
-    Slink        *link = &net->Link[k];
 
-    switch (link->Type)
+    for (k = 1; k <= net->Nlinks; k++)
     {
-    case EN_CVPIPE:
+        switch (net->Link[k].Type)
-    case EN_PIPE:
+        {
-        pipecoeff(pr, k);
+        case EN_CVPIPE:
-        break;
+        case EN_PIPE:
-    case EN_PUMP:
+            pipecoeff(pr, k);
-        pumpcoeff(pr, k);
+            break;
-        break;
+        case EN_PUMP:
-    case EN_PBV:
+            pumpcoeff(pr, k);
-        pbvcoeff(pr, k);
+            break;
-        break;
+        case EN_PBV:
-    case EN_TCV:
+            pbvcoeff(pr, k);
-        tcvcoeff(pr, k);
+            break;
-        break;
+        case EN_TCV:
-    case EN_GPV:
+            tcvcoeff(pr, k);
-        gpvcoeff(pr, k);
+            break;
-        break;
+        case EN_GPV:
-    case EN_FCV:
+            gpvcoeff(pr, k);
-    case EN_PRV:
+            break;
-    case EN_PSV:
+        case EN_FCV:
-        if (hyd->LinkSetting[k] == MISSING) {
+        case EN_PRV:
-            valvecoeff(pr, k);
+        case EN_PSV:
+            if (hyd->LinkSetting[k] == MISSING) valvecoeff(pr, k);
+            else hyd->solver.P[k] = 0.0;
         }
-        else sol->P[k] = 0.0;
     }
 }
 
@@ -162,13 +169,23 @@ void   matrixcoeffs(EN_Project *pr)
     solver_t     *sol = &hyd->solver;
     EN_Network   *net = &pr->network;
 
+    // Reset values of all diagonal coeffs. (Aii), off-diagonal
+    // coeffs. (Aij), r.h.s. coeffs. (F) and node flow balance (X_tmp)
     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));
+
+    // Compute matrix coeffs. from links, emitters, and nodal demands
     linkcoeffs(pr);
     emittercoeffs(pr);
+    demandcoeffs(pr);
+
+    // Update nodal flow balances with demands and add onto r.h.s. coeffs.
     nodecoeffs(pr);
+
+    // Finally, find coeffs. for PRV/PSV/FCV control valves whose
+    // status is not fixed to OPEN/CLOSED
     valvecoeffs(pr);
 }
 
@@ -189,7 +206,7 @@ void  linkcoeffs(EN_Project *pr)
     solver_t     *sol = &hyd->solver;
     Slink *link;
 
-    // Examine each link of network */
+    // Examine each link of network
     for (k = 1; k <= net->Nlinks; k++)
     {
         if (sol->P[k] == 0.0) continue;
@@ -197,7 +214,7 @@ void  linkcoeffs(EN_Project *pr)
         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)
+        // Update net nodal inflows (X_tmp), 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];
@@ -212,21 +229,18 @@ void  linkcoeffs(EN_Project *pr)
             sol->F[sol->Row[n1]] += sol->Y[k];     // RHS coeff.
         }
 
-        // Node n1 is a tank
+        // Node n1 is a tank/reservoir
-        else {
+        else sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); 
-            sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]); 
-        }
 
         // Node n2 is junction
-        if (n2 <= net->Njuncs) {
+        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
+        // Node n2 is a tank/reservoir
-        else {
+        else sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); 
-            sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]); 
-        }
     }
 }
 
@@ -236,7 +250,7 @@ void  nodecoeffs(EN_Project *pr)
 **----------------------------------------------------------------
 **  Input:   none
 **  Output:  none
-**  Purpose: completes calculation of nodal flow imbalance (X)
+**  Purpose: completes calculation of nodal flow imbalance (X_tmp)
 **           & flow correction (F) arrays
 **----------------------------------------------------------------
 */
@@ -250,7 +264,7 @@ void  nodecoeffs(EN_Project *pr)
     // flow imbalance & add imbalance to RHS array F.
     for (i = 1; i <= net->Njuncs; i++)
     {
-        hyd->X_tmp[i] -= hyd->NodeDemand[i];
+        hyd->X_tmp[i] -= hyd->DemandFlows[i];
         sol->F[sol->Row[i]] += hyd->X_tmp[i];
     }
 }
@@ -276,15 +290,20 @@ void  valvecoeffs(EN_Project *pr)
     // Examine each valve
     for (i = 1; i <= net->Nvalves; i++)
     {
+        // Find valve's link index
         valve = &net->Valve[i];
-        k = valve->Link;                         // Link index of valve
+        k = valve->Link;
+
+        // Coeffs. for fixed status valves have already been computed
+        if (hyd->LinkSetting[k] == MISSING) continue;
+
+        // Start & end nodes of valve's link        
         link = &net->Link[k];
-        if (hyd->LinkSetting[k] == MISSING) {
+        n1 = link->N1;
-            continue;                            // Valve status fixed
-        }
-        n1 = link->N1;                           // Start & end nodes
         n2 = link->N2;
-        switch (link->Type)                      // Call valve-specific function 
+
+        // Call valve-specific function
+        switch (link->Type) 
         {
         case EN_PRV:
             prvcoeff(pr, k, n1, n2);
@@ -330,17 +349,17 @@ void  emittercoeffs(EN_Project *pr)
     for (i = 1; i <= net->Njuncs; i++)
     {
         node = &net->Node[i];
-        if (node->Ke == 0.0) {
+        if (node->Ke == 0.0) continue;
-            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) {
+        if (p < hyd->RQtol)
+        {
             p = 1.0 / hyd->RQtol;
         }
-        else {
+        else
+        {
             p = 1.0 / p;                         // inverse head loss gradient
         }
         y = SGN(q)*z*p;                          // head loss / gradient
@@ -351,6 +370,173 @@ void  emittercoeffs(EN_Project *pr)
 }
 
 
+double emitflowchange(EN_Project *pr, int i)
+/*
+**--------------------------------------------------------------
+**   Input:   i = node index
+**   Output:  returns change in flow at an emitter node
+**   Purpose: computes flow change at an emitter node
+**--------------------------------------------------------------
+*/
+{
+    double ke, p;
+    hydraulics_t *hyd = &pr->hydraulics;
+    Snode *node = &pr->network.Node[i];
+
+    ke = MAX(CSMALL, node->Ke);
+    p = hyd->Qexp * ke * pow(ABS(hyd->EmitterFlows[i]), (hyd->Qexp - 1.0));
+    if (p < hyd->RQtol) p = 1 / hyd->RQtol;
+    else                p = 1.0 / p;
+    return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El));
+}
+
+
+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
+**            demand function.
+**--------------------------------------------------------------
+*/
+{
+    hydraulics_t *hyd = &pr->hydraulics;
+
+    // If required pressure equals minimum pressure, use a linear demand
+    // curve with a 0.01 PSI pressure range to approximate an all or
+    // nothing demand solution
+    if (hyd->Preq == hyd->Pmin)
+    {
+        *dp = 0.01 / PSIperFT;
+        *n = 1.0;
+    }
+
+    // Otherwise use the user-supplied demand curve parameters
+    else
+    {
+        *dp = hyd->Preq - hyd->Pmin;
+        *n = 1.0 / hyd->Pexp;
+    }
+}
+
+
+void  demandcoeffs(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+**   Input:   none
+**   Output:  none
+**   Purpose: computes matrix coeffs. for pressure dependent demands
+**
+**   Note: Pressure dependent demands are modelled like emitters
+**         with Hloss = Pserv * (D / Dfull)^(1/Pexp)
+**         where D (actual demand) is zero for negative pressure
+**         and is Dfull above pressure Pserv.
+**--------------------------------------------------------------
+*/
+{
+    int i, row;
+    double  dp,         // pressure range over which demand can vary (ft)
+            n,          // exponent in head loss v. demand function
+            hloss,      // head loss in supplying demand (ft)
+            hgrad;      // gradient of demand head loss (ft/cfs)    
+
+    hydraulics_t *hyd = &pr->hydraulics;
+    solver_t     *sol = &hyd->solver;
+    EN_Network   *net = &pr->network;
+
+    // Get demand function parameters
+    if (hyd->DemandModel == DDA) return;
+    demandparams(pr, &dp, &n);
+
+    // Examine each junction node
+    for (i = 1; i <= net->Njuncs; i++)
+    {
+        // Skip junctions with non-positive demands
+        if (hyd->NodeDemand[i] <= 0.0) continue;
+
+        // Find head loss for demand outflow at node's elevation
+        demandheadloss(hyd->DemandFlows[i], hyd->NodeDemand[i], dp, n,
+                    &hloss, &hgrad);
+
+        // Update row of solution matrix A & its r.h.s. F
+        row = sol->Row[i];
+        sol->Aii[row] += 1.0 / hgrad;
+        sol->F[row] += (hloss + net->Node[i].El + hyd->Pmin) / hgrad;
+    }
+}
+
+
+double demandflowchange(EN_Project *pr, int i, double dp, double n)
+/*
+**--------------------------------------------------------------
+**   Input:   i  = node index
+**            dp = pressure range fro demand funtion (ft)
+**            n  = exponent in head v. demand function
+**   Output:  returns change in pressure dependent demand flow
+**   Purpose: computes change in outflow at at a node subject to
+**            pressure dependent demands
+**--------------------------------------------------------------
+*/
+{
+    double hloss, hgrad;
+    hydraulics_t *hyd = &pr->hydraulics;
+
+    demandheadloss(hyd->DemandFlows[i], hyd->NodeDemand[i], dp, n, &hloss, &hgrad);
+    return (hloss - hyd->NodeHead[i] + pr->network.Node[i].El + hyd->Pmin) / hgrad;
+}
+
+
+void demandheadloss(double d, double dfull, double dp, double n,
+    double *hloss, double *hgrad)
+    /*
+    **--------------------------------------------------------------
+    **   Input:   d     = actual junction demand (cfs)
+    **            dfull = full junction demand required (cfs)
+    **            dp    = pressure range for demand function (ft)
+    **            n     = exponent in head v. demand function
+    **   Output:  hloss = head loss delivering demand d (ft)
+    **            hgrad = gradient of head loss (ft/cfs)
+    **  Purpose:  computes head loss and its gradient for delivering
+    **            a pressure dependent demand flow.
+    **--------------------------------------------------------------
+    */
+{
+    const double RB = 1.0e9;
+    const double EPS = 0.001;
+    double r = d / dfull;
+
+    // Use upper barrier function for excess demand above full value
+    if (r > 1.0)
+    {
+        *hgrad = RB;
+        *hloss = dp + RB * (d - dfull);
+    }
+
+    // Use lower barrier function for negative demand
+    else if (r < 0)
+    {
+        *hgrad = RB;
+        *hloss = RB * d;
+    }
+
+    // Use linear head loss function for near zero demand
+    else if (r < EPS)
+    {
+        *hgrad = dp * pow(EPS, n - 1.0) / dfull;
+        *hloss = (*hgrad) * d;
+    }
+
+    // Otherwise use power head loss function
+    else
+    {
+        *hgrad = n * dp * pow(r, n - 1.0) / dfull;
+        *hloss = (*hgrad) * d / n;
+    }
+}
+
+
 void  pipecoeff(EN_Project *pr, int k)
 /*
 **--------------------------------------------------------------
@@ -382,8 +568,9 @@ void  pipecoeff(EN_Project *pr, int k)
         return;
     }
 
-    // ... head loss formula is Darcy-Weisbach
+    // Use custom function for Darcy-Weisbach formula
-    if (hyd->Formflag == DW) {
+    if (hyd->Formflag == DW)
+    {
         DWpipecoeff(pr, k);
         return;
     }
@@ -550,7 +737,8 @@ void  pumpcoeff(EN_Project *pr, int k)
     Spump *pump;
 
     // Use high resistance pipe if pump closed or cannot deliver head
-    if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0) {
+    if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0)
+    {
         sol->P[k] = 1.0 / CBIG;
         sol->Y[k] = hyd->LinkFlows[k];
         return;
@@ -580,9 +768,7 @@ void  pumpcoeff(EN_Project *pr, int k)
     h0 = SQR(setting) * pump->H0;
     n = pump->N;
     r = pump->R * pow(setting, 2.0 - n);
-    if (n != 1.0) {
+    if (n != 1.0) r = n * r * pow(q, 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);
@@ -615,15 +801,9 @@ void  curvecoeff(EN_Project *pr, int i, double q, double *h0, double *r)
 
     // Find linear segment of curve that brackets flow q
     k2 = 0;
-    while (k2 < npts && x[k2] < q)
+    while (k2 < npts && x[k2] < q) k2++;
-        k2++;
+    if (k2 == 0) k2++;
-
+    else if (k2 == npts)  k2--;
-    if (k2 == 0)
-        k2++;
-
-    else if (k2 == npts)
-        k2--;
-
     k1 = k2 - 1;
 
     // Compute slope and intercept of this segment
@@ -653,13 +833,12 @@ void  gpvcoeff(EN_Project *pr, int k)
     solver_t     *sol = &hyd->solver;
 
     // Treat as a pipe if valve closed
-    if (hyd->LinkStatus[k] == CLOSED) {
+    if (hyd->LinkStatus[k] == CLOSED) valvecoeff(pr, k);                          
-        valvecoeff(pr, k);                          
-    }
 
     // Otherwise utilize headloss curve
     // whose index is stored in K
-    else {
+    else
+    {
         // Find slope & intercept of headloss curve.
         q = ABS(hyd->LinkFlows[k]);
         q = MAX(q, TINY);
@@ -687,18 +866,22 @@ void  pbvcoeff(EN_Project *pr, int k)
     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) {
+    if (hyd->LinkSetting[k] == MISSING || hyd->LinkSetting[k] == 0.0)
+    {
         valvecoeff(pr, k);         
     }
 
     // If valve is active
-    else {
+    else
+    {
         // Treat as a pipe if minor loss > valve setting
-        if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k]) {
+        if (link->Km * SQR(hyd->LinkFlows[k]) > hyd->LinkSetting[k])
+        {
             valvecoeff(pr, k);        
         }
         // Otherwise force headloss across valve to be equal to setting
-        else {
+        else
+        {
             sol->P[k] = CBIG;
             sol->Y[k] = hyd->LinkSetting[k] * CBIG;
         }
@@ -723,7 +906,8 @@ void  tcvcoeff(EN_Project *pr, int k)
     km = link->Km;
 
     // If valve not fixed OPEN or CLOSED, compute its loss coeff.
-    if (hyd->LinkSetting[k] != MISSING) {
+    if (hyd->LinkSetting[k] != MISSING)
+    {
         link->Km = 0.02517 * hyd->LinkSetting[k] / (SQR(link->Diam)*SQR(link->Diam));
     }
 
@@ -768,7 +952,8 @@ void  prvcoeff(EN_Project *pr, int k, int n1, int n2)
         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) {
+        if (hyd->X_tmp[n2] < 0.0)
+        {
             sol->F[i] += hyd->X_tmp[n2];
         }
         return;
@@ -818,7 +1003,8 @@ void  psvcoeff(EN_Project *pr, int k, int n1, int n2)
         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) {
+        if (hyd->X_tmp[n1] > 0.0)
+        {
             sol->F[j] += hyd->X_tmp[n1];
         }
         return;
@@ -919,9 +1105,7 @@ void valvecoeff(EN_Project *pr, int k)
     if (link->Km > 0.0)
     {
         p = 2.0 * link->Km * fabs(flow);
-        if (p < hyd->RQtol) {
+        if (p < hyd->RQtol) p = hyd->RQtol;
-            p = hyd->RQtol;
-        }
         sol->P[k] = 1.0 / p;
         sol->Y[k] = flow / 2.0;
     }

src/hydraul.c

@@ -60,6 +60,20 @@ AUTHOR:     L. Rossman
 
 #define   QZERO  1.e-6  /* Equivalent to zero flow */
 
+// Local functions
+int     allocmatrix(EN_Project *pr);
+void    freematrix(EN_Project *pr);
+void    initlinkflow(EN_Project *pr, int, char, double);
+void    setlinkflow(EN_Project *pr, int, double);
+void    demands(EN_Project *pr);
+int     controls(EN_Project *pr);
+long    timestep(EN_Project *pr);
+void    controltimestep(EN_Project *pr, long *);
+void    ruletimestep(EN_Project *pr, long *);
+void    addenergy(EN_Project *pr, long);
+void    tanklevels(EN_Project *pr, long);
+
+
 int  openhyd(EN_Project *pr)
 /*
  *--------------------------------------------------------------
@@ -194,15 +208,12 @@ int   runhyd(EN_Project *pr, long *t)
 
     /* Solve network hydraulic equations */
     errcode = hydsolve(pr,&iter,&relerr);
+
     if (!errcode) {
         /* Report new status & save results */
         if (rep->Statflag) {
             writehydstat(pr,iter,relerr);
         }
-
-        /* solution info */
-        hyd->relativeError = relerr;
-        hyd->iterations = iter;
      
 /*** Updated 3/1/01 ***/
         /* If system unbalanced and no extra trials */
@@ -316,6 +327,7 @@ int  allocmatrix(EN_Project *pr)
    s->Aii = (double *) calloc(net->Nnodes+1,sizeof(double));
    s->Aij = (double *) calloc(hyd->Ncoeffs+1,sizeof(double));
    s->F   = (double *) calloc(net->Nnodes+1,sizeof(double));
+   hyd->DemandFlows = (double *)calloc(net->Nnodes + 1, sizeof(double));
    hyd->EmitterFlows   = (double *) calloc(net->Nnodes+1,sizeof(double));
    s->P   = (double *) calloc(net->Nlinks+1,sizeof(double));
    s->Y   = (double *) calloc(net->Nlinks+1,sizeof(double));
@@ -324,6 +336,7 @@ int  allocmatrix(EN_Project *pr)
    ERRCODE(MEMCHECK(s->Aii));
    ERRCODE(MEMCHECK(s->Aij));
    ERRCODE(MEMCHECK(s->F));
+   ERRCODE(MEMCHECK(hyd->DemandFlows));
    ERRCODE(MEMCHECK(hyd->EmitterFlows));
    ERRCODE(MEMCHECK(s->P));
    ERRCODE(MEMCHECK(s->Y));
@@ -348,6 +361,7 @@ void  freematrix(EN_Project *pr)
    free(s->Aii);
    free(s->Aij);
    free(s->F);
+   free(hyd->DemandFlows);
    free(hyd->EmitterFlows);
    free(s->P);
    free(s->Y);
@@ -551,7 +565,6 @@ void  setlinksetting(EN_Project *pr, int index, double value, StatType *s, doubl
 } 
 
 
-
 void  demands(EN_Project *pr)
 /*
 **--------------------------------------------------------------------
@@ -591,6 +604,9 @@ void  demands(EN_Project *pr)
          sum += djunc;
       }
       hyd->NodeDemand[i] = sum;
+
+      // Initialize pressure dependent demand
+      hyd->DemandFlows[i] = sum;
    }
 
    /* Update head at fixed grade nodes with time patterns. */

src/hydstatus.c

@@ -0,0 +1,462 @@
+/*
+*********************************************************************
+
+HYDSTATUS.C --  Hydraulic status updating for the EPANET Program
+
+*******************************************************************
+*/
+
+#include <stdio.h>
+#include "types.h"
+#include "funcs.h"
+
+// External functions
+int  valvestatus(EN_Project *pr);
+int  linkstatus(EN_Project *pr);
+
+// Local functions
+static StatType cvstatus(EN_Project *pr, StatType, double, double);
+static StatType pumpstatus(EN_Project *pr, int, double);
+static StatType prvstatus(EN_Project *pr, int, StatType, double, double, double);
+static StatType psvstatus(EN_Project *pr, int, StatType, double, double, double);
+static StatType fcvstatus(EN_Project *pr, int, StatType, double, double);
+static void     tankstatus(EN_Project *pr, int, int, int);
+
+
+int  valvestatus(EN_Project *pr)
+/*
+**-----------------------------------------------------------------
+**  Input:   none
+**  Output:  returns 1 if any pressure or flow control valve
+**           changes status, 0 otherwise
+**  Purpose: updates status for PRVs & PSVs whose status
+**           is not fixed to OPEN/CLOSED
+**-----------------------------------------------------------------
+*/
+{
+    int   change = FALSE,            // Status change flag
+          i, k,                      // Valve & link indexes
+          n1, n2;                    // Start & end nodes
+    StatType  status;                // Valve status settings
+    double hset;                     // Valve head setting
+    Slink *link;
+
+    EN_Network       *net = &pr->network;
+    hydraulics_t     *hyd = &pr->hydraulics;
+    report_options_t *rep = &pr->report;
+
+    // Examine each valve
+    for (i = 1; i <= net->Nvalves; i++)
+    {
+        // Get valve's link and its index
+        k = net->Valve[i].Link;
+        link = &net->Link[k];
+
+        // Ignore valve if its status is fixed to OPEN/CLOSED
+        if (hyd->LinkSetting[k] == MISSING) continue;
+        
+        // Get start/end node indexes & save current status
+        n1 = link->N1;
+        n2 = link->N2;
+        status = hyd->LinkStatus[k];
+
+        // Evaluate valve's new status
+        switch (link->Type)
+        {
+        case EN_PRV:
+            hset = net->Node[n2].El + hyd->LinkSetting[k];
+            hyd->LinkStatus[k] = prvstatus(pr, k, status, hset,
+                                     hyd->NodeHead[n1], hyd->NodeHead[n2]);
+            break;
+        case EN_PSV:
+            hset = net->Node[n1].El + hyd->LinkSetting[k];
+            hyd->LinkStatus[k] = psvstatus(pr, k, status, hset,
+                                     hyd->NodeHead[n1], hyd->NodeHead[n2]);
+            break;
+        default:
+            continue;
+        }
+
+        // Check for a status change
+        if (status != hyd->LinkStatus[k])
+        {
+            if (rep->Statflag == FULL)
+            {
+                writestatchange(pr, k, status, hyd->LinkStatus[k]);
+            }
+            change = TRUE;
+        }
+    }
+    return change;
+}                       /* End of valvestatus() */
+
+
+int  linkstatus(EN_Project *pr)
+/*
+**--------------------------------------------------------------
+**  Input:   none
+**  Output:  returns 1 if any link changes status, 0 otherwise
+**  Purpose: determines new status for pumps, CVs, FCVs & pipes
+**           to tanks.
+**--------------------------------------------------------------
+*/
+{
+    int change = FALSE,             // Status change flag
+        k,                          // Link index
+        n1,                         // Start node index
+        n2;                         // End node index
+    double dh;                      // Head difference across link
+    StatType  status;               // Current status
+
+    EN_Network       *net = &pr->network;
+    hydraulics_t     *hyd = &pr->hydraulics;
+    report_options_t *rep = &pr->report;
+    Slink *link;
+
+    // Examine each link
+    for (k = 1; k <= net->Nlinks; k++)
+    {
+        link = &net->Link[k];
+        n1 = link->N1;
+        n2 = link->N2;
+        dh = hyd->NodeHead[n1] - hyd->NodeHead[n2];
+
+        // Re-open temporarily closed links (status = XHEAD or TEMPCLOSED)
+        status = hyd->LinkStatus[k];
+        if (status == XHEAD || status == TEMPCLOSED)
+        {
+            hyd->LinkStatus[k] = OPEN;
+        }
+
+        // Check for status changes in CVs and pumps
+        if (link->Type == EN_CVPIPE)
+        {
+            hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh,
+                                          hyd->LinkFlows[k]);
+        }
+        if (link->Type == EN_PUMP && hyd->LinkStatus[k] >= OPEN &&
+            hyd->LinkSetting[k] > 0.0)
+        {
+            hyd->LinkStatus[k] = pumpstatus(pr, k, -dh);
+        }
+
+        // Check for status changes in non-fixed FCVs
+        if (link->Type == EN_FCV && hyd->LinkSetting[k] != MISSING)
+        {
+            hyd->LinkStatus[k] = fcvstatus(pr, k, status, hyd->NodeHead[n1],
+                                           hyd->NodeHead[n2]);
+        }
+
+        // Check for flow into (out of) full (empty) tanks
+        if (n1 > net->Njuncs || n2 > net->Njuncs)
+        {
+            tankstatus(pr, k, n1, n2);
+        }
+
+        // Note any change in link status; do not revise link flow
+        if (status != hyd->LinkStatus[k])
+        {
+            change = TRUE;
+            if (rep->Statflag == FULL)
+            {
+                writestatchange(pr, k, status, hyd->LinkStatus[k]);
+            }
+        }
+    }
+    return change;
+}
+
+
+StatType  cvstatus(EN_Project *pr, StatType s, double dh, double q)
+/*
+**--------------------------------------------------
+**  Input:   s  = current link status
+**           dh = head loss across link
+**           q  = link flow
+**  Output:  returns new link status
+**  Purpose: updates status of a check valve link.
+**--------------------------------------------------
+*/
+{
+    hydraulics_t *hyd = &pr->hydraulics;
+
+    // Prevent reverse flow through CVs
+    if (ABS(dh) > hyd->Htol)
+    {
+        if (dh < -hyd->Htol)     return CLOSED;
+        else if (q < -hyd->Qtol) return CLOSED;
+        else                     return OPEN;
+    }
+    else
+    {
+        if (q < -hyd->Qtol) return CLOSED;
+        else                return s;
+    }
+}
+
+
+StatType  pumpstatus(EN_Project *pr, int k, double dh)
+/*
+**--------------------------------------------------
+**  Input:   k  = link index
+**           dh = head gain across link
+**  Output:  returns new pump status
+**  Purpose: updates status of an open pump.
+**--------------------------------------------------
+*/
+{
+    int   p;
+    double hmax;
+    hydraulics_t *hyd = &pr->hydraulics;
+    EN_Network *net = &pr->network;
+
+    // Find maximum head (hmax) pump can deliver
+    p = findpump(net, k);
+    if (net->Pump[p].Ptype == CONST_HP)
+    {
+        // Use huge value for constant HP pump
+        hmax = BIG;
+    }
+    else
+    {
+        // Use speed-adjusted shut-off head for other pumps
+        hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax;
+    }
+
+    // Check if currrent head gain exceeds pump's max. head
+    if (dh > hmax + hyd->Htol) return XHEAD;
+
+    // No check is made to see if flow exceeds pump's max. flow
+    return OPEN;
+}
+
+
+StatType  prvstatus(EN_Project *pr, int k, StatType s, double hset,
+                    double h1, double h2)
+/*
+**-----------------------------------------------------------
+**  Input:   k    = link index
+**           s    = current status
+**           hset = valve head setting
+**           h1   = head at upstream node
+**           h2   = head at downstream node
+**  Output:  returns new valve status
+**  Purpose: updates status of a pressure reducing valve.
+**-----------------------------------------------------------
+*/
+{
+    StatType  status;             // Valve's new status
+    double hml;                   // Head loss when fully opened
+    hydraulics_t *hyd = &pr->hydraulics;
+
+    double htol = hyd->Htol;
+    Slink *link = &pr->network.Link[k];
+
+    // Head loss when fully open
+    hml = link->Km * SQR(hyd->LinkFlows[k]);
+
+    // Rules for updating valve's status from current value s
+    status = s;
+    switch (s)
+    {
+    case ACTIVE:
+        if (hyd->LinkFlows[k] < -hyd->Qtol)  status = CLOSED;
+        else if (h1 - hml < hset - htol)     status = OPEN;
+        else                                 status = ACTIVE;
+        break;
+
+    case OPEN:
+        if (hyd->LinkFlows[k] < -hyd->Qtol)  status = CLOSED;
+        else if (h2 >= hset + htol)          status = ACTIVE;
+        else                                 status = OPEN;
+        break;
+
+    case CLOSED:
+        if (h1 >= hset + htol && h2 < hset - htol)   status = ACTIVE;
+        else if (h1 < hset - htol && h1 > h2 + htol) status = OPEN;
+        else                                         status = CLOSED;
+        break;
+
+    case XPRESSURE:
+        if (hyd->LinkFlows[k] < -hyd->Qtol) status = CLOSED;
+        break;
+
+    default:
+        break;
+    }
+    return status;
+}
+
+
+StatType  psvstatus(EN_Project *pr, int k, StatType s, double hset,
+                    double h1, double h2)
+/*
+**-----------------------------------------------------------
+**  Input:   k    = link index
+**           s    = current status
+**           hset = valve head setting
+**           h1   = head at upstream node
+**           h2   = head at downstream node
+**  Output:  returns new valve status
+**  Purpose: updates status of a pressure sustaining valve.
+**-----------------------------------------------------------
+*/
+{
+    StatType  status;             // Valve's new status
+    double hml;                   // Head loss when fully opened
+    hydraulics_t *hyd = &pr->hydraulics;
+
+    double htol = hyd->Htol;
+    Slink *link = &pr->network.Link[k];
+
+    // Head loss when fully open
+    hml = link->Km * SQR(hyd->LinkFlows[k]);
+
+    // Rules for updating valve's status from current value s
+    status = s;
+    switch (s)
+    {
+    case ACTIVE:
+        if (hyd->LinkFlows[k] < -hyd->Qtol) status = CLOSED;
+        else if (h2 + hml > hset + htol)    status = OPEN;
+        else                                status = ACTIVE;
+        break;
+
+    case OPEN:
+        if (hyd->LinkFlows[k] < -hyd->Qtol) status = CLOSED;
+        else if (h1 < hset - htol)          status = ACTIVE;
+        else                                status = OPEN;
+        break;
+
+    case CLOSED:
+        if (h2 > hset + htol && h1 > h2 + htol)       status = OPEN;
+        else if (h1 >= hset + htol && h1 > h2 + htol) status = ACTIVE;
+        else                                          status = CLOSED;
+        break;
+
+    case XPRESSURE:
+        if (hyd->LinkFlows[k] < -hyd->Qtol) status = CLOSED;
+        break;
+
+    default:
+        break;
+    }
+    return status;
+}
+
+
+StatType  fcvstatus(EN_Project *pr, int k, StatType s, double h1, double h2)
+/*
+**-----------------------------------------------------------
+**  Input:   k    = link index
+**           s    = current status
+**           h1   = head at upstream node
+**           h2   = head at downstream node
+**  Output:  returns new valve status
+**  Purpose: updates status of a flow control valve.
+**
+**    Valve status changes to XFCV if flow reversal.
+**    If current status is XFCV and current flow is
+**    above setting, then valve becomes active.
+**    If current status is XFCV, and current flow
+**    positive but still below valve setting, then
+**    status remains same.
+**-----------------------------------------------------------
+*/
+{
+    StatType  status;             // New valve status
+    hydraulics_t *hyd = &pr->hydraulics;
+
+    status = s;
+    if (h1 - h2 < -hyd->Htol)
+    {
+        status = XFCV;
+    }
+    else if (hyd->LinkFlows[k] < -hyd->Qtol)
+    {
+        status = XFCV;
+    }
+    else if (s == XFCV && hyd->LinkFlows[k] >= hyd->LinkSetting[k])
+    {
+        status = ACTIVE;
+    }
+    return status;
+}
+
+
+void  tankstatus(EN_Project *pr, int k, int n1, int n2)
+/*
+**----------------------------------------------------------------
+**  Input:   k  = link index
+**           n1 = start node of link
+**           n2 = end node of link
+**  Output:  none
+**  Purpose: closes link flowing into full or out of empty tank
+**----------------------------------------------------------------
+*/
+{
+    int   i, n;
+    double h, q;
+    Stank *tank;
+
+    hydraulics_t *hyd = &pr->hydraulics;
+    EN_Network *net = &pr->network;
+    Slink *link = &net->Link[k];
+
+    // Return if link is closed
+    if (hyd->LinkStatus[k] <= CLOSED) return;
+
+    // Make node n1 be the tank, reversing flow (q) if need be
+    q = hyd->LinkFlows[k];
+    i = n1 - net->Njuncs;
+    if (i <= 0)
+    {
+        i = n2 - net->Njuncs;
+        if (i <= 0) return;
+        n = n1;
+        n1 = n2;
+        n2 = n;
+        q = -q;
+    }
+
+    // Ignore reservoirs
+    tank = &net->Tank[i];
+    if (tank->A == 0.0) return;
+
+    // Find head difference across link
+    h = hyd->NodeHead[n1] - hyd->NodeHead[n2];
+
+    // If tank is full, then prevent flow into it
+    if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol)
+    {
+        // Case 1: Link is a pump discharging into tank
+        if (link->Type == EN_PUMP)
+        {
+            if (link->N2 == n1) hyd->LinkStatus[k] = TEMPCLOSED;
+        }
+
+        // Case 2: Downstream head > tank head
+        // (e.g., an open outflow check valve would close)
+        else if (cvstatus(pr, OPEN, h, q) == CLOSED) 
+        {
+            hyd->LinkStatus[k] = TEMPCLOSED;
+        }
+    }
+
+    // If tank is empty, then prevent flow out of it
+    if (hyd->NodeHead[n1] <= tank->Hmin + hyd->Htol)
+    {
+        // Case 1: Link is a pump discharging from tank
+        if (link->Type == EN_PUMP)
+        {
+            if (link->N1 == n1) hyd->LinkStatus[k] = TEMPCLOSED;
+        }
+
+        // Case 2: Tank head > downstream head
+        // (e.g., a closed outflow check valve would open)
+        else if (cvstatus(pr, CLOSED, h, q) == OPEN)
+        {
+            hyd->LinkStatus[k] = TEMPCLOSED;
+        }
+    }
+}

src/inpfile.c

@@ -26,14 +26,14 @@ data describing a piping network to a file in EPANET's text format.
 #else
 #include <stdlib.h>
 #endif
+
 #include "hash.h"
 #include "text.h"
 #include "types.h"
-#include "epanet2.h"
 #include "funcs.h"
 #include <math.h>
-#define EXTERN extern
+//#define EXTERN extern
-#include "vars.h"
+//#include "vars.h"
 
 /* Defined in enumstxt.h in EPANET.C */
 extern char *LinkTxt[];
@@ -618,11 +618,15 @@ 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]);
-      fprintf(f, "\n HEADERROR           %-.8f", hyd->HeadErrorLimit * pr->Ucf[HEAD]);
+  fprintf(f, "\n FLOWCHANGE          %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]);
-  }
+
-  if (hyd->FlowChangeLimit > 0.0) {
+  if (hyd->DemandModel == PDA)
-      fprintf(f, "\n FLOWCHANGE          %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]);
+  {
+    fprintf(f, "\n DEMAND MODEL        PDA");
+    fprintf(f, "\n MINIMUM PRESSURE    %-.4f", hyd->Pmin * pr->Ucf[PRESSURE]);
+    fprintf(f, "\n REQUIRED PRESSURE   %-.4f", hyd->Preq * pr->Ucf[PRESSURE]);
+    fprintf(f, "\n PRESSURE EXPONENT   %-.4f", hyd->Pexp);
   }
 
   /* Write [REPORT] section */

src/input1.c

@@ -33,8 +33,6 @@ AUTHOR:     L. Rossman
 #include "epanet2.h"
 #include "funcs.h"
 #include <math.h>
-#define EXTERN extern
-#include "vars.h"
 
 /*
   --------------------- Module Global Variables  ----------------------
@@ -128,8 +126,12 @@ void setdefaults(EN_Project *pr)
   hyd->Qtol = QTOL;        /* Default flow tolerance         */
   hyd->Hacc = HACC;        /* Default hydraulic accuracy     */
 
-  hyd->FlowChangeLimit = 0.0; /* Default flow change limit   */
+  hyd->FlowChangeLimit = 0.0; // Default flow change limit
-  hyd->HeadErrorLimit = 0.0;  /* Default head error limit    */
+  hyd->HeadErrorLimit = 0.0;  // Default head error limit
+  hyd->DemandModel = DDA;     // Demand driven analysis
+  hyd->Pmin = 0.0;            // Minimum demand pressure (ft)
+  hyd->Preq = 0.0;            // Required demand pressure (ft)
+  hyd->Pexp = 0.5;            // Pressure function exponent
 
   qu->Ctol = MISSING;      /* No pre-set quality tolerance   */
   hyd->MaxIter = MAXITER;  /* Default max. hydraulic trials  */
@@ -581,6 +583,8 @@ void convertunits(EN_Project *pr)
       demand->Base /= pr->Ucf[DEMAND];
     }
   }
+  hyd->Pmin /= pr->Ucf[PRESSURE];
+  hyd->Preq /= pr->Ucf[PRESSURE];
 
   /* Convert emitter discharge coeffs. to head loss coeff. */
   ucf = pow(pr->Ucf[FLOW], hyd->Qexp) / pr->Ucf[PRESSURE];

src/input2.c

@@ -37,8 +37,6 @@ The following utility functions are all called from INPUT3.C
 #include "epanet2.h"
 #include "funcs.h"
 #include <math.h>
-#define EXTERN extern
-#include "vars.h"
 
 #define MAXERRS 10 /* Max. input errors reported        */
 

src/input3.c

@@ -32,12 +32,11 @@ All functions in this module are called from newline() in INPUT2.C.
 #include "text.h"
 #include "types.h"
 #include <math.h>
-#define EXTERN extern
-#include "vars.h"
 
 /* Defined in enumstxt.h in EPANET.C */
 extern char *MixTxt[];
 extern char *Fldname[];
+extern char *DemandModelTxt[];
 
 /* Defined in INPUT2.C */
 
@@ -1763,6 +1762,7 @@ int optionchoice(EN_Project *pr, int n)
 **    VERIFY              filename
 **    UNBALANCED          STOP/CONTINUE {Niter}
 **    PATTERN             id
+**    DEMAND MODEL        DDA/PDA/PPA
 **--------------------------------------------------------------
 */
 {
@@ -1771,131 +1771,125 @@ int optionchoice(EN_Project *pr, int n)
   quality_t *qu = &pr->quality;
   parser_data_t *par = &pr->parser;
   out_file_t *out = &pr->out_files;
-
+  int choice;
   
   /* Check if 1st token matches a parameter name and */
   /* process the input for the matched parameter     */
-  if (n < 0)
+  if (n < 0) return (201);
-    return (201);
+  if (match(par->Tok[0], w_UNITS))
-  if (match(par->Tok[0], w_UNITS)) {
+  {
-    if (n < 1)
+    if (n < 1) return (0);
-      return (0);
+    else if (match(par->Tok[1], w_CFS))   par->Flowflag = CFS;
-    else if (match(par->Tok[1], w_CFS))
+    else if (match(par->Tok[1], w_GPM))   par->Flowflag = GPM;
-      par->Flowflag = CFS;
+    else if (match(par->Tok[1], w_AFD))   par->Flowflag = AFD;
-    else if (match(par->Tok[1], w_GPM))
+    else if (match(par->Tok[1], w_MGD))   par->Flowflag = MGD;
-      par->Flowflag = GPM;
+    else if (match(par->Tok[1], w_IMGD))  par->Flowflag = IMGD;
-    else if (match(par->Tok[1], w_AFD))
+    else if (match(par->Tok[1], w_LPS))   par->Flowflag = LPS;
-      par->Flowflag = AFD;
+    else if (match(par->Tok[1], w_LPM))   par->Flowflag = LPM;
-    else if (match(par->Tok[1], w_MGD))
+    else if (match(par->Tok[1], w_CMH))   par->Flowflag = CMH;
-      par->Flowflag = MGD;
+    else if (match(par->Tok[1], w_CMD))   par->Flowflag = CMD;
-    else if (match(par->Tok[1], w_IMGD))
+    else if (match(par->Tok[1], w_MLD))   par->Flowflag = MLD;
-      par->Flowflag = IMGD;
+    else if (match(par->Tok[1], w_SI))    par->Flowflag = LPS;
-    else if (match(par->Tok[1], w_LPS))
+    else return (201);
-      par->Flowflag = LPS;
+  }
-    else if (match(par->Tok[1], w_LPM))
+
-      par->Flowflag = LPM;
+  else if (match(par->Tok[0], w_PRESSURE))
-    else if (match(par->Tok[1], w_CMH))
+  {
-      par->Flowflag = CMH;
+    if (n < 1) return (0);
-    else if (match(par->Tok[1], w_CMD))
+    else if (match(par->Tok[1], w_EXPONENT)) return -1;
-      par->Flowflag = CMD;
+    else if (match(par->Tok[1], w_PSI))    par->Pressflag = PSI;
-    else if (match(par->Tok[1], w_MLD))
+    else if (match(par->Tok[1], w_KPA))    par->Pressflag = KPA;
-      par->Flowflag = MLD;
+    else if (match(par->Tok[1], w_METERS)) par->Pressflag = METERS;
-    else if (match(par->Tok[1], w_SI))
+    else return (201);
-      par->Flowflag = LPS;
+  }
-    else
+  
-      return (201);
+  else if (match(par->Tok[0], w_HEADLOSS))
-  } else if (match(par->Tok[0], w_PRESSURE)) {
+  {
-    if (n < 1)
+    if (n < 1)  return (0);
-      return (0);
+    else if (match(par->Tok[1], w_HW))   hyd->Formflag = HW;
-    else if (match(par->Tok[1], w_PSI))
+    else if (match(par->Tok[1], w_DW))   hyd->Formflag = DW;
-      par->Pressflag = PSI;
+    else if (match(par->Tok[1], w_CM))   hyd->Formflag = CM;
-    else if (match(par->Tok[1], w_KPA))
+    else return (201);
-      par->Pressflag = KPA;
+  }
-    else if (match(par->Tok[1], w_METERS))
+  
-      par->Pressflag = METERS;
+  else if (match(par->Tok[0], w_HYDRAULIC))
-    else
+  {
-      return (201);
+    if (n < 2) return (0);
-  } else if (match(par->Tok[0], w_HEADLOSS)) {
+    else if (match(par->Tok[1], w_USE))  out->Hydflag = USE;
-    if (n < 1)
+    else if (match(par->Tok[1], w_SAVE)) out->Hydflag = SAVE;
-      return (0);
+    else return (201);
-    else if (match(par->Tok[1], w_HW))
-      hyd->Formflag = HW;
-    else if (match(par->Tok[1], w_DW))
-      hyd->Formflag = DW;
-    else if (match(par->Tok[1], w_CM))
-      hyd->Formflag = CM;
-    else
-      return (201);
-  } else if (match(par->Tok[0], w_HYDRAULIC)) {
-    if (n < 2)
-      return (0);
-    else if (match(par->Tok[1], w_USE))
-      out->Hydflag = USE;
-    else if (match(par->Tok[1], w_SAVE))
-      out->Hydflag = SAVE;
-    else
-      return (201);
     strncpy(out->HydFname, par->Tok[2], MAXFNAME);
-  } else if (match(par->Tok[0], w_QUALITY)) {
+  }
-    if (n < 1)
+
-      return (0);
+  else if (match(par->Tok[0], w_QUALITY))
-    else if (match(par->Tok[1], w_NONE))
+  {
-      qu->Qualflag = NONE;
+    if (n < 1) return (0);
-    else if (match(par->Tok[1], w_CHEM))
+    else if (match(par->Tok[1], w_NONE))  qu->Qualflag = NONE;
-      qu->Qualflag = CHEM;
+    else if (match(par->Tok[1], w_CHEM))  qu->Qualflag = CHEM;
-    else if (match(par->Tok[1], w_AGE))
+    else if (match(par->Tok[1], w_AGE))   qu->Qualflag = AGE;
-      qu->Qualflag = AGE;
+    else if (match(par->Tok[1], w_TRACE)) qu->Qualflag = TRACE;
-    else if (match(par->Tok[1], w_TRACE))
+    else
-      qu->Qualflag = TRACE;
+    {
-    else {
+        qu->Qualflag = CHEM;
-      qu->Qualflag = CHEM;
+        strncpy(qu->ChemName, par->Tok[1], MAXID);
-      strncpy(qu->ChemName, par->Tok[1], MAXID);
+        if (n >= 2) strncpy(qu->ChemUnits, par->Tok[2], MAXID);
-      if (n >= 2)
-        strncpy(qu->ChemUnits, par->Tok[2], MAXID);
     }
     if (qu->Qualflag == TRACE) /* Source tracing option */
     {
       /* Copy Trace Node ID to par->Tok[0] for error reporting */
       strcpy(par->Tok[0], "");
-      if (n < 2)
+      if (n < 2) return (212);
-        return (212);
       strcpy(par->Tok[0], par->Tok[2]);
       qu->TraceNode = findnode(net,par->Tok[2]);
-      if (qu->TraceNode == 0)
+      if (qu->TraceNode == 0)  return (212);
-        return (212);
       strncpy(qu->ChemName, u_PERCENT, MAXID);
       strncpy(qu->ChemUnits, par->Tok[2], MAXID);
     }
-    if (qu->Qualflag == AGE) {
+    if (qu->Qualflag == AGE)
+    {
       strncpy(qu->ChemName, w_AGE, MAXID);
       strncpy(qu->ChemUnits, u_HOURS, MAXID);
     }
-  } else if (match(par->Tok[0], w_MAP)) {
+  }
-    if (n < 1)
+  
-      return (0);
+  else if (match(par->Tok[0], w_MAP))
+  {
+    if (n < 1) return (0);
     strncpy(pr->MapFname, par->Tok[1], MAXFNAME); /* Map file name */
-  } else if (match(par->Tok[0], w_VERIFY)) {
+  }
+  
+  else if (match(par->Tok[0], w_VERIFY))
+  {
     /* Backward compatibility for verification file */
-  } else if (match(par->Tok[0], w_UNBALANCED)) /* Unbalanced option */
+  }
+  
+  else if (match(par->Tok[0], w_UNBALANCED)) /* Unbalanced option */
   {
-    if (n < 1)
+    if (n < 1) return (0);
-      return (0);
+    if (match(par->Tok[1], w_STOP)) hyd->ExtraIter = -1;
-    if (match(par->Tok[1], w_STOP))
+    else if (match(par->Tok[1], w_CONTINUE))
-      hyd->ExtraIter = -1;
+    {
-    else if (match(par->Tok[1], w_CONTINUE)) {
+      if (n >= 2)  hyd->ExtraIter = atoi(par->Tok[2]);
-      if (n >= 2)
+      else         hyd->ExtraIter = 0;
-        hyd->ExtraIter = atoi(par->Tok[2]);
+    }
-      else
+    else return (201);
-        hyd->ExtraIter = 0;
+  }
-    } else
+  
-      return (201);
+  else if (match(par->Tok[0], w_PATTERN)) /* Pattern option */
-  } else if (match(par->Tok[0], w_PATTERN)) /* Pattern option */
   {
-    if (n < 1)
+    if (n < 1) return (0);
-      return (0);
     strncpy(par->DefPatID, par->Tok[1], MAXID);
-  } else
+  }
-    return (-1);
+
+  else if (match(par->Tok[0], w_DEMAND))
+  {
+      if (n < 2) return 0;
+      if (!match(par->Tok[1], w_MODEL)) return -1;
+      choice = findmatch(par->Tok[2], DemandModelTxt);
+      if (choice < 0) return 201;
+      hyd->DemandModel = choice;
+  }
+  else return (-1);
   return (0);
 } /* end of optionchoice */
 
@@ -1916,6 +1910,9 @@ int optionvalue(EN_Project *pr, int n)
 
 **    HEADLIMIT           value
 **    FLOWLIMIT           value
+**    MINIMUM PRESSURE    value
+**    REQUIRED PRESSURE   value
+**    PRESSURE EXPONENT   value
 
 **    TOLERANCE           value
 **    SEGMENTS            value  (not used)
@@ -1939,39 +1936,40 @@ int optionvalue(EN_Project *pr, int n)
   double y;
 
   /* Check for obsolete SEGMENTS keyword */
-  //if (match(par->Tok[0], w_SEGMENTS))
+  if (match(tok0, w_SEGMENTS)) return (0);
-  if (match(tok0, w_SEGMENTS))
-          return (0);
 
   /* Check for missing value (which is permissible) */
   if (match(tok0, w_SPECGRAV) || match(tok0, w_EMITTER) ||
-      match(tok0, w_DEMAND))
+      match(tok0, w_DEMAND)   || match(tok0, w_MINIMUM) ||
+      match(tok0, w_REQUIRED) || match(tok0, w_PRESSURE) ||
+      match(tok0, w_PRECISION))
+  {
     nvalue = 2;
-  if (n < nvalue)
+  }
-    return (0);
+  if (n < nvalue) return (0);
 
   /* Check for valid numerical input */
-  if (!getfloat(par->Tok[nvalue], &y))
+  if (!getfloat(par->Tok[nvalue], &y)) return (213);
-    return (213);
 
   /* Check for WQ tolerance option (which can be 0) */
-  if (match(tok0, w_TOLERANCE)) {
+  if (match(tok0, w_TOLERANCE))
-    if (y < 0.0)
+  {
-      return (213);
+    if (y < 0.0) return (213);
     qu->Ctol = y; /* Quality tolerance*/
     return (0);
   }
 
   /* Check for Diffusivity option */
-  if (match(tok0, w_DIFFUSIVITY)) {
+  if (match(tok0, w_DIFFUSIVITY))
-    if (y < 0.0)
+  {
-      return (213);
+    if (y < 0.0) return (213);
     qu->Diffus = y;
     return (0);
   }
 
   /* Check for Damping Limit option */ 
-  if (match(tok0, w_DAMPLIMIT)) {
+  if (match(tok0, w_DAMPLIMIT))
+  {
     hyd->DampLimit = y;
     return (0);
   }
@@ -1992,41 +1990,51 @@ int optionvalue(EN_Project *pr, int n)
       return 0;
   }
 
+  /* Check for pressure dependent demand params */
+  else if (match(tok0, w_MINIMUM))
+  {
+      if (y < 0.0) return 213;
+      hyd->Pmin = y;
+      return 0;
+  }
+  else if (match(tok0, w_REQUIRED))
+  {
+      if (y < 0.0) return 213;
+      hyd->Preq = y;
+      return 0;
+  }
+  else if (match(tok0, w_PRESSURE))
+  {
+      if (y < 0.0) return 213;
+      hyd->Pexp = y;
+      return 0;
+  }
+
   /* All other options must be > 0 */
-  if (y <= 0.0)
+  if (y <= 0.0) return (213);
-    return (213);
 
   /* Assign value to specified option */
-  if (match(tok0, w_VISCOSITY))
+  if (match(tok0, w_VISCOSITY))     hyd->Viscos = y; /* Viscosity */
-    hyd->Viscos = y; /* Viscosity */
+  else if (match(tok0, w_SPECGRAV)) hyd->SpGrav = y; /* Spec. gravity */
-  else if (match(tok0, w_SPECGRAV))
+  else if (match(tok0, w_TRIALS))   hyd->MaxIter = (int)y;  /* Max. trials */
-    hyd->SpGrav = y; /* Spec. gravity */
-  else if (match(tok0, w_TRIALS))
-    hyd->MaxIter = (int)y;                   /* Max. trials */
   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(tok0, w_HTOL))
+  else if (match(tok0, w_HTOL))  hyd->Htol = y;
-    hyd->Htol = y;
+  else if (match(tok0, w_QTOL))  hyd->Qtol = y;
-  else if (match(tok0, w_QTOL))
+  else if (match(tok0, w_RQTOL))
-    hyd->Qtol = y;
+  {
-  else if (match(tok0, w_RQTOL)) {
+    if (y >= 1.0) return (213);
-    if (y >= 1.0)
-      return (213);
     hyd->RQtol = y;
-  } else if (match(tok0, w_CHECKFREQ))
+  }
-    hyd->CheckFreq = (int)y;
+  else if (match(tok0, w_CHECKFREQ))  hyd->CheckFreq = (int)y;
-  else if (match(tok0, w_MAXCHECK))
+  else if (match(tok0, w_MAXCHECK))   hyd->MaxCheck = (int)y;
-    hyd->MaxCheck = (int)y;
+  else if (match(tok0, w_EMITTER))    hyd->Qexp = 1.0 / y;
-  else if (match(tok0, w_EMITTER))
+  else if (match(tok0, w_DEMAND))     hyd->Dmult = y;
-    hyd->Qexp = 1.0 / y;
+  else return (201);
-  else if (match(tok0, w_DEMAND))
-    hyd->Dmult = y;
-  else
-    return (201);
   return (0);
 } /* end of optionvalue */
 
@@ -2091,7 +2099,8 @@ int getpumpcurve(EN_Project *pr, int n)
   return (0);
 }
 
-int powercurve(double h0, double h1, double h2, double q1, double q2, double *a, double *b, double *c)
+int powercurve(double h0, double h1, double h2, double q1, double q2,
+               double *a, double *b, double *c)
 /*
 **---------------------------------------------------------
 **  Input:   h0 = shutoff head

src/output.c

@@ -19,14 +19,12 @@ AUTHOR:     L. Rossman
 #else
 #include <stdlib.h>
 #endif
-#include "epanet2.h"
+
 #include "funcs.h"
 #include "text.h"
 #include "types.h"
 #include <math.h>
-#define EXTERN extern
 #include "hash.h"
-#include "vars.h"
 
 /* write x[1] to x[n] to file */
 size_t f_save(REAL4 *x, int n, FILE *file) {

src/quality.c

@@ -54,15 +54,13 @@ AUTHOR:     L. Rossman
 #else
 #include <stdlib.h>
 #endif
+
 #include "hash.h"
 #include "text.h"
 #include "types.h"
-#include "epanet2.h"
 #include "funcs.h"
 #include <math.h>
-#define EXTERN extern
 #include "mempool.h"
-#include "vars.h"
 
 /*
 ** Macros to identify upstream & downstream nodes of a link

src/report.c

@@ -32,15 +32,13 @@ formatted string S to the report file.
 #else
 #include <stdlib.h>
 #endif
-#include "epanet2.h"
+
 #include "funcs.h"
 #include "hash.h"
 #include "text.h"
 #include "types.h"
 #include <math.h>
 #include <time.h>
-#define EXTERN extern
-#include "vars.h"
 
 #undef WINDOWS
 #ifdef _WIN32
@@ -50,12 +48,13 @@ formatted string S to the report file.
 
 #define MAXCOUNT 10 /* Max. # of disconnected nodes listed */
 
-/* Defined in enumstxt.h in EPANET.C */
+/* Defined in enumstxt.h */
 extern char *NodeTxt[];
 extern char *LinkTxt[];
 extern char *StatTxt[];
 extern char *TstatTxt[];
 extern char *RptFormTxt[];
+extern char *DemandModelTxt[];
 
 typedef REAL4 *Pfloat;
 void writenodetable(EN_Project *pr, Pfloat *);
@@ -224,6 +223,8 @@ void writesummary(EN_Project *pr)
   writeline(pr, s);
   sprintf(s, FMT25, RptFormTxt[hyd->Formflag]);
   writeline(pr, s);
+  sprintf(s, FMT25a, DemandModelTxt[hyd->DemandModel]);
+  writeline(pr, s);
   sprintf(s, FMT26, time->Hstep * pr->Ucf[TIME], rep->Field[TIME].Units);
   writeline(pr, s);
   sprintf(s, FMT27, hyd->Hacc);

src/rules.c

@@ -29,13 +29,11 @@ AUTHOR:     L. Rossman
 #else
 #include <stdlib.h>
 #endif
-#include "epanet2.h"
+
 #include "funcs.h"
 #include "hash.h"
 #include "text.h"
 #include "types.h"
-#define EXTERN extern
-#include "vars.h"
 
 enum Rulewords {
   r_RULE,

src/text.h

@@ -145,6 +145,12 @@ AUTHOR:     L. Rossman
 #define   w_FLOWCHANGE  "FLOWCHANGE"
 #define   w_HEADERROR   "HEADERROR"
 
+#define   w_MODEL       "MODEL"
+#define   w_DDA         "DDA"
+#define   w_PDA         "PDA"
+#define   w_REQUIRED    "REQ"
+#define   w_EXPONENT    "EXP"
+
 #define   w_SECONDS     "SEC"
 #define   w_MINUTES     "MIN"
 #define   w_HOURS       "HOU"
@@ -323,6 +329,9 @@ AUTHOR:     L. Rossman
 #define   t_perM3       "  /m3"
 #define   t_perMGAL     "/Mgal"
 #define   t_DIFFER      "DIFFERENTIAL"
+#define   t_FIXED       "Fixed Demands"
+#define   t_POWER       "Power Function"
+#define   t_ORIFICE     "Orifice Flow"
 
 
 /* ------------------ Format Messages ------------------*/
@@ -362,6 +371,7 @@ AUTHOR:     L. Rossman
 #define FMT23  "    Number of Pumps ................... %-d"
 #define FMT24  "    Number of Valves .................. %-d"
 #define FMT25  "    Headloss Formula .................. %s"
+#define FMT25a "    Nodal Demand Model ................ %s"
 #define FMT26  "    Hydraulic Timestep ................ %-.2f %s"
 #define FMT27  "    Hydraulic Accuracy ................ %-.6f"
 

src/types.h

@@ -298,6 +298,11 @@ typedef enum {
     MAX_ENERGY_STATS
 } EnergyStats;
 
+typedef enum {
+    DDA,        // Demand Driven Analysis    
+    PDA         // Pressure Driven Analysis
+} DemandModelType;
+
 /*
 ------------------------------------------------------
    Global Data Structures                             
@@ -759,7 +764,8 @@ typedef struct {
 
 typedef struct {
   double  
-  *NodeDemand,           /* Node actual demand           */
+  *NodeDemand,           // Node actual total outflow
+  *DemandFlows,          // Demand outflows
   *EmitterFlows,         /* Emitter flows                */
   *LinkSetting,          /* Link settings                */
   *LinkFlows,            /* Link flows                   */
@@ -768,10 +774,13 @@ typedef struct {
   Qtol,                  /* Flow rate tolerance          */
   RQtol,                 /* Flow resistance tolerance    */
   Hexp,                  /* Exponent in headloss formula */
-  Qexp,                  /* Exponent in orifice formula  */
+  Qexp,                  /* Exponent in emitter formula  */
+  Pexp,                  // Exponent in demand formula
+  Pmin,                  // Pressure needed for any demand
+  Preq,                  // Pressure needed for full demand
   Dmult,                 /* Demand multiplier            */
-  Hacc,                  /* Hydraulics solution accuracy */
 
+  Hacc,                  /* Hydraulics solution accuracy */
   FlowChangeLimit,       /* Hydraulics flow change limit */
   HeadErrorLimit,        /* Hydraulics head error limit  */
 
@@ -787,7 +796,8 @@ typedef struct {
   
   int
   DefPat,                /* Default demand pattern       */
-  Epat;                  /* Energy cost time pattern     */
+  Epat,                  /* Energy cost time pattern     */
+  DemandModel;           // Fixed or pressure dependent
 
   StatType  
   *LinkStatus,           /* Link status                  */
@@ -805,8 +815,10 @@ typedef struct {
   Formflag;              /* Hydraulic formula flag       */
   
   /* Info about hydraulic solution */
-  double relativeError;
+  double RelativeError;
-  int iterations;
+  double MaxHeadError;
+  double MaxFlowChange;
+  int    Iterations;
   
   /* Flag used to halt taking further time steps */
   int Haltflag;