diff --git a/ReleaseNotes2_2.md b/ReleaseNotes2_2.md new file mode 100644 index 0000000..d09fbf8 --- /dev/null +++ b/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 *Df* is the full demand required, *Pmin* is the pressure below which demand is zero, *Preq* is the pressure required to deliver the full required demand and *Pexp* is an exponent. When *P < Pmin* demand is 0 and when *P > Preq* demand equals *Df*. + +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 *Pmin* | 0 +| REQUIRED PRESSURE | value for *Preq* | 0 +| PRESSURE EXPONENT | value for *Pexp* | 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. + - *P0* is allowed to equal to *Preq*. 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 *Pmin*. + +## 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. diff --git a/include/epanet2.h b/include/epanet2.h index 81b469e..bfc1905 100644 --- a/include/epanet2.h +++ b/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_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); @@ -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 diff --git a/src/enumstxt.h b/src/enumstxt.h index c3abbc3..6688a1a 100755 --- a/src/enumstxt.h +++ b/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, diff --git a/src/epanet.c b/src/epanet.c index f0f0ff8..d89327d 100644 --- a/src/epanet.c +++ b/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" -#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); } @@ -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 -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; } @@ -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 - HeadlossFormula = headloss formula flag - + ** UnitsType = flow units 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; - *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; } @@ -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); } diff --git a/src/funcs.h b/src/funcs.h index c8dfb05..e274ded 100755 --- a/src/funcs.h +++ b/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 */ diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 9b830d9..00b0b82 100644 --- a/src/hydcoeffs.c +++ b/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 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; } diff --git a/src/hydraul.c b/src/hydraul.c index 545ed37..3984ec0 100755 --- a/src/hydraul.c +++ b/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. */ diff --git a/src/hydsolver.c b/src/hydsolver.c index beeff1d..d08076a 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -22,7 +22,7 @@ The solver implements Todini's Global Gradient Algorithm. #include "types.h" #include "funcs.h" -// Error in network being hydraulically balanced +// Hydraulic balance error for network being analyzed typedef struct { double maxheaderror; double maxflowerror; @@ -33,22 +33,25 @@ typedef struct { } Hydbalance; // External functions -//int hydsolve(EN_Project *pr, int *iter, double *relerr); +//int hydsolve(EN_Project *pr, int *iter, double *relerr); +//void headlosscoeffs(EN_Project *pr); +//void matrixcoeffs(EN_Project *pr); + +extern int valvestatus(EN_Project *pr); //(see HYDSTATUS.C) +extern int linkstatus(EN_Project *pr); //(see HYDSTATUS.C) // Local functions static int badvalve(EN_Project *pr, int); -static int valvestatus(EN_Project *pr); -static int linkstatus(EN_Project *pr); -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); static int pswitch(EN_Project *pr); static double newflows(EN_Project *pr, Hydbalance *hbal); -static double emitflowchange(EN_Project *pr, int i); +static void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum, + double *dqsum); +static void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum, + double *dqsum); +static void newdemandflows(EN_Project *pr, Hydbalance *hbal, double *qsum, + double *dqsum); + static void checkhydbalance(EN_Project *pr, Hydbalance *hbal); static int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal); static void reporthydbal(EN_Project *pr, Hydbalance *hbal); @@ -75,122 +78,117 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr) ** is > 0 then future computed flow changes are only 60% of ** their full value. A complete status check on all links ** is made when convergence is achieved. If convergence is -** not achieved in hyd->MaxIter trials and hyd->ExtraIter > 0 then -** another hyd->ExtraIter trials are made with no status changes +** not achieved in MaxIter trials and ExtraIter > 0 then +** another ExtraIter trials are made with no status changes ** made to any links and a warning message is generated. ** ** This procedure calls linsolve() which appears in SMATRIX.C. **------------------------------------------------------------------- */ { - int i; /* Node index */ - int errcode = 0; /* Node causing solution error */ - int nextcheck; /* Next status check trial */ - int maxtrials; /* Max. trials for convergence */ - double newerr; /* New convergence error */ - int valveChange; /* Valve status change flag */ - int statChange; /* Non-valve status change flag */ - Hydbalance hydbal; /* Hydraulic balance errors */ + int i; // Node index + int errcode = 0; // Node causing solution error + int nextcheck; // Next status check trial + int maxtrials; // Max. trials for convergence + double newerr; // New convergence error + int valveChange; // Valve status change flag + int statChange; // Non-valve status change flag + Hydbalance hydbal; // Hydraulic balance errors EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; solver_t *sol = &hyd->solver; report_options_t *rep = &pr->report; - /* Initialize status checking & relaxation factor */ + // Initialize status checking & relaxation factor nextcheck = hyd->CheckFreq; hyd->RelaxFactor = 1.0; - /* Repeat iterations until convergence or trial limit is exceeded. */ - /* (hyd->ExtraIter used to increase trials in case of status cycling.) */ - if (pr->report.Statflag == FULL) { + // Repeat iterations until convergence or trial limit is exceeded. + // (ExtraIter used to increase trials in case of status cycling.) + if (pr->report.Statflag == FULL) + { writerelerr(pr, 0, 0); } maxtrials = hyd->MaxIter; - if (hyd->ExtraIter > 0) { + if (hyd->ExtraIter > 0) + { maxtrials += hyd->ExtraIter; } *iter = 1; - while (*iter <= maxtrials) { + while (*iter <= maxtrials) + { /* ** Compute coefficient matrices A & F and solve A*H = F ** where H = heads, A = Jacobian coeffs. derived from ** head loss gradients, & F = flow correction terms. ** Solution for H is returned in F from call to linsolve(). */ - for (i = 1; i <= net->Nlinks; i++) { - hlosscoeff(pr, i); - } + headlosscoeffs(pr); matrixcoeffs(pr); errcode = linsolve(pr, net->Njuncs); - /* Take action depending on error code */ - if (errcode < 0) { - break; /* Memory allocation problem */ - } - if (errcode > 0) { /* Ill-conditioning problem */ - /* If control valve causing problem, fix its status & continue, */ - /* otherwise end the iterations with no solution. */ - if (badvalve(pr, sol->Order[errcode])) { - continue; - } + // Quit if memory allocation error + if (errcode < 0) break; + + // Matrix ill-conditioning problem - if control valve causing problem, + // fix its status & continue, otherwise quit with no solution. + if (errcode > 0) + { + if (badvalve(pr, sol->Order[errcode])) continue; else break; } - /* Update current solution. */ - /* (Row[i] = row of solution matrix corresponding to node i). */ - for (i = 1; i <= net->Njuncs; i++) { - hyd->NodeHead[i] = sol->F[sol->Row[i]]; /* Update heads */ + // Update current solution. + // (Row[i] = row of solution matrix corresponding to node i) + for (i = 1; i <= net->Njuncs; i++) + { + hyd->NodeHead[i] = sol->F[sol->Row[i]]; // Update heads } - newerr = newflows(pr, &hydbal); /* Update flows */ + newerr = newflows(pr, &hydbal); // Update flows *relerr = newerr; - /* Write convergence error to status report if called for */ - if (rep->Statflag == FULL) { + // Write convergence error to status report if called for + if (rep->Statflag == FULL) + { writerelerr(pr, *iter, *relerr); } - /* Apply solution damping & check for change in valve status */ + // Apply solution damping & check for change in valve status hyd->RelaxFactor = 1.0; valveChange = FALSE; - if (hyd->DampLimit > 0.0) { - if (*relerr <= hyd->DampLimit) { + if (hyd->DampLimit > 0.0) + { + if (*relerr <= hyd->DampLimit) + { hyd->RelaxFactor = 0.6; valveChange = valvestatus(pr); } } - else { + else + { valveChange = valvestatus(pr); } - /* Check for convergence */ - if (hasconverged(pr, relerr, &hydbal)) { - /* We have convergence. Quit if we are into extra iterations. */ - if (*iter > hyd->MaxIter) { - break; - } + // Check for convergence + if (hasconverged(pr, relerr, &hydbal)) + { + // We have convergence - quit if we are into extra iterations + if (*iter > hyd->MaxIter) break; - /* Quit if no status changes occur. */ + // Quit if no status changes occur statChange = FALSE; - if (valveChange) { - statChange = TRUE; - } - if (linkstatus(pr)) { - statChange = TRUE; - } - if (pswitch(pr)) { - statChange = TRUE; - } - if (!statChange) { - break; - } + if (valveChange) statChange = TRUE; + if (linkstatus(pr)) statChange = TRUE; + if (pswitch(pr)) statChange = TRUE; + if (!statChange) break; - /* We have a status change so continue the iterations */ + // We have a status change so continue the iterations nextcheck = *iter + hyd->CheckFreq; } - /* No convergence yet. See if its time for a periodic status */ - /* check on pumps, CV's, and pipes connected to tanks. */ + // No convergence yet - see if its time for a periodic status + // check on pumps, CV's, and pipes connected to tank else if (*iter <= hyd->MaxCheck && *iter == nextcheck) { linkstatus(pr); @@ -199,20 +197,28 @@ int hydsolve(EN_Project *pr, int *iter, double *relerr) (*iter)++; } - /* Iterations ended. Report any errors. */ - if (errcode < 0) errcode = 101; /* Memory allocation error */ + // Iterations ended - report any errors. + if (errcode < 0) errcode = 101; // Memory allocation error else if (errcode > 0) { - writehyderr(pr, sol->Order[errcode]); /* Ill-conditioned eqns. error */ + writehyderr(pr, sol->Order[errcode]); // Ill-conditioned matrix error errcode = 110; } - /* Add any emitter flows to junction demands */ - for (i = 1; i <= net->Njuncs; i++) { - hyd->NodeDemand[i] += hyd->EmitterFlows[i]; + // Replace junction demands with total outflow values + for (i = 1; i <= net->Njuncs; i++) + { + hyd->NodeDemand[i] = hyd->DemandFlows[i] + hyd->EmitterFlows[i]; } + + // Save convergence info + hyd->RelativeError = *relerr; + hyd->MaxHeadError = hydbal.maxheaderror; + hyd->MaxFlowChange = hydbal.maxflowchange; + hyd->Iterations = *iter; + return(errcode); -} /* End of netsolve */ +} int badvalve(EN_Project *pr, int n) @@ -235,513 +241,44 @@ int badvalve(EN_Project *pr, int n) report_options_t *rep = &pr->report; time_options_t *time = &pr->time_options; Slink *link; + EN_LinkType t; - for (i = 1; i <= net->Nvalves; i++) { + for (i = 1; i <= net->Nvalves; i++) + { k = net->Valve[i].Link; link = &net->Link[k]; n1 = link->N1; n2 = link->N2; - if (n == n1 || n == n2) { - EN_LinkType t = link->Type; - if (t == EN_PRV || t == EN_PSV || t == EN_FCV) { - if (hyd->LinkStatus[k] == ACTIVE) { - if (rep->Statflag == FULL) { + if (n == n1 || n == n2) + { + t = link->Type; + if (t == EN_PRV || t == EN_PSV || t == EN_FCV) + { + if (hyd->LinkStatus[k] == ACTIVE) + { + if (rep->Statflag == FULL) + { sprintf(pr->Msg, FMT61, clocktime(rep->Atime, time->Htime), link->ID); writeline(pr, pr->Msg); } - if (link->Type == EN_FCV) { + if (link->Type == EN_FCV) + { hyd->LinkStatus[k] = XFCV; } - else { + else + { hyd->LinkStatus[k] = XPRESSURE; } - return(1); + return 1; } } - return(0); + return 0; } } - return(0); + return 0; } -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; - - for (i = 1; i <= net->Nvalves; i++) /* Examine each valve */ - { - k = net->Valve[i].Link; /* Link index of valve */ - link = &net->Link[k]; - if (hyd->LinkSetting[k] == MISSING) { - continue; /* Valve status fixed */ - } - n1 = link->N1; /* Start & end nodes */ - n2 = link->N2; - status = hyd->LinkStatus[k]; /* Save current status */ - -// if (s != CLOSED /* No change if flow is */ -// && ABS(hyd->LinkFlows[k]) < hyd->Qtol) continue; /* negligible. */ - - switch (link->Type) /* Evaluate new status: */ - { - 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; - -//// FCV status checks moved back into the linkstatus() function //// -// case FCV: S[k] = fcvstatus(k,s,hyd->NodeHead[n1],hyd->NodeHead[n2]); -// break; - - default: - continue; - } - - /*** Updated 9/7/00 ***/ - /* Do not reset flow in valve if its status changes. */ - /* This strategy improves convergence. */ - - /* Check for 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() */ - - -/*** Updated 9/7/00 ***/ -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 */ - StatType status; /* Current status */ - - EN_Network *net = &pr->network; - hydraulics_t *hyd = &pr->hydraulics; - report_options_t *rep = &pr->report; - Slink *link; - - /* Examine each Slink */ - 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 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]); - } - - //if (S[k] <= CLOSED) hyd->LinkFlows[k] = QZERO; - //else setlinkflow(k, dh); - } - } - return(change); -} /* End of linkstatus */ - - -StatType cvstatus(EN_Project *pr, StatType s, double dh, double q) -/* -**-------------------------------------------------- -** Input: s = current status -** dh = headloss -** q = flow -** Output: returns new link status -** Purpose: updates status of a check valve. -**-------------------------------------------------- -*/ -{ - 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 -** 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; - - /* Prevent reverse flow through pump */ - p = findpump(net, k); - if (net->Pump[p].Ptype == CONST_HP) { - hmax = BIG; - } - else { - hmax = SQR(hyd->LinkSetting[k]) * net->Pump[p].Hmax; - } - if (dh > hmax + hyd->Htol) { - return(XHEAD); - } - - /*** Flow higher than pump curve no longer results in a status change ***/ - /* Check if pump cannot deliver flow */ - //if (hyd->LinkFlows[k] > K[k]*Pump[p].Qmax + hyd->Qtol) return(XFLOW); - 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; /* New valve status */ - double hml; /* Minor headloss */ - hydraulics_t *hyd = &pr->hydraulics; - - double htol = hyd->Htol; - Slink *link = &pr->network.Link[k]; - - status = s; - if (hyd->LinkSetting[k] == MISSING) - return(status); /* Status fixed by user */ - hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ - - /*** Status rules below have changed. ***/ - - 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; /* New valve status */ - double hml; /* Minor headloss */ - hydraulics_t *hyd = &pr->hydraulics; - - double htol = hyd->Htol; - - Slink *link = &pr->network.Link[k]; - - status = s; - if (hyd->LinkSetting[k] == MISSING) { - return(status); /* Status fixed by user */ - } - hml = link->Km*SQR(hyd->LinkFlows[k]); /* Head loss when open */ - -/*** Status rules below have changed. ***/ - - 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); -} - - -/*** Updated 9/7/00 ***/ -/*** Updated 11/19/01 ***/ -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]; - - /* Make node n1 be the tank */ - 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; - } - h = hyd->NodeHead[n1] - hyd->NodeHead[n2]; - tank = &net->Tank[i]; - /* Skip reservoirs & closed links */ - if (tank->A == 0.0 || hyd->LinkStatus[k] <= CLOSED) { - return; - } - - /* If tank 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 */ - /* (i.e., an open outflow check valve would close) */ - else if (cvstatus(pr, OPEN, h, q) == CLOSED) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } - - /* If tank 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 */ - /* (i.e., a closed outflow check valve would open) */ - else if (cvstatus(pr, CLOSED, h, q) == OPEN) { - hyd->LinkStatus[k] = TEMPCLOSED; - } - } -} /* End of tankstatus */ - - int pswitch(EN_Project *pr) /* **-------------------------------------------------------------- @@ -752,98 +289,133 @@ int pswitch(EN_Project *pr) **-------------------------------------------------------------- */ { - int i, /* Control statement index */ - k, /* Link being controlled */ - n, /* Node controlling Slink */ - reset, /* Flag on control conditions */ - change, /* Flag for status or setting change */ - anychange = 0; /* Flag for 1 or more changes */ - char s; /* Current link status */ + int i, // Control statement index + k, // Index of link being controlled + n, // Node controlling link k + reset, // Flag on control conditions + change, // Flag for status or setting change + anychange = 0; // Flag for 1 or more control actions + char s; // Current link status EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; report_options_t *rep = &pr->report; Slink *link; - /* Check each control statement */ + // Check each control statement for (i = 1; i <= net->Ncontrols; i++) { reset = 0; - if ((k = net->Control[i].Link) <= 0) { - continue; - } + k = net->Control[i].Link; + if (k <= 0) continue; - /* Determine if control based on a junction, not a tank */ - if ((n = net->Control[i].Node) > 0 && n <= net->Njuncs) { - /* Determine if control conditions are satisfied */ - if (net->Control[i].Type == LOWLEVEL - && hyd->NodeHead[n] <= net->Control[i].Grade + hyd->Htol) { + // Determine if control based on a junction, not a tank + n = net->Control[i].Node; + if (n > 0 && n <= net->Njuncs) + { + // Determine if control conditions are satisfied + if (net->Control[i].Type == LOWLEVEL && + hyd->NodeHead[n] <= net->Control[i].Grade + hyd->Htol) + { reset = 1; } - if (net->Control[i].Type == HILEVEL - && hyd->NodeHead[n] >= net->Control[i].Grade - hyd->Htol) { + if (net->Control[i].Type == HILEVEL && + hyd->NodeHead[n] >= net->Control[i].Grade - hyd->Htol) + { reset = 1; } } - /* Determine if control forces a status or setting change */ + // Determine if control forces a status or setting change if (reset == 1) { link = &net->Link[k]; change = 0; s = hyd->LinkStatus[k]; - if (link->Type == EN_PIPE) { - if (s != net->Control[i].Status) { - change = 1; - } + if (link->Type == EN_PIPE) + { + if (s != net->Control[i].Status) change = 1; } - if (link->Type == EN_PUMP) { - if (hyd->LinkSetting[k] != net->Control[i].Setting) { - change = 1; - } + if (link->Type == EN_PUMP) + { + if (hyd->LinkSetting[k] != net->Control[i].Setting) change = 1; } - if (link->Type >= EN_PRV) { - if (hyd->LinkSetting[k] != net->Control[i].Setting) { - change = 1; - } - else if (hyd->LinkSetting[k] == MISSING && s != net->Control[i].Status) { + if (link->Type >= EN_PRV) + { + if (hyd->LinkSetting[k] != net->Control[i].Setting) change = 1; + else if (hyd->LinkSetting[k] == MISSING && s != net->Control[i].Status) + { change = 1; } } - /* If a change occurs, update status & setting */ - if (change) { + // If a change occurs, update status & setting + if (change) + { hyd->LinkStatus[k] = net->Control[i].Status; - if (link->Type > EN_PIPE) { + if (link->Type > EN_PIPE) + { hyd->LinkSetting[k] = net->Control[i].Setting; } - if (rep->Statflag == FULL) { + if (rep->Statflag == FULL) + { writestatchange(pr, k, s, hyd->LinkStatus[k]); } - - /* Re-set flow if status has changed */ - // if (S[k] != s) initlinkflow(k, S[k], K[k]); anychange = 1; } } } return(anychange); -} /* End of pswitch */ +} double newflows(EN_Project *pr, Hydbalance *hbal) /* **---------------------------------------------------------------- -** Input: none +** Input: hbal = ptr. to hydraulic balance information ** Output: returns solution convergence error +** Purpose: updates link, emitter & demand flows after new +** nodal heads are computed. +**---------------------------------------------------------------- +*/ +{ + double dqsum, // Network flow change + qsum; // Network total flow + + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + + // Initialize sum of flows & corrections + qsum = 0.0; + dqsum = 0.0; + hbal->maxflowchange = 0.0; + hbal->maxflowlink = 1; + hbal->maxflownode = -1; + + // Update flows in all real and virtual links + newlinkflows(pr, hbal, &qsum, &dqsum); + newemitterflows(pr, hbal, &qsum, &dqsum); + newdemandflows(pr, hbal, &qsum, &dqsum); + + // Return ratio of total flow corrections to total flow + if (qsum > hyd->Hacc) return (dqsum / qsum); + else return dqsum; +} + + +void newlinkflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) +/* +**---------------------------------------------------------------- +** Input: hbal = ptr. to hydraulic balance information +** qsum = sum of current system flows +** dqsum = sum of system flow changes +** Output: updates hbal, qsum and dqsum ** Purpose: updates link flows after new nodal heads computed **---------------------------------------------------------------- */ { double dh, /* Link head loss */ dq; /* Link flow change */ - double dqsum, /* Network flow change */ - qsum; /* Network total flow */ int k, n, n1, n2; EN_Network *net = &pr->network; @@ -851,114 +423,152 @@ double newflows(EN_Project *pr, Hydbalance *hbal) solver_t *sol = &hyd->solver; Slink *link; - /* Initialize net inflows (i.e., demands) at tanks */ - for (n = net->Njuncs + 1; n <= net->Nnodes; n++) { + // Initialize net inflows (i.e., demands) at fixed grade nodes + for (n = net->Njuncs + 1; n <= net->Nnodes; n++) + { hyd->NodeDemand[n] = 0.0; } - /* Initialize sum of flows & corrections */ - qsum = 0.0; - dqsum = 0.0; - - hbal->maxflowchange = 0.0; - hbal->maxflowlink = 1; - - /* Update flows in all links */ + // Examine each link for (k = 1; k <= net->Nlinks; k++) { + // Get link and its end nodes link = &net->Link[k]; - /* - ** Apply flow update formula: - ** dq = Y - P*(new head loss) - ** P = 1/(dh/dq) - ** Y = P*(head loss based on current flow) - ** where P & Y were computed in newcoeffs(). - */ - n1 = link->N1; n2 = link->N2; + + // Apply flow update formula: + // dq = Y - P * (new head loss) + // P = 1 / (previous head loss gradient) + // Y = P * (previous head loss) + // where P & Y were computed in hlosscoeff() in hydcoeffs.c + dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; dq = sol->Y[k] - sol->P[k] * dh; - /* Adjust flow change by the relaxation factor */ + // Adjust flow change by the relaxation factor dq *= hyd->RelaxFactor; - /* Prevent flow in constant HP pumps from going negative */ - if (link->Type == EN_PUMP) { + // Prevent flow in constant HP pumps from going negative + if (link->Type == EN_PUMP) + { n = findpump(net, k); - if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlows[k]) { + if (net->Pump[n].Ptype == CONST_HP && dq > hyd->LinkFlows[k]) + { dq = hyd->LinkFlows[k] / 2.0; } } + + // Update link flow and system flow summation hyd->LinkFlows[k] -= dq; + *qsum += ABS(hyd->LinkFlows[k]); + *dqsum += ABS(dq); - /* Update sum of absolute flows & flow corrections */ - qsum += ABS(hyd->LinkFlows[k]); - dqsum += ABS(dq); - - /* Update identity of link with max. flow change */ - if (ABS(dq) > hbal->maxflowchange) { + // Update identity of element with max. flow change + if (ABS(dq) > hbal->maxflowchange) + { hbal->maxflowchange = ABS(dq); hbal->maxflowlink = k; - } + hbal->maxflownode = -1; + } - /* Update net flows to tanks */ + // Update net flows to fixed grade nodes if (hyd->LinkStatus[k] > CLOSED) { - if (n1 > net->Njuncs) { - hyd->NodeDemand[n1] -= hyd->LinkFlows[k]; - } - if (n2 > net->Njuncs) { - hyd->NodeDemand[n2] += hyd->LinkFlows[k]; - } + if (n1 > net->Njuncs) hyd->NodeDemand[n1] -= hyd->LinkFlows[k]; + if (n2 > net->Njuncs) hyd->NodeDemand[n2] += hyd->LinkFlows[k]; } } - - /* Update emitter flows */ - for (k = 1; k <= net->Njuncs; k++) - { - if (net->Node[k].Ke == 0.0) { - continue; - } - dq = emitflowchange(pr, k); - hyd->EmitterFlows[k] -= dq; - qsum += ABS(hyd->EmitterFlows[k]); - dqsum += ABS(dq); - } - - /* Return ratio of total flow corrections to total flow */ - if (qsum > hyd->Hacc) { - return(dqsum / qsum); - } - else { - return(dqsum); - } -} /* End of newflows */ +} -double emitflowchange(EN_Project *pr, int i) +void newemitterflows(EN_Project *pr, Hydbalance *hbal, double *qsum, + double *dqsum) /* -**-------------------------------------------------------------- -** Input: i = node index -** Output: returns change in flow at an emitter node -** Purpose: computes flow change at an emitter node -**-------------------------------------------------------------- +**---------------------------------------------------------------- +** Input: hbal = ptr. to hydraulic balance information +** qsum = sum of current system flows +** dqsum = sum of system flow changes +** Output: updates hbal, qsum and dqsum +** Purpose: updates nodal emitter flows after new nodal heads computed +**---------------------------------------------------------------- */ { - double ke, p; + double dq; + int k; + EN_Network *net = &pr->network; hydraulics_t *hyd = &pr->hydraulics; - EN_Network *n = &pr->network; - Snode *node = &n->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; + // Examine each network junction + for (k = 1; k <= net->Njuncs; k++) + { + // Skip junction if it does not have an emitter + if (net->Node[k].Ke == 0.0) continue; + + // Find emitter flow change (see hydcoeffs.c) + dq = emitflowchange(pr, k); + hyd->EmitterFlows[k] -= dq; + + // Update system flow summation + *qsum += ABS(hyd->EmitterFlows[k]); + *dqsum += ABS(dq); + + // Update identity of element with max. flow change + if (ABS(dq) > hbal->maxflowchange) + { + hbal->maxflowchange = ABS(dq); + hbal->maxflownode = k; + hbal->maxflowlink = -1; + } } - else { - p = 1.0 / p; +} + + +void newdemandflows(EN_Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) +/* +**---------------------------------------------------------------- +** Input: hbal = ptr. to hydraulic balance information +** qsum = sum of current system flows +** dqsum = sum of system flow changes +** Output: updates hbal, qsum and dqsum +** Purpose: updates nodal pressure dependent demand flows after +** new nodal heads computed +**---------------------------------------------------------------- +*/ +{ + double dp, // pressure range over which demand can vary (ft) + dq, // change in demand flow (cfs) + n; // exponent in head loss v. demand function + int k; + EN_Network *net = &pr->network; + hydraulics_t *hyd = &pr->hydraulics; + + // Get demand function parameters + if (hyd->DemandModel == DDA) return; + demandparams(pr, &dp, &n); + + // Examine each junction + for (k = 1; k <= net->Njuncs; k++) + { + // Skip junctions with no positive demand + if (hyd->NodeDemand[k] <= 0.0) continue; + + // Find change in demand flow (see hydcoeffs.c) + dq = demandflowchange(pr, k, dp, n); + hyd->DemandFlows[k] -= dq; + + // Update system flow summation + *qsum += ABS(hyd->DemandFlows[k]); + *dqsum += ABS(dq); + + // Update identity of element with max. flow change + if (ABS(dq) > hbal->maxflowchange) + { + hbal->maxflowchange = ABS(dq); + hbal->maxflownode = k; + hbal->maxflowlink = -1; + } } - return(hyd->EmitterFlows[i] / hyd->Qexp - p * (hyd->NodeHead[i] - node->El)); } @@ -980,9 +590,10 @@ void checkhydbalance(EN_Project *pr, Hydbalance *hbal) Slink *link; hbal->maxheaderror = 0.0; hbal->maxheadlink = 1; - for (k = 1; k <= net->Nlinks; k++) { + headlosscoeffs(pr); + for (k = 1; k <= net->Nlinks; k++) + { if (hyd->LinkStatus[k] <= CLOSED) continue; - hlosscoeff(pr, k); if (sol->P[k] == 0.0) continue; link = &net->Link[k]; n1 = link->N1; @@ -990,7 +601,8 @@ void checkhydbalance(EN_Project *pr, Hydbalance *hbal) dh = hyd->NodeHead[n1] - hyd->NodeHead[n2]; headloss = sol->Y[k] / sol->P[k]; headerror = ABS(dh - headloss); - if (headerror > hbal->maxheaderror) { + if (headerror > hbal->maxheaderror) + { hbal->maxheaderror = headerror; hbal->maxheadlink = k; } @@ -1013,13 +625,14 @@ int hasconverged(EN_Project *pr, double *relerr, Hydbalance *hbal) if (*relerr > hyd->Hacc) return 0; checkhydbalance(pr, hbal); - if (pr->report.Statflag == FULL) { + if (pr->report.Statflag == FULL) + { reporthydbal(pr, hbal); } if (hyd->HeadErrorLimit > 0.0 && - hbal->maxheaderror > hyd->HeadErrorLimit) return 0; + hbal->maxheaderror > hyd->HeadErrorLimit) return 0; if (hyd->FlowChangeLimit > 0.0 && - hbal->maxflowchange > hyd->FlowChangeLimit) return 0; + hbal->maxflowchange > hyd->FlowChangeLimit) return 0; return 1; } @@ -1038,11 +651,13 @@ void reporthydbal(EN_Project *pr, Hydbalance *hbal) double herror = hbal->maxheaderror * pr->Ucf[HEAD]; int qlink = hbal->maxflowlink; int hlink = hbal->maxheadlink; - if (qlink >= 1) { + if (qlink >= 1) + { sprintf(pr->Msg, FMT66, qchange, pr->network.Link[qlink].ID); writeline(pr, pr->Msg); } - if (hlink >= 1) { + if (hlink >= 1) + { sprintf(pr->Msg, FMT67, herror, pr->network.Link[hlink].ID); writeline(pr, pr->Msg); } diff --git a/src/hydstatus.c b/src/hydstatus.c new file mode 100644 index 0000000..f5041c6 --- /dev/null +++ b/src/hydstatus.c @@ -0,0 +1,462 @@ +/* +********************************************************************* + +HYDSTATUS.C -- Hydraulic status updating for the EPANET Program + +******************************************************************* +*/ + +#include +#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; + } + } +} diff --git a/src/inpfile.c b/src/inpfile.c index 8bf894d..f5ef6b8 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -26,14 +26,14 @@ data describing a piping network to a file in EPANET's text format. #else #include #endif + #include "hash.h" #include "text.h" #include "types.h" -#include "epanet2.h" #include "funcs.h" #include -#define EXTERN extern -#include "vars.h" +//#define EXTERN extern +//#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]); - } - if (hyd->FlowChangeLimit > 0.0) { - fprintf(f, "\n FLOWCHANGE %-.8f", hyd->FlowChangeLimit * pr->Ucf[FLOW]); + fprintf(f, "\n HEADERROR %-.8f", hyd->HeadErrorLimit * pr->Ucf[HEAD]); + 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 */ diff --git a/src/input1.c b/src/input1.c index b59096a..46da66a 100644 --- a/src/input1.c +++ b/src/input1.c @@ -33,8 +33,6 @@ AUTHOR: L. Rossman #include "epanet2.h" #include "funcs.h" #include -#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]; diff --git a/src/input2.c b/src/input2.c index b0a1b13..e4836b5 100644 --- a/src/input2.c +++ b/src/input2.c @@ -37,8 +37,6 @@ The following utility functions are all called from INPUT3.C #include "epanet2.h" #include "funcs.h" #include -#define EXTERN extern -#include "vars.h" #define MAXERRS 10 /* Max. input errors reported */ diff --git a/src/input3.c b/src/input3.c index 9c66f63..58f9ce2 100644 --- a/src/input3.c +++ b/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 -#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 diff --git a/src/output.c b/src/output.c index 39687cf..d9ec660 100644 --- a/src/output.c +++ b/src/output.c @@ -19,14 +19,12 @@ AUTHOR: L. Rossman #else #include #endif -#include "epanet2.h" + #include "funcs.h" #include "text.h" #include "types.h" #include -#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) { diff --git a/src/quality.c b/src/quality.c index 03e29b6..a7286db 100644 --- a/src/quality.c +++ b/src/quality.c @@ -54,15 +54,13 @@ AUTHOR: L. Rossman #else #include #endif + #include "hash.h" #include "text.h" #include "types.h" -#include "epanet2.h" #include "funcs.h" #include -#define EXTERN extern #include "mempool.h" -#include "vars.h" /* ** Macros to identify upstream & downstream nodes of a link diff --git a/src/report.c b/src/report.c index 8d2101e..80eb45a 100644 --- a/src/report.c +++ b/src/report.c @@ -32,15 +32,13 @@ formatted string S to the report file. #else #include #endif -#include "epanet2.h" + #include "funcs.h" #include "hash.h" #include "text.h" #include "types.h" #include #include -#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); diff --git a/src/rules.c b/src/rules.c index 39d6a6a..d5da392 100644 --- a/src/rules.c +++ b/src/rules.c @@ -29,13 +29,11 @@ AUTHOR: L. Rossman #else #include #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, diff --git a/src/text.h b/src/text.h index 85159a9..4684fc8 100755 --- a/src/text.h +++ b/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" diff --git a/src/types.h b/src/types.h index a039d8c..f3ac921 100755 --- a/src/types.h +++ b/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; - int iterations; + double RelativeError; + double MaxHeadError; + double MaxFlowChange; + int Iterations; /* Flag used to halt taking further time steps */ int Haltflag;