From 6a0f56947687d6fdfa10866e1fd5fd57e1dff666 Mon Sep 17 00:00:00 2001
From: Paul Stapor <paul.stapor@helmholtz-muenchen.de>
Date: Mon, 30 Jul 2018 15:36:49 +0200
Subject: [PATCH] Fix state positivity in Newton solver (#378)

* Fixed the handling of state positivity in the Newton solver (Improving convergence).
---
 src/steadystateproblem.cpp | 24 +++++++++++++-----------
 1 file changed, 13 insertions(+), 11 deletions(-)

diff --git a/src/steadystateproblem.cpp b/src/steadystateproblem.cpp
index 55c822b36f..02318c9500 100644
--- a/src/steadystateproblem.cpp
+++ b/src/steadystateproblem.cpp
@@ -122,11 +122,6 @@ void SteadystateProblem::applyNewtonsMethod(ReturnData *rdata,
         Ensure positivity of the state */
     x_newton = *x;
     N_VAbs(x_newton.getNVector(), x_newton.getNVector());
-    for (ix = 0; ix < model->nx; ix++) {
-        if (x_newton[ix] < newtonSolver->atol) {
-            x_newton[ix] = newtonSolver->atol;
-        }
-    }
     N_VDiv(xdot.getNVector(), x_newton.getNVector(), rel_x_newton.getNVector());
     double res_rel = sqrt(N_VDotProd(rel_x_newton.getNVector(), rel_x_newton.getNVector()));
     x_old = *x;
@@ -149,10 +144,6 @@ void SteadystateProblem::applyNewtonsMethod(ReturnData *rdata,
         
         /* Try a full, undamped Newton step */
         N_VLinearSum(1.0, x_old.getNVector(), gamma, delta.getNVector(), x->getNVector());
-        /* Ensure positivity of the state */
-        for (ix = 0; ix < model->nx; ix++)
-            if ((*x)[ix] < newtonSolver->atol)
-                (*x)[ix] = newtonSolver->atol;
         
         /* Compute new xdot and residuals */
         model->fxdot(*t, x, &dx, &xdot);
@@ -169,8 +160,19 @@ void SteadystateProblem::applyNewtonsMethod(ReturnData *rdata,
             compNewStep = TRUE;
             /* Check residuals vs tolerances */
             converged = (res_abs < newtonSolver->atol) || (res_rel < newtonSolver->rtol);
-            /* increase dampening factor (superfluous, if converged) */
-            gamma = fmin(1.0, 2.0 * gamma);
+            
+            if (converged) {
+                /* Ensure positivity of the found state */
+                for (ix = 0; ix < model->nx; ix++) {
+                    if ((*x)[ix] < 0.0) {
+                        (*x)[ix] = 0.0;
+                        converged = FALSE;
+                    }
+                }
+            } else {
+                /* increase dampening factor (superfluous, if converged) */
+                gamma = fmin(1.0, 2.0 * gamma);
+            }
         } else {
             /* Reduce dampening factor */
             gamma = gamma / 4.0;