PDA changes to improve convergence

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

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