[NL] added conv crit implementations

NumLib/ODESolver/ConvergenceCriterionDeltaX.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 "ConvergenceCriterionDeltaX.h"
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "MathLib/LinAlg/LinAlg.h"
+
+namespace NumLib
+{
+ConvergenceCriterionDeltaX::ConvergenceCriterionDeltaX(
+    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 ConvergenceCriterionDeltaX::checkDeltaX(const GlobalVector& minus_delta_x,
+                                             GlobalVector const& x)
+{
+    auto error_dx = MathLib::LinAlg::norm(minus_delta_x, _norm_type);
+    auto norm_x = MathLib::LinAlg::norm(x, _norm_type);
+
+    INFO("Convergence criterion: |dx|=%.4e, |x|=%.4e, |dx|/|x|=%.4e", error_dx,
+         norm_x, error_dx / norm_x);
+
+    bool satisfied_abs = false;
+    bool satisfied_rel = false;
+
+    if (_abstol) {
+        satisfied_abs = error_dx < *_abstol;
+    }
+    if (_reltol) {
+        satisfied_rel = error_dx < *_reltol * norm_x;
+    }
+
+    _satisfied = _satisfied && (satisfied_abs || satisfied_rel);
+}
+
+std::unique_ptr<ConvergenceCriterionDeltaX> createConvergenceCriterionDeltaX(
+    const BaseLib::ConfigTree& config)
+{
+    config.checkConfigParameter("type", "DeltaX");
+
+    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<ConvergenceCriterionDeltaX>(
+        new ConvergenceCriterionDeltaX(std::move(abstol), std::move(reltol),
+                                       norm_type));
+}
+
+}  // NumLib

NumLib/ODESolver/ConvergenceCriterionDeltaX.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_CONVERGENCECRITERIONDELTAX_H
+#define NUMLIB_CONVERGENCECRITERIONDELTAX_H
+
+#include <boost/optional.hpp>
+#include "ConvergenceCriterion.h"
+#include "MathLib/LinAlg/LinAlgEnums.h"
+
+namespace NumLib
+{
+class ConvergenceCriterionDeltaX final : public ConvergenceCriterion
+{
+public:
+    explicit ConvergenceCriterionDeltaX(
+        boost::optional<double>&& absolute_tolerance,
+        boost::optional<double>&& relative_tolerance,
+        MathLib::VecNormType norm_type);
+
+    bool hasDeltaXCheck() const override { return true; }
+    bool hasResidualCheck() const override { return false; }
+
+    void checkDeltaX(const GlobalVector& minus_delta_x,
+                     GlobalVector const& x) override;
+    void checkResidual(const GlobalVector& /*residual*/) override {}
+
+    void reset() override { _satisfied = true; }
+    bool isSatisfied() const { return _satisfied; }
+private:
+    const boost::optional<double> _abstol;
+    const boost::optional<double> _reltol;
+    const MathLib::VecNormType _norm_type;
+    bool _satisfied = true;
+};
+
+std::unique_ptr<ConvergenceCriterionDeltaX> createConvergenceCriterionDeltaX(
+    BaseLib::ConfigTree const& config);
+
+}  // NumLib
+
+#endif  // NUMLIB_CONVERGENCECRITERIONDELTAX_H

NumLib/ODESolver/ConvergenceCriterionPerComponent.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_CONVERGENCECRITERIONPERCOMPONENT_H
+#define NUMLIB_CONVERGENCECRITERIONPERCOMPONENT_H
+
+#include "ConvergenceCriterion.h"
+
+namespace MeshLib
+{
+class Mesh;
+}  // namespace MeshLib
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+
+class ConvergenceCriterionPerComponent : public ConvergenceCriterion
+{
+public:
+    virtual void setDOFTable(NumLib::LocalToGlobalIndexMap const& dof_table,
+                             MeshLib::Mesh const& mesh) = 0;
+};
+
+}  // namespace NumLib
+
+#endif  // NUMLIB_CONVERGENCECRITERIONPERCOMPONENT_H

NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.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 "ConvergenceCriterionPerComponentDeltaX.h"
+#include <logog/include/logog.hpp>
+
+#include "BaseLib/ConfigTree.h"
+#include "NumLib/DOF/LocalToGlobalIndexMap.h"
+#include "NumLib/DOF/DOFTableUtil.h"
+
+namespace NumLib
+{
+ConvergenceCriterionPerComponentDeltaX::ConvergenceCriterionPerComponentDeltaX(
+    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)
+{
+    if (_abstols.size() != _reltols.size())
+        OGS_FATAL(
+            "The number of absolute and relative tolerances given must be the "
+            "same.");
+}
+
+void ConvergenceCriterionPerComponentDeltaX::checkDeltaX(
+    const GlobalVector& minus_delta_x, GlobalVector const& x)
+{
+    if ((!_dof_table) || (!_mesh))
+        OGS_FATAL("D.o.f. table or mesh have not been set.");
+
+    bool satisfied_abs = true;
+    bool satisfied_rel = true;
+
+    for (unsigned global_component = 0; global_component < _abstols.size();
+         ++global_component)
+    {
+        // TODO short cut if tol <= 0.0
+        auto error_dx = norm(minus_delta_x, global_component, _norm_type,
+                             *_dof_table, *_mesh);
+        auto norm_x =
+            norm(x, global_component, _norm_type, *_dof_table, *_mesh);
+
+        INFO(
+            "Convergence criterion, component %u: |dx|=%.4e, |x|=%.4e, "
+            "|dx|/|x|=%.4e",
+            error_dx, global_component, norm_x, error_dx / norm_x);
+
+        satisfied_abs = satisfied_abs && error_dx < _abstols[global_component];
+        satisfied_rel =
+            satisfied_rel && error_dx < _reltols[global_component] * norm_x;
+    }
+
+    _satisfied = _satisfied && (satisfied_abs || satisfied_rel);
+}
+
+void ConvergenceCriterionPerComponentDeltaX::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<ConvergenceCriterionPerComponentDeltaX>
+createConvergenceCriterionPerComponentDeltaX(const BaseLib::ConfigTree& config)
+{
+    config.checkConfigParameter("type", "PerComponentDeltaX");
+
+    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<ConvergenceCriterionPerComponentDeltaX>(
+        new ConvergenceCriterionPerComponentDeltaX(
+            std::move(*abstols), std::move(*reltols), norm_type));
+}
+
+}  // NumLib

NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.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_CONVERGENCECRITERIONPERCOMPONENTDELTAX_H
+#define NUMLIB_CONVERGENCECRITERIONPERCOMPONENTDELTAX_H
+
+#include "MathLib/LinAlg/LinAlgEnums.h"
+#include "ConvergenceCriterionPerComponent.h"
+
+namespace NumLib
+{
+class LocalToGlobalIndexMap;
+
+class ConvergenceCriterionPerComponentDeltaX
+    : public ConvergenceCriterionPerComponent
+{
+public:
+    explicit ConvergenceCriterionPerComponentDeltaX(
+        std::vector<double>&& absolute_tolerances,
+        std::vector<double>&& relative_tolerances,
+        MathLib::VecNormType norm_type);
+
+    bool hasDeltaXCheck() const override { return true; }
+    bool hasResidualCheck() const override { return false; }
+
+    void checkDeltaX(const GlobalVector& minus_delta_x,
+                     GlobalVector const& x) override;
+    void checkResidual(const GlobalVector& /*residual*/) override {}
+
+    void reset() override { _satisfied = true; }
+    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;
+};
+
+std::unique_ptr<ConvergenceCriterionPerComponentDeltaX>
+createConvergenceCriterionPerComponentDeltaX(BaseLib::ConfigTree const& config);
+
+}  // namespace NumLib
+
+#endif  // NUMLIB_CONVERGENCECRITERIONPERCOMPONENTDELTAX_H