From 4c5cca1a501638d7669e7fcdd037f1308c6bd8e4 Mon Sep 17 00:00:00 2001
From: Dmitri Naumov <github@naumov.de>
Date: Wed, 27 Jan 2021 00:19:57 +0100
Subject: [PATCH] [NL] Local Newton. Add increment tolerance.

The increment tolerance is used as an additional convergence
criterion.

Introduction of the new tolerance requires renaming of the
residuum tolerance.
---
 .../nonlinear_solver/t_error_tolerance.md     |  5 +++-
 .../nonlinear_solver/t_increment_tolerance.md |  1 +
 .../nonlinear_solver/t_residuum_tolerance.md  |  1 +
 .../CreateNewtonRaphsonSolverParameters.cpp   | 30 +++++++++++++++----
 NumLib/NewtonRaphson.h                        | 21 +++++++++----
 ...TwoPhaseFlowWithPrhoMaterialProperties.cpp |  5 ++--
 6 files changed, 50 insertions(+), 13 deletions(-)
 create mode 100644 Documentation/ProjectFile/nonlinear_solver/t_increment_tolerance.md
 create mode 100644 Documentation/ProjectFile/nonlinear_solver/t_residuum_tolerance.md

diff --git a/Documentation/ProjectFile/nonlinear_solver/t_error_tolerance.md b/Documentation/ProjectFile/nonlinear_solver/t_error_tolerance.md
index ccae39ea3ad..d6d98d1f91f 100644
--- a/Documentation/ProjectFile/nonlinear_solver/t_error_tolerance.md
+++ b/Documentation/ProjectFile/nonlinear_solver/t_error_tolerance.md
@@ -1 +1,4 @@
-\copydoc NumLib::NewtonRaphson::_tolerance_squared
+Deprecated flag for residuum tolerance. Use `residuum_tolerance` and
+`increment_tolerance` in the future.
+
+If set, increment tolerance will be set to zero.
diff --git a/Documentation/ProjectFile/nonlinear_solver/t_increment_tolerance.md b/Documentation/ProjectFile/nonlinear_solver/t_increment_tolerance.md
new file mode 100644
index 00000000000..9a26ebd7ac4
--- /dev/null
+++ b/Documentation/ProjectFile/nonlinear_solver/t_increment_tolerance.md
@@ -0,0 +1 @@
+\copydoc NumLib::NewtonRaphson::_increment_tolerance_squared
diff --git a/Documentation/ProjectFile/nonlinear_solver/t_residuum_tolerance.md b/Documentation/ProjectFile/nonlinear_solver/t_residuum_tolerance.md
new file mode 100644
index 00000000000..80989e1f68c
--- /dev/null
+++ b/Documentation/ProjectFile/nonlinear_solver/t_residuum_tolerance.md
@@ -0,0 +1 @@
+\copydoc NumLib::NewtonRaphson::_residuum_tolerance_squared
diff --git a/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.cpp b/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.cpp
index c037578659a..35fa51b8c0c 100644
--- a/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.cpp
+++ b/MaterialLib/SolidModels/CreateNewtonRaphsonSolverParameters.cpp
@@ -29,10 +29,30 @@ NumLib::NewtonRaphsonSolverParameters createNewtonRaphsonSolverParameters(
 
     auto const error_tolerance =
         //! \ogs_file_param{nonlinear_solver__error_tolerance}
-        config.getConfigParameter<double>("error_tolerance");
-
-    DBUG("\terror_tolerance: {:g}.", error_tolerance);
-
-    return {maximum_iterations, error_tolerance};
+        config.getConfigParameterOptional<double>("error_tolerance");
+    if (error_tolerance)
+    {
+        WARN(
+            "The 'error_tolerance' tag for the Newton-Raphson solver is "
+            "deprecated.\n"
+            "Use new tags 'residuum_tolerance' and 'increment_tolerance'.\n"
+            "For now we use residuum_tolerance {} and increment_tolerance 0.",
+            *error_tolerance);
+        return {maximum_iterations, *error_tolerance, 0};
+    }
+
+    auto const residuum_tolerance =
+        //! \ogs_file_param{nonlinear_solver__residuum_tolerance}
+        config.getConfigParameter<double>("residuum_tolerance");
+
+    DBUG("\tresiduum_tolerance: {:g}.", residuum_tolerance);
+
+    auto const increment_tolerance =
+        //! \ogs_file_param{nonlinear_solver__increment_tolerance}
+        config.getConfigParameter<double>("increment_tolerance");
+
+    DBUG("\tincrement_tolerance: {:g}.", increment_tolerance);
+
+    return {maximum_iterations, residuum_tolerance, increment_tolerance};
 }
 }  // namespace MaterialLib
diff --git a/NumLib/NewtonRaphson.h b/NumLib/NewtonRaphson.h
index 633063cf2f9..b78af5dfafb 100644
--- a/NumLib/NewtonRaphson.h
+++ b/NumLib/NewtonRaphson.h
@@ -21,7 +21,8 @@ namespace NumLib
 struct NewtonRaphsonSolverParameters
 {
     int const maximum_iterations;
-    double const error_tolerance;
+    double const residuum_tolerance;
+    double const increment_tolerance;
 };
 
 /// Newton-Raphson solver for system of equations using an Eigen linear solvers
@@ -44,8 +45,10 @@ public:
           _residual_update(residual_update),
           _solution_update(solution_update),
           _maximum_iterations(solver_parameters.maximum_iterations),
-          _tolerance_squared(solver_parameters.error_tolerance *
-                             solver_parameters.error_tolerance)
+          _residuum_tolerance_squared(solver_parameters.residuum_tolerance *
+                                      solver_parameters.residuum_tolerance),
+          _increment_tolerance_squared(solver_parameters.increment_tolerance *
+                                       solver_parameters.increment_tolerance)
     {
     }
 
@@ -64,7 +67,7 @@ public:
             _jacobian_update(jacobian);
             _residual_update(residual);
 
-            if (residual.squaredNorm() < _tolerance_squared)
+            if (residual.squaredNorm() < _residuum_tolerance_squared)
             {
                 break;  // convergence criteria fulfilled.
             }
@@ -76,6 +79,11 @@ public:
 
             _solution_update(increment);
 
+            if (increment.squaredNorm() < _increment_tolerance_squared)
+            {
+                break;  // increment to small.
+            }
+
             // DBUG("Local Newton: Iteration #{:d} |dx| = {:g}, |r| = {:g}",
             //      iteration, increment.norm(), residual.norm());
         } while (iteration++ < _maximum_iterations);
@@ -99,6 +107,9 @@ private:
     ResidualUpdate _residual_update;
     SolutionUpdate _solution_update;
     const int _maximum_iterations;  ///< Maximum number of iterations.
-    const double _tolerance_squared;    ///< Error tolerance for the residual.
+    const double
+        _residuum_tolerance_squared;  ///< Error tolerance for the residuum.
+    const double
+        _increment_tolerance_squared;  ///< Error tolerance for the increment.
 };
 }  // namespace NumLib
diff --git a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
index 126341f582a..b15d56a9cf0 100644
--- a/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
+++ b/ProcessLib/TwoPhaseFlowWithPrho/TwoPhaseFlowWithPrhoMaterialProperties.cpp
@@ -205,14 +205,15 @@ bool TwoPhaseFlowWithPrhoMaterialProperties::computeConstitutiveRelation(
         // criteria available from the input file configuration. See Ehlers
         // material model implementation for the example.
         const int maximum_iterations(20);
-        const double tolerance(1.e-14);
+        const double residuum_tolerance(1.e-14);
+        const double increment_tolerance(0);
 
         auto newton_solver = NumLib::NewtonRaphson<
             decltype(linear_solver), LocalJacobianMatrix,
             decltype(update_jacobian), LocalResidualVector,
             decltype(update_residual), decltype(update_solution)>(
             linear_solver, update_jacobian, update_residual, update_solution,
-            {maximum_iterations, tolerance});
+            {maximum_iterations, residuum_tolerance, increment_tolerance});
 
         auto const success_iterations = newton_solver.solve(J_loc);
 
-- 
GitLab