[NL] added function for relative tolerance check

NumLib/ODESolver/ConvergenceCriterion.cpp

@@ -37,4 +37,12 @@ std::unique_ptr<ConvergenceCriterion> createConvergenceCriterion(
     OGS_FATAL("There is no convergence criterion of type `%s'.", type.c_str());
 }
 
+bool checkRelativeTolerance(const double reltol, const double numerator,
+                            const double denominator)
+{
+    auto const eps = std::numeric_limits<double>::epsilon();
+    return std::abs(numerator) <
+           std::abs(reltol) * (std::abs(denominator) + eps);
+}
+
 }  // NumLib

NumLib/ODESolver/ConvergenceCriterion.h

@@ -77,6 +77,13 @@ public:
 std::unique_ptr<ConvergenceCriterion> createConvergenceCriterion(
     BaseLib::ConfigTree const& config);
 
+//! Returns if |numerator/denominator| < |reltol|.
+//! This method copes with the case that denominator = 0 by always adding
+//! epsilon to the denominator.
+bool checkRelativeTolerance(double const reltol,
+                            double const numerator,
+                            double const denominator);
+
 } // namespace NumLib
 
 #endif  // NUMLIB_CONVERGENCECRITERION_H

NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp

@@ -45,7 +45,7 @@ void ConvergenceCriterionDeltaX::checkDeltaX(const GlobalVector& minus_delta_x,
         satisfied_abs = error_dx < *_abstol;
     }
     if (_reltol) {
-        satisfied_rel = error_dx < *_reltol * norm_x;
+        satisfied_rel = checkRelativeTolerance(*_reltol, error_dx, norm_x);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);

NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp

@@ -58,7 +58,8 @@ void ConvergenceCriterionPerComponentDeltaX::checkDeltaX(
 
         satisfied_abs = satisfied_abs && error_dx < _abstols[global_component];
         satisfied_rel =
-            satisfied_rel && error_dx < _reltols[global_component] * norm_x;
+            satisfied_rel && checkRelativeTolerance(_reltols[global_component],
+                                                    error_dx, norm_x);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);

NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp

@@ -65,9 +65,10 @@ void ConvergenceCriterionPerComponentResidual::checkResidual(
         }
 
         satisfied_abs = satisfied_abs && norm_res < _abstols[global_component];
-        satisfied_rel = satisfied_rel &&
-                        norm_res < _reltols[global_component] *
-                                       _residual_norms_0[global_component];
+        satisfied_rel =
+            satisfied_rel &&
+            checkRelativeTolerance(_reltols[global_component], norm_res,
+                                   _residual_norms_0[global_component]);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);

NumLib/ODESolver/ConvergenceCriterionResidual.cpp

@@ -48,7 +48,8 @@ void ConvergenceCriterionResidual::checkResidual(const GlobalVector& residual)
         satisfied_abs = norm_res < *_abstol;
     }
     if (_reltol && !_is_first_iteration) {
-        satisfied_rel = norm_res < *_reltol * _residual_norm_0;
+        satisfied_rel =
+            checkRelativeTolerance(*_reltol, norm_res, _residual_norm_0);
     }
 
     _satisfied = _satisfied && (satisfied_abs || satisfied_rel);