From 533466b796418e2796b61a75068cdec7e0ace140 Mon Sep 17 00:00:00 2001
From: Lew Rossman <LRossman@cinci.rr.com>
Date: Thu, 19 Nov 2020 16:23:24 -0500
Subject: [PATCH] PDA changes to improve convergence

---
 src/hydcoeffs.c | 2 +-
 src/hydsolver.c | 4 ++++
 2 files changed, 5 insertions(+), 1 deletion(-)

diff --git a/src/hydcoeffs.c b/src/hydcoeffs.c
index 9f4ac3d..7b2dc47 100644
--- a/src/hydcoeffs.c
+++ b/src/hydcoeffs.c
@@ -486,7 +486,7 @@ void demandheadloss(Project *pr, int i, double dp, double n,
     }
 
     // Use power head loss function for demand less than full
-    else if (r < 1.0)
+    else if (r <= 1.0)                                                         //(2.2.1)
     {
         *hgrad = n * dp * pow(r, n - 1.0) / dfull;
         // ... use linear function for very small gradient
diff --git a/src/hydsolver.c b/src/hydsolver.c
index fcb8258..7f80a29 100644
--- a/src/hydsolver.c
+++ b/src/hydsolver.c
@@ -553,6 +553,10 @@ void newdemandflows(Project *pr, Hydbalance *hbal, double *qsum, double *dqsum)
         dh = hyd->NodeHead[i] - net->Node[i].El - hyd->Pmin;
         dq = (hloss - dh) / hgrad;
         dq *= hyd->RelaxFactor;
+
+        // Prevent a flow change greater than full demand                      //(2.2.1)
+        if (fabs(dq) > hyd->NodeDemand[i])
+            dq = 0.5 * SGN(dq) * hyd->NodeDemand[i];
         hyd->DemandFlow[i] -= dq;
 
         // Update system flow summation