Update epanet2_enums.h

include/epanet2_enums.h

@@ -9,7 +9,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 07/17/2023
+ Last Updated: 06/25/2024
  ******************************************************************************
 */
 
@@ -65,7 +65,10 @@ typedef enum {
   EN_CANOVERFLOW  = 26, //!< Tank can overflow (= 1) or not (= 0)
   EN_DEMANDDEFICIT = 27,//!< Amount that full demand is reduced under PDA (read only)
   EN_NODE_INCONTROL = 28, //!< Is present in any simple or rule-based control (= 1) or not (= 0)
-  EN_EMITTERFLOW    = 29  //!< Current emitter flow (read only)
+  EN_EMITTERFLOW    = 29, //!< Current emitter flow (read only)
+  EN_LEAKAGEFLOW    = 30, //!< Current leakage flow (read only)
+  EN_DEMANDFLOW     = 31, //!< Current consumer demand delivered (read only)
+  EN_FULLDEMAND     = 32  //!< Current consumer demand requested (read only)
 } EN_NodeProperty;
 
 /// Link properties
@@ -99,7 +102,10 @@ typedef enum {
   EN_PUMP_EPAT    = 22, //!< Pump energy price time pattern index
   EN_LINK_INCONTROL = 23,  //!< Is present in any simple or rule-based control (= 1) or not (= 0)
   EN_GPV_CURVE    = 24, //!< GPV head loss v. flow curve index
-  EN_PCV_CURVE    = 25  //!< PCV loss coeff. curve index
+  EN_PCV_CURVE    = 25, //!< PCV loss coeff. curve index
+  EN_LEAK_AREA    = 26, //!< Pipe leak area (sq mm per 100 length units)
+  EN_LEAK_EXPAN   = 27, //!< Leak expansion rate (sq mm per unit of pressure head)
+  EN_LINK_LEAKAGE = 28  //!< Current leakage rate (read only)
 } EN_LinkProperty;
 
 /// Time parameters
@@ -152,7 +158,8 @@ typedef enum {
   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_DEMANDREDUCTION = 6, //!< % demand reduction at pressure deficient nodes
+  EN_LEAKAGELOSS     = 7  //!< % flow lost to system leakage
 } EN_AnalysisStatistic;
 
 /// Types of network objects

src/enumstxt.h

@@ -7,7 +7,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 02/05/2023
+ Last Updated: 06/24/2024
  ******************************************************************************
 */
 
@@ -127,7 +127,7 @@ char *SectTxt[]         = {s_TITLE,     s_JUNCTIONS, s_RESERVOIRS,
                            s_REACTIONS, s_MIXING,    s_REPORT,
                            s_TIMES,     s_OPTIONS,   s_COORDS,
                            s_VERTICES,  s_LABELS,    s_BACKDROP,
-                           s_TAGS,      s_END,
+                           s_TAGS,      s_LEAKAGE,   s_END,
                            NULL};
 
 char *Fldname[]         = {t_ELEV,      t_DEMAND,    t_HEAD,

src/epanet.c

@@ -1044,6 +1044,9 @@ int DLLEXPORT EN_getstatistic(EN_Project p, int type, double *value)
     case EN_DEMANDREDUCTION:
         *value = p->hydraul.DemandReduction;
         break;
+    case EN_LEAKAGELOSS:    
+        *value = p->hydraul.LeakageLoss;
+        break;
     case EN_MASSBALANCE:
         *value = p->quality.MassBalance.ratio;
         break;
@@ -1864,8 +1867,10 @@ int DLLEXPORT EN_addnode(EN_Project p, const char *id, int nodeType, int *index)
     hyd->NodeDemand = (double *)realloc(hyd->NodeDemand, size);
     qual->NodeQual = (double *)realloc(qual->NodeQual, size);
     hyd->NodeHead = (double *)realloc(hyd->NodeHead, size);
-    hyd->DemandFlow = (double *)realloc(hyd->DemandFlow, size);
+    hyd->FullDemand = (double *)realloc(hyd->FullDemand, size);
     hyd->EmitterFlow = (double *)realloc(hyd->EmitterFlow, size);
+    hyd->LeakageFlow = (double *)realloc(hyd->LeakageFlow, size);
+    hyd->DemandFlow = (double *)realloc(hyd->DemandFlow, size);
 
     // Actions taken when a new Junction is added
     if (nodeType == EN_JUNCTION)
@@ -2256,7 +2261,7 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val
                 Ucf[VOLUME];
         break;
 
-    case EN_DEMAND:
+    case EN_DEMAND: // Consumer Demand + Emitter Flow + Leakage Flow
         v = hyd->NodeDemand[index] * Ucf[FLOW];
         break;
 
@@ -2336,11 +2341,13 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val
 
     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];
+        // FullDemand contains node's required consumer demand
+        // while DemandFlow contains delivered consumer demand
+        if (hyd->FullDemand[index] <= 0.0) return 0;
+        v = (hyd->FullDemand[index] - hyd->DemandFlow[index]) /
+             hyd->FullDemand[index];
+        if (v < TINY) v = 0.0;
+        v *= 100.0;
         break;
 
     case EN_NODE_INCONTROL:
@@ -2350,6 +2357,18 @@ int DLLEXPORT EN_getnodevalue(EN_Project p, int index, int property, double *val
     case EN_EMITTERFLOW:
         v = hyd->EmitterFlow[index] * Ucf[FLOW];
         break;
+
+    case EN_LEAKAGEFLOW:
+        v =  hyd->LeakageFlow[index] * Ucf[FLOW];
+        break;
+
+    case EN_DEMANDFLOW:  // Consumer demand delivered
+        v = hyd->DemandFlow[index] * Ucf[FLOW];
+        break;
+
+    case EN_FULLDEMAND:  // Consumer demand requested
+        v = hyd->FullDemand[index] * Ucf[FLOW];
+        break;
 
     default:
         return 251;
@@ -2367,7 +2386,7 @@ int DLLEXPORT EN_getnodevalues(EN_Project p, int property, double *values)
 **----------------------------------------------------------------
 */
 {
-    int status = 0, i = 0;
+    int i, status = 0;
 
     for (i = 1; i <= p->network.Nnodes; i++)
     {
@@ -3352,6 +3371,8 @@ int DLLEXPORT EN_addlink(EN_Project p, const char *id, int linkType,
     }
     link->Kb = 0;
     link->Kw = 0;
+    link->LeakArea = 0;
+    link->LeakExpan = 0;
     link->R = 0;
     link->Rc = 0;
     link->Rpt = 0;
@@ -3923,6 +3944,18 @@ int DLLEXPORT EN_getlinkvalue(EN_Project p, int index, int property, double *val
         v = (double)incontrols(p, LINK, index);
         break;
 
+    case EN_LEAK_AREA:
+        v = Link[index].LeakArea * Ucf[LENGTH];
+        break;
+
+    case EN_LEAK_EXPAN:
+        v = Link[index].LeakExpan * Ucf[LENGTH];
+        break;
+
+    case EN_LINK_LEAKAGE:
+        v = leakage_getlinkleakage(p, index) * Ucf[FLOW];
+        break;
+
     default:
         return 251;
     }
@@ -3939,7 +3972,7 @@ int DLLEXPORT EN_getlinkvalues(EN_Project p, int property, double *values)
 **----------------------------------------------------------------
 */
 {
-    int status = 0, i = 0;
+    int i, status = 0;
     for(i = 1; i <= p->network.Nlinks; i++)
     {
         status = EN_getlinkvalue(p, i, property, &values[i-1]);
@@ -4163,6 +4196,16 @@ int DLLEXPORT EN_setlinkvalue(EN_Project p, int index, int property, double valu
         }
         break;
 
+    case EN_LEAK_AREA:  // leak area in sq mm per 100 pipe length units
+        if (value < 0.0) return 211;
+        Link[index].LeakArea = value / Ucf[LENGTH];
+        break;
+
+    case EN_LEAK_EXPAN:  // leak area expansion slope (sq mm per unit of head)
+        if (value < 0.0) return 211;
+        Link[index].LeakExpan = value / Ucf[LENGTH];
+        break;
+
     default:
         return 251;
     }

src/flowbalance.c

@@ -0,0 +1,191 @@
+/*
+ ******************************************************************************
+ Project:      OWA EPANET
+ Version:      2.3
+ Module:       flowbalance.c
+ Description:  computes components of network's flow balance
+ Authors:      see AUTHORS
+ Copyright:    see AUTHORS
+ License:      see LICENSE
+ Last Updated: 06/14/2024
+ ******************************************************************************
+*/
+
+#include "types.h"
+
+// Exported functions (declared in funcs.h)
+//void flowbalance_start(Project *);
+//void flowbalance_update(Project *, long);
+//void flowbalance_end(Project *);
+
+void flowbalance_start(Project *pr)
+/*
+**-------------------------------------------------------------------
+**  Input:   none
+**  Output:  none
+**  Purpose: initializes components of the network's flow balance.
+**-------------------------------------------------------------------
+*/
+{
+    Hydraul *hyd = &pr->hydraul;
+    hyd->FlowBalance.totalInflow = 0.0;
+    hyd->FlowBalance.totalOutflow = 0.0;
+    hyd->FlowBalance.consumerDemand = 0.0;
+    hyd->FlowBalance.emitterDemand = 0.0;
+    hyd->FlowBalance.leakageDemand = 0.0;
+    hyd->FlowBalance.deficitDemand = 0.0;
+    hyd->FlowBalance.storageDemand = 0.0;
+    hyd->FlowBalance.ratio = 0.0;
+}    
+
+void flowbalance_update(Project *pr, long hstep)
+/*
+**-------------------------------------------------------------------
+**  Input:   hstep = time step (sec)
+**  Output:  none
+**  Purpose: updates components of the system flow balance.
+**-------------------------------------------------------------------
+*/
+{
+    Network *net = &pr->network;
+    Hydraul *hyd = &pr->hydraul;
+    Times   *time = &pr->times;
+
+    int i, j;
+    double v, dt, deficit, fullDemand;
+    SflowBalance flowBalance;
+
+    // Determine current time interval in seconds
+    if (time->Dur == 0) dt = 1.0;
+    else if (time->Htime < time->Dur)
+    {
+        dt = (double) hstep;
+    }
+    else return;
+
+    // Initialize local flow balance
+    flowBalance.totalInflow = 0.0;
+    flowBalance.totalOutflow = 0.0;
+    flowBalance.consumerDemand = 0.0;
+    flowBalance.emitterDemand = 0.0;
+    flowBalance.leakageDemand = 0.0;
+    flowBalance.deficitDemand = 0.0;
+    flowBalance.storageDemand = 0.0;
+    fullDemand = 0.0;
+
+    // Initialize demand deficiency & leakage loss
+    hyd->DeficientNodes = 0;
+    hyd->DemandReduction = 0.0;
+    hyd->LeakageLoss = 0.0;
+
+    // Examine each junction node
+    for (i = 1; i <= net->Njuncs; i++) 
+    {
+        // Accumulate consumer demand flow
+        v = hyd->DemandFlow[i];
+        if (v < 0.0)
+            flowBalance.totalInflow += (-v);
+        else
+        {
+            fullDemand += hyd->FullDemand[i];
+            flowBalance.consumerDemand += v;
+            flowBalance.totalOutflow += v;
+        }
+
+        // Accumulate emitter and leakage flow
+        v = hyd->EmitterFlow[i];
+        flowBalance.emitterDemand += v;
+        flowBalance.totalOutflow += v;
+        v = hyd->LeakageFlow[i];
+        flowBalance.leakageDemand += v;
+        flowBalance.totalOutflow += v;
+
+        // Accumulate demand deficit flow
+        if (hyd->DemandModel == PDA && hyd->FullDemand[i] > 0.0)
+        {
+            deficit = hyd->FullDemand[i] - hyd->DemandFlow[i];
+            if (deficit > TINY)
+            {
+                hyd->DeficientNodes++;
+                flowBalance.deficitDemand += deficit;
+            }
+        }
+    }
+
+    // Examine each tank/reservoir node
+    for (j = 1; j <= net->Ntanks; j++)
+    {
+        i = net->Tank[j].Node;
+        v = hyd->NodeDemand[i];
+
+        // For a snapshot analysis or a reservoir node
+        if (time->Dur == 0 || net->Tank[j].A == 0.0)
+        {
+            if (v >= 0.0)
+                flowBalance.totalOutflow += v;
+            else
+                flowBalance.totalInflow += (-v);
+        }
+
+        // For tank under extended period analysis
+        else
+            flowBalance.storageDemand += v;
+    }
+
+    // Find % demand reduction & % leakage for current period
+    if (fullDemand > 0.0)
+        hyd->DemandReduction = flowBalance.deficitDemand / fullDemand * 100.0;
+    if (flowBalance.totalInflow > 0.0)
+        hyd->LeakageLoss = flowBalance.leakageDemand / flowBalance.totalInflow * 100.0;
+
+    // Update flow balance for entire run
+    hyd->FlowBalance.totalInflow += flowBalance.totalInflow * dt;
+    hyd->FlowBalance.totalOutflow += flowBalance.totalOutflow * dt;
+    hyd->FlowBalance.consumerDemand += flowBalance.consumerDemand * dt;
+    hyd->FlowBalance.emitterDemand += flowBalance.emitterDemand * dt;
+    hyd->FlowBalance.leakageDemand += flowBalance.leakageDemand * dt;
+    hyd->FlowBalance.deficitDemand += flowBalance.deficitDemand * dt;
+    hyd->FlowBalance.storageDemand += flowBalance.storageDemand * dt;
+}
+
+void flowbalance_end(Project *pr)
+/*
+**-------------------------------------------------------------------
+**  Input:   none
+**  Output:  none
+**  Purpose: finalizes components of the system flow balance.
+**-------------------------------------------------------------------
+*/
+{
+    Hydraul *hyd = &pr->hydraul;
+    Times   *time = &pr->times;
+
+    double seconds, qin, qout, qstor, r;
+
+    if (time->Htime > 0)
+        seconds = time->Htime;
+    else
+        seconds = 1.0;        
+    hyd->FlowBalance.totalInflow /= seconds;
+    hyd->FlowBalance.totalOutflow /= seconds;
+    hyd->FlowBalance.consumerDemand /= seconds;
+    hyd->FlowBalance.emitterDemand /= seconds;
+    hyd->FlowBalance.leakageDemand /= seconds;
+    hyd->FlowBalance.deficitDemand /= seconds;
+    hyd->FlowBalance.storageDemand /= seconds;
+
+    qin = hyd->FlowBalance.totalInflow;
+    qout = hyd->FlowBalance.totalOutflow;
+    qstor = hyd->FlowBalance.storageDemand;
+    if (qstor > 0.0)
+        qout += qstor;
+    else
+        qin -= qstor;
+    if (qin == qout)
+        r = 1.0;
+    else if (qin > 0.0)
+        r = qout / qin;
+    else
+        r = 0.0;
+    hyd->FlowBalance.ratio = r;
+}