From a89f339525bd0774b35fe43b9ff4e873723c17bc Mon Sep 17 00:00:00 2001 From: Lew Rossman Date: Mon, 22 Jul 2019 09:50:41 -0400 Subject: [PATCH] PDA fixes --- ReleaseNotes2_2.md | 47 +++-- include/epanet2.bas | 7 +- include/epanet2.pas | 417 ++++++++++++++++++++++++++++++++++++++++ include/epanet2.vb | 7 +- include/epanet2_2.h | 36 ++-- include/epanet2_enums.h | 17 +- src/epanet.c | 46 +++-- src/errors.dat | 1 + src/funcs.h | 7 +- src/hydcoeffs.c | 128 ++++-------- src/hydraul.c | 3 +- src/hydsolver.c | 99 ++++++++-- src/input1.c | 15 +- src/input3.c | 8 + src/report.c | 93 +++++---- src/text.h | 4 +- src/types.h | 10 +- tests/test_pda.cpp | 90 +++++++++ 18 files changed, 814 insertions(+), 221 deletions(-) create mode 100644 include/epanet2.pas create mode 100644 tests/test_pda.cpp diff --git a/ReleaseNotes2_2.md b/ReleaseNotes2_2.md index 0e6dbad..189bbda 100644 --- a/ReleaseNotes2_2.md +++ b/ReleaseNotes2_2.md @@ -12,7 +12,7 @@ one would use: `EN_getnodevalue(ph, nodeIndex, EN_ELEVATION, &elev)` -where `ph` is the handle assigned to the project. +where `ph` is the handle assigned to the project. Two new functions have been added to the API to manage the creation and deletion of project handles. `EN_createproject` creates a new project along with its handle, while `EN_deleteproject` deletes a project. An example of using the thread-safe version of the API is shown below: ``` @@ -47,32 +47,32 @@ int buildandrunEpanet(char *rptfile) err = EN_createproject(&ph); if (err) return err; EN_init(ph, rptfile, "", EN_GPM, EN_HW); - + //Add a junction node with 710 ft elevation and 500 gpm demand EN_addnode(ph, "J1", EN_JUNCTION, &index); EN_setjuncdata(ph, index, 710, 500, ""); - + // Add a reservoir node at 800 ft elevation EN_addnode(ph, "R1", EN_RESERVOIR, &index); EN_setnodevalue(ph, index, EN_ELEVATION, 800); - + // Add a 5280 ft long, 14-inch pipe with C-factor of 100 // from the reservoir to the demand node EN_addlink(ph, "P1", EN_PIPE, "R1", "J1", &index); EN_setpipedata(ph, index, 5280, 14, 100, 0); - + // Solve for hydraulics and report nodal results EN_setreport(ph, "NODES ALL"); err = EN_solveH(ph); if (!err) err = EN_report(ph); - + // Close and delete the project EN_close(ph); EN_deleteproject(ph); return err; } ``` -Instead of using `EN_open` to load data from a file, `EN_init` is used to initialize a project with the names of its report and binary files, and its flow units and head loss formula. The legacy-style API has a complementary set of functions for building networks from code. +Instead of using `EN_open` to load data from a file, `EN_init` is used to initialize a project with the names of its report and binary files, and its flow units and head loss formula. The legacy-style API has a complementary set of functions for building networks from code. ## Additional Convergence Parameters @@ -82,13 +82,13 @@ Two new analysis options have been added to provide more rigorous convergence cr `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. +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 `EN_getstatistic` (or `ENgetstatistic`) function to retrieve their computed values (even when their option values are 0) after a hydraulic solution has been completed. ## More Efficient Node Re-ordering 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 efficient **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, where most of the work is for node-reordering, by an average of 55%. Since MMD did not reduce the time needed to solve the linear equations generated at each iteration of the hydraulic solver longer extended period simulations will not exhibit a similar speedup. +EPANET's original node re-ordering scheme has been replaced by the more efficient **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, where most of the work is for node-reordering, by an average of 55%. Since MMD did not reduce the time needed to solve the linear equations generated at each iteration of the hydraulic solver longer extended period simulations will not exhibit a similar speedup. ## Improved Handling of Near-Zero Flows @@ -98,7 +98,7 @@ The hydraulic solver has been modified to use a linear head loss v. flow relatio ## 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. +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: @@ -110,10 +110,10 @@ To implement pressure driven analysis four new parameters have been added to the | 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 | +| `DEMAND MODEL` | either DDA or PDA | DDA | +| `MINIMUM PRESSURE` | value for Pmin | 0 +| `REQUIRED PRESSURE` | value for Preq | 0.1 +| `PRESSURE EXPONENT` | value for Pexp | 0.5 | These parameters can also be set and retrieved in code using the following API functions ``` @@ -126,10 +126,12 @@ int EN_setdemandmodel(EN_Project ph, int modelType, double pMin, double pReq, do int EN_getdemandmodel(EN_Project ph, 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 left as a future feature to implement. - - Pmin 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. + - Preq must be at least 0.1 (either psi or m) higher than Pmin to avoid numerical issues caused by having too steep a demand curve. + - Using `EN_DEFICIENTNODES` as the argument to `EN_getstatistic` (or `ENgetstatistic`) will retrieve the number of nodes that are pressure deficient. These are nodes with positive required demand whose pressure is below 0 under **DDA** or below Preq under **PDA**. + - Using `EN_DEMANDREDUCTION` as an argument will retrieve the total percent reduction of demands at pressure deficient nodes under **PDA**. + - Using `EN_DEMANDDEFICIT` with the `EN_getnodevalue` (or `ENgetnodevalue`) function will return the amount of demand reduction produced by a **PDA** at any particular node. ## Tank Overflows EPANET has always prevented tanks from overflowing by closing any links that supply inflow to a full tank. A new option `EN_CANOVERFLOW`, has been added to the list of Tank node properties. When set to 1 it will allow its tank to overflow when it becomes full. The spillage rate is returned in the tank's EN_DEMAND property. The default value for `EN_CANOVERFLOW` is 0 indicating that the tank cannot overflow. @@ -163,7 +165,7 @@ With this change EPANET 2.2 now produces perfect mass balances when tested again |`EN_deleterule`|Deletes a rule-based control from the project| |`EN_setnodeid`|Changes the ID name for a node| |`EN_setjuncdata` |Sets values for a junction's parameters | -|`EN_settankdata` |Sets values for a tank's parameters| +|`EN_settankdata` |Sets values for a tank's parameters| |`EN_setlinkid`|Changes the ID name for a link| |`EN_setlinknodes`|Sets a link's start- and end-nodes| |`EN_setlinktype`|Changes the type of a specific link| @@ -205,7 +207,7 @@ In addition to these new functions, a tank's volume curve `EN_VOLCURVE` can be s - `EN_PUMP_ECOST` (average energy price) - `EN_PUMP_EPAT` (energy pricing pattern) - `EN_LINKPATTERN` (speed setting pattern) - + Access to the following global energy options have been added to `EN_getoption` and `EN_setoption`: - `EN_GLOBALEFFIC` (global pump efficiency) - `EN_GLOBALPRICE` (global average energy price per kW-Hour) @@ -215,6 +217,7 @@ Access to the following global energy options have been added to `EN_getoption` ## New API Constants ### Node value types: - `EN_CANOVERFLOW` +- `EN_DEMANDDEFICIT` ### Link value types: - `EN_PUMP_STATE` @@ -254,11 +257,13 @@ Access to the following global energy options have been added to `EN_getoption` - `EN_WALLORDER` - `EN_TANKORDER` - `EN_CONCENLIMIT` - + ### Simulation statistic types: - `EN_MAXHEADERROR` - `EN_MAXFLOWCHANGE` - `EN_MASSBALANCE` + - `EN_DEFICIENTNODES` + - `EN_DEMANDREDUCTION` ### Action code types: - `EN_UNCONDITIONAL` @@ -276,7 +281,7 @@ Access to the following global energy options have been added to `EN_getoption` - `EN_PDA` ## Documentation -Doxygen files have been created to generate a complete Users Guide for version 2.2's API. The guide's format is similar to the original EPANET Programmer's Toolkit help file and can be produced as a set of HTML pages, a Windows help file or a PDF document. - +Doxygen files have been created to generate a complete Users Guide for version 2.2's API. The guide's format is similar to the original EPANET Programmer's Toolkit help file and can be produced as a set of HTML pages, a Windows help file or a PDF document. + ## Authors contributing to this release: - List item diff --git a/include/epanet2.bas b/include/epanet2.bas index 391b0e9..37f2804 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 03/17/2019 +'Last updated on 07/18/2019 ' These are codes used by the DLL functions Public Const EN_ELEVATION = 0 ' Node parameters @@ -35,7 +35,8 @@ Public Const EN_MIXFRACTION = 22 Public Const EN_TANK_KBULK = 23 Public Const EN_TANKVOLUME = 24 Public Const EN_MAXVOLUME = 25 -Public Const EN_CANOVERFLOW = 26 +Public Const EN_CANOVERFLOW = 26 +Public Const EN_DEMANDDEFICIT = 27 Public Const EN_DIAMETER = 0 ' Link parameters Public Const EN_LENGTH = 1 @@ -83,6 +84,8 @@ Public Const EN_RELATIVEERROR = 1 Public Const EN_MAXHEADERROR = 2 Public Const EN_MAXFLOWCHANGE = 3 Public Const EN_MASSBALANCE = 4 +Public Const EN_DEFICIENTNODES = 5 +Public Const EN_DEMANDREDUCTION = 6 Public Const EN_NODE = 0 ' Component types Public Const EN_LINK = 1 diff --git a/include/epanet2.pas b/include/epanet2.pas new file mode 100644 index 0000000..830a29c --- /dev/null +++ b/include/epanet2.pas @@ -0,0 +1,417 @@ +unit epanet2; + +{ Declarations of imported procedures from the EPANET PROGRAMMERs TOOLKIT } +{ (EPANET2.DLL) } + +{Last updated on 7/18/19} + +interface + +const + +{ These are codes used by the DLL functions } + EN_MAXID = 31; { Max. # characters in ID name } + EN_MAXMSG = 255; { Max. # characters in strings } + + EN_ELEVATION = 0; { Node parameters } + EN_BASEDEMAND = 1; + EN_PATTERN = 2; + EN_EMITTER = 3; + EN_INITQUAL = 4; + EN_SOURCEQUAL = 5; + EN_SOURCEPAT = 6; + EN_SOURCETYPE = 7; + EN_TANKLEVEL = 8; + EN_DEMAND = 9; + EN_HEAD = 10; + EN_PRESSURE = 11; + EN_QUALITY = 12; + EN_SOURCEMASS = 13; + EN_INITVOLUME = 14; + EN_MIXMODEL = 15; + EN_MIXZONEVOL = 16; + + EN_TANKDIAM = 17; + EN_MINVOLUME = 18; + EN_VOLCURVE = 19; + EN_MINLEVEL = 20; + EN_MAXLEVEL = 21; + EN_MIXFRACTION = 22; + EN_TANK_KBULK = 23; + EN_TANKVOLUME = 24; + EN_MAXVOLUME = 25; + EN_CANOVERFLOW = 26; + EN_DEMANDDEFICIT = 27; + + EN_DIAMETER = 0; { Link parameters } + EN_LENGTH = 1; + EN_ROUGHNESS = 2; + EN_MINORLOSS = 3; + EN_INITSTATUS = 4; + EN_INITSETTING = 5; + EN_KBULK = 6; + EN_KWALL = 7; + EN_FLOW = 8; + EN_VELOCITY = 9; + EN_HEADLOSS = 10; + EN_STATUS = 11; + EN_SETTING = 12; + EN_ENERGY = 13; + EN_LINKQUAL = 14; + EN_LINKPATTERN = 15; + EN_PUMP_STATE = 16; + EN_PUMP_EFFIC = 17; + EN_PUMP_POWER = 18; + EN_PUMP_HCURVE = 19; + EN_PUMP_ECURVE = 20; + EN_PUMP_ECOST = 21; + EN_PUMP_EPAT = 22; + + EN_DURATION = 0; { Time parameters } + EN_HYDSTEP = 1; + EN_QUALSTEP = 2; + EN_PATTERNSTEP = 3; + EN_PATTERNSTART = 4; + EN_REPORTSTEP = 5; + EN_REPORTSTART = 6; + EN_RULESTEP = 7; + EN_STATISTIC = 8; + EN_PERIODS = 9; + EN_STARTTIME = 10; + EN_HTIME = 11; + EN_QTIME = 12; + EN_HALTFLAG = 13; + EN_NEXTEVENT = 14; + EN_NEXTEVENTTANK = 15; + + EN_ITERATIONS = 0; { Analysis statistics } + EN_RELATIVEERROR = 1; + EN_MAXHEADERROR = 2; + EN_MAXFLOWCHANGE = 3; + EN_MASSBALANCE = 4; + EN_DEFICIENTNODES = 5; + EN_DEMANDREDUCTION = 6; + + EN_NODE = 0; { Component Types } + EN_LINK = 1; + EN_TIMEPAT = 2; + EN_CURVE = 3; + EN_CONTROL = 4; + EN_RULE = 5; + + EN_NODECOUNT = 0; { Component counts } + EN_TANKCOUNT = 1; + EN_LINKCOUNT = 2; + EN_PATCOUNT = 3; + EN_CURVECOUNT = 4; + EN_CONTROLCOUNT = 5; + EN_RULECOUNT = 6; + + EN_JUNCTION = 0; { Node types } + EN_RESERVOIR = 1; + EN_TANK = 2; + + EN_CVPIPE = 0; { Link types } + EN_PIPE = 1; + EN_PUMP = 2; + EN_PRV = 3; + EN_PSV = 4; + EN_PBV = 5; + EN_FCV = 6; + EN_TCV = 7; + EN_GPV = 8; + + EN_CLOSED = 0; { Link status types } + EN_OPEN = 1; + + EN_PUMP_XHEAD = 0; { Pump state types } + EN_PUMP_CLOSED = 2; + EN_PUMP_OPEN = 3; + EN_PUMP_XFLOW = 5; + + EN_NONE = 0; { Quality analysis types } + EN_CHEM = 1; + EN_AGE = 2; + EN_TRACE = 3; + + EN_CONCEN = 0; { Source quality types } + EN_MASS = 1; + EN_SETPOINT = 2; + EN_FLOWPACED = 3; + + EN_HW = 0; { Head loss formulas } + EN_DW = 1; + EN_CM = 2; + + EN_CFS = 0; { Flow units types } + EN_GPM = 1; + EN_MGD = 2; + EN_IMGD = 3; + EN_AFD = 4; + EN_LPS = 5; + EN_LPM = 6; + EN_MLD = 7; + EN_CMH = 8; + EN_CMD = 9; + + EN_DDA = 0; { Demand model types } + EN_PDA = 1; + + EN_TRIALS = 0; { Option types } + EN_ACCURACY = 1; + EN_TOLERANCE = 2; + EN_EMITEXPON = 3; + EN_DEMANDMULT = 4; + EN_HEADERROR = 5; + EN_FLOWCHANGE = 6; + EN_HEADLOSSFORM = 7; + EN_GLOBALEFFIC = 8; + EN_GLOBALPRICE = 9; + EN_GLOBALPATTERN = 10; + EN_DEMANDCHARGE = 11; + EN_SP_GRAVITY = 12; + EN_SP_VISCOS = 13; + EN_EXTRA_ITER = 14; + EN_CHECKFREQ = 15; + EN_MAXCHECK = 16; + EN_DAMPLIMIT = 17; + EN_SP_DIFFUS = 18; + EN_BULKORDER = 19; + EN_WALLORDER = 20; + EN_TANKORDER = 21; + EN_CONCENLIMIT = 22; + + EN_LOWLEVEL = 0; { Control types } + EN_HILEVEL = 1; + EN_TIMER = 2; + EN_TIMEOFDAY = 3; + + EN_SERIES = 0; { Report statistic types } + EN_AVERAGE = 1; + EN_MINIMUM = 2; + EN_MAXIMUM = 3; + EN_RANGE = 4; + + EN_MIX1 = 0; { Tank mixing models } + EN_MIX2 = 1; + EN_FIFO = 2; + EN_LIFO = 3; + + EN_NOSAVE = 0; { Hydraulics flags } + EN_SAVE = 1; + EN_INITFLOW = 10; + EN_SAVE_AND_INIT = 11; + + EN_CONST_HP = 0; { Pump curve types } + EN_POWER_FUNC = 1; + EN_CUSTOM = 2; + EN_NOCURVE = 3; + + EN_VOLUME_CURVE = 0; { Curve types } + EN_PUMP_CURVE = 1; + EN_EFFIC_CURVE = 2; + EN_HLOSS_CURVE = 3; + EN_GENERIC_CURVE = 4; + + EN_UNCONDITIONAL = 0; { Deletion action codes } + EN_CONDITIONAL = 1; + + EN_NO_REPORT = 0; { Status reporting levels } + EN_NORMAL_REPORT = 1; + EN_FULL_REPORT = 2; + + EN_R_NODE = 6; { Rule-based control objects } + EN_R_LINK = 7; + EN_R_SYSTEM = 8; + + EN_R_DEMAND = 0; { Rule-based control variables } + EN_R_HEAD = 1; + EN_R_GRADE = 2; + EN_R_LEVEL = 3; + EN_R_PRESSURE = 4; + EN_R_FLOW = 5; + EN_R_STATUS = 6; + EN_R_SETTING = 7; + EN_R_POWER = 8; + EN_R_TIME = 9; + EN_R_CLOCKTIME = 10; + EN_R_FILLTIME = 11; + EN_R_DRAINTIME = 12; + + EN_R_EQ = 0; { Rule-based control operators } + EN_R_NE = 1; + EN_R_LE = 2; + EN_R_GE = 3; + EN_R_LT = 4; + EN_R_GT = 5; + EN_R_IS = 6; + EN_R_NOT = 7; + EN_R_BELOW = 8; + EN_R_ABOVE = 9; + + EN_R_IS_OPEN = 1; { Rule-based control link status } + EN_R_IS_CLOSED = 2; + EN_R_IS_ACTIVE = 3; + + EpanetLib = 'epanet2.dll'; + +{Project Functions} + function ENepanet(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar; F4: Pointer): Integer; stdcall; external EpanetLib; + function ENinit(F2: PAnsiChar; F3: PAnsiChar; UnitsType: Integer; HeadlossType: Integer): Integer; stdcall; external EpanetLib; + function ENopen(F1: PAnsiChar; F2: PAnsiChar; F3: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetcount(Code: Integer; var Count: Integer): Integer; stdcall; external EpanetLib; + function ENgettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsettitle(Line1: PAnsiChar; Line2: PAnsiChar; Line3: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsetcomment(ObjType: Integer; Index: Integer; Comment: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsaveinpfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENclose: Integer; stdcall; external EpanetLib; + +{Hydraulic Functions} + function ENsolveH: Integer; stdcall; external EpanetLib; + 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 ENcloseH: Integer; stdcall; external EpanetLib; + function ENsavehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENusehydfile(F: PAnsiChar): Integer; stdcall; external EpanetLib; + +{Quality Functions} + 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 ENcloseQ: Integer; stdcall; external EpanetLib; + +{Reporting Functions} + function ENwriteline(S: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENreport: Integer; stdcall; external EpanetLib; + function ENcopyreport(F: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENclearreport: Integer; stdcall; external EpanetLib; + function ENresetreport: Integer; stdcall; external EpanetLib; + function ENsetreport(S: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsetstatusreport(Code: Integer): Integer; stdcall; external EpanetLib; + function ENgetversion(var Version: Integer): Integer; stdcall; external EpanetLib; + function ENgeterror(Errcode: Integer; Errmsg: PAnsiChar; MaxLen: Integer): Integer; stdcall; external EpanetLib; + function ENgetstatistic(StatType: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + +{Analysis Options Functions} + function ENgetoption(Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + 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 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; + +{Node Functions} + function ENaddnode(ID: PAnsiChar; NodeType: Integer): Integer; stdcall; external EpanetLib; + function ENdeletenode(Index: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; + function ENgetnodeindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetnodeid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsetnodeid(Index: Integer; NewID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetnodetype(Index: Integer; var Code: Integer): Integer; stdcall; external EpanetLib; + function ENgetnodevalue(Index: Integer; Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + function ENsetnodevalue(Index: Integer; Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; + function ENsetjuncdata(Index: Integer; Elev: Single; Dmnd: Single; DmndPat: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsettankdata(Index: Integer; Elev, InitLvl, MinLvl, MaxLvl, Diam, MinVol: Single; VolCurve: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetcoord(Index: Integer; var X: Double; var Y: Double): Integer; stdcall; external EpanetLib; + function ENsetcoord(Index: Integer; X: Double; Y: Double): Integer; stdcall; external EpanetLib; + +{Demand Functions} + function ENgetdemandmodel(var Model: Integer; var Pmin: Single; var Preq: Single; var Pexp: Single): Integer; stdcall; external EpanetLib; + function ENsetdemandmodel(Model: Integer; Pmin: Single; Preq: Single; Pexp: Single): Integer; stdcall; external EpanetLib; + function ENgetnumdemands(NodeIndex: Integer; var NumDemands: Integer): Integer; stdcall; external EpanetLib; + function ENadddemand(NodeIndex: Integer; BaseDemand: Single; PatIndex: Integer; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENdeletedemand(NodeIndex: Integer; DemandIndex: Integer): Integer; stdcall; external EpanetLib; + function ENgetdemandindex(NodeIndex: Integer; DemandName: PAnsiString; var DemandIndex: Integer): Integer; stdcall; external EpanetLib; + function ENgetbasedemand(NodeIndex: Integer; DemandIndex: Integer; var BaseDemand: Single): Integer; stdcall; external EpanetLib; + function ENsetbasedemand(NodeIndex: Integer; DemandIndex: Integer; BaseDemand: Single): Integer; stdcall; external EpanetLib; + function ENgetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; var PatIndex: Integer): Integer; stdcall; external EpanetLib; + function ENsetdemandpattern(NodeIndex: Integer; DemandIndex: Integer; PatIndex: Integer): Integer; stdcall; external EpanetLib; + function ENgetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsetdemandname(NodeIndex: Integer; DemandIndex: Integer; DemandName: PAnsiChar): Integer; stdcall; external EpanetLib; + +{Link Functions} + function ENaddlink(ID: PAnsiChar; LinkType: Integer; FromNode: PAnsiChar; ToNode: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENdeletelink(Index: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; + function ENgetlinkindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetlinkid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsetlinkid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetlinktype(Index: Integer; var Code: Integer): Integer; stdcall; external EpanetLib; + function ENsetlinktype(var Index: Integer; LinkType: Integer; ActionCode: Integer): Integer; stdcall; external EpanetLib; + function ENgetlinknodes(Index: Integer; var Node1: Integer; var Node2: Integer): Integer; stdcall; external EpanetLib; + function ENsetlinknodes(Index: Integer; Node1: Integer; Node2: Integer): Integer; stdcall; external EpanetLib; + function ENgetlinkvalue(Index: Integer; Code: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + function ENsetlinkvalue(Index: Integer; Code: Integer; Value: Single): Integer; stdcall; external EpanetLib; + function ENsetpipedata(Index: Integer; Length: Single; Diam: Single; Rough: Single; Mloss:Single): Integer; stdcall; external EpanetLib; + +{Pump Functions} + function ENgetpumptype(LinkIndex: Integer; var PumpType: Integer): Integer; stdcall; external EpanetLib; + function ENgetheadcurveindex(LinkIndex: Integer; var CurveIndex: Integer): Integer; stdcall; external EpanetLib; + function ENsetheadcurveindex(LinkIndex: Integer; CurveIndex: Integer): Integer; stdcall; external EpanetLib; + +{Pattern Functions} + function ENaddpattern(ID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENdeletepattern(Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetpatternindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetpatternid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENsetpatternid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENgetpatternlen(Index: Integer; var Len: Integer): Integer; stdcall; external EpanetLib; + function ENgetpatternvalue(Index: Integer; Period: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + function ENsetpatternvalue(Index: Integer; Period: Integer; Value: Single): Integer; stdcall; external EpanetLib; + function ENgetaveragepatternvalue(Index: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + function ENsetpattern(Index: Integer; F: array of Single; N: Integer): Integer; stdcall; external EpanetLib; + +{Curve Functions} + function ENaddcurve(ID: PAnsiChar): Integer; stdcall; external EpanetLib; + function ENdeletecurve(Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetcurveindex(ID: PAnsiChar; var Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetcurveid(Index: Integer; ID: PAnsiChar): Integer; stdcall; external EpanetLib; + 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 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: array of Single; var Y: array of Single): Integer; stdcall; external EpanetLib; + function ENsetcurve(Index: Integer; X: array of Single; Y: array of Single; N: Integer): Integer; stdcall; external EpanetLib; + +{Control Functions} + function ENaddcontrol(Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single; var Index: Integer): Integer; stdcall; external EpanetLib; + function ENdeletecontrol(Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetcontrol(Index: Integer; var Ctype: Integer; var Link: Integer; var Setting: Single; var Node: Integer; var Level: Single): Integer; stdcall; external EpanetLib; + function ENsetcontrol(Index: Integer; Ctype: Integer; Link: Integer; Setting: Single; Node: Integer; Level: Single): Integer; stdcall; external EpanetLib; + + {Rule-Based Control Functions} + function ENaddrule(Rule: PAnsiString): Integer; stdcall; external EpanetLib; + function ENdeleterule(Index: Integer): Integer; stdcall; external EpanetLib; + function ENgetrule(Index: Integer; var Npremises: Integer; var NthenActions: Integer; + var NelseActions: Integer; var Priority: Single): Integer; stdcall; external EpanetLib; + function ENgetruleID(Index: Integer; ID: PAnsiString): Integer; stdcall; external EpanetLib; + function ENsetrulepriority(Index: Integer; Priority: Single): Integer; stdcall; external EpanetLib; + function ENgetpremise(RuleIndex: Integer; PremiseIndex: Integer; var LogOp: Integer; + var ObjType: Integer; var ObjIndex: Integer; var Param: Integer; var RelOp: Integer; + var Status: Integer; var Value: Single): Integer; stdcall; external EpanetLib; + function ENsetpremise(RuleIndex: Integer; PremiseIndex: Integer; LogOp: Integer; ObjType: Integer; + ObjIndex: Integer; Param: Integer; RelOp: Integer; Status: Integer; Value: Single): Integer; stdcall; external EpanetLib; + function ENsetpremiseindex(RuleIndex: Integer; PremiseIndex: Integer; ObjIndex: Integer): Integer; stdcall; external EpanetLib; + function ENsetpremisestatus(RuleIndex: Integer; PremiseIndex: Integer; Status: Integer): Integer; stdcall; external EpanetLib; + function ENsetpremisevalue(RuleIndex: Integer; PremiseIndex: Integer; Value: Single): Integer; stdcall; external EpanetLib; + function ENgetthenaction(RuleIndex: Integer; ActionIndex: Integer; var LinkIndex: Integer; + var Status: Integer; var Setting: Single): Integer; stdcall; external EpanetLib; + function ENsetthenaction(RuleIndex: Integer; ActionIndex: Integer; LinkIndex: Integer; + Status: Integer; Setting: Single): Integer; stdcall; external EpanetLib; + function ENgetelseaction(RuleIndex: Integer; ActionIndex: Integer; var LinkIndex: Integer; + var Status: Integer; var Setting: Single): Integer; stdcall; external EpanetLib; + function ENsetelseaction(RuleIndex: Integer; ActionIndex: Integer; LinkIndex: Integer; + Status: Integer; Setting: Single): Integer; stdcall; external EpanetLib; + +implementation + +end. diff --git a/include/epanet2.vb b/include/epanet2.vb index 347233f..92622e3 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 03/17/2019 +'Last updated on 07/18/2019 Imports System.Runtime.InteropServices Imports System.Text @@ -40,7 +40,8 @@ Public Const EN_TANK_KBULK = 23 Public Const EN_TANKVOLUME = 24 Public Const EN_MAXVOLUME = 25 -Public Const EN_CANOVERFLOW = 26 +Public Const EN_CANOVERFLOW = 26 +Public Const EN_DEMANDDEFICIT = 27 Public Const EN_DIAMETER = 0 ' Link parameters Public Const EN_LENGTH = 1 @@ -88,6 +89,8 @@ Public Const EN_RELATIVEERROR = 1 Public Const EN_MAXHEADERROR = 2 Public Const EN_MAXFLOWCHANGE = 3 Public Const EN_MASSBALANCE = 4 +Public Const EN_DEFICIENTNODES = 5 +Public Const EN_DEMANDREDUCTION = 6 Public Const EN_NODE = 0 ' Component types Public Const EN_LINK = 1 diff --git a/include/epanet2_2.h b/include/epanet2_2.h index a63fbe5..9ed1b5c 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: 05/30/2019 + Last Updated: 07/20/2019 ****************************************************************************** */ @@ -94,7 +94,7 @@ typedef struct Project *EN_Project; the pviewprog argument should be `NULL`. */ int DLLEXPORT EN_runproject(EN_Project ph, const char *inpFile, const char *rptFile, - const char *outputFile, void (*pviewprog)(char *)); + const char *outFile, void (*pviewprog)(char *)); /** @brief Initializes an EPANET project. @@ -129,9 +129,9 @@ typedef struct Project *EN_Project; /** @brief Retrieves the title lines of the project @param ph an EPANET project handle. - @param[out] line1 first title line - @param[out] line2 second title line - @param[out] line3 third title line + @param[out] out_line1 first title line + @param[out] out_line2 second title line + @param[out] out_line3 third title line @return an error code */ int DLLEXPORT EN_gettitle(EN_Project ph, char *out_line1, char *out_line2, char *out_line3); @@ -151,7 +151,7 @@ typedef struct Project *EN_Project; @param ph an EPANET project handle. @param object a type of object (either EN_NODE, EN_LINK, EN_TIMEPAT or EN_CURVE) @param index the object's index starting from 1 - @param[out] comment the comment string assigned to the object + @param[out] out_comment the comment string assigned to the object @return an error code */ int DLLEXPORT EN_getcomment(EN_Project ph, int object, int index, char *out_comment); @@ -223,7 +223,7 @@ typedef struct Project *EN_Project; EN_solveH(ph); EN_solveQ(ph); EN_report(ph); - EN_deleteproject(&ph); + EN_deleteproject(ph); \endcode */ int DLLEXPORT EN_solveH(EN_Project ph); @@ -621,7 +621,7 @@ typedef struct Project *EN_Project; /** @brief Returns the text of an error message generated by an error code. @param errcode an error code. - @param[out] errmsg the error message generated by the error code + @param[out] out_errmsg the error message generated by the error code @param maxLen maximum number of characters that errmsg can hold @return an error code @@ -707,8 +707,8 @@ typedef struct Project *EN_Project; @brief Gets information about the type of water quality analysis requested. @param ph an EPANET project handle. @param[out] qualType type of analysis to run (see @ref EN_QualityType). - @param[out] chemName name of chemical constituent. - @param[out] chemUnits concentration units of the constituent. + @param[out] out_chemName name of chemical constituent. + @param[out] out_chemUnits concentration units of the constituent. @param[out] traceNode index of the node being traced (if applicable). @return an error code. */ @@ -787,7 +787,7 @@ typedef struct Project *EN_Project; @brief Gets the ID name of a node given its index. @param ph an EPANET project handle. @param index a node's index (starting from 1). - @param[out] id the node's ID name. + @param[out] out_id the node's ID name. @return an error code The ID name must be sized to hold at least @ref EN_SizeLimits "EN_MAXID" characters. @@ -1032,7 +1032,7 @@ typedef struct Project *EN_Project; @param ph an EPANET project handle. @param nodeIndex a node's index (starting from 1). @param demandIndex the index of one of the node's demand categories (starting from 1). - @param[out] demandName The name of the selected category. + @param[out] out_demandName The name of the selected category. @return an error code. \b demandName must be sized to contain at least @ref EN_SizeLimits "EN_MAXID" characters. @@ -1111,7 +1111,7 @@ typedef struct Project *EN_Project; @brief Gets the ID name of a link given its index. @param ph an EPANET project handle. @param index a link's index (starting from 1). - @param[out] id The link's ID name. + @param[out] out_id The link's ID name. @return an error code. The ID name must be sized to hold at least @ref EN_SizeLimits "EN_MAXID" characters. @@ -1141,7 +1141,7 @@ typedef struct Project *EN_Project; /** @brief Changes a link's type. @param ph an EPANET project handle. - @param[in,out] index the link's index before [in] and after [out] the type change. + @param[in,out] inout_index the link's index before [in] and after [out] the type change. @param linkType the new type to change the link to (see @ref EN_LinkType). @param actionCode the action taken if any controls contain the link. @return an error code. @@ -1283,7 +1283,7 @@ typedef struct Project *EN_Project; @brief Retrieves the ID name of a time pattern given its index. @param ph an EPANET project handle. @param index a time pattern index (starting from 1). - @param[out] id the time pattern's ID name. + @param[out] out_id the time pattern's ID name. @return an error code. The ID name must be sized to hold at least @ref EN_SizeLimits "EN_MAXID" characters. @@ -1391,7 +1391,7 @@ typedef struct Project *EN_Project; @brief Retrieves the ID name of a curve given its index. @param ph an EPANET project handle. @param index a curve's index (starting from 1). - @param[out] id the curve's ID name. + @param[out] out_id the curve's ID name. @return an error code. The ID name must be sized to hold at least @ref EN_SizeLimits "EN_MAXID" characters. @@ -1455,7 +1455,7 @@ typedef struct Project *EN_Project; @brief Retrieves all of a curve's data. @param ph an EPANET project handle. @param index a curve's index (starting from 1). - @param[out] id the curve's ID name. + @param[out] out_id the curve's ID name. @param[out] nPoints the number of data points on the curve. @param[out] xValues the curve's x-values. @param[out] yValues the curve's y-values. @@ -1590,7 +1590,7 @@ typedef struct Project *EN_Project; @brief Gets the ID name of a rule-based control given its index. @param ph an EPANET project handle. @param index the rule's index (starting from 1). - @param[out] id the rule's ID name. + @param[out] out_id the rule's ID name. @return Error code. The ID name must be sized to hold at least @ref EN_SizeLimits "EN_MAXID" characters. diff --git a/include/epanet2_enums.h b/include/epanet2_enums.h index f7b5cb2..71bb797 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: 05/30/2019 + Last Updated: 07/18/2019 ****************************************************************************** */ @@ -62,7 +62,8 @@ typedef enum { EN_TANK_KBULK = 23, //!< Tank bulk decay coefficient EN_TANKVOLUME = 24, //!< Current computed tank volume (read only) EN_MAXVOLUME = 25, //!< Tank maximum volume (read only) - EN_CANOVERFLOW = 26 //!< Tank can overflow (= 1) or not (= 0) + EN_CANOVERFLOW = 26, //!< Tank can overflow (= 1) or not (= 0) + EN_DEMANDDEFICIT = 27 //!< Amount that full demand is reduced under PDA (read only) } EN_NodeProperty; /// Link properties @@ -128,11 +129,13 @@ and the cumulative water quality mass balance error at the current simulation ti can be retrieved with @ref EN_getstatistic. */ typedef enum { - EN_ITERATIONS = 0, //!< Number of hydraulic iterations taken - EN_RELATIVEERROR = 1, //!< Sum of link flow changes / sum of link flows - EN_MAXHEADERROR = 2, //!< Largest head loss error for links - EN_MAXFLOWCHANGE = 3, //!< Largest flow change in links - EN_MASSBALANCE = 4 //!< Cumulative water quality mass balance ratio + EN_ITERATIONS = 0, //!< Number of hydraulic iterations taken + EN_RELATIVEERROR = 1, //!< Sum of link flow changes / sum of link flows + EN_MAXHEADERROR = 2, //!< Largest head loss error for links + EN_MAXFLOWCHANGE = 3, //!< Largest flow change in links + EN_MASSBALANCE = 4, //!< Cumulative water quality mass balance ratio + EN_DEFICIENTNODES = 5, //!< Number of pressure deficient nodes + EN_DEMANDREDUCTION = 6 //!< % demand reduction at pressure deficient nodes } EN_AnalysisStatistic; /// Types of network objects diff --git a/src/epanet.c b/src/epanet.c index 71d67f7..eb6b5cf 100644 --- a/src/epanet.c +++ b/src/epanet.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 05/15/2019 + Last Updated: 07/18/2019 ****************************************************************************** */ @@ -68,7 +68,6 @@ int DLLEXPORT EN_deleteproject(EN_Project p) remove(p->TmpOutFname); remove(p->TmpStatFname); free(p); - p = NULL; return 0; } @@ -1040,19 +1039,25 @@ int DLLEXPORT EN_getstatistic(EN_Project p, int type, double *value) switch (type) { case EN_ITERATIONS: - *value = (double)p->hydraul.Iterations; + *value = p->hydraul.Iterations; break; case EN_RELATIVEERROR: - *value = (double)p->hydraul.RelativeError; + *value = p->hydraul.RelativeError; break; case EN_MAXHEADERROR: - *value = (double)(p->hydraul.MaxHeadError * p->Ucf[HEAD]); + *value = p->hydraul.MaxHeadError * p->Ucf[HEAD]; break; case EN_MAXFLOWCHANGE: - *value = (double)(p->hydraul.MaxFlowChange * p->Ucf[FLOW]); + *value = p->hydraul.MaxFlowChange * p->Ucf[FLOW]; + break; + case EN_DEFICIENTNODES: + *value = p->hydraul.DeficientNodes; + break; + case EN_DEMANDREDUCTION: + *value = p->hydraul.DemandReduction; break; case EN_MASSBALANCE: - *value = (double)(p->quality.MassBalance.ratio); + *value = p->quality.MassBalance.ratio; break; default: *value = 0.0; @@ -2050,7 +2055,6 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val int nJuncs = net->Njuncs; double *Ucf = p->Ucf; - double *NodeDemand = hyd->NodeDemand; double *NodeHead = hyd->NodeHead; double *NodeQual = qual->NodeQual; @@ -2128,7 +2132,7 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val break; case EN_DEMAND: - v = NodeDemand[index] * Ucf[FLOW]; + v = hyd->NodeDemand[index] * Ucf[FLOW]; break; case EN_HEAD: @@ -2204,7 +2208,16 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val if (Node[index].Type != TANK) return 0; v = Tank[index - nJuncs].CanOverflow; break; - + + case EN_DEMANDDEFICIT: + if (index > nJuncs) return 0; + // After an analysis, DemandFlow contains node's required demand + // while NodeDemand contains delivered demand + emitter flow + if (hyd->DemandFlow[index] < 0.0) return 0; + v = (hyd->DemandFlow[index] - + (hyd->NodeDemand[index] - hyd->EmitterFlow[index])) * Ucf[FLOW]; + break; + default: return 251; } @@ -2702,9 +2715,9 @@ int DLLEXPORT EN_getdemandmodel(EN_Project p, int *model, double *pmin, */ { *model = p->hydraul.DemandModel; - *pmin = (double)(p->hydraul.Pmin * p->Ucf[PRESSURE]); - *preq = (double)(p->hydraul.Preq * p->Ucf[PRESSURE]); - *pexp = (double)(p->hydraul.Pexp); + *pmin = p->hydraul.Pmin * p->Ucf[PRESSURE]; + *preq = p->hydraul.Preq * p->Ucf[PRESSURE]; + *pexp = p->hydraul.Pexp; return 0; } @@ -2722,7 +2735,12 @@ int DLLEXPORT EN_setdemandmodel(EN_Project p, int model, double pmin, */ { if (model < 0 || model > EN_PDA) return 251; - if (pmin > preq || pexp <= 0.0) return 209; + if (model == EN_PDA) + { + if (pexp <= 0.0) return 208; + if (pmin < 0.0) return 208; + if (preq - pmin < MINPDIFF) return 208; + } p->hydraul.DemandModel = model; p->hydraul.Pmin = pmin / p->Ucf[PRESSURE]; p->hydraul.Preq = preq / p->Ucf[PRESSURE]; diff --git a/src/errors.dat b/src/errors.dat index cf808a8..6d3a142 100644 --- a/src/errors.dat +++ b/src/errors.dat @@ -22,6 +22,7 @@ DAT(204,"undefined link") DAT(205,"undefined time pattern") DAT(206,"undefined curve") DAT(207,"attempt to control CV/GPV link") +DAT(208,"illegal PDA pressure limits") DAT(209,"illegal node property value") DAT(211,"illegal link property value") DAT(212,"undefined trace node") diff --git a/src/funcs.h b/src/funcs.h index 687d3d3..e520b29 100755 --- a/src/funcs.h +++ b/src/funcs.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 05/15/2019 + Last Updated: 07/08/2019 ****************************************************************************** */ #ifndef FUNCS_H @@ -157,9 +157,8 @@ double tankgrade(Project *, int, double); void resistcoeff(Project *, int); void headlosscoeffs(Project *); void matrixcoeffs(Project *); -void emitheadloss(Project *, int, double *, double *); -double demandflowchange(Project *, int, double, double); -void demandparams(Project *, double *, double *); +void emitterheadloss(Project *, int, double *, double *); +void demandheadloss(Project *, int, double, double, double *, double *); // ------- QUALITY.C -------------------- diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c index d97bd93..07da0dc 100644 --- a/src/hydcoeffs.c +++ b/src/hydcoeffs.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 05/15/2019 + Last Updated: 07/08/2019 ****************************************************************************** */ @@ -38,9 +38,8 @@ const double CBIG = 1.e8; //void resistcoeff(Project *, int ); //void headlosscoeffs(Project *); //void matrixcoeffs(Project *); -//void emitheadloss(Project *, int, double *, double *); -//double demandflowchange(Project *, int, double, double); -//void demandparams(Project *, double *, double *); +//void emitterheadloss(Project *, int, double *, double *); +//void demandheadloss(Project *, int, double, double, double *, double *); // Local functions static void linkcoeffs(Project *pr); @@ -48,8 +47,6 @@ static void nodecoeffs(Project *pr); static void valvecoeffs(Project *pr); static void emittercoeffs(Project *pr); static void demandcoeffs(Project *pr); -static void demandheadloss(double d, double dfull, double dp, - double n, double *hloss, double *hgrad); static void pipecoeff(Project *pr, int k); static void DWpipecoeff(Project *pr, int k); @@ -370,7 +367,7 @@ void emittercoeffs(Project *pr) if (node->Ke == 0.0) continue; // Find emitter head loss and gradient - emitheadloss(pr, i, &hloss, &hgrad); + emitterheadloss(pr, i, &hloss, &hgrad); // Row of solution matrix row = sm->Row[i]; @@ -385,7 +382,7 @@ void emittercoeffs(Project *pr) } -void emitheadloss(Project *pr, int i, double *hloss, double *hgrad) +void emitterheadloss(Project *pr, int i, double *hloss, double *hgrad) /* **------------------------------------------------------------- ** Input: i = node index @@ -424,37 +421,6 @@ void emitheadloss(Project *pr, int i, double *hloss, double *hgrad) } -void demandparams(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. -**-------------------------------------------------------------- -*/ -{ - Hydraul *hyd = &pr->hydraul; - - // 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(Project *pr) /* **-------------------------------------------------------------- @@ -479,81 +445,57 @@ void demandcoeffs(Project *pr) n, // exponent in head loss v. demand function hloss, // head loss in supplying demand (ft) hgrad; // gradient of demand head loss (ft/cfs) - + // Get demand function parameters if (hyd->DemandModel == DDA) return; - demandparams(pr, &dp, &n); + dp = hyd->Preq - hyd->Pmin; + n = 1.0 / hyd->Pexp; // 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->DemandFlow[i], hyd->NodeDemand[i], dp, n, - &hloss, &hgrad); - + demandheadloss(pr, i, dp, n, &hloss, &hgrad); + // Update row of solution matrix A & its r.h.s. F - row = sm->Row[i]; - sm->Aii[row] += 1.0 / hgrad; - sm->F[row] += (hloss + net->Node[i].El + hyd->Pmin) / hgrad; + if (hgrad > 0.0) + { + row = sm->Row[i]; + sm->Aii[row] += 1.0 / hgrad; + sm->F[row] += (hloss + net->Node[i].El + hyd->Pmin) / hgrad; + } } } - -double demandflowchange(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 -**-------------------------------------------------------------- -*/ -{ - Hydraul *hyd = &pr->hydraul; - - double hloss, hgrad; - - demandheadloss(hyd->DemandFlow[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, +void demandheadloss(Project *pr, int i, 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) +** Input: i = junction index +** dp = pressure range for demand function (ft) +** n = exponent in head v. demand function +** Output: hloss = pressure dependent demand head loss (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; + Hydraul *hyd = &pr->hydraul; + + const double EPS = 0.01; + double d = hyd->DemandFlow[i]; + double dfull = hyd->NodeDemand[i]; 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) + if (r <= 0) { - *hgrad = RB; - *hloss = RB * d; + *hgrad = CBIG; + *hloss = CBIG * d; } // Use linear head loss function for near zero demand @@ -568,7 +510,15 @@ void demandheadloss(double d, double dfull, double dp, double n, { *hgrad = n * dp * pow(r, n - 1.0) / dfull; *hloss = (*hgrad) * d / n; + + // Add on barrier function for any demand above full value + if (r >= 1.0) + { + *hgrad += CBIG; + *hloss += CBIG * (d - dfull); + } } + } diff --git a/src/hydraul.c b/src/hydraul.c index 0918a92..d3cd761 100755 --- a/src/hydraul.c +++ b/src/hydraul.c @@ -195,7 +195,7 @@ int runhyd(Project *pr, long *t) int iter; // Iteration count int errcode; // Error code double relerr; // Solution accuracy - + // Find new demands & control actions *t = time->Htime; demands(pr); @@ -203,7 +203,6 @@ int runhyd(Project *pr, long *t) // Solve network hydraulic equations errcode = hydsolve(pr,&iter,&relerr); - if (!errcode) { // Report new status & save results diff --git a/src/hydsolver.c b/src/hydsolver.c index 7811c06..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: 05/15/2019 + Last Updated: 07/15/2019 ****************************************************************************** */ @@ -50,6 +50,7 @@ static void newdemandflows(Project *, Hydbalance *, double *, double *); static void checkhydbalance(Project *, Hydbalance *); static int hasconverged(Project *, double *, Hydbalance *); +static int pdaconverged(Project *); static void reporthydbal(Project *, Hydbalance *); @@ -92,12 +93,17 @@ int hydsolve(Project *pr, int *iter, double *relerr) int valveChange; // Valve status change flag int statChange; // Non-valve status change flag Hydbalance hydbal; // Hydraulic balance errors + double fullDemand; // Full demand for a node (cfs) // Initialize status checking & relaxation factor nextcheck = hyd->CheckFreq; hyd->RelaxFactor = 1.0; + + // Initialize convergence criteria and PDA results hydbal.maxheaderror = 0.0; hydbal.maxflowchange = 0.0; + hyd->DeficientNodes = 0; + hyd->DemandReduction = 0.0; // Repeat iterations until convergence or trial limit is exceeded. // (ExtraIter used to increase trials in case of status cycling.) @@ -189,10 +195,12 @@ int hydsolve(Project *pr, int *iter, double *relerr) errcode = 110; } - // Replace junction demands with total outflow values + // Store actual junction outflow in NodeDemand & full demand in DemandFlow for (i = 1; i <= net->Njuncs; i++) { + fullDemand = hyd->NodeDemand[i]; hyd->NodeDemand[i] = hyd->DemandFlow[i] + hyd->EmitterFlow[i]; + hyd->DemandFlow[i] = fullDemand; } // Save convergence info @@ -200,7 +208,6 @@ int hydsolve(Project *pr, int *iter, double *relerr) hyd->MaxHeadError = hydbal.maxheaderror; hyd->MaxFlowChange = hydbal.maxflowchange; hyd->Iterations = *iter; - return errcode; } @@ -484,11 +491,12 @@ void newemitterflows(Project *pr, Hydbalance *hbal, double *qsum, if (net->Node[i].Ke == 0.0) continue; // Find emitter head loss and gradient - emitheadloss(pr, i, &hloss, &hgrad); + emitterheadloss(pr, i, &hloss, &hgrad); // Find emitter flow change dh = hyd->NodeHead[i] - net->Node[i].El; dq = (hloss - dh) / hgrad; + dq *= hyd->RelaxFactor; hyd->EmitterFlow[i] -= dq; // Update system flow summation @@ -523,32 +531,39 @@ void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum) 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; - + n, // exponent in head loss v. demand function + hloss, // current head loss through outflow junction (ft) + hgrad, // head loss gradient with respect to flow (ft/cfs) + dh; // new head loss through outflow junction (ft) + int i; + // Get demand function parameters if (hyd->DemandModel == DDA) return; - demandparams(pr, &dp, &n); + dp = MAX((hyd->Preq - hyd->Pmin), MINPDIFF); + n = 1.0 / hyd->Pexp; // Examine each junction - for (k = 1; k <= net->Njuncs; k++) + for (i = 1; i <= net->Njuncs; i++) { // Skip junctions with no positive demand - if (hyd->NodeDemand[k] <= 0.0) continue; - + if (hyd->NodeDemand[i] <= 0.0) continue; + // Find change in demand flow (see hydcoeffs.c) - dq = demandflowchange(pr, k, dp, n); - hyd->DemandFlow[k] -= dq; + demandheadloss(pr, i, dp, n, &hloss, &hgrad); + dh = hyd->NodeHead[i] - net->Node[i].El - hyd->Pmin; + dq = (hloss - dh) / hgrad; + dq *= hyd->RelaxFactor; + hyd->DemandFlow[i] -= dq; // Update system flow summation - *qsum += ABS(hyd->DemandFlow[k]); + *qsum += ABS(hyd->DemandFlow[i]); *dqsum += ABS(dq); // Update identity of element with max. flow change if (ABS(dq) > hbal->maxflowchange) { hbal->maxflowchange = ABS(dq); - hbal->maxflownode = k; + hbal->maxflownode = i; hbal->maxflowlink = -1; } } @@ -605,20 +620,72 @@ int hasconverged(Project *pr, double *relerr, Hydbalance *hbal) */ { Hydraul *hyd = &pr->hydraul; - + + // Check that total relative flow change is small enough if (*relerr > hyd->Hacc) return 0; + + // Find largest head loss error and absolute flow change checkhydbalance(pr, hbal); if (pr->report.Statflag == FULL) { reporthydbal(pr, hbal); } + + // Check that head loss error and flow change criteria are met if (hyd->HeadErrorLimit > 0.0 && hbal->maxheaderror > hyd->HeadErrorLimit) return 0; if (hyd->FlowChangeLimit > 0.0 && hbal->maxflowchange > hyd->FlowChangeLimit) return 0; + + // Check for pressure driven analysis convergence + if (hyd->DemandModel == PDA) return pdaconverged(pr); return 1; } +int pdaconverged(Project *pr) +/* +**-------------------------------------------------------------- +** Input: none +** Output: returns 1 if PDA converged, 0 if not +** Purpose: checks if pressure driven analysis has converged +** and updates total demand deficit +**-------------------------------------------------------------- +*/ +{ + Hydraul *hyd = &pr->hydraul; + + const double TOL = 0.001; + int i, converged = 1; + double totalDemand = 0.0, totalReduction = 0.0; + + hyd->DeficientNodes = 0; + hyd->DemandReduction = 0.0; + + // Add up number of junctions with demand deficits + for (i = 1; i <= pr->network.Njuncs; i++) + { + // Skip nodes whose required demand is non-positive + if (hyd->NodeDemand[i] <= 0.0) continue; + + // Check for negative demand flow or positive demand flow at negative pressure + if (hyd->DemandFlow[i] < -TOL) converged = 0; + if (hyd->DemandFlow[i] > TOL && + hyd->NodeHead[i] - pr->network.Node[i].El - hyd->Pmin < -TOL) + converged = 0; + + // Accumulate total required demand and demand deficit + if (hyd->DemandFlow[i] + 0.0001 < hyd->NodeDemand[i]) + { + hyd->DeficientNodes++; + totalDemand += hyd->NodeDemand[i]; + totalReduction += hyd->NodeDemand[i] - hyd->DemandFlow[i]; + } + } + if (totalDemand > 0.0) + hyd->DemandReduction = totalReduction / totalDemand * 100.0; + return converged; +} + void reporthydbal(Project *pr, Hydbalance *hbal) /* diff --git a/src/input1.c b/src/input1.c index fd972e1..10a6f95 100644 --- a/src/input1.c +++ b/src/input1.c @@ -7,7 +7,7 @@ Description: retrieves network data from an EPANET input file Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE -Last Updated: 05/15/2019 +Last Updated: 07/08/2019 ****************************************************************************** */ @@ -109,7 +109,7 @@ void setdefaults(Project *pr) 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->Preq = MINPDIFF; // Required demand pressure (ft) hyd->Pexp = 0.5; // Pressure function exponent hyd->MaxIter = MAXITER; // Default max. hydraulic trials hyd->ExtraIter = -1; // Stop if network unbalanced @@ -267,7 +267,7 @@ void adjustdata(Project *pr) // Revise pressure units depending on flow units if (parser->Unitsflag != SI) parser->Pressflag = PSI; else if (parser->Pressflag == PSI) parser->Pressflag = METERS; - + // Store value of viscosity & diffusivity ucf = 1.0; if (parser->Unitsflag == SI) ucf = SQR(MperFT); @@ -322,9 +322,9 @@ void adjustdata(Project *pr) if (tank->Kb == MISSING) tank->Kb = qual->Kbulk; } - // Use default pattern if none assigned to a demand - parser->DefPat = findpattern(net, parser->DefPatID); - if (parser->DefPat > 0) for (i = 1; i <= net->Nnodes; i++) + // Use default pattern if none assigned to a demand + parser->DefPat = findpattern(net, parser->DefPatID); + if (parser->DefPat > 0) for (i = 1; i <= net->Nnodes; i++) { node = &net->Node[i]; for (demand = node->D; demand != NULL; demand = demand->next) @@ -559,7 +559,8 @@ void convertunits(Project *pr) demand->Base /= pr->Ucf[DEMAND]; } } - + + // Convert PDA pressure limits hyd->Pmin /= pr->Ucf[PRESSURE]; hyd->Preq /= pr->Ucf[PRESSURE]; diff --git a/src/input3.c b/src/input3.c index 0770a83..940e572 100644 --- a/src/input3.c +++ b/src/input3.c @@ -1924,12 +1924,20 @@ int optionvalue(Project *pr, int n) else if (match(tok0, w_MINIMUM)) { if (y < 0.0) return setError(parser, nvalue, 213); + // Required pressure still at default value + if (hyd->Preq == MINPDIFF) + hyd->Preq = y + MINPDIFF; + // Required pressure already entered + else if (hyd->Preq - y < MINPDIFF) + return setError(parser, nvalue, 208); hyd->Pmin = y; return 0; } else if (match(tok0, w_REQUIRED)) { if (y < 0.0) return setError(parser, nvalue, 213); + if (y - hyd->Pmin < MINPDIFF) + return setError(parser, nvalue, 208); hyd->Preq = y; return 0; } diff --git a/src/report.c b/src/report.c index d9d67bb..10a9774 100644 --- a/src/report.c +++ b/src/report.c @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 06/20/2019 + Last Updated: 07/15/2019 ****************************************************************************** */ @@ -358,7 +358,15 @@ void writehydstat(Project *pr, int iter, double relerr) { if (relerr <= hyd->Hacc) sprintf(s1, FMT58, atime, iter); else sprintf(s1, FMT59, atime, iter, relerr); - writeline(pr, s1); + writeline(pr, s1); + if (hyd->DemandModel == PDA && hyd->DeficientNodes > 0) + { + if (hyd->DeficientNodes == 1) + sprintf(s1, FMT69a, hyd->DemandReduction); + else + sprintf(s1, FMT69b, hyd->DeficientNodes, hyd->DemandReduction); + writeline(pr, s1); + } } // Display status changes for tanks: @@ -1064,14 +1072,9 @@ int writehydwarn(Project *pr, int iter, double relerr) int i, j; char flag = 0; int s; - Snode *Node = net->Node; - Slink *Link = net->Link; - Spump *Pump = net->Pump; - Svalve *Valve = net->Valve; - const int Njuncs = net->Njuncs; - double *NodeDemand = hyd->NodeDemand; - double *LinkFlow = hyd->LinkFlow; - double *LinkSetting = hyd->LinkSetting; + Snode *node; + Slink *link; + Spump *pump; // Check if system unstable if (iter > hyd->MaxIter && relerr <= hyd->Hacc) @@ -1081,29 +1084,41 @@ int writehydwarn(Project *pr, int iter, double relerr) flag = 2; } - // Check for negative pressures - for (i = 1; i <= Njuncs; i++) + // Check for pressure deficient nodes + if (hyd->DemandModel == DDA) { - Snode *node = &Node[i]; - if (hyd->NodeHead[i] < node->El && NodeDemand[i] > 0.0) + hyd->DeficientNodes = 0; + for (i = 1; i <= net->Njuncs; i++) { - sprintf(pr->Msg, WARN06, clocktime(rpt->Atime, time->Htime)); - if (rpt->Messageflag) writeline(pr, pr->Msg); + node = &net->Node[i]; + if (hyd->NodeHead[i] < node->El && hyd->NodeDemand[i] > 0.0) + hyd->DeficientNodes++; + } + if (hyd->DeficientNodes > 0) + { + if (rpt->Messageflag) + { + sprintf(pr->Msg, WARN06, clocktime(rpt->Atime, time->Htime)); + writeline(pr, pr->Msg); + } flag = 6; - break; } } // Check for abnormal valve condition for (i = 1; i <= net->Nvalves; i++) { - j = Valve[i].Link; + j = net->Valve[i].Link; + link = &net->Link[j]; if (hyd->LinkStatus[j] >= XFCV) { - sprintf(pr->Msg, WARN05, LinkTxt[Link[j].Type], Link[j].ID, - StatTxt[hyd->LinkStatus[j]], - clocktime(rpt->Atime, time->Htime)); - if (rpt->Messageflag) writeline(pr, pr->Msg); + if (rpt->Messageflag) + { + sprintf(pr->Msg, WARN05, LinkTxt[link->Type], link->ID, + StatTxt[hyd->LinkStatus[j]], + clocktime(rpt->Atime, time->Htime)); + writeline(pr, pr->Msg); + } flag = 5; } } @@ -1111,18 +1126,22 @@ int writehydwarn(Project *pr, int iter, double relerr) // Check for abnormal pump condition for (i = 1; i <= net->Npumps; i++) { - j = Pump[i].Link; + pump = &net->Pump[i]; + j = pump->Link; s = hyd->LinkStatus[j]; if (hyd->LinkStatus[j] >= OPEN) { - if (LinkFlow[j] > LinkSetting[j] * Pump[i].Qmax) s = XFLOW; - if (LinkFlow[j] < 0.0) s = XHEAD; + if (hyd->LinkFlow[j] > hyd->LinkSetting[j] * pump->Qmax) s = XFLOW; + if (hyd->LinkFlow[j] < 0.0) s = XHEAD; } if (s == XHEAD || s == XFLOW) { - sprintf(pr->Msg, WARN04, Link[j].ID, StatTxt[s], - clocktime(rpt->Atime, time->Htime)); - if (rpt->Messageflag) writeline(pr, pr->Msg); + if (rpt->Messageflag) + { + sprintf(pr->Msg, WARN04, net->Link[j].ID, StatTxt[s], + clocktime(rpt->Atime, time->Htime)); + writeline(pr, pr->Msg); + } flag = 4; } } @@ -1130,9 +1149,12 @@ int writehydwarn(Project *pr, int iter, double relerr) // Check if system is unbalanced if (iter > hyd->MaxIter && relerr > hyd->Hacc) { - sprintf(pr->Msg, WARN01, clocktime(rpt->Atime, time->Htime)); - if (hyd->ExtraIter == -1) strcat(pr->Msg, t_HALTED); - if (rpt->Messageflag) writeline(pr, pr->Msg); + if (rpt->Messageflag) + { + sprintf(pr->Msg, WARN01, clocktime(rpt->Atime, time->Htime)); + if (hyd->ExtraIter == -1) strcat(pr->Msg, t_HALTED); + writeline(pr, pr->Msg); + } flag = 1; } @@ -1162,9 +1184,12 @@ void writehyderr(Project *pr, int errnode) Snode *Node = net->Node; - sprintf(pr->Msg, FMT62, clocktime(rpt->Atime, time->Htime), - Node[errnode].ID); - if (rpt->Messageflag) writeline(pr, pr->Msg); + if (rpt->Messageflag) + { + sprintf(pr->Msg, FMT62, clocktime(rpt->Atime, time->Htime), + Node[errnode].ID); + writeline(pr, pr->Msg); + } writehydstat(pr, 0, 0); disconnected(pr); } diff --git a/src/text.h b/src/text.h index 2ce73da..a4ac84b 100755 --- a/src/text.h +++ b/src/text.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 06/20/2019 + Last Updated: 07/15/2019 ****************************************************************************** */ @@ -428,6 +428,8 @@ #define FMT66 " maximum flow change = %.4f for Link %s" #define FMT67 " maximum flow change = %.4f for Node %s" #define FMT68 " maximum head error = %.4f for Link %s\n" +#define FMT69a " 1 node had its demand reduced by a total of %.2f%%" +#define FMT69b " %-d nodes had demands reduced by a total of %.2f%%" //----- Energy Report Table ------------------------------- diff --git a/src/types.h b/src/types.h index c07d4aa..7e89a3f 100755 --- a/src/types.h +++ b/src/types.h @@ -7,7 +7,7 @@ Authors: see AUTHORS Copyright: see AUTHORS License: see LICENSE - Last Updated: 06/20/2019 + Last Updated: 07/08/2019 ****************************************************************************** */ @@ -42,7 +42,6 @@ typedef int INT4; #define MAXLINE 1024 // Max. # characters read from input line #define MAXFNAME 259 // Max. # characters in file name #define MAXTOKS 40 // Max. items per line of input -#define TZERO 1.E-4 // Zero time tolerance #define TRUE 1 #define FALSE 0 #define FULL 2 @@ -53,6 +52,7 @@ typedef int INT4; // @ 20 deg C (sq ft/sec) #define VISCOS 1.1E-5 // Kinematic viscosity of water // @ 20 deg C (sq ft/sec) +#define MINPDIFF 0.1 // PDA min. pressure difference (psi or m) #define SEPSTR " \t\n\r" // Token separator characters #ifdef M_PI #define PI M_PI @@ -689,7 +689,7 @@ typedef struct { double *NodeHead, // Node hydraulic heads *NodeDemand, // Node demand + emitter flows - *DemandFlow, // Demand outflows + *DemandFlow, // Work array of demand flows *EmitterFlow, // Emitter outflows *LinkFlow, // Link flows *LinkSetting, // Link settings @@ -716,6 +716,7 @@ typedef struct { RelativeError, // Total flow change / total flow MaxHeadError, // Max. error for link head loss MaxFlowChange, // Max. change in link flow + DemandReduction, // % demand reduction at pressure deficient nodes RelaxFactor, // Relaxation factor for flow updating *P, // Inverse of head loss derivatives *Y, // Flow correction factors @@ -731,7 +732,8 @@ typedef struct { CheckFreq, // Hydraulic trials between status checks MaxCheck, // Hydraulic trials limit on status checks OpenHflag, // Hydraulic system opened flag - Haltflag; // Flag to halt simulation + Haltflag, // Flag to halt simulation + DeficientNodes; // Number of pressure deficient nodes StatusType *LinkStatus, // Link status diff --git a/tests/test_pda.cpp b/tests/test_pda.cpp new file mode 100644 index 0000000..7811daf --- /dev/null +++ b/tests/test_pda.cpp @@ -0,0 +1,90 @@ +/* + ****************************************************************************** + Project: OWA EPANET + Version: 2.2 + Module: test_pda.cpp + Description: Tests EPANET toolkit api functions + Authors: see AUTHORS + Copyright: see AUTHORS + License: see LICENSE + Last Updated: 07/20/2019 + ****************************************************************************** +*/ + +/* + Tests the Pressure Driven Analysis option +*/ + +#include + +#include "test_toolkit.hpp" + +BOOST_AUTO_TEST_SUITE (test_pda) + +BOOST_AUTO_TEST_CASE(test_pda) + +{ + int error = 0; + int index; + double count, reduction; + + EN_Project ph = NULL; + error = EN_createproject(&ph); + error = EN_open(ph, DATA_PATH_NET1, DATA_PATH_RPT, ""); + + // Set Demand Multiplier to 10 to cause negative pressures + error = EN_setoption(ph, EN_DEMANDMULT, 10); + BOOST_REQUIRE(error == 0); + + // Run single period analysis + error = EN_settimeparam(ph, EN_DURATION, 0); + BOOST_REQUIRE(error == 0); + + // Solve hydraulics with default DDA option + // which will return with neg. pressure warning code + error = EN_solveH(ph); + BOOST_REQUIRE(error == 6); + + // Check that 4 demand nodes have negative pressures + error = EN_getstatistic(ph, EN_DEFICIENTNODES, &count); + BOOST_REQUIRE(error == 0); + BOOST_REQUIRE(count == 4); + + // Switch to PDA with pressure limits of 20 - 100 psi + error = EN_setdemandmodel(ph, EN_PDA, 20, 100, 0.5); + BOOST_REQUIRE(error == 0); + + // Solve hydraulics again + error = EN_solveH(ph); + BOOST_REQUIRE(error == 0); + + // Check that 6 nodes had demand reductions totaling 32.66% + error = EN_getstatistic(ph, EN_DEFICIENTNODES, &count); + BOOST_REQUIRE(error == 0); + BOOST_REQUIRE(count == 6); + error = EN_getstatistic(ph, EN_DEMANDREDUCTION, &reduction); + BOOST_REQUIRE(error == 0); + BOOST_REQUIRE(abs(reduction - 32.66) < 0.01); + + // Check that Junction 12 had full demand + error = EN_getnodeindex(ph, (char *)"12", &index); + BOOST_REQUIRE(error == 0); + error = EN_getnodevalue(ph, index, EN_DEMANDDEFICIT, &reduction); + BOOST_REQUIRE(error == 0); + BOOST_REQUIRE(abs(reduction) < 0.01); + + // Check that Junction 21 had a demand deficit of 413.67 + error = EN_getnodeindex(ph, (char *)"21", &index); + BOOST_REQUIRE(error == 0); + error = EN_getnodevalue(ph, index, EN_DEMANDDEFICIT, &reduction); + BOOST_REQUIRE(error == 0); + BOOST_REQUIRE(abs(reduction - 413.67) < 0.01); + + // Clean up + error = EN_close(ph); + BOOST_REQUIRE(error == 0); + error = EN_deleteproject(ph); + BOOST_REQUIRE(error == 0); +} + +BOOST_AUTO_TEST_SUITE_END()