From 4ff7f49ed8fc8a2fb9e4efde25d68da7a32370f4 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Tue, 23 Aug 2016 09:53:15 +0200
Subject: [PATCH] [NL] added function for relative tolerance check
---
NumLib/ODESolver/ConvergenceCriterion.cpp | 8 ++++++++
NumLib/ODESolver/ConvergenceCriterion.h | 7 +++++++
NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp | 2 +-
.../ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp | 3 ++-
.../ConvergenceCriterionPerComponentResidual.cpp | 7 ++++---
NumLib/ODESolver/ConvergenceCriterionResidual.cpp | 3 ++-
6 files changed, 24 insertions(+), 6 deletions(-)
diff --git a/NumLib/ODESolver/ConvergenceCriterion.cpp b/NumLib/ODESolver/ConvergenceCriterion.cpp
index e89a854f690..42191dc81cf 100644
--- a/NumLib/ODESolver/ConvergenceCriterion.cpp
+++ b/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
diff --git a/NumLib/ODESolver/ConvergenceCriterion.h b/NumLib/ODESolver/ConvergenceCriterion.h
index aaf07e8c15f..45d473d8823 100644
--- a/NumLib/ODESolver/ConvergenceCriterion.h
+++ b/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
diff --git a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
index 5ab4ab6edbd..2c2249c1dbd 100644
--- a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp
+++ b/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);
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
index fe5f73b4528..9615ab813d1 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp
+++ b/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);
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
index dc94fb9b6e9..155857c71fc 100644
--- a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
+++ b/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);
diff --git a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
index ff0d719c078..141245a4492 100644
--- a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp
+++ b/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);
--
GitLab