[NL] added residual convergence criterion

NumLib/ODESolver/ConvergenceCriterion.cpp

@@ -12,6 +12,7 @@
 #include "BaseLib/Error.h"
 
 #include "ConvergenceCriterionDeltaX.h"
+#include "ConvergenceCriterionResidual.h"
 #include "ConvergenceCriterionPerComponentDeltaX.h"
 
 namespace NumLib
@@ -24,6 +25,8 @@ std::unique_ptr<ConvergenceCriterion> createConvergenceCriterion(
 
     if (type == "DeltaX") {
         return createConvergenceCriterionDeltaX(config);
+    } else if (type == "Residual") {
+        return createConvergenceCriterionResidual(config);
     } else if (type == "PerComponentDeltaX") {
         return createConvergenceCriterionPerComponentDeltaX(config);
     }

NumLib/ODESolver/ConvergenceCriterionResidual.cpp

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#include "ConvergenceCriterionResidual.h"
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "MathLib/LinAlg/LinAlg.h"
+
+namespace NumLib
+{
+ConvergenceCriterionResidual::ConvergenceCriterionResidual(
+    boost::optional<double>&& absolute_tolerance,
+    boost::optional<double>&& relative_tolerance,
+    MathLib::VecNormType norm_type)
+    : _abstol(std::move(absolute_tolerance)),
+      _reltol(std::move(relative_tolerance)),
+      _norm_type(norm_type)
+{
+    if ((!_abstol) && (!_reltol))
+        OGS_FATAL(
+            "At least one of absolute or relative tolerance has to be "
+            "specified.");
+}
+
+void ConvergenceCriterionResidual::checkResidual(const GlobalVector& residual)
+{
+    auto norm_res = MathLib::LinAlg::norm(residual, _norm_type);
+
+    if (_is_first_iteration) {
+        INFO("Convergence criterion: |r0|=%.4e", norm_res);
+        _residual_norm_0 = norm_res;
+    } else {
+        INFO("Convergence criterion: |r|=%.4e |r0|=%.4e |r|/|r0|=%.4e",
+             norm_res, _residual_norm_0, norm_res / _residual_norm_0);
+    }
+
+    bool satisfied_abs = false;
+    bool satisfied_rel = false;
+
+    if (_abstol) {
+        satisfied_abs = norm_res < *_abstol;
+    }
+    if (_reltol && !_is_first_iteration) {
+        satisfied_rel = norm_res < *_reltol * _residual_norm_0;
+    }
+
+    _satisfied = _satisfied && (satisfied_abs || satisfied_rel);
+}
+
+std::unique_ptr<ConvergenceCriterionResidual>
+createConvergenceCriterionResidual(const BaseLib::ConfigTree& config)
+{
+    config.checkConfigParameter("type", "Residual");
+
+    auto abstol = config.getConfigParameterOptional<double>("abstol");
+    auto reltol = config.getConfigParameterOptional<double>("reltol");
+    auto const norm_type_str =
+        config.getConfigParameter<std::string>("norm_type");
+    auto const norm_type = MathLib::convertStringToVecNormType(norm_type_str);
+
+    if (norm_type == MathLib::VecNormType::INVALID)
+        OGS_FATAL("Unknown vector norm type `%s'.", norm_type_str.c_str());
+
+    return std::unique_ptr<ConvergenceCriterionResidual>(
+        new ConvergenceCriterionResidual(std::move(abstol), std::move(reltol),
+                                         norm_type));
+}
+
+}  // NumLib

NumLib/ODESolver/ConvergenceCriterionResidual.h

+/**
+ * \copyright
+ * Copyright (c) 2012-2016, OpenGeoSys Community (http://www.opengeosys.org)
+ *            Distributed under a Modified BSD License.
+ *              See accompanying file LICENSE.txt or
+ *              http://www.opengeosys.org/project/license
+ *
+ */
+
+#ifndef NUMLIB_CONVERGENCECRITERIONRESIDUAL_H
+#define NUMLIB_CONVERGENCECRITERIONRESIDUAL_H
+
+#include <boost/optional.hpp>
+#include "ConvergenceCriterion.h"
+#include "MathLib/LinAlg/LinAlgEnums.h"
+
+namespace NumLib
+{
+class ConvergenceCriterionResidual final : public ConvergenceCriterion
+{
+public:
+    explicit ConvergenceCriterionResidual(
+        boost::optional<double>&& absolute_tolerance,
+        boost::optional<double>&& relative_tolerance,
+        MathLib::VecNormType norm_type);
+
+    bool hasDeltaXCheck() const override { return false; }
+    bool hasResidualCheck() const override { return true; }
+
+    void checkDeltaX(const GlobalVector& /*minus_delta_x*/,
+                     GlobalVector const& /*x*/) override {}
+    void checkResidual(const GlobalVector& residual) override;
+
+    void preFirstIteration() override { _is_first_iteration = true; }
+    void reset() override { _satisfied = true; _is_first_iteration = false; }
+    bool isSatisfied() const { return _satisfied; }
+private:
+    const boost::optional<double> _abstol;
+    const boost::optional<double> _reltol;
+    const MathLib::VecNormType _norm_type;
+    bool _satisfied = true;
+    bool _is_first_iteration = true;
+    double _residual_norm_0;
+};
+
+std::unique_ptr<ConvergenceCriterionResidual> createConvergenceCriterionResidual(
+    BaseLib::ConfigTree const& config);
+
+}  // NumLib
+
+#endif  // NUMLIB_CONVERGENCECRITERIONRESIDUAL_H