[NL] added per component residual check

NumLib/ODESolver/ConvergenceCriterion.cpp

@@ -14,6 +14,7 @@
 #include "ConvergenceCriterionDeltaX.h"
 #include "ConvergenceCriterionResidual.h"
 #include "ConvergenceCriterionPerComponentDeltaX.h"
+#include "ConvergenceCriterionPerComponentResidual.h"
 
 namespace NumLib
 {
@@ -29,6 +30,8 @@ std::unique_ptr<ConvergenceCriterion> createConvergenceCriterion(
         return createConvergenceCriterionResidual(config);
     } else if (type == "PerComponentDeltaX") {
         return createConvergenceCriterionPerComponentDeltaX(config);
+    } else if (type == "PerComponentResidual") {
+        return createConvergenceCriterionPerComponentResidual(config);
     }
 
     OGS_FATAL("There is no convergence criterion of type `%s'.", type.c_str());

NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.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 "ConvergenceCriterionPerComponentResidual.h"
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+
+namespace NumLib
+{
+ConvergenceCriterionPerComponentResidual::
+    ConvergenceCriterionPerComponentResidual(
+        std::vector<double>&& absolute_tolerances,
+        std::vector<double>&& relative_tolerances,
+        MathLib::VecNormType norm_type)
+    : _abstols(std::move(absolute_tolerances)),
+      _reltols(std::move(relative_tolerances)),
+      _norm_type(norm_type),
+      _residual_norms_0(_abstols.size())
+{
+    if (_abstols.size() != _reltols.size())
+        OGS_FATAL(
+            "The number of absolute and relative tolerances given must be the "
+            "same.");
+
+    if (_abstols.empty())
+        OGS_FATAL("The given tolerances vector is empty.");
+}
+
+void ConvergenceCriterionPerComponentResidual::checkResidual(
+    const GlobalVector& residual)
+{
+    if ((!_dof_table) || (!_mesh))
+        OGS_FATAL("D.o.f. table or mesh have not been set.");
+
+    bool satisfied_abs = true;
+    // Make sure that in the first iteration the relative residual tolerance is
+    // not satisfied.
+    bool satisfied_rel = _is_first_iteration ? false : true;
+
+    for (unsigned global_component = 0; global_component < _abstols.size();
+         ++global_component)
+    {
+        // TODO short cut if tol <= 0.0
+        auto norm_res = norm(residual, global_component, _norm_type,
+                             *_dof_table, *_mesh);
+
+        if (_is_first_iteration) {
+            INFO("Convergence criterion: |r0|=%.4e", norm_res);
+            _residual_norms_0[global_component] = norm_res;
+        } else {
+            auto const norm_res0 = _residual_norms_0[global_component];
+            INFO(
+                "Convergence criterion, component %u: |r|=%.4e, |r0|=%.4e, "
+                "|r|/|r0|=%.4e",
+                norm_res, global_component, norm_res0, norm_res / norm_res0);
+        }
+
+        satisfied_abs = satisfied_abs && norm_res < _abstols[global_component];
+        satisfied_rel = satisfied_rel &&
+                        norm_res < _reltols[global_component] *
+                                       _residual_norms_0[global_component];
+    }
+
+    _satisfied = _satisfied && (satisfied_abs || satisfied_rel);
+}
+
+void ConvergenceCriterionPerComponentResidual::setDOFTable(
+    const LocalToGlobalIndexMap& dof_table, MeshLib::Mesh const& mesh)
+{
+    _dof_table = &dof_table;
+    _mesh = &mesh;
+
+    if (_dof_table->getNumberOfComponents() != _abstols.size())
+        OGS_FATAL(
+            "The number of components in the DOF table and the number of "
+            "tolerances given do not match.");
+}
+
+std::unique_ptr<ConvergenceCriterionPerComponentResidual>
+createConvergenceCriterionPerComponentResidual(
+    const BaseLib::ConfigTree& config)
+{
+    config.checkConfigParameter("type", "PerComponentResidual");
+
+    auto abstols =
+        config.getConfigParameterOptional<std::vector<double>>("abstols");
+    auto reltols =
+        config.getConfigParameterOptional<std::vector<double>>("reltols");
+    auto const norm_type_str =
+        config.getConfigParameter<std::string>("norm_type");
+
+    if ((!abstols) && (!reltols))
+        OGS_FATAL(
+            "At least one of absolute or relative tolerance has to be "
+            "specified.");
+    if (!abstols) {
+        abstols = std::vector<double>(reltols->size());
+    } else if (!reltols) {
+        reltols = std::vector<double>(abstols->size());
+    }
+
+    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<ConvergenceCriterionPerComponentResidual>(
+        new ConvergenceCriterionPerComponentResidual(
+            std::move(*abstols), std::move(*reltols), norm_type));
+}
+
+}  // NumLib

NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.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_CONVERGENCECRITERIONPERCOMPONENTRESIDUAL_H
+#define NUMLIB_CONVERGENCECRITERIONPERCOMPONENTRESIDUAL_H
+
+#include <vector>
+#include "MathLib/LinAlg/LinAlgEnums.h"
+#include "ConvergenceCriterionPerComponent.h"
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+
+class ConvergenceCriterionPerComponentResidual
+    : public ConvergenceCriterionPerComponent
+{
+public:
+    explicit ConvergenceCriterionPerComponentResidual(
+        std::vector<double>&& absolute_tolerances,
+        std::vector<double>&& relative_tolerances,
+        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 override { return _satisfied; }
+
+    void setDOFTable(const LocalToGlobalIndexMap& dof_table,
+                     MeshLib::Mesh const& mesh) override;
+
+private:
+    const std::vector<double> _abstols;
+    const std::vector<double> _reltols;
+    const MathLib::VecNormType _norm_type;
+    bool _satisfied = true;
+    LocalToGlobalIndexMap const* _dof_table = nullptr;
+    MeshLib::Mesh const* _mesh = nullptr;
+    bool _is_first_iteration = true;
+    std::vector<double> _residual_norms_0;
+};
+
+std::unique_ptr<ConvergenceCriterionPerComponentResidual>
+createConvergenceCriterionPerComponentResidual(
+    BaseLib::ConfigTree const& config);
+
+}  // namespace NumLib
+
+#endif  // NUMLIB_CONVERGENCECRITERIONPERCOMPONENTRESIDUAL_H