JIANG/EPANET
1
0
Fork 0
Code Issues Pull Requests Actions 4 Packages Projects Releases Wiki Activity

Merge pull request #210 from LRossman/contributor-lr

Browse Source 
Pressure Dependent Demands added to address issue 163
This commit is contained in:
Michael Tryby
2018-08-24 17:30:38 -04:00
committed by GitHub
parent 761199b08e 93061191c5
commit c4b7c90634
25 changed files with 1524 additions and 1028 deletions
Download Patch File Download Diff File Expand all files Collapse all files

3
README.md
View File

@@ -13,8 +13,9 @@ The EPANET Library is a pressurized pipe network hydraulic and water quality ana
Please see the [`version 2.1` Release Notes](https://github.com/OpenWaterAnalytics/EPANET/blob/master/ReleaseNotes2_1.md) for information relevant to users of the previous official version (2.00.12). If you would like to contribute by addressing any of the outstanding [Issues](https://github.com/OpenWaterAnalytics/EPANET/issues), then please comment on the Issue, then Fork this repo to your own account and base your commits on the [`dev` branch](https://github.com/OpenWaterAnalytics/EPANET/tree/dev). Once you are finished, you can open a Pull Request to test the code and discuss merging your changes back into the community respository.
A step-by-step tutorial on how to contribute to EPANET using GitHub is also [available] (http://www.slideshare.net/demetriseliades/contributing-to-epanet-using-github-in-windows).
A step-by-step tutorial on how to contribute to EPANET using GitHub is also [available](http://www.slideshare.net/demetriseliades/contributing-to-epanet-using-github-in-windows).
__Note:__ This repository is not affiliated with, or endorsed by, the USEPA. For the last "official" release of EPANET (2.00.12 UI and Toolkit) please go to the [EPA's GitHub repo](https://github.com/USEPA/Water-Distribution-Network-Model) or [the USEPA website](http://www2.epa.gov/water-research/epanet). It is also not the graphical user interface version. This is the hydraulic and water quality solver engine.
However, if you are interested in extending EPANET for academic, personal, or commercial use, then you've come to the right place. For community discussion, FAQ, and roadmapping of the project, go to the [Community Forum](http://community.wateranalytics.org/category/epanet).

96
ReleaseNotes2_2.md Normal file
View File

@@ -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.

8
include/epanet2.bas
View File

@@ -78,6 +78,8 @@ Public Const EN_NEXTEVENT = 14
Public Const EN_ITERATIONS = 0
Public Const EN_RELATIVEERROR = 1
Public Const EN_MAXHEADERROR = 2
Public Const EN_MAXFLOWCHANGE = 3
Public Const EN_NODECOUNT = 0 'Component counts
Public Const EN_TANKCOUNT = 1
@@ -126,6 +128,9 @@ Public Const EN_MLD = 7
Public Const EN_CMH = 8
Public Const EN_CMD = 9
Public Const EN_DDA = 0 ' Demand driven analysis
Public Const EN_PDA = 1 ' Pressure driven analysis
Public Const EN_TRIALS = 0 ' Misc. options
Public Const EN_ACCURACY = 1
Public Const EN_TOLERANCE = 2
@@ -233,6 +238,9 @@ Public Const EN_G_CURVE = 4 ' General\default curve
Declare Function ENgetversion Lib "epanet2.dll" (value As Long) As Long
Declare Function ENgetdemandmodel Lib "epanet2.dll" (type as long, pmin as Single, preq as Single, pexp as Single) As Long
Declare Function ENsetdemandmodel Lib "epanet2.dll" (ByVal type as long, ByVal pmin as Single, ByVal preq as Single, ByVal pexp as Single) As Long
Declare Function ENsetflowunits Lib "epanet2.dll" (ByVal code As Long) As Long
Declare Function ENsetcontrol Lib "epanet2.dll" (ByVal Cindex As Long, ByVal Ctype As Long, ByVal Lindex As Long, ByVal setting As Single, ByVal Nindex As Long, ByVal Level As Single) As Long
Declare Function ENsetnodevalue Lib "epanet2.dll" (ByVal index As Long, ByVal code As Long, ByVal value As Single) As Long

55
include/epanet2.h
View File

@@ -141,10 +141,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 {
@@ -153,7 +154,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;
@@ -208,6 +209,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 {
@@ -272,8 +277,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
@@ -524,7 +527,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
@@ -806,7 +831,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.
@@ -1131,15 +1156,14 @@ extern "C" {
int DLLEXPORT ENdeletelink(int linkIndex);
/***************************************************
Threadsafe versions of all epanet functions
***************************************************/
int DLLEXPORT EN_alloc(EN_Project **p);
int DLLEXPORT EN_free(EN_Project *p);
int DLLEXPORT EN_createproject(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);
@@ -1171,6 +1195,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);
@@ -1236,11 +1264,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

4
src/enumstxt.h
View File

@@ -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,

102
src/epanet.c
View File

@@ -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"
#include "vars.h"
// This single global variable is used only when the library is called
// 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);
}
@@ -536,23 +550,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
int DLLEXPORT EN_alloc(EN_Project **p)
/// Create an EPANET project
int DLLEXPORT EN_createproject(EN_Project **p)
{
EN_Project *project = calloc(1, sizeof(EN_Project));
*p = project;
return 0;
EN_Project *project = calloc(1, sizeof(EN_Project));
if (project == NULL) return 101;
*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;
}
@@ -562,9 +585,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
HeadlossFormula = headloss formula flag
** UnitsType = flow units flag
** HeadlossFormula = headloss formula flag
** Output: none
** Returns: error code
** Purpose: opens EPANET
@@ -1508,6 +1530,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;
@@ -1661,11 +1707,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;
}
@@ -1729,8 +1781,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;
*y = p->network.Coord[index].Y;
*x = (EN_API_FLOAT_TYPE)p->network.Coord[index].X;
*y = (EN_API_FLOAT_TYPE)p->network.Coord[index].Y;
return 0;
}
@@ -4258,7 +4310,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);
@@ -4729,7 +4781,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);
}
@@ -4863,7 +4915,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);
}
@@ -4911,7 +4963,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);
}

22
src/funcs.h
View File

@@ -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 */

360
src/hydcoeffs.c
View File

@@ -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 hlosscoeff(EN_Project *pr, int k);
//void matrixcoeffs(EN_Project *pr);
//void resistcoeff(EN_Project *pr, int k);
//void headlosscoeffs(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:
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);
switch (net->Link[k].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 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
else {
sol->F[sol->Row[n2]] += (sol->P[k] * hyd->NodeHead[n1]);
}
// Node n1 is a tank/reservoir
else 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
else {
sol->F[sol->Row[n1]] += (sol->P[k] * hyd->NodeHead[n2]);
}
// Node n2 is a tank/reservoir
else 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) {
continue; // Valve status fixed
}
n1 = link->N1; // Start & end nodes
n1 = link->N1;
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) {
continue;
}
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) {
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
if (hyd->Formflag == DW) {
// Use custom function for Darcy-Weisbach formula
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) {
r = n * r * pow(q, n - 1.0);
}
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);
@@ -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)
k2++;
if (k2 == 0)
k2++;
else if (k2 == npts)
k2--;
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
@@ -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) {
valvecoeff(pr, k);
}
if (hyd->LinkStatus[k] == CLOSED) 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) {
p = hyd->RQtol;
}
if (p < hyd->RQtol) p = hyd->RQtol;
sol->P[k] = 1.0 / p;
sol->Y[k] = flow / 2.0;
}

26
src/hydraul.c
View File

@@ -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. */

1007
src/hydsolver.c
View File

File diff suppressed because it is too large Load Diff

462
src/hydstatus.c Normal file
View File

@@ -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;
}
}
}

13
src/inpfile.c
View File

@@ -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
#include "vars.h"
//#define EXTERN extern
//#include "vars.h"
/* Defined in enumstxt.h in EPANET.C */
extern char *LinkTxt[];
@@ -624,6 +624,13 @@ int saveinpfile(EN_Project *pr, char *fname)
if (hyd->FlowChangeLimit > 0.0) {
fprintf(f, "\n FLOWCHANGE %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]);
}
if (hyd->DemandModel == PDA)
{
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 */

12
src/input1.c
View File

@@ -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->HeadErrorLimit = 0.0; /* Default head error limit */
hyd->FlowChangeLimit = 0.0; // Default flow change 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];

2
src/input2.c
View File

@@ -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 */

301
src/input3.c
View File

@@ -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)
return (201);
if (match(par->Tok[0], w_UNITS)) {
if (n < 1)
return (0);
else if (match(par->Tok[1], w_CFS))
par->Flowflag = CFS;
else if (match(par->Tok[1], w_GPM))
par->Flowflag = GPM;
else if (match(par->Tok[1], w_AFD))
par->Flowflag = AFD;
else if (match(par->Tok[1], w_MGD))
par->Flowflag = MGD;
else if (match(par->Tok[1], w_IMGD))
par->Flowflag = IMGD;
else if (match(par->Tok[1], w_LPS))
par->Flowflag = LPS;
else if (match(par->Tok[1], w_LPM))
par->Flowflag = LPM;
else if (match(par->Tok[1], w_CMH))
par->Flowflag = CMH;
else if (match(par->Tok[1], w_CMD))
par->Flowflag = CMD;
else if (match(par->Tok[1], w_MLD))
par->Flowflag = MLD;
else if (match(par->Tok[1], w_SI))
par->Flowflag = LPS;
else
return (201);
} else if (match(par->Tok[0], w_PRESSURE)) {
if (n < 1)
return (0);
else if (match(par->Tok[1], w_PSI))
par->Pressflag = PSI;
else if (match(par->Tok[1], w_KPA))
par->Pressflag = KPA;
else if (match(par->Tok[1], w_METERS))
par->Pressflag = METERS;
else
return (201);
} else if (match(par->Tok[0], w_HEADLOSS)) {
if (n < 1)
return (0);
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);
if (n < 0) return (201);
if (match(par->Tok[0], w_UNITS))
{
if (n < 1) return (0);
else if (match(par->Tok[1], w_CFS)) par->Flowflag = CFS;
else if (match(par->Tok[1], w_GPM)) par->Flowflag = GPM;
else if (match(par->Tok[1], w_AFD)) par->Flowflag = AFD;
else if (match(par->Tok[1], w_MGD)) par->Flowflag = MGD;
else if (match(par->Tok[1], w_IMGD)) par->Flowflag = IMGD;
else if (match(par->Tok[1], w_LPS)) par->Flowflag = LPS;
else if (match(par->Tok[1], w_LPM)) par->Flowflag = LPM;
else if (match(par->Tok[1], w_CMH)) par->Flowflag = CMH;
else if (match(par->Tok[1], w_CMD)) par->Flowflag = CMD;
else if (match(par->Tok[1], w_MLD)) par->Flowflag = MLD;
else if (match(par->Tok[1], w_SI)) par->Flowflag = LPS;
else return (201);
}
else if (match(par->Tok[0], w_PRESSURE))
{
if (n < 1) return (0);
else if (match(par->Tok[1], w_EXPONENT)) return -1;
else if (match(par->Tok[1], w_PSI)) par->Pressflag = PSI;
else if (match(par->Tok[1], w_KPA)) par->Pressflag = KPA;
else if (match(par->Tok[1], w_METERS)) par->Pressflag = METERS;
else return (201);
}
else if (match(par->Tok[0], w_HEADLOSS))
{
if (n < 1) return (0);
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[1], w_NONE))
qu->Qualflag = NONE;
else if (match(par->Tok[1], w_CHEM))
qu->Qualflag = CHEM;
else if (match(par->Tok[1], w_AGE))
qu->Qualflag = AGE;
else if (match(par->Tok[1], w_TRACE))
qu->Qualflag = TRACE;
else {
qu->Qualflag = CHEM;
strncpy(qu->ChemName, par->Tok[1], MAXID);
if (n >= 2)
strncpy(qu->ChemUnits, par->Tok[2], MAXID);
}
else if (match(par->Tok[0], w_QUALITY))
{
if (n < 1) return (0);
else if (match(par->Tok[1], w_NONE)) qu->Qualflag = NONE;
else if (match(par->Tok[1], w_CHEM)) qu->Qualflag = CHEM;
else if (match(par->Tok[1], w_AGE)) qu->Qualflag = AGE;
else if (match(par->Tok[1], w_TRACE)) qu->Qualflag = TRACE;
else
{
qu->Qualflag = CHEM;
strncpy(qu->ChemName, par->Tok[1], 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)
return (212);
if (n < 2) return (212);
strcpy(par->Tok[0], par->Tok[2]);
qu->TraceNode = findnode(net,par->Tok[2]);
if (qu->TraceNode == 0)
return (212);
if (qu->TraceNode == 0) 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)
return (0);
if (match(par->Tok[1], w_STOP))
hyd->ExtraIter = -1;
else if (match(par->Tok[1], w_CONTINUE)) {
if (n >= 2)
hyd->ExtraIter = atoi(par->Tok[2]);
else
hyd->ExtraIter = 0;
} else
return (201);
} else if (match(par->Tok[0], w_PATTERN)) /* Pattern option */
if (n < 1) return (0);
if (match(par->Tok[1], w_STOP)) hyd->ExtraIter = -1;
else if (match(par->Tok[1], w_CONTINUE))
{
if (n >= 2) hyd->ExtraIter = atoi(par->Tok[2]);
else hyd->ExtraIter = 0;
}
else return (201);
}
else if (match(par->Tok[0], w_PATTERN)) /* Pattern option */
{
if (n < 1)
return (0);
if (n < 1) 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))
return (213);
if (!getfloat(par->Tok[nvalue], &y)) return (213);
/* Check for WQ tolerance option (which can be 0) */
if (match(tok0, w_TOLERANCE)) {
if (y < 0.0)
return (213);
if (match(tok0, w_TOLERANCE))
{
if (y < 0.0) return (213);
qu->Ctol = y; /* Quality tolerance*/
return (0);
}
/* Check for Diffusivity option */
if (match(tok0, w_DIFFUSIVITY)) {
if (y < 0.0)
return (213);
if (match(tok0, w_DIFFUSIVITY))
{
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)
return (213);
if (y <= 0.0) return (213);
/* Assign value to specified option */
if (match(tok0, w_VISCOSITY))
hyd->Viscos = y; /* Viscosity */
else if (match(tok0, w_SPECGRAV))
hyd->SpGrav = y; /* Spec. gravity */
else if (match(tok0, w_TRIALS))
hyd->MaxIter = (int)y; /* Max. trials */
if (match(tok0, w_VISCOSITY)) hyd->Viscos = y; /* Viscosity */
else if (match(tok0, w_SPECGRAV)) 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))
hyd->Htol = y;
else if (match(tok0, w_QTOL))
hyd->Qtol = y;
else if (match(tok0, w_RQTOL)) {
if (y >= 1.0)
return (213);
else if (match(tok0, w_HTOL)) hyd->Htol = y;
else if (match(tok0, w_QTOL)) hyd->Qtol = y;
else if (match(tok0, w_RQTOL))
{
if (y >= 1.0) return (213);
hyd->RQtol = y;
} else if (match(tok0, w_CHECKFREQ))
hyd->CheckFreq = (int)y;
else if (match(tok0, w_MAXCHECK))
hyd->MaxCheck = (int)y;
else if (match(tok0, w_EMITTER))
hyd->Qexp = 1.0 / y;
else if (match(tok0, w_DEMAND))
hyd->Dmult = y;
else
return (201);
}
else if (match(tok0, w_CHECKFREQ)) hyd->CheckFreq = (int)y;
else if (match(tok0, w_MAXCHECK)) hyd->MaxCheck = (int)y;
else if (match(tok0, w_EMITTER)) hyd->Qexp = 1.0 / y;
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

4
src/output.c
View File

@@ -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) {

4
src/quality.c
View File

@@ -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

9
src/report.c
View File

@@ -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);

4
src/rules.c
View File

@@ -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,

13
src/text.h
View File

@@ -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"
@@ -420,7 +430,8 @@ AUTHOR: L. Rossman
/*** End of update ***/
#define FMT66 " maximum flow change = %.4f for Link %s"
#define FMT67 " maximum head error = %.4f for Link %s\n"
#define FMT67 " maximum flow change = %.4f for Node %s"
#define FMT68 " maximum head error = %.4f for Link %s\n"
/* -------------------- Energy Report Table ------------------- */
#define FMT71 "Energy Usage:"

24
src/types.h
View File

@@ -299,6 +299,11 @@ typedef enum {
MAX_ENERGY_STATS
} EnergyStats;
typedef enum {
DDA, // Demand Driven Analysis
PDA // Pressure Driven Analysis
} DemandModelType;
/*
------------------------------------------------------
Global Data Structures
@@ -760,7 +765,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 */
@@ -769,10 +775,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 */
@@ -788,7 +797,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 */
@@ -806,8 +816,10 @@ typedef struct {
Formflag; /* Hydraulic formula flag */
/* Info about hydraulic solution */
double relativeError;
int iterations;
double RelativeError;
double MaxHeadError;
double MaxFlowChange;
int Iterations;
/* Flag used to halt taking further time steps */
int Haltflag;

4
tools/before-test.cmd
View File

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

4
tools/run-nrtest.cmd
View File

@@ -18,12 +18,12 @@ setlocal
set NRTEST_SCRIPT_PATH=%~1
set TEST_SUITE_PATH=%~2
set BENCHMARK_VER=220dev1
set BENCHMARK_VER=220dev2
set NRTEST_EXECUTE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest execute
set TEST_APP_PATH=apps\epanet-%3.json
set TESTS=tests\examples tests\exeter tests\large tests\network_one tests\small tests\tanks tests\valves
set TESTS=tests\examples tests\exeter tests\large tests\network_one tests\press_depend tests\small tests\tanks tests\valves
set TEST_OUTPUT_PATH=benchmark\epanet-%3
set NRTEST_COMPARE_CMD=python %NRTEST_SCRIPT_PATH%\nrtest compare

8
win_build/WinSDK/Makefile.bat
View File

@@ -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 hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /link /DLL
cl -o epanet2.dll epanet.c hash.c hydraul.c hydcoeffs.c hydsolver.c hydstatus.c genmmd.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 hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.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 hydstatus.c genmmd.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 hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.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 hydstatus.c genmmd.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 hydcoeffs.c hydsolver.c inpfile.c input1.c input2.c input3.c mempool.c output.c quality.c report.c rules.c smatrix.c genmmd.c /I ..\include /I ..\run /I ..\src /link
cl -o epanet2.exe epanet.c ..\run\main.c hash.c hydraul.c hydcoeffs.c hydstatus.c genmmd.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

5
win_build/WinSDK/epanet2.def
View File

@@ -93,4 +93,7 @@ EXPORTS
ENdeletelink = _ENdeletelink@4
ENdeletenode = _ENdeletenode@4
ENsetlinktype = _ENsetlinktype@8
ENgetcurvetype = _ENgetcurvetype@8
ENgetdemandmodel = _ENgetdemandmodel@16
ENsetdemandmodel = _ENsetdemandmodel@16
ENgetcurvetype = _ENgetcurvetype@8