From 80f9acfe4d6626b6d8ae8acf21187c752485967c Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Tue, 4 Feb 2020 10:01:23 -0500 Subject: [PATCH 1/8] Network building enhancements --- ReleaseNotes2_3.md | 11 ++++++++ include/epanet2.bas | 5 +++- include/epanet2.def | 6 ++-- include/epanet2.h | 6 +++- include/epanet2.pas | 25 +++++++++++------ include/epanet2.vb | 5 +++- include/epanet2_2.h | 24 ++++++++++++++-- include/epanet2_enums.h | 5 ++-- src/epanet.c | 62 ++++++++++++++++++++++++++++++++++++++++- src/epanet2.c | 12 +++++++- src/errors.dat | 2 +- src/funcs.h | 4 ++- src/hydraul.c | 11 ++------ src/input2.c | 59 +-------------------------------------- src/project.c | 34 ++++++++++++++++++++-- src/quality.c | 10 ++++++- 16 files changed, 188 insertions(+), 93 deletions(-) create mode 100644 ReleaseNotes2_3.md diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md new file mode 100644 index 0000000..1f922fc --- /dev/null +++ b/ReleaseNotes2_3.md @@ -0,0 +1,11 @@ +> +## Release Notes for EPANET 2.3 + +This document describes the changes and updates that have been made in version 2.3 of EPANET. + + - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files to be opened by the toolkit. + - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). + - A `EN_setvertex` function was added to allow API clients to change the coordinates of a link's vertex. + - The index of a General Purpose Valve's (GPV's) head loss curve was added to the list of editable Link Properties using the symbolic constant name `EN_GPV_CURVE`. + - The `EN_getlinkvalue` and `EN_setlinkvalue` functions were updated to get and set the value of `EN_GPV_CURVE`. + diff --git a/include/epanet2.bas b/include/epanet2.bas index 4b9df87..7b40f36 100644 --- a/include/epanet2.bas +++ b/include/epanet2.bas @@ -5,7 +5,7 @@ Attribute VB_Name = "Module1" 'Declarations of functions in the EPANET PROGRAMMERs TOOLKIT '(EPANET2.DLL) -'Last updated on 11/04/2019 +'Last updated on 02/01/2020 ' These are codes used by the DLL functions Public Const EN_ELEVATION = 0 ' Node parameters @@ -62,6 +62,7 @@ Public Const EN_PUMP_HCURVE = 19 Public Const EN_PUMP_ECURVE = 20 Public Const EN_PUMP_ECOST = 21 Public Const EN_PUMP_EPAT = 22 +Public Const EN_GPV_CURVE = 23 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 @@ -340,6 +341,7 @@ Public Const EN_MISSING As Double = -1.0E10 Declare Function ENsetpipedata Lib "epanet2.dll" (ByVal index As Long, ByVal length As Single, ByVal diam As Single, ByVal rough As Single, ByVal mloss As Single) As Long Declare Function ENgetvertexcount Lib "epanet2.dll" (ByVal index As Long, count As Long) As Long Declare Function ENgetvertex Lib "epanet2.dll" (ByVal index As Long, ByVal vertex As Long, x As Double, y As Double) As Long + Declare Function ENsetvertex Lib "epanet2.dll" (ByVal index As Long, ByVal vertex As Long, ByVal x As Double, ByVal y As Double) As Long Declare Function ENsetvertices Lib "epanet2.dll" (ByVal index As Long, xCoords As Any, yCoords As Any, ByVal count As Long) As Long 'Pump Functions @@ -367,6 +369,7 @@ Public Const EN_MISSING As Double = -1.0E10 Declare Function ENsetcurveid Lib "epanet2.dll" (ByVal index As Long, ByVal newid As String) As Long Declare Function ENgetcurvelen Lib "epanet2.dll" (ByVal index As Long, len_ As Long) As Long Declare Function ENgetcurvetype Lib "epanet2.dll" (ByVal index As Long, type_ As Long) As Long + Declare Function ENsetcurvetype Lib "epanet2.dll" (ByVal index As Long, ByVal type_ As Long) As Long Declare Function ENgetcurvevalue Lib "epanet2.dll" (ByVal curveIndex As Long, ByVal pointIndex As Long, x As Single, y As Single) As Long Declare Function ENsetcurvevalue Lib "epanet2.dll" (ByVal curveIndex As Long, ByVal pointIndex As Long, ByVal x As Single, ByVal y As Single) As Long Declare Function ENgetcurve Lib "epanet2.dll" (ByVal index As Long, ByVal id As String, nPoints As Long, xValues As Any, yValues As Any) As Long diff --git a/include/epanet2.def b/include/epanet2.def index 9bf41b4..66b87ed 100644 --- a/include/epanet2.def +++ b/include/epanet2.def @@ -91,7 +91,8 @@ EXPORTS ENsetcontrol = _ENsetcontrol@24 ENsetcoord = _ENsetcoord@20 ENsetcurve = _ENsetcurve@16 - ENsetcurveid = _ENsetcurveid@8 + ENsetcurveid = _ENsetcurveid@8 + ENsetcurvetype = _ENsetcurvetype@8 ENsetcurvevalue = _ENsetcurvevalue@16 ENsetdemandmodel = _ENsetdemandmodel@16 ENsetdemandname = _ENsetdemandname@12 @@ -122,7 +123,8 @@ EXPORTS ENsettankdata = _ENsettankdata@32 ENsetthenaction = _ENsetthenaction@20 ENsettimeparam = _ENsettimeparam@8 - ENsettitle = _ENsettitle@12 + ENsettitle = _ENsettitle@12 + ENsetvertex = _ENsetvertex@24 ENsetvertices = _ENsetvertices@16 ENsolveH = _ENsolveH@0 ENsolveQ = _ENsolveQ@0 diff --git a/include/epanet2.h b/include/epanet2.h index 19c9c67..3779d66 100644 --- a/include/epanet2.h +++ b/include/epanet2.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 10/29/2019 + Last Updated: 02/01/2020 ****************************************************************************** */ @@ -290,6 +290,8 @@ extern "C" { int DLLEXPORT ENgetvertex(int index, int vertex, double *x, double *y); + int DLLEXPORT ENsetvertex(int index, int vertex, double x, double y); + int DLLEXPORT ENsetvertices(int index, double *x, double *y, int count); /******************************************************************** @@ -349,6 +351,8 @@ extern "C" { int DLLEXPORT ENgetcurvelen(int index, int *len); int DLLEXPORT ENgetcurvetype(int index, int *type); + + int DLLEXPORT ENsetcurvetype(int index, int type); int DLLEXPORT ENgetcurvevalue(int curveIndex, int pointIndex, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y); diff --git a/include/epanet2.pas b/include/epanet2.pas index a532134..6eaf7bb 100644 --- a/include/epanet2.pas +++ b/include/epanet2.pas @@ -3,7 +3,7 @@ unit epanet2; { Declarations of imported procedures from the EPANET PROGRAMMERs TOOLKIT } { (EPANET2.DLL) } -{Last updated on 11/12/19} +{Last updated on 02/01/2020} interface @@ -67,6 +67,7 @@ const EN_PUMP_ECURVE = 20; EN_PUMP_ECOST = 21; EN_PUMP_EPAT = 22; + EN_GPV_CURVE = 23; EN_DURATION = 0; { Time parameters } EN_HYDSTEP = 1; @@ -253,8 +254,12 @@ const EN_R_IS_OPEN = 1; { Rule-based control link status } EN_R_IS_CLOSED = 2; EN_R_IS_ACTIVE = 3; - + +{$ifdef WINDOWS} EpanetLib = 'epanet2.dll'; +{$else} + EpanetLib = 'libepanet2.so'; +{$endif} {Project Functions} function ENepanet(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar; F4: Pointer): Integer; stdcall; external EpanetLib; @@ -273,8 +278,8 @@ const function ENsaveH: Integer; stdcall; external EpanetLib; function ENopenH: Integer; stdcall; external EpanetLib; function ENinitH(SaveFlag: Integer): Integer; stdcall; external EpanetLib; - function ENrunH(var T: LongInt): Integer; stdcall; external EpanetLib; - function ENnextH(var Tstep: LongInt): Integer; stdcall; external EpanetLib; + function ENrunH(var T: Integer): Integer; stdcall; external EpanetLib; + function ENnextH(var Tstep: Integer): Integer; stdcall; external EpanetLib; function ENcloseH: Integer; stdcall; external EpanetLib; function ENsavehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; function ENusehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; @@ -283,9 +288,9 @@ const function ENsolveQ: Integer; stdcall; external EpanetLib; function ENopenQ: Integer; stdcall; external EpanetLib; function ENinitQ(SaveFlag: Integer): Integer; stdcall; external EpanetLib; - function ENrunQ(var T: LongInt): Integer; stdcall; external EpanetLib; - function ENnextQ(var Tstep: LongInt): Integer; stdcall; external EpanetLib; - function ENstepQ(var Tleft: LongInt): Integer; stdcall; external EpanetLib; + function ENrunQ(var T: Integer): Integer; stdcall; external EpanetLib; + function ENnextQ(var Tstep: Integer): Integer; stdcall; external EpanetLib; + function ENstepQ(var Tleft: Integer): Integer; stdcall; external EpanetLib; function ENcloseQ: Integer; stdcall; external EpanetLib; {Reporting Functions} @@ -306,8 +311,8 @@ const function ENsetoption(Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; function ENgetflowunits(var Code: Integer): Integer; stdcall; external EpanetLib; function ENsetflowunits(Code: Integer): Integer; stdcall; external EpanetLib; - function ENgettimeparam(Code: Integer; var Value: LongInt): Integer; stdcall; external EpanetLib; - function ENsettimeparam(Code: Integer; Value: LongInt): Integer; stdcall; external EpanetLib; + function ENgettimeparam(Code: Integer; var Value: Integer): Integer; stdcall; external EpanetLib; + function ENsettimeparam(Code: Integer; Value: Integer): Integer; stdcall; external EpanetLib; function ENgetqualinfo(var QualType: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; var TraceNode: Integer): Integer; stdcall; external EpanetLib; function ENgetqualtype(var QualCode: Integer; var TraceNode: Integer): Integer; stdcall; external EpanetLib; function ENsetqualtype(QualCode: Integer; ChemName: PAnsiChar; ChemUnits: PAnsiChar; TraceNodeID: PAnsiChar): Integer; stdcall; external EpanetLib; @@ -356,6 +361,7 @@ const function ENgetvertexcount(Index: Integer; var Count: Integer): Integer; stdcall; external EpanetLib; function ENgetvertex(Index: Integer; Vertex: Integer; var X: Double; var Y: Double): Integer; stdcall; external EpanetLib; + function ENsetvertex(Index: Integer; Vertex: Integer; X: Double; Y: Double): Integer; stdcall; external EpanetLib; function ENsetvertices(Index: Integer; var X: Double; var Y: Double; Count: Integer): Integer; stdcall; external EpanetLib; {Pump Functions} @@ -383,6 +389,7 @@ const function ENsetcurveid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; function ENgetcurvelen(Index: Integer; var Len: Integer): Integer; stdcall; external EpanetLib; function ENgetcurvetype(Index: Integer; var CurveType: Integer): Integer; stdcall; external EpanetLib; + function ENsetcurvetype(Index: Integer; CurveType: Integer): Integer; stdcall; external EpanetLib; function ENgetcurvevalue(CurveIndex: Integer; PointIndex: Integer; var X: Single; var Y: Single): Integer; stdcall; external EpanetLib; function ENsetcurvevalue(CurveIndex: Integer; PointIndex: Integer; X: Single; Y: Single): Integer; stdcall; external EpanetLib; function ENgetcurve(Index: Integer; ID: PAnsiChar; var N: Integer; var X: Single; var Y: Single): Integer; stdcall; external EpanetLib; diff --git a/include/epanet2.vb b/include/epanet2.vb index bfe65c2..5d721bc 100644 --- a/include/epanet2.vb +++ b/include/epanet2.vb @@ -4,7 +4,7 @@ 'Declarations of functions in the EPANET PROGRAMMERs TOOLKIT '(EPANET2.DLL) for use with VB.Net. -'Last updated on 11/04/2019 +'Last updated on 02/01/2020 Imports System.Runtime.InteropServices Imports System.Text @@ -67,6 +67,7 @@ Public Const EN_PUMP_HCURVE = 19 Public Const EN_PUMP_ECURVE = 20 Public Const EN_PUMP_ECOST = 21 Public Const EN_PUMP_EPAT = 22 +Public Const EN_GPV_CURVE = 23 Public Const EN_DURATION = 0 ' Time parameters Public Const EN_HYDSTEP = 1 @@ -345,6 +346,7 @@ Public Const EN_MISSING As Double = -1.0E10 Declare Function ENsetpipedata Lib "epanet2.dll" (ByVal index As Int32, ByVal length As Single, ByVal diam As Single, ByVal rough As Single, ByVal mloss As Single) As Int32 Declare Function ENgetvertexcount Lib "epanet2.dll" (ByVal index As Int32, count As Int32) As Int32 Declare Function ENgetvertex Lib "epanet2.dll" (ByVal index As Int32, ByVal vertex As Int32, x As Double, y As Double) As Int32 + Declare Function ENsetvertex Lib "epanet2.dll" (ByVal index As Int32, ByVal vertex As Int32, ByVal x As Double, ByVal y As Double) As Int32 Declare Function ENsetvertices Lib "epanet2.dll" (ByVal index As Int32, xCoords As Any, yCoords As Any, ByVal count As Int32) As Int32 'Pump Functions @@ -372,6 +374,7 @@ Public Const EN_MISSING As Double = -1.0E10 Declare Function ENsetcurveid Lib "epanet2.dll" (ByVal index As Int32, ByVal newid As String) As Int32 Declare Function ENgetcurvelen Lib "epanet2.dll" (ByVal index As Int32, len_ As Int32) As Int32 Declare Function ENgetcurvetype Lib "epanet2.dll" (ByVal index As Int32, type_ As Int32) As Int32 + Declare Function ENsetcurvetype Lib "epanet2.dll" (ByVal index As Int32, ByVal type_ As Int32) As Int32 Declare Function ENgetcurvevalue Lib "epanet2.dll" (ByVal curveIndex As Int32, ByVal pointIndex As Int32, x As Single, y As Single) As Int32 Declare Function ENsetcurvevalue Lib "epanet2.dll" (ByVal curveIndex As Int32, ByVal pointIndex As Int32, ByVal x As Single, ByVal y As Single) As Int32 Declare Function ENgetcurve Lib "epanet2.dll" (ByVal index As Int32, ByVal id As String, nPoints As Int32, xValues As Any, yValues As Any) As Int32 diff --git a/include/epanet2_2.h b/include/epanet2_2.h index cbfece3..4760986 100644 --- a/include/epanet2_2.h +++ b/include/epanet2_2.h @@ -11,7 +11,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 10/29/2019 + Last Updated: 02/01/2020 ****************************************************************************** */ @@ -1238,7 +1238,7 @@ typedef struct Project *EN_Project; int DLLEXPORT EN_getvertexcount(EN_Project ph, int index, int *count); /** - @brief Retrieves the coordinate's of a vertex point assigned to a link. + @brief Retrieves the coordinates of a vertex point assigned to a link. @param ph an EPANET project handle. @param index a link's index (starting from 1). @param vertex a vertex point index (starting from 1). @@ -1248,6 +1248,17 @@ typedef struct Project *EN_Project; */ int DLLEXPORT EN_getvertex(EN_Project ph, int index, int vertex, double *x, double *y); + /** + @brief Sets the coordinates of a vertex point assigned to a link. + @param ph an EPANET project handle. + @param index a link's index (starting from 1). + @param vertex a vertex point index (starting from 1). + @param x the vertex's X-coordinate value. + @param y the vertex's Y-coordinate value. + @return an error code. + */ + int DLLEXPORT EN_setvertex(EN_Project ph, int index, int vertex, double x, double y); + /** @brief Assigns a set of internal vertex points to a link. @param ph an EPANET project handle. @@ -1475,6 +1486,15 @@ typedef struct Project *EN_Project; */ int DLLEXPORT EN_getcurvetype(EN_Project ph, int index, int *type); + /** + @brief Sets a curve's type. + @param ph an EPANET project handle. + @param index a curve's index (starting from 1). + @param type the curve's type (see @ref EN_CurveType). + @return an error code. + */ + int DLLEXPORT EN_setcurvetype(EN_Project ph, int index, int type); + /** @brief Retrieves the value of a single data point for a curve. @param ph an EPANET project handle. diff --git a/include/epanet2_enums.h b/include/epanet2_enums.h index 32d8171..7a665b6 100644 --- a/include/epanet2_enums.h +++ b/include/epanet2_enums.h @@ -9,7 +9,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 11/06/2019 + Last Updated: 02/01/2020 ****************************************************************************** */ @@ -94,7 +94,8 @@ typedef enum { EN_PUMP_HCURVE = 19, //!< Pump head v. flow curve index EN_PUMP_ECURVE = 20, //!< Pump efficiency v. flow curve index EN_PUMP_ECOST = 21, //!< Pump average energy price - EN_PUMP_EPAT = 22 //!< Pump energy price time pattern index + EN_PUMP_EPAT = 22, //!< Pump energy price time pattern index + EN_GPV_CURVE = 23 //!< GPV head loss v. flow curve index } EN_LinkProperty; /// Time parameters diff --git a/src/epanet.c b/src/epanet.c index efc0535..82ece34 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 11/15/2019 + Last Updated: 02/01/2020 ****************************************************************************** */ @@ -3778,6 +3778,12 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val v = (double)Pump[findpump(&p->network, index)].Epat; } break; + + case EN_GPV_CURVE: + if (Link[index].Type == GPV) + { + v = Link[index].Kc; + } default: return 251; @@ -3985,6 +3991,14 @@ int DLLEXPORT EN_setlinkvalue(EN_Project p, int index, int property, double valu net->Pump[pumpIndex].Epat = patIndex; } break; + + case EN_GPV_CURVE: + if (Link[index].Type == GPV) + { + curveIndex = ROUND(value); + if (curveIndex < 0 || curveIndex > net->Ncurves) return 206; + Link[index].Kc = curveIndex; + } default: return 251; @@ -4088,6 +4102,35 @@ int DLLEXPORT EN_getvertex(EN_Project p, int index, int vertex, double *x, doubl *y = vertices->Y[vertex - 1]; return 0; } + +int DLLEXPORT EN_setvertex(EN_Project p, int index, int vertex, double x, double y) +/*---------------------------------------------------------------- +** Input: index = link index +** vertex = index of a link vertex point +** x = vertex point's X-coordinate +** y = vertex point's Y-coordinate +** Returns: error code +** Purpose: sets the coordinates of a vertex point in a link +**---------------------------------------------------------------- +*/ +{ + Network *net = &p->network; + + Slink *Link = net->Link; + Pvertices vertices; + + // Check that link exists + if (!p->Openflag) return 102; + if (index <= 0 || index > net->Nlinks) return 204; + + // Check that vertex exists + vertices = Link[index].Vertices; + if (vertices == NULL) return 255; + if (vertex <= 0 || vertex > vertices->Npts) return 255; + vertices->X[vertex - 1] = x; + vertices->Y[vertex - 1] = y; + return 0; +} int DLLEXPORT EN_setvertices(EN_Project p, int index, double *x, double *y, int count) /*---------------------------------------------------------------- @@ -4699,6 +4742,23 @@ int DLLEXPORT EN_getcurvetype(EN_Project p, int index, int *type) return 0; } +int DLLEXPORT EN_setcurvetype(EN_Project p, int index, int type) +/*---------------------------------------------------------------- +** Input: index = data curve index +** type = type of data curve (see EN_CurveType) +** Returns: error code +** Purpose: sets the type assigned to a data curve +**---------------------------------------------------------------- +*/ +{ + Network *net = &p->network; + if (!p->Openflag) return 102; + if (index < 1 || index > net->Ncurves) return 206; + if (type < 0 || type > EN_GENERIC_CURVE) return 251; + net->Curve[index].Type = type; + return 0; +} + int DLLEXPORT EN_getcurvevalue(EN_Project p, int curveIndex, int pointIndex, double *x, double *y) /*---------------------------------------------------------------- diff --git a/src/epanet2.c b/src/epanet2.c index f5940cd..1706936 100644 --- a/src/epanet2.c +++ b/src/epanet2.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 11/02/2019 + Last Updated: 02/01/2020 ****************************************************************************** */ @@ -522,6 +522,11 @@ int DLLEXPORT ENgetvertex(int index, int vertex, double *x, double *y) return EN_getvertex(_defaultProject, index, vertex, x, y); } +int DLLEXPORT ENsetvertex(int index, int vertex, double x, double y) +{ + return EN_setvertex(_defaultProject, index, vertex, x, y); +} + int DLLEXPORT ENsetvertices(int index, double *x, double *y, int count) { return EN_setvertices(_defaultProject, index, x, y, count); @@ -662,6 +667,11 @@ int DLLEXPORT ENgetcurvetype(int index, int *type) return EN_getcurvetype(_defaultProject, index, type); } +int DLLEXPORT ENsetcurvetype(int index, int type) +{ + return EN_setcurvetype(_defaultProject, index, type); +} + int DLLEXPORT ENgetcurvevalue(int curveIndex, int pointIndex, EN_API_FLOAT_TYPE *x, EN_API_FLOAT_TYPE *y) { diff --git a/src/errors.dat b/src/errors.dat index 0930b52..4768cc3 100644 --- a/src/errors.dat +++ b/src/errors.dat @@ -43,7 +43,7 @@ DAT(225,"invalid lower/upper levels for tank") DAT(226,"no head curve or power rating for pump") DAT(227,"invalid head curve for pump") DAT(230,"nonincreasing x-values for curve") -DAT(233,"network has unconnected node") +DAT(233,"network has unconnected nodes") // These errors apply only to API functions DAT(240,"nonexistent source") diff --git a/src/funcs.h b/src/funcs.h index d2bf39e..366e5c7 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 11/15/2019 + Last Updated: 02/03/2020 ****************************************************************************** */ #ifndef FUNCS_H @@ -29,6 +29,8 @@ void freeadjlists(Network *); int incontrols(Project *, int, int); int valvecheck(Project *, int, int, int, int); +int unlinked(Project *); + int findnode(Network *, char *); int findlink(Network *, char *); int findtank(Network *, int); diff --git a/src/hydraul.c b/src/hydraul.c index 732108c..ad24f0e 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 12/05/2019 + Last Updated: 02/03/2020 ****************************************************************************** */ @@ -64,14 +64,7 @@ int openhyd(Project *pr) ERRCODE(allocmatrix(pr)); // Check for unconnected nodes - if (!errcode) for (i = 1; i <= pr->network.Njuncs; i++) - { - if (pr->network.Adjlist[i] == NULL) - { - errcode = 233; - break; - } - } + ERRCODE(unlinked(pr)); // Initialize link flows if (!errcode) for (i = 1; i <= pr->network.Nlinks; i++) diff --git a/src/input2.c b/src/input2.c index 766ce24..a4d7fc8 100644 --- a/src/input2.c +++ b/src/input2.c @@ -7,7 +7,7 @@ Description: reads and interprets network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 10/29/2019 +Last Updated: 02/03/2020 ****************************************************************************** */ @@ -37,7 +37,6 @@ extern int powercurve(double, double, double, double, double, double *, static int newline(Project *, int, char *); static int addpattern(Network *, char *); static int addcurve(Network *, char *); -static int unlinked(Project *); static int getpumpparams(Project *); static void inperrmsg(Project *, int, int, char *); @@ -130,11 +129,6 @@ int netsize(Project *pr) parser->MaxNodes = parser->MaxJuncs + parser->MaxTanks; parser->MaxLinks = parser->MaxPipes + parser->MaxPumps + parser->MaxValves; if (parser->MaxPats < 1) parser->MaxPats = 1; - if (!errcode) - { - if (parser->MaxJuncs < 1) errcode = 223; // Not enough nodes - else if (parser->MaxTanks == 0) errcode = 224; // No tanks - } return errcode; } @@ -263,9 +257,6 @@ int readdata(Project *pr) // Check for errors if (errsum > 0) errcode = 200; - // Check for unlinked nodes - if (!errcode) errcode = unlinked(pr); - // Determine pump curve parameters if (!errcode) errcode = getpumpparams(pr); @@ -572,54 +563,6 @@ int addcurve(Network *network, char *id) return 0; } -int unlinked(Project *pr) -/* -**-------------------------------------------------------------- -** Input: none -** Output: returns error code if any unlinked junctions found -** Purpose: checks for unlinked junctions in network -** -** NOTE: unlinked tanks have no effect on computations. -**-------------------------------------------------------------- -*/ -{ - Network *net = &pr->network; - int *marked; - int i, err, errcode; - - errcode = 0; - err = 0; - - // Create an array to record number of links incident on each node - marked = (int *)calloc(net->Nnodes + 1, sizeof(int)); - ERRCODE(MEMCHECK(marked)); - if (errcode) return errcode; - memset(marked, 0, (net->Nnodes + 1) * sizeof(int)); - - // Mark end nodes of each link - for (i = 1; i <= net->Nlinks; i++) - { - marked[net->Link[i].N1]++; - marked[net->Link[i].N2]++; - } - - // Check each junction - for (i = 1; i <= net->Njuncs; i++) - { - // If not marked then error - if (marked[i] == 0) - { - err++; - sprintf(pr->Msg, "Error 233: %s %s", geterrmsg(233, pr->Msg), net->Node[i].ID); - writeline(pr, pr->Msg); - } - if (err >= MAXERRS) break; - } - if (err > 0) errcode = 200; - free(marked); - return errcode; -} - int findmatch(char *line, char *keyword[]) /* **-------------------------------------------------------------- diff --git a/src/project.c b/src/project.c index 2d19729..03f0d9a 100644 --- a/src/project.c +++ b/src/project.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 11/15/2019 + Last Updated: 02/03/2020 ****************************************************************************** */ @@ -795,6 +795,34 @@ int valvecheck(Project *pr, int index, int type, int j1, int j2) return 0; } +int unlinked(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns error code if any unlinked junctions found +** Purpose: checks for unlinked junctions in network +** +** NOTE: unlinked tanks have no effect on computations. +**-------------------------------------------------------------- +*/ +{ + Network *net = &pr->network; + int i, count = 0; + + for (i = 1; i <= net->Njuncs; i++) + { + if (pr->network.Adjlist[i] == NULL) + { + count++; + sprintf(pr->Msg, "Error 233: %s %s", geterrmsg(233, pr->Msg), net->Node[i].ID); + writeline(pr, pr->Msg); + } + if (count >= 10) break; + } + if (count > 0) return 233; + return 0; +} + int findnode(Network *network, char *id) /*---------------------------------------------------------------- ** Input: id = node ID @@ -912,8 +940,8 @@ void adjustpattern(int *pat, int index) **---------------------------------------------------------------- */ { - if (*pat == index) *pat = 0; - else if (*pat > index) (*pat)--; + if (*pat == index) *pat = 0; + else if (*pat > index) (*pat)--; } void adjustpatterns(Network *network, int index) diff --git a/src/quality.c b/src/quality.c index 05d6519..41869fe 100644 --- a/src/quality.c +++ b/src/quality.c @@ -7,7 +7,7 @@ Description: implements EPANET's water quality engine Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 05/15/2019 +Last Updated: 02/03/2020 ****************************************************************************** */ @@ -63,8 +63,16 @@ int openqual(Project *pr) // Build nodal adjacency lists if they don't already exist if (net->Adjlist == NULL) { + // Check for too few nodes & no fixed grade nodes + if (net->Nnodes < 2) return 223; + if (net->Ntanks == 0) return 224; + + // Build adjacency lists errcode = buildadjlists(net); if (errcode ) return errcode; + + // Check for unconnected nodes + if (errcode = unlinked(pr)) return errcode; } // Create a memory pool for water quality segments From 3ee30ce01951522646a4cdc1e6950ac5bccd4b05 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Fri, 7 Feb 2020 10:44:52 -0500 Subject: [PATCH 2/8] Replaced status checking for pumps & FCVs See ReleaseNotes2_3.md. --- ReleaseNotes2_3.md | 7 +++- src/hydcoeffs.c | 94 ++++++++++++++++++++-------------------------- src/hydsolver.c | 12 +++--- src/hydstatus.c | 14 +------ src/input2.c | 13 +++++-- 5 files changed, 64 insertions(+), 76 deletions(-) diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index 1f922fc..f6b2f68 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -3,9 +3,14 @@ This document describes the changes and updates that have been made in version 2.3 of EPANET. - - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files to be opened by the toolkit. + - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files can be opened by the toolkit. - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). - A `EN_setvertex` function was added to allow API clients to change the coordinates of a link's vertex. - The index of a General Purpose Valve's (GPV's) head loss curve was added to the list of editable Link Properties using the symbolic constant name `EN_GPV_CURVE`. - The `EN_getlinkvalue` and `EN_setlinkvalue` functions were updated to get and set the value of `EN_GPV_CURVE`. + - For `EN_CUSTOM` type pump curves the maximum head value is now extrapolated to the y-axis intercept instead of being based on the first curve data point. Similarly, the maximum flow value is extrapolated to the x-axis intercept. + - Status checking for a pump not able to deliver enough head has been replaced by adding a penalty term to the pump's operating curve that prevents it from having negative flow (i.e., from crossing the y-axis). + - Status checking for Flow Control Valves has been eliminated by using a continuous head v. flow function. If the current flow is below the valve setting then the normal open head loss relation is used; otherwise a linear penalty function is applied to any flow in excess of the setting. Warnings are no longer issued when the valve operates fully opened at flows below the setting. + + diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index 9f4ac3d..b461946 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 10/04/2019 + Last Updated: 02/07/2020 ****************************************************************************** */ @@ -59,9 +59,9 @@ static void valvecoeff(Project *pr, int k); static void gpvcoeff(Project *pr, int k); static void pbvcoeff(Project *pr, int k); static void tcvcoeff(Project *pr, int k); +static void fcvcoeff(Project *pr, int k); static void prvcoeff(Project *pr, int k, int n1, int n2); static void psvcoeff(Project *pr, int k, int n1, int n2); -static void fcvcoeff(Project *pr, int k, int n1, int n2); void resistcoeff(Project *pr, int k) @@ -152,6 +152,8 @@ void headlosscoeffs(Project *pr) gpvcoeff(pr, k); break; case FCV: + fcvcoeff(pr, k); + break; case PRV: case PSV: if (hyd->LinkSetting[k] == MISSING) valvecoeff(pr, k); @@ -285,8 +287,8 @@ void valvecoeffs(Project *pr) ** Input: none ** Output: none ** Purpose: computes coeffs. of the linearized hydraulic eqns. -** contributed by PRVs, PSVs & FCVs whose status is -** not fixed to OPEN/CLOSED +** contributed by PRVs & PSVs whose status is not +** fixed to OPEN/CLOSED **-------------------------------------------------------------- */ { @@ -313,19 +315,8 @@ void valvecoeffs(Project *pr) n2 = link->N2; // Call valve-specific function - switch (link->Type) - { - case PRV: - prvcoeff(pr, k, n1, n2); - break; - case PSV: - psvcoeff(pr, k, n1, n2); - break; - case FCV: - fcvcoeff(pr, k, n1, n2); - break; - default: continue; - } + if (link->Type == PRV) prvcoeff(pr, k, n1, n2); + if (link->Type == PSV) psvcoeff(pr, k, n1, n2); } } @@ -701,20 +692,31 @@ void pumpcoeff(Project *pr, int k) } // Obtain reference to pump object - q = ABS(hyd->LinkFlow[k]); p = findpump(&pr->network, k); pump = &pr->network.Pump[p]; + + // Prevent negative flow + q = hyd->LinkFlow[k]; + if (q < 0.0) + { + hloss = -(SQR(setting) * pump->Hmax) + CBIG * q; + hgrad = CBIG; + hyd->P[k] = 1.0 / hgrad; + hyd->Y[k] = hloss / hgrad; + return; + } // If no pump curve treat pump as an open valve if (pump->Ptype == NOCURVE) { hyd->P[k] = 1.0 / CSMALL; - hyd->Y[k] = hyd->LinkFlow[k]; + hyd->Y[k] = q; return; } // Get pump curve coefficients for custom pump curve // (Other pump types have pre-determined coeffs.) + q = ABS(q); if (pump->Ptype == CUSTOM) { // Find intercept (h0) & slope (r) of pump curve @@ -1044,12 +1046,10 @@ void psvcoeff(Project *pr, int k, int n1, int n2) } -void fcvcoeff(Project *pr, int k, int n1, int n2) +void fcvcoeff(Project *pr, int k) /* **-------------------------------------------------------------- ** Input: k = link index -** n1 = upstream node of valve -** n2 = downstream node of valve ** Output: none ** Purpose: computes solution matrix coeffs. for flow control ** valve @@ -1059,40 +1059,28 @@ void fcvcoeff(Project *pr, int k, int n1, int n2) Hydraul *hyd = &pr->hydraul; Smatrix *sm = &hyd->smatrix; - int i, j; // Rows in solution matrix - double q; // Valve flow setting + double qset; // Valve flow setting + double flow; // Current valve flow + double hloss, hgrad; // Head loss & gradient - q = hyd->LinkSetting[k]; - i = sm->Row[n1]; - j = sm->Row[n2]; - - // If valve active, break network at valve and treat - // flow setting as external demand at upstream node - // and external supply at downstream node. - - if (hyd->LinkStatus[k] == ACTIVE) - { - hyd->Xflow[n1] -= q; - hyd->Xflow[n2] += q; - hyd->Y[k] = hyd->LinkFlow[k] - q; - sm->F[i] -= q; - sm->F[j] += q; - hyd->P[k] = 1.0 / CBIG; - sm->Aij[sm->Ndx[k]] -= hyd->P[k]; - sm->Aii[i] += hyd->P[k]; - sm->Aii[j] += hyd->P[k]; - } - - // Otherwise treat valve as an open pipe - - else + // Treat as a regular valve if status fixed or flow below setting + qset = hyd->LinkSetting[k]; + flow = hyd->LinkFlow[k]; + if (qset == MISSING || hyd->LinkStatus[k] <= CLOSED || flow < qset) { valvecoeff(pr, k); - sm->Aij[sm->Ndx[k]] -= hyd->P[k]; - sm->Aii[i] += hyd->P[k]; - sm->Aii[j] += hyd->P[k]; - sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]); - sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]); + } + + // Otherwise prevent flow from exceeding the setting + else + { + hyd->LinkFlow[k] = qset; + valvecoeff(pr, k); + hloss = hyd->Y[k] / hyd->P[k] + CBIG * (flow - qset); + hgrad = CBIG; + hyd->P[k] = 1.0 / hgrad; + hyd->Y[k] = hloss / hgrad; + hyd->LinkFlow[k] = flow; } } diff --git a/src/hydsolver.c b/src/hydsolver.c index fcb8258..53057c2 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -8,7 +8,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 07/15/2019 + Last Updated: 02/07/2020 ****************************************************************************** */ @@ -111,6 +111,7 @@ int hydsolve(Project *pr, int *iter, double *relerr) maxtrials = hyd->MaxIter; if (hyd->ExtraIter > 0) maxtrials += hyd->ExtraIter; *iter = 1; + headlosscoeffs(pr); while (*iter <= maxtrials) { // Compute coefficient matrices A & F and solve A*H = F @@ -118,7 +119,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) // head loss gradients, & F = flow correction terms. // Solution for H is returned in F from call to linsolve(). - headlosscoeffs(pr); matrixcoeffs(pr); errcode = linsolve(sm, net->Njuncs); @@ -138,6 +138,9 @@ int hydsolve(Project *pr, int *iter, double *relerr) } newerr = newflows(pr, &hydbal); // Update flows *relerr = newerr; + + // Compute head loss coeffs. for new flows + headlosscoeffs(pr); // Write convergence error to status report if called for if (rpt->Statflag == FULL) @@ -243,7 +246,7 @@ int badvalve(Project *pr, int n) if (n == n1 || n == n2) { t = link->Type; - if (t == PRV || t == PSV || t == FCV) + if (t == PRV || t == PSV) { if (hyd->LinkStatus[k] == ACTIVE) { @@ -253,8 +256,7 @@ int badvalve(Project *pr, int n) clocktime(rpt->Atime, time->Htime), link->ID); writeline(pr, pr->Msg); } - if (link->Type == FCV) hyd->LinkStatus[k] = XFCV; - else hyd->LinkStatus[k] = XPRESSURE; + hyd->LinkStatus[k] = XPRESSURE; return 1; } } diff --git a/src/hydstatus.c b/src/hydstatus.c index 87d45c6..76ca823 100644 --- a/src/hydstatus.c +++ b/src/hydstatus.c @@ -7,7 +7,7 @@ Description: updates hydraulic status of network elements Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 05/15/2019 +Last Updated: 02/07/2020 ****************************************************************************** */ @@ -141,18 +141,6 @@ int linkstatus(Project *pr) hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh, hyd->LinkFlow[k]); } - if (link->Type == 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 == 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) diff --git a/src/input2.c b/src/input2.c index a4d7fc8..136d297 100644 --- a/src/input2.c +++ b/src/input2.c @@ -7,7 +7,7 @@ Description: reads and interprets network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 02/03/2020 +Last Updated: 02/07/2020 ****************************************************************************** */ @@ -429,9 +429,14 @@ int updatepumpparams(Project *pr, int pumpindex) { if (curve->Y[m] >= curve->Y[m - 1]) return 227; } - pump->Qmax = curve->X[npts - 1]; - pump->Q0 = (curve->X[0] + pump->Qmax) / 2.0; - pump->Hmax = curve->Y[0]; + pump->Q0 = (curve->X[0] + curve->X[npts-1]) / 2.0; + + // Extend curve to find Hmax (at 0 flow) and Qmax (at 0 head) + b = (curve->Y[1] - curve->Y[0]) / (curve->X[1] - curve->X[0]); + pump->Hmax = curve->Y[0] + b * curve->X[0]; + b = (curve->Y[npts-1] - curve->Y[npts-2]) / + (curve->X[npts-1] - curve->X[npts-2]); + pump->Qmax = curve->X[npts-1] - curve->Y[npts-1] / b; } // Compute shape factors & limits of power function curves From 7e3d5e1fe07202463842342f7b165fe855cc1b92 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Fri, 7 Feb 2020 13:49:54 -0500 Subject: [PATCH 3/8] Changed benchmark results for EN_getstatistic test --- ReleaseNotes2_3.md | 1 + tests/test_report.cpp | 7 +++++-- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index f6b2f68..8d47c19 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -11,6 +11,7 @@ This document describes the changes and updates that have been made in version 2 - For `EN_CUSTOM` type pump curves the maximum head value is now extrapolated to the y-axis intercept instead of being based on the first curve data point. Similarly, the maximum flow value is extrapolated to the x-axis intercept. - Status checking for a pump not able to deliver enough head has been replaced by adding a penalty term to the pump's operating curve that prevents it from having negative flow (i.e., from crossing the y-axis). - Status checking for Flow Control Valves has been eliminated by using a continuous head v. flow function. If the current flow is below the valve setting then the normal open head loss relation is used; otherwise a linear penalty function is applied to any flow in excess of the setting. Warnings are no longer issued when the valve operates fully opened at flows below the setting. + - The maximum link head loss error convergence criterion is now evaluated using the most recently computed link flows instead of flows from the previous iteration. diff --git a/tests/test_report.cpp b/tests/test_report.cpp index 5b5850c..c2c2a0d 100644 --- a/tests/test_report.cpp +++ b/tests/test_report.cpp @@ -25,8 +25,11 @@ BOOST_FIXTURE_TEST_CASE(test_rprt_anlysstats, FixtureOpenClose) std::vector test(5); double *array = test.data(); - std::vector ref = {3.0, 7.0799498320679432e-06, 1.6680242187483429e-08, - 0.0089173150106518495, 0.99999998187144024}; + std::vector ref = +// {3.0, 7.0799498320679432e-06, 1.6680242187483429e-08, // v2.2 +// 0.0089173150106518495, 0.99999998187144024}; + {3.0, 8.3792202148e-6, 2.63983750e-8, // v2.3_dev + 0.0112012133155924, 0.9999999807954413}; error = EN_solveH(ph); BOOST_REQUIRE(error == 0); From a3537b767aabb2bb721438993becd1380561cfe0 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Sun, 9 Feb 2020 16:58:37 -0500 Subject: [PATCH 4/8] Revert "Changed benchmark results for EN_getstatistic test" This reverts commit 7e3d5e1fe07202463842342f7b165fe855cc1b92. --- ReleaseNotes2_3.md | 1 - tests/test_report.cpp | 7 ++----- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index 8d47c19..f6b2f68 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -11,7 +11,6 @@ This document describes the changes and updates that have been made in version 2 - For `EN_CUSTOM` type pump curves the maximum head value is now extrapolated to the y-axis intercept instead of being based on the first curve data point. Similarly, the maximum flow value is extrapolated to the x-axis intercept. - Status checking for a pump not able to deliver enough head has been replaced by adding a penalty term to the pump's operating curve that prevents it from having negative flow (i.e., from crossing the y-axis). - Status checking for Flow Control Valves has been eliminated by using a continuous head v. flow function. If the current flow is below the valve setting then the normal open head loss relation is used; otherwise a linear penalty function is applied to any flow in excess of the setting. Warnings are no longer issued when the valve operates fully opened at flows below the setting. - - The maximum link head loss error convergence criterion is now evaluated using the most recently computed link flows instead of flows from the previous iteration. diff --git a/tests/test_report.cpp b/tests/test_report.cpp index c2c2a0d..5b5850c 100644 --- a/tests/test_report.cpp +++ b/tests/test_report.cpp @@ -25,11 +25,8 @@ BOOST_FIXTURE_TEST_CASE(test_rprt_anlysstats, FixtureOpenClose) std::vector test(5); double *array = test.data(); - std::vector ref = -// {3.0, 7.0799498320679432e-06, 1.6680242187483429e-08, // v2.2 -// 0.0089173150106518495, 0.99999998187144024}; - {3.0, 8.3792202148e-6, 2.63983750e-8, // v2.3_dev - 0.0112012133155924, 0.9999999807954413}; + std::vector ref = {3.0, 7.0799498320679432e-06, 1.6680242187483429e-08, + 0.0089173150106518495, 0.99999998187144024}; error = EN_solveH(ph); BOOST_REQUIRE(error == 0); From a9ab376aa766911143b108cce82936473f6a538e Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Sun, 9 Feb 2020 17:25:33 -0500 Subject: [PATCH 5/8] Revert "Replaced status checking for pumps & FCVs" This reverts commit 3ee30ce01951522646a4cdc1e6950ac5bccd4b05. --- ReleaseNotes2_3.md | 7 +--- src/hydcoeffs.c | 88 ++++++++++++++++++++++++++-------------------- src/hydsolver.c | 12 +++---- src/hydstatus.c | 14 +++++++- src/input2.c | 13 +++---- 5 files changed, 73 insertions(+), 61 deletions(-) diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index f6b2f68..1f922fc 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -3,14 +3,9 @@ This document describes the changes and updates that have been made in version 2.3 of EPANET. - - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files can be opened by the toolkit. + - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files to be opened by the toolkit. - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). - A `EN_setvertex` function was added to allow API clients to change the coordinates of a link's vertex. - The index of a General Purpose Valve's (GPV's) head loss curve was added to the list of editable Link Properties using the symbolic constant name `EN_GPV_CURVE`. - The `EN_getlinkvalue` and `EN_setlinkvalue` functions were updated to get and set the value of `EN_GPV_CURVE`. - - For `EN_CUSTOM` type pump curves the maximum head value is now extrapolated to the y-axis intercept instead of being based on the first curve data point. Similarly, the maximum flow value is extrapolated to the x-axis intercept. - - Status checking for a pump not able to deliver enough head has been replaced by adding a penalty term to the pump's operating curve that prevents it from having negative flow (i.e., from crossing the y-axis). - - Status checking for Flow Control Valves has been eliminated by using a continuous head v. flow function. If the current flow is below the valve setting then the normal open head loss relation is used; otherwise a linear penalty function is applied to any flow in excess of the setting. Warnings are no longer issued when the valve operates fully opened at flows below the setting. - - diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index b461946..9f4ac3d 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/07/2020 + Last Updated: 10/04/2019 ****************************************************************************** */ @@ -59,9 +59,9 @@ static void valvecoeff(Project *pr, int k); static void gpvcoeff(Project *pr, int k); static void pbvcoeff(Project *pr, int k); static void tcvcoeff(Project *pr, int k); -static void fcvcoeff(Project *pr, int k); static void prvcoeff(Project *pr, int k, int n1, int n2); static void psvcoeff(Project *pr, int k, int n1, int n2); +static void fcvcoeff(Project *pr, int k, int n1, int n2); void resistcoeff(Project *pr, int k) @@ -152,8 +152,6 @@ void headlosscoeffs(Project *pr) gpvcoeff(pr, k); break; case FCV: - fcvcoeff(pr, k); - break; case PRV: case PSV: if (hyd->LinkSetting[k] == MISSING) valvecoeff(pr, k); @@ -287,8 +285,8 @@ void valvecoeffs(Project *pr) ** Input: none ** Output: none ** Purpose: computes coeffs. of the linearized hydraulic eqns. -** contributed by PRVs & PSVs whose status is not -** fixed to OPEN/CLOSED +** contributed by PRVs, PSVs & FCVs whose status is +** not fixed to OPEN/CLOSED **-------------------------------------------------------------- */ { @@ -315,8 +313,19 @@ void valvecoeffs(Project *pr) n2 = link->N2; // Call valve-specific function - if (link->Type == PRV) prvcoeff(pr, k, n1, n2); - if (link->Type == PSV) psvcoeff(pr, k, n1, n2); + switch (link->Type) + { + case PRV: + prvcoeff(pr, k, n1, n2); + break; + case PSV: + psvcoeff(pr, k, n1, n2); + break; + case FCV: + fcvcoeff(pr, k, n1, n2); + break; + default: continue; + } } } @@ -692,31 +701,20 @@ void pumpcoeff(Project *pr, int k) } // Obtain reference to pump object + q = ABS(hyd->LinkFlow[k]); p = findpump(&pr->network, k); pump = &pr->network.Pump[p]; - - // Prevent negative flow - q = hyd->LinkFlow[k]; - if (q < 0.0) - { - hloss = -(SQR(setting) * pump->Hmax) + CBIG * q; - hgrad = CBIG; - hyd->P[k] = 1.0 / hgrad; - hyd->Y[k] = hloss / hgrad; - return; - } // If no pump curve treat pump as an open valve if (pump->Ptype == NOCURVE) { hyd->P[k] = 1.0 / CSMALL; - hyd->Y[k] = q; + hyd->Y[k] = hyd->LinkFlow[k]; return; } // Get pump curve coefficients for custom pump curve // (Other pump types have pre-determined coeffs.) - q = ABS(q); if (pump->Ptype == CUSTOM) { // Find intercept (h0) & slope (r) of pump curve @@ -1046,10 +1044,12 @@ void psvcoeff(Project *pr, int k, int n1, int n2) } -void fcvcoeff(Project *pr, int k) +void fcvcoeff(Project *pr, int k, int n1, int n2) /* **-------------------------------------------------------------- ** Input: k = link index +** n1 = upstream node of valve +** n2 = downstream node of valve ** Output: none ** Purpose: computes solution matrix coeffs. for flow control ** valve @@ -1059,28 +1059,40 @@ void fcvcoeff(Project *pr, int k) Hydraul *hyd = &pr->hydraul; Smatrix *sm = &hyd->smatrix; - double qset; // Valve flow setting - double flow; // Current valve flow - double hloss, hgrad; // Head loss & gradient + int i, j; // Rows in solution matrix + double q; // Valve flow setting - // Treat as a regular valve if status fixed or flow below setting - qset = hyd->LinkSetting[k]; - flow = hyd->LinkFlow[k]; - if (qset == MISSING || hyd->LinkStatus[k] <= CLOSED || flow < qset) + q = hyd->LinkSetting[k]; + i = sm->Row[n1]; + j = sm->Row[n2]; + + // If valve active, break network at valve and treat + // flow setting as external demand at upstream node + // and external supply at downstream node. + + if (hyd->LinkStatus[k] == ACTIVE) { - valvecoeff(pr, k); + hyd->Xflow[n1] -= q; + hyd->Xflow[n2] += q; + hyd->Y[k] = hyd->LinkFlow[k] - q; + sm->F[i] -= q; + sm->F[j] += q; + hyd->P[k] = 1.0 / CBIG; + sm->Aij[sm->Ndx[k]] -= hyd->P[k]; + sm->Aii[i] += hyd->P[k]; + sm->Aii[j] += hyd->P[k]; } - - // Otherwise prevent flow from exceeding the setting + + // Otherwise treat valve as an open pipe + else { - hyd->LinkFlow[k] = qset; valvecoeff(pr, k); - hloss = hyd->Y[k] / hyd->P[k] + CBIG * (flow - qset); - hgrad = CBIG; - hyd->P[k] = 1.0 / hgrad; - hyd->Y[k] = hloss / hgrad; - hyd->LinkFlow[k] = flow; + sm->Aij[sm->Ndx[k]] -= hyd->P[k]; + sm->Aii[i] += hyd->P[k]; + sm->Aii[j] += hyd->P[k]; + sm->F[i] += (hyd->Y[k] - hyd->LinkFlow[k]); + sm->F[j] -= (hyd->Y[k] - hyd->LinkFlow[k]); } } diff --git a/src/hydsolver.c b/src/hydsolver.c index 53057c2..fcb8258 100644 --- a/src/hydsolver.c +++ b/src/hydsolver.c @@ -8,7 +8,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 02/07/2020 + Last Updated: 07/15/2019 ****************************************************************************** */ @@ -111,7 +111,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) maxtrials = hyd->MaxIter; if (hyd->ExtraIter > 0) maxtrials += hyd->ExtraIter; *iter = 1; - headlosscoeffs(pr); while (*iter <= maxtrials) { // Compute coefficient matrices A & F and solve A*H = F @@ -119,6 +118,7 @@ int hydsolve(Project *pr, int *iter, double *relerr) // head loss gradients, & F = flow correction terms. // Solution for H is returned in F from call to linsolve(). + headlosscoeffs(pr); matrixcoeffs(pr); errcode = linsolve(sm, net->Njuncs); @@ -138,9 +138,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) } newerr = newflows(pr, &hydbal); // Update flows *relerr = newerr; - - // Compute head loss coeffs. for new flows - headlosscoeffs(pr); // Write convergence error to status report if called for if (rpt->Statflag == FULL) @@ -246,7 +243,7 @@ int badvalve(Project *pr, int n) if (n == n1 || n == n2) { t = link->Type; - if (t == PRV || t == PSV) + if (t == PRV || t == PSV || t == FCV) { if (hyd->LinkStatus[k] == ACTIVE) { @@ -256,7 +253,8 @@ int badvalve(Project *pr, int n) clocktime(rpt->Atime, time->Htime), link->ID); writeline(pr, pr->Msg); } - hyd->LinkStatus[k] = XPRESSURE; + if (link->Type == FCV) hyd->LinkStatus[k] = XFCV; + else hyd->LinkStatus[k] = XPRESSURE; return 1; } } diff --git a/src/hydstatus.c b/src/hydstatus.c index 76ca823..87d45c6 100644 --- a/src/hydstatus.c +++ b/src/hydstatus.c @@ -7,7 +7,7 @@ Description: updates hydraulic status of network elements Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 02/07/2020 +Last Updated: 05/15/2019 ****************************************************************************** */ @@ -141,6 +141,18 @@ int linkstatus(Project *pr) hyd->LinkStatus[k] = cvstatus(pr, hyd->LinkStatus[k], dh, hyd->LinkFlow[k]); } + if (link->Type == 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 == 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) diff --git a/src/input2.c b/src/input2.c index 136d297..a4d7fc8 100644 --- a/src/input2.c +++ b/src/input2.c @@ -7,7 +7,7 @@ Description: reads and interprets network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 02/07/2020 +Last Updated: 02/03/2020 ****************************************************************************** */ @@ -429,14 +429,9 @@ int updatepumpparams(Project *pr, int pumpindex) { if (curve->Y[m] >= curve->Y[m - 1]) return 227; } - pump->Q0 = (curve->X[0] + curve->X[npts-1]) / 2.0; - - // Extend curve to find Hmax (at 0 flow) and Qmax (at 0 head) - b = (curve->Y[1] - curve->Y[0]) / (curve->X[1] - curve->X[0]); - pump->Hmax = curve->Y[0] + b * curve->X[0]; - b = (curve->Y[npts-1] - curve->Y[npts-2]) / - (curve->X[npts-1] - curve->X[npts-2]); - pump->Qmax = curve->X[npts-1] - curve->Y[npts-1] / b; + pump->Qmax = curve->X[npts - 1]; + pump->Q0 = (curve->X[0] + pump->Qmax) / 2.0; + pump->Hmax = curve->Y[0]; } // Compute shape factors & limits of power function curves From 71a733ed85cc89e77267cb578c07394a7464c8a3 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Fri, 10 Jul 2020 15:42:34 -0400 Subject: [PATCH 6/8] Add missing break statement after update from dev --- ReleaseNotes2_3.md | 4 ++-- src/epanet.c | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/ReleaseNotes2_3.md b/ReleaseNotes2_3.md index 1f922fc..83d49c2 100644 --- a/ReleaseNotes2_3.md +++ b/ReleaseNotes2_3.md @@ -3,9 +3,9 @@ This document describes the changes and updates that have been made in version 2.3 of EPANET. - - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files to be opened by the toolkit. + - The check for at least two nodes, one tank/reservoir and no unconnected junction nodes was moved from `EN_open` to `EN_openH` and `EN_openQ` so that partial network data files could be opened by the toolkit. - A `EN_setcurvetype` function was added to allow API clients to set a curve's type (e.g., `EN_PUMP_CURVE,` `EN_VOLUME_CURVE,` etc.). - - A `EN_setvertex` function was added to allow API clients to change the coordinates of a link's vertex. + - A `EN_setvertex` function was added to allow API clients to change the coordinates of a single link vertex. - The index of a General Purpose Valve's (GPV's) head loss curve was added to the list of editable Link Properties using the symbolic constant name `EN_GPV_CURVE`. - The `EN_getlinkvalue` and `EN_setlinkvalue` functions were updated to get and set the value of `EN_GPV_CURVE`. diff --git a/src/epanet.c b/src/epanet.c index 00c0f41..47882ba 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -3796,6 +3796,7 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val { v = Link[index].Kc; } + break; case EN_LINK_INCONTROL: v = (double)incontrols(p, LINK, index); From 951c7ce727208278a327b0dcab720055c2ff1fb7 Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Sat, 11 Jul 2020 17:12:49 -0400 Subject: [PATCH 7/8] Refactor tank mixing zone volume Addresses issue #597. --- src/epanet.c | 21 +++++++-------------- src/inpfile.c | 2 +- src/input1.c | 1 - src/input3.c | 4 ++-- src/qualreact.c | 2 +- src/qualroute.c | 2 +- src/types.h | 4 ++-- 7 files changed, 14 insertions(+), 22 deletions(-) diff --git a/src/epanet.c b/src/epanet.c index 47882ba..be8c83d 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -1847,7 +1847,7 @@ int DLLEXPORT EN_addnode(EN_Project p, char *id, int nodeType, int *index) tank->Pat = 0; tank->Vcurve = 0; tank->MixModel = 0; - tank->V1max = 10000; + tank->V1frac = 1; tank->CanOverflow = FALSE; } net->Nnodes++; @@ -2164,7 +2164,9 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val case EN_MIXZONEVOL: v = 0.0; - if (index > nJuncs) v = Tank[index - nJuncs].V1max * Ucf[VOLUME]; + if (index > nJuncs) + v = Tank[index - nJuncs].V1frac * Tank[index - nJuncs].Vmax * + Ucf[VOLUME]; break; case EN_DEMAND: @@ -2224,9 +2226,9 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val case EN_MIXFRACTION: v = 1.0; - if (index > nJuncs && Tank[index - nJuncs].Vmax > 0.0) + if (index > nJuncs) { - v = Tank[index - nJuncs].V1max / Tank[index - nJuncs].Vmax; + v = Tank[index - nJuncs].V1frac; } break; @@ -2293,7 +2295,6 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu int i, j, n; Psource source; double hTmp; - double vTmp; if (!p->Openflag) return 102; if (index <= 0 || index > nNodes) return 203; @@ -2418,9 +2419,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu // shape below Hmin. Vmin can always be changed by setting // EN_MINVOLUME in a subsequent function call. Tank[j].V0 = tankvolume(p, j, Tank[j].H0); // new init. volume - vTmp = Tank[j].Vmax; // old max. volume Tank[j].Vmax = tankvolume(p, j, Tank[j].Hmax); // new max. volume - Tank[j].V1max *= Tank[j].Vmax / vTmp; // new mix zone volume break; case EN_MINVOLUME: @@ -2450,9 +2449,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu // Since Vmin changes the other volumes need updating Tank[j].V0 = tankvolume(p, j, Tank[j].H0); // new init. volume - vTmp = Tank[j].Vmax; // old max. volume Tank[j].Vmax = tankvolume(p, j, Tank[j].Hmax); // new max. volume - Tank[j].V1max *= Tank[j].Vmax / vTmp; // new mix zone volume } break; @@ -2477,9 +2474,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu Tank[j].Vcurve = i; // assign curve to tank Tank[j].Vmin = tankvolume(p, j, Tank[j].Hmin); // new min. volume Tank[j].V0 = tankvolume(p, j, Tank[j].H0); // new init. volume - vTmp = Tank[j].Vmax; // old max. volume Tank[j].Vmax = tankvolume(p, j, Tank[j].Hmax); // new max. volume - Tank[j].V1max *= Tank[j].Vmax / vTmp; // new mix zone volume Tank[j].A = (curve->Y[n] - curve->Y[0]) / // nominal area (curve->X[n] - curve->X[0]); break; @@ -2521,9 +2516,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu if (value > curve->X[n]) return 225; // new level is off curve } Tank[j].Hmax = hTmp; // new max. head - vTmp = Tank[j].Vmax; // old max. volume Tank[j].Vmax = tankvolume(p, j, hTmp); // new max. volume - Tank[j].V1max *= Tank[j].Vmax / vTmp; // new mix zone volume break; case EN_MIXMODEL: @@ -2542,7 +2535,7 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu j = index - nJuncs; if (Tank[j].A > 0.0) { - Tank[j].V1max = value * Tank[j].Vmax; + Tank[j].V1frac = value; } break; diff --git a/src/inpfile.c b/src/inpfile.c index 39e361f..c8eb1a9 100644 --- a/src/inpfile.c +++ b/src/inpfile.c @@ -523,7 +523,7 @@ int saveinpfile(Project *pr, const char *fname) tank = &net->Tank[i]; if (tank->A == 0.0) continue; fprintf(f, "\n %-31s %-8s %12.4f", net->Node[tank->Node].ID, - MixTxt[tank->MixModel], (tank->V1max / tank->Vmax)); + MixTxt[tank->MixModel], tank->V1frac); } // Write [REACTIONS] section diff --git a/src/input1.c b/src/input1.c index 10a6f95..e79da83 100644 --- a/src/input1.c +++ b/src/input1.c @@ -588,7 +588,6 @@ void convertunits(Project *pr) tank->Kb /= SECperDAY; tank->V = tank->V0; tank->C = node->C0; - tank->V1max *= tank->Vmax; } // Convert hydraulic convergence criteria diff --git a/src/input3.c b/src/input3.c index 28570fc..e89d6e1 100644 --- a/src/input3.c +++ b/src/input3.c @@ -252,7 +252,7 @@ int tankdata(Project *pr) tank->Vcurve = curve; tank->MixModel = MIX1; // Completely mixed - tank->V1max = 1.0; // Mixing compartment size fraction + tank->V1frac = 1.0; // Mixing compartment size fraction return 0; } @@ -1288,7 +1288,7 @@ int mixingdata(Project *pr) i = j - net->Njuncs; if (net->Tank[i].A == 0.0) return 0; net->Tank[i].MixModel = (char)m; - net->Tank[i].V1max = v; + net->Tank[i].V1frac = v; return 0; } diff --git a/src/qualreact.c b/src/qualreact.c index 1042527..1fc5710 100644 --- a/src/qualreact.c +++ b/src/qualreact.c @@ -525,7 +525,7 @@ void tankmix2(Project *pr, int i, double vin, double win, double vnet) if (mixzone == NULL || stagzone == NULL) return; // Full mixing zone volume - vmz = tank->V1max; + vmz = tank->V1frac * tank->Vmax; // Tank is filling vt = 0.0; diff --git a/src/qualroute.c b/src/qualroute.c index 592f92f..b71f218 100644 --- a/src/qualroute.c +++ b/src/qualroute.c @@ -612,7 +612,7 @@ void initsegs(Project *pr) if (net->Tank[j].MixModel == MIX2) { // ... mixing zone segment - v1 = MAX(0, v - net->Tank[j].V1max); + v1 = MAX(0, v - net->Tank[j].V1frac * net->Tank[j].Vmax); qual->FirstSeg[k]->v = v1; // ... stagnant zone segment diff --git a/src/types.h b/src/types.h index 50dc50a..902a22d 100755 --- a/src/types.h +++ b/src/types.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 10/29/2019 + Last Updated: 07/11/2020 ****************************************************************************** */ @@ -428,7 +428,7 @@ typedef struct // Tank Object int Pat; // fixed grade time pattern int Vcurve; // volume v. elev. curve index MixType MixModel; // type of mixing model - double V1max; // mixing compartment size + double V1frac; // mixing compartment fraction int CanOverflow; // tank can overflow or not } Stank; From 87d1b2ed7943083bb3665393829cd4a2f9107faf Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Wed, 15 Jul 2020 14:58:58 -0400 Subject: [PATCH 8/8] Update modules.dox --- doc/modules.dox | 2 ++ 1 file changed, 2 insertions(+) diff --git a/doc/modules.dox b/doc/modules.dox index ca2ba4f..249c0ec 100644 --- a/doc/modules.dox +++ b/doc/modules.dox @@ -197,6 +197,7 @@ These are the toolkit's enumerated types whose members are used as function argu @fn int EN_setheadcurveindex(EN_Project ph, int pumpIndex, int curveIndex) @fn int EN_getvertexcount(EN_Project ph, int index, int *count) @fn int EN_getvertex(EN_Project ph, int index, int vertex, double *x, double *y) +@fn int EN_setvertex(EN_Project ph, int index, int vertex, double x, double y) @fn int EN_setvertices(EN_Project ph, int index, double *x, double *y, int count) @} */ @@ -227,6 +228,7 @@ These are the toolkit's enumerated types whose members are used as function argu @fn int EN_setcurveid(EN_Project ph, int index, char *id) @fn int EN_getcurvelen(EN_Project ph, int index, int *len) @fn int EN_getcurvetype(EN_Project ph, int index, int *type) +@fn int EN_setcurvetype(EN_Project ph, int index, int type) @fn int EN_getcurvevalue(EN_Project ph, int curveIndex, int pointIndex, double *x, double *y) @fn int EN_setcurvevalue(EN_Project ph, int curveIndex, int pointIndex, double x, double y) @fn int EN_getcurve(EN_Project ph, int curveIndex, char* id, int *nPoints, double **xValues, double **yValues)