Merge pull request #505 from OpenWaterAnalytics/feature-tankoverflow

Adds tank overflow feature

ReleaseNotes2_2.md

@@ -131,6 +131,11 @@ for the thread-safe API. Some additional points regarding the new **PDA** option
  - 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.
 
+## 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.
+
+For the input file, a new field has been appended to the data supplied for each tank in the `[TANKS]` section of the file. A value of **YES** indicates that the tank can over flow while **NO** (the default) indicates that it cannot. For the volume curve field that precedes it, an asterisk (*) can be used if the tank does not utilize a volume curve.
+
 ## Improved Water Quality Mass Balance
 
 As described by [Davis et al. (2018)](https://www.drink-water-eng-sci.net/11/25/2018/dwes-11-25-2018.pdf) EPANET's water quality simulations can result in some significant mass balance errors when modeling short term mass injections (errors are much smaller for continuous source flows). The entire water quality engine has been re-written to eliminate these errors. It still uses the Lagrangian Time Driven transport method but now analyzes each network node in topologically sorted order rather than in arbitrary order.
@@ -208,6 +213,8 @@ Access to the following global energy options have been added to  `EN_getoption`
  - `EN_DEMANDCHARGE` (price per maximum kW of energy consumption)
 
 ## New API Constants
+### Node value types:
+- `EN_CANOVERFLOW`
 
 ### Link value types:
 - `EN_PUMP_STATE`

doc/input-file.dox

@@ -912,21 +912,24 @@ One line for each junction containing:
 - Nominal diameter, ft (m)
 - Minimum volume, cubic ft (cubic meters)
 - Volume curve ID (optional)
+- Overflow indicator (<b>YES / NO</b>) (optional)
 
 __Remarks:__
 1. Water surface elevation equals bottom elevation plus water level.
 2. Non-cylindrical tanks can be modeled by specifying a curve of volume versus water depth in the @ref CurvesPage section.
 3. If a volume curve is supplied the diameter value can be any non-zero number
 4. Minimum volume (tank volume at minimum water level) can be zero for a cylindrical tank or if a volume curve is supplied.
-5. A network must contain at least one tank or reservoir.
+5. If the overflow indicator is \b YES then the tank is allowed to overflow once it reaches it maximum water level. The default is no overflow.
+6. If the tank does not use a volume curve then an asterisk (*) can be used as a placeholder for it if an overflow indicator is specified.
+7. A network must contain at least one tank or reservoir.
 
 __Example:__
 @code
 [TANKS] 
-;ID   Elev.  InitLvl  MinLvl  MaxLvl  Diam  MinVol  VolCurve 
+;ID   Elev.  InitLvl  MinLvl  MaxLvl  Diam  MinVol  VolCurve  Overflow 
-;----------------------------------------------------------- 
+;--------------------------------------------------------------------- 
-;Cylindrical tank 
+;Cylindrical tank that can overflow
-T1    100     15       5       25     120   0 
+T1    100     15       5       25     120   0       *          YES 
 
 ;Non-cylindrical tank with arbitrary diameter
 T2   100     15       5       25     1     0        VC1    

include/epanet2.bas

@@ -35,6 +35,7 @@ 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_DIAMETER = 0      ' Link parameters
 Public Const EN_LENGTH = 1

include/epanet2.vb

@@ -38,8 +38,9 @@ Public Const EN_MAXLEVEL = 21
 Public Const EN_MIXFRACTION = 22
 Public Const EN_TANK_KBULK = 23
 
-Public Const EN_TANKVOLUME = 24 'ES
+Public Const EN_TANKVOLUME = 24
 Public Const EN_MAXVOLUME = 25
+Public Const EN_CANOVERFLOW  = 26
 
 Public Const EN_DIAMETER = 0      ' Link parameters
 Public Const EN_LENGTH = 1

include/epanet2_enums.h

@@ -61,7 +61,8 @@ typedef enum {
   EN_MIXFRACTION  = 22, //!< Tank mixing fraction
   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_MAXVOLUME    = 25, //!< Tank maximum volume (read only)
+  EN_CANOVERFLOW  = 26  //!< Tank can overflow (= 1) or not (= 0)
 } EN_NodeProperty;
 
 /// Link properties

src/enumstxt.h

@@ -7,7 +7,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 11/27/2018
+ Last Updated: 06/20/2019
  ******************************************************************************
 */
 
@@ -38,7 +38,8 @@ char *StatTxt[]         = {t_XHEAD,
                            t_XFCV,
                            t_XPRESSURE,
                            t_FILLING,
-                           t_EMPTYING};
+                           t_EMPTYING,
+                           t_OVERFLOWING};
 
 char *FormTxt[]         = {w_HW,
                            w_DW,

src/epanet.c

@@ -1802,6 +1802,7 @@ int DLLEXPORT EN_addnode(EN_Project p, char *id, int nodeType, int *index)
         tank->Vcurve = 0;
         tank->MixModel = 0;
         tank->V1max = 10000;
+        tank->CanOverflow = FALSE;
     }
     net->Nnodes++;
     p->parser.MaxNodes = net->Nnodes;
@@ -2197,6 +2198,11 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val
         v = tankvolume(p, index - nJuncs, NodeHead[index]) * Ucf[VOLUME];
         break;
 
+    case EN_CANOVERFLOW:
+        if (Node[index].Type != TANK) return 0;
+        v = Tank[index - nJuncs].CanOverflow;
+        break;
+
     default:
         return 251;
     }
@@ -2495,6 +2501,11 @@ int DLLEXPORT EN_setnodevalue(EN_Project p, int index, int property, double valu
         }
         break;
 
+    case EN_CANOVERFLOW:
+        if (Node[index].Type != TANK) return 0;
+        Tank[index - nJuncs].CanOverflow = (value != 0.0);
+        break;
+
     default:
         return 251;
     }

src/hydstatus.c

@@ -441,7 +441,7 @@ void  tankstatus(Project *pr, int k, int n1, int n2)
     h = hyd->NodeHead[n1] - hyd->NodeHead[n2];
 
     // If tank is full, then prevent flow into it
-    if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol)
+    if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol && !tank->CanOverflow)
     {
         // Case 1: Link is a pump discharging into tank
         if (link->Type == PUMP)

src/inpfile.c

@@ -195,8 +195,10 @@ int saveinpfile(Project *pr, const char *fname)
                     sqrt(4.0 * tank->A / PI) * pr->Ucf[ELEV],
                     tank->Vmin * SQR(pr->Ucf[ELEV]) * pr->Ucf[ELEV]);
             if ((j = tank->Vcurve) > 0) sprintf(s1, "%s", net->Curve[j].ID);
+            else if (tank->CanOverflow) strcpy(s1, "*");
             else strcpy(s1, " ");
             fprintf(f, "\n%s %-31s", s, s1);
+            if (tank->CanOverflow) fprintf(f, "  YES  ");
             if (node->Comment) fprintf(f, " ;%s", node->Comment);
         }
     }

src/input3.c

@@ -7,7 +7,7 @@ Description:  parses network data from a line of an EPANET input file
 Authors:      see AUTHORS
 Copyright:    see AUTHORS
 License:      see LICENSE
-Last Updated: 05/24/2019
+Last Updated: 06/19/2019
 ******************************************************************************
 */
 
@@ -138,7 +138,8 @@ int tankdata(Project *pr)
     int    i,               // Node index
            n,               // # data items
            pattern = 0,     // Time pattern index
-           curve = 0;       // Curve index
+           curve = 0,       // Curve index
+           overflow = FALSE;// Overflow indicator    
     double el = 0.0,        // Elevation
            initlevel = 0.0, // Initial level
            minlevel = 0.0,  // Minimum level
@@ -185,12 +186,24 @@ int tankdata(Project *pr)
         if (n >= 7 && !getfloat(parser->Tok[6], &minvol)) return setError(parser, 6, 202);
 
         // If volume curve supplied check it exists
-        if (n == 8)
+        if (n >= 8)
+        {
+            if (strlen(parser->Tok[7]) > 0 && *(parser->Tok[7]) != '*')
             {
                 curve = findcurve(net, parser->Tok[7]);
                 if (curve == 0) return setError(parser, 7, 206);
                 net->Curve[curve].Type = VOLUME_CURVE;
             }
+        }
+
+        // Parse overflow indicator if present
+        if (n >= 9)
+        {
+            if (match(parser->Tok[8], w_YES)) overflow = TRUE;
+            else if (match(parser->Tok[8], w_NO)) overflow = FALSE;
+            else return setError(parser, 8, 213);
+        }
+
         if (initlevel < 0.0) return setError(parser, 2, 209);
         if (minlevel  < 0.0) return setError(parser, 3, 209);
         if (maxlevel  < 0.0) return setError(parser, 4, 209);
@@ -216,6 +229,7 @@ int tankdata(Project *pr)
     tank->A = diam;
     tank->Pat = pattern;
     tank->Kb = MISSING;
+    tank->CanOverflow = overflow;
 
     //*******************************************************************
     // NOTE: The min, max, & initial volumes set here are based on a

src/report.c

@@ -7,7 +7,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 05/15/2019
+ Last Updated: 06/20/2019
  ******************************************************************************
 */
 
@@ -370,8 +370,13 @@ void writehydstat(Project *pr, int iter, double relerr)
     n = net->Tank[i].Node;
     NodeDemand = hyd->NodeDemand;
     if (ABS(NodeDemand[n]) < 0.001) newstat = CLOSED;
-    else if (NodeDemand[n] > 0.0)   newstat = FILLING;
     else if (NodeDemand[n] < 0.0)   newstat = EMPTYING;
+    else if (NodeDemand[n] > 0.0)
+    {
+        if (Tank[i].A > 0.0 && ABS(hyd->NodeHead[n] - Tank[i].Hmax) < 0.001)
+            newstat = OVERFLOWING;
+        else newstat = FILLING;
+    }
     else newstat = hyd->OldStatus[net->Nlinks + i];
     if (newstat != hyd->OldStatus[net->Nlinks + i])
     {

src/text.h

@@ -7,7 +7,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 11/27/2018
+ Last Updated: 06/20/2019
  ******************************************************************************
 */
 
@@ -282,6 +282,7 @@
 #define   t_XPRESSURE   "open but cannot deliver pressure"
 #define   t_FILLING     "filling"
 #define   t_EMPTYING    "emptying"
+#define   t_OVERFLOWING "overflowing"
 
 #define   t_ELEV        "Elevation"
 #define   t_DEMAND      "Demand"

src/types.h

@@ -7,7 +7,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 05/15/2019
+ Last Updated: 06/20/2019
  ******************************************************************************
 */
 
@@ -200,7 +200,8 @@ typedef enum {
   XFCV,          // FCV cannot supply flow
   XPRESSURE,     // valve cannot supply pressure
   FILLING,       // tank filling
-  EMPTYING       // tank emptying
+  EMPTYING,      // tank emptying
+  OVERFLOWING    // tank overflowing
 } StatusType;
 
 typedef enum {
@@ -417,6 +418,7 @@ typedef struct             // Tank Object
   int     Vcurve;          // volume v. elev. curve index
   MixType MixModel;        // type of mixing model
   double  V1max;           // mixing compartment size
+  int     CanOverflow;     // tank can overflow or not
 } Stank;
 
 typedef struct             // Pump Object

tests/test_overflow.cpp

@@ -0,0 +1,119 @@
+/*
+ ******************************************************************************
+ Project:      OWA EPANET
+ Version:      2.2
+ Module:       test_overflow.cpp
+ Description:  Tests EPANET toolkit api functions
+ Authors:      see AUTHORS
+ Copyright:    see AUTHORS
+ License:      see LICENSE
+ Last Updated: 06/16/2019
+ ******************************************************************************
+*/
+
+/*
+   Tests the EN_CANOVERFLOW option for Tank nodes
+*/
+
+#include <boost/test/unit_test.hpp>
+
+#include "test_toolkit.hpp"
+
+BOOST_AUTO_TEST_SUITE (test_overflow)
+
+BOOST_AUTO_TEST_CASE(test_tank_overflow)
+
+{
+    int error = 0;
+    int Nindex, Lindex;
+    double level, spillage, spillage2, inflow;
+    char testFile[] = "test_overflow.inp";
+
+    EN_Project ph = NULL;
+
+    error = EN_createproject(&ph);
+    error = EN_open(ph, DATA_PATH_NET1, DATA_PATH_RPT, "");
+
+    // Get index of the tank and its inlet/outlet pipe
+    error = EN_getnodeindex(ph, (char *)"2", &Nindex);
+    BOOST_REQUIRE(error == 0);
+    error = EN_getlinkindex(ph, (char *)"110", &Lindex);
+    BOOST_REQUIRE(error == 0);
+
+    // Set initial & maximum level to 130
+    error = EN_setnodevalue(ph, Nindex, EN_TANKLEVEL, 130);
+    BOOST_REQUIRE(error == 0);
+    error = EN_setnodevalue(ph, Nindex, EN_MAXLEVEL, 130);
+    BOOST_REQUIRE(error == 0);
+
+    // Set duration to 1 hr
+    error = EN_settimeparam(ph, EN_DURATION, 3600);
+    BOOST_REQUIRE(error == 0);
+
+    // Solve hydraulics with default of no tank spillage allowed
+    error = EN_solveH(ph);
+    BOOST_REQUIRE(error == 0);
+
+    // Check that tank remains full
+    error = EN_getnodevalue(ph, Nindex, EN_TANKLEVEL, &level);
+    BOOST_REQUIRE(error == 0);
+    BOOST_REQUIRE(abs(level - 130.0) < 0.0001);
+
+    // Check that there is no spillage
+    error = EN_getnodevalue(ph, Nindex, EN_DEMAND, &spillage);
+    BOOST_REQUIRE(error == 0);
+    BOOST_REQUIRE(abs(spillage) < 0.0001);
+
+    // Check that inflow link is closed
+    error = EN_getlinkvalue(ph, Lindex, EN_FLOW, &inflow);
+    BOOST_REQUIRE(error == 0);
+    BOOST_REQUIRE(abs(inflow) < 0.0001);
+
+    // Turn tank overflow option on
+    error = EN_setnodevalue(ph, Nindex, EN_CANOVERFLOW, 1);
+    BOOST_REQUIRE(error == 0);
+
+    // Solve hydraulics again
+    error = EN_solveH(ph);
+    BOOST_REQUIRE(error == 0);
+
+    // Check that tank remains full
+    error = EN_getnodevalue(ph, Nindex, EN_TANKLEVEL, &level);
+    BOOST_REQUIRE(error == 0);
+    BOOST_REQUIRE(abs(level - 130.0) < 0.0001);
+
+    // Check that there is spillage equal to tank inflow
+    // (inflow has neg. sign since tank is start node of inflow pipe)
+    error = EN_getnodevalue(ph, Nindex, EN_DEMAND, &spillage);
+    BOOST_REQUIRE(error == 0);
+    BOOST_REQUIRE(spillage > 0.0001);
+    error = EN_getlinkvalue(ph, Lindex, EN_FLOW, &inflow);
+    BOOST_REQUIRE(error == 0);
+    BOOST_REQUIRE(abs(-inflow - spillage) < 0.0001);    
+    
+    // Save project to file and then close it
+    error = EN_saveinpfile(ph, testFile);
+    BOOST_REQUIRE(error == 0);
+    error = EN_close(ph);
+    BOOST_REQUIRE(error == 0);
+    
+    // Re-open saved file & run it
+    error = EN_open(ph, testFile, DATA_PATH_RPT, "");
+    BOOST_REQUIRE(error == 0);
+    error = EN_solveH(ph);
+    BOOST_REQUIRE(error == 0);
+
+    // Check that tank spillage has same value as before
+    error = EN_getnodevalue(ph, Nindex, EN_DEMAND, &spillage2);
+    BOOST_REQUIRE(error == 0);
+    BOOST_REQUIRE(abs(spillage - spillage2) < 0.0001);
+    
+    // Clean up
+    error = EN_close(ph);
+    BOOST_REQUIRE(error == 0);
+    error = EN_deleteproject(&ph);
+    BOOST_REQUIRE(error == 0);
+    
+}
+
+BOOST_AUTO_TEST_SUITE_END()