From 409a6004556d4796e9b340052e15d10d76a8718c Mon Sep 17 00:00:00 2001
From: Lew Rossman <LRossman@cinci.rr.com>
Date: Fri, 16 Apr 2021 17:01:54 -0400
Subject: [PATCH] Improved flow checks for full and empty tanks

Improves  the way link flow into full and out of empty tanks is checked. Also corrects the head loss and gradient calculation for constant horsepower pumps at very low and very high flow rates (this was affecting one of the examples used to test the full/empty tank fix).
---
 src/hydcoeffs.c |  8 ++++++--
 src/hydstatus.c | 17 +++++++++++++----
 2 files changed, 19 insertions(+), 6 deletions(-)

diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c
index 9f4ac3d..e044dff 100644
--- a/src/hydcoeffs.c
+++ b/src/hydcoeffs.c
@@ -691,6 +691,8 @@ void  pumpcoeff(Project *pr, int k)
            hgrad;            // Head loss gradient
     Spump  *pump;
 
+    double qstar;
+
     // Use high resistance pipe if pump closed or cannot deliver head
     setting = hyd->LinkSetting[k];
     if (hyd->LinkStatus[k] <= CLOSED || setting == 0.0)
@@ -748,12 +750,14 @@ void  pumpcoeff(Project *pr, int k)
             if (hgrad > CBIG)
             {
                 hgrad = CBIG;
-                hloss = -hgrad * hyd->LinkFlow[k];
+                qstar = sqrt(-r / hgrad);
+                hloss = (r / qstar) - hgrad * (qstar - q);
             }
             else if (hgrad < hyd->RQtol)
             {
                 hgrad = hyd->RQtol;
-                hloss = -hgrad * hyd->LinkFlow[k];
+                qstar = sqrt(-r / hgrad);
+                hloss = (r / qstar) - hgrad * (qstar - q);
             }
             // ... otherwise compute head loss from pump curve
             else
diff --git a/src/hydstatus.c b/src/hydstatus.c
index 3d7c106..4013ec9 100644
--- a/src/hydstatus.c
+++ b/src/hydstatus.c
@@ -437,12 +437,21 @@ void  tankstatus(Project *pr, int k, int n1, int n2)
     tank = &net->Tank[i];
     if (tank->A == 0.0) return;
 
-    if (hyd->NodeHead[n1] >= tank->Hmax - hyd->Htol && !tank->CanOverflow)
+    // Can't add flow to a full tank
+    if (hyd->NodeHead[n1] >= tank->Hmax && !tank->CanOverflow)
     {
-        if (q < 0.0) hyd->LinkStatus[k] = TEMPCLOSED;
+        if (link->Type == PUMP && link->N2 == n1)
+            hyd->LinkStatus[k] = TEMPCLOSED;
+        else if (q < 0.0)
+            hyd->LinkStatus[k] = TEMPCLOSED;
     }
-    if (hyd->NodeHead[n1] <= tank->Hmin + hyd->Htol)
+
+    // Can't remove flow from an empty tank
+    if (hyd->NodeHead[n1] <= tank->Hmin)
     {
-        if (q > 0.0) hyd->LinkStatus[k] = TEMPCLOSED;
+        if (link->Type == PUMP && link->N1 == n1)
+            hyd->LinkStatus[k] = TEMPCLOSED;
+        else if (q > 0.0)
+            hyd->LinkStatus[k] = TEMPCLOSED;
     }
 }