From 68b73a14f1a20daa292da8d40554ef2c40b74baf Mon Sep 17 00:00:00 2001
From: Lew Rossman <LRossman@outlook.com>
Date: Wed, 26 Jun 2024 17:50:24 -0400
Subject: [PATCH] Fixes for PDA

---
 src/flowbalance.c  | 25 ++++++++++---------------
 src/hydsolver.c    | 26 +++++++++++++++++++++-----
 tests/test_pda.cpp |  3 ---
 3 files changed, 31 insertions(+), 23 deletions(-)

diff --git a/src/flowbalance.c b/src/flowbalance.c
index 2c0efa66..aea6f199 100644
--- a/src/flowbalance.c
+++ b/src/flowbalance.c
@@ -73,9 +73,7 @@ void updateflowbalance(Project *pr, long hstep)
     flowBalance.storageDemand = 0.0;
     fullDemand = 0.0;
     
-    // Initialize demand deficiency & leakage loss
-    hyd->DeficientNodes = 0;
-    hyd->DemandReduction = 0.0;
+    // Initialize leakage loss
     hyd->LeakageLoss = 0.0;
 
     // Examine each junction node
@@ -104,11 +102,8 @@ void updateflowbalance(Project *pr, long hstep)
         if (hyd->DemandModel == PDA && hyd->FullDemand[i] > 0.0)
         {
             deficit = hyd->FullDemand[i] - hyd->DemandFlow[i];
-            if (deficit > TINY)
-            {
-                hyd->DeficientNodes++;
+            if (deficit > 0.0)
                 flowBalance.deficitDemand += deficit;
-            }
         }
     }
     
@@ -118,8 +113,8 @@ void updateflowbalance(Project *pr, long hstep)
         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)
+        // For a reservoir node
+        if (net->Tank[j].A == 0.0)
         {
             if (v >= 0.0)
                 flowBalance.totalOutflow += v;
@@ -127,16 +122,16 @@ void updateflowbalance(Project *pr, long hstep)
                 flowBalance.totalInflow += (-v);
         }
         
-        // For tank under extended period analysis
+        // For tank
         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;
+    // Find % leakage for current period
+    v = flowBalance.totalInflow;
+    if (flowBalance.storageDemand < 0.0) v += (-flowBalance.storageDemand);
+    if (v > 0.0)
+        hyd->LeakageLoss = flowBalance.leakageDemand / v * 100.0;
 
     // Update flow balance for entire run
     hyd->FlowBalance.totalInflow += flowBalance.totalInflow * dt;
diff --git a/src/hydsolver.c b/src/hydsolver.c
index d88823fb..5d76bc60 100644
--- a/src/hydsolver.c
+++ b/src/hydsolver.c
@@ -8,7 +8,7 @@
  Authors:      see AUTHORS
  Copyright:    see AUTHORS
  License:      see LICENSE
- Last Updated: 06/15/2024
+ Last Updated: 06/26/2024
  ******************************************************************************
 */
 
@@ -703,10 +703,15 @@ int pdaconverged(Project *pr)
     Hydraul *hyd = &pr->hydraul;
 
     const double QTOL = 0.0001;  // 0.0001 cfs ~= 0.005 gpm ~= 0.2 lpm)
-    int i;
+    int i, converged = 1;
+
+    double totalDemand = 0.0, totalReduction = 0.0;
     double dp = hyd->Preq - hyd->Pmin;
     double p, q, r;
-    
+
+    hyd->DeficientNodes = 0;
+    hyd->DemandReduction = 0.0;
+       
     // Examine each network junction
     for (i = 1; i <= pr->network.Njuncs; i++)
     {
@@ -726,9 +731,20 @@ int pdaconverged(Project *pr)
         }
         
         // Check if demand has not converged
-        if (fabs(q - hyd->DemandFlow[i]) > QTOL) return 0;
+        if (fabs(q - hyd->DemandFlow[i]) > QTOL)
+            converged = 0;
+
+        // Accumulate demand deficient node count and demand deficit
+        if (hyd->DemandFlow[i] + QTOL < hyd->FullDemand[i])
+        {
+            hyd->DeficientNodes++;
+            totalDemand += hyd->FullDemand[i];
+            totalReduction += hyd->FullDemand[i] - hyd->DemandFlow[i];
+        }
     }            
-    return 1;
+    if (totalDemand > 0.0)
+        hyd->DemandReduction = totalReduction / totalDemand * 100.0;
+    return converged;
 }
 
 
diff --git a/tests/test_pda.cpp b/tests/test_pda.cpp
index 7e951c0d..6091def7 100644
--- a/tests/test_pda.cpp
+++ b/tests/test_pda.cpp
@@ -79,9 +79,6 @@ BOOST_AUTO_TEST_CASE(test_pda_model)
     error = EN_getnodeindex(ph, (char *)"21", &index);
     BOOST_REQUIRE(error == 0);
     error = EN_getnodevalue(ph, index, EN_DEMANDDEFICIT, &reduction);
-
-    printf("\nreduction = %f", reduction);
-
     BOOST_REQUIRE(error == 0);
     BOOST_REQUIRE(abs(reduction - 413.67) < 0.01);