From dd30c1c033e4513b33f895aa98f20c7d6fc235e0 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Mon, 22 Aug 2016 16:04:31 +0200
Subject: [PATCH] [NL] added per component residual check
---
NumLib/ODESolver/ConvergenceCriterion.cpp | 3 +
...nvergenceCriterionPerComponentResidual.cpp | 121 ++++++++++++++++++
...ConvergenceCriterionPerComponentResidual.h | 61 +++++++++
3 files changed, 185 insertions(+)
create mode 100644 NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
create mode 100644 NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h
diff --git a/NumLib/ODESolver/ConvergenceCriterion.cpp b/NumLib/ODESolver/ConvergenceCriterion.cpp
index 7e9fba73c6d..e89a854f690 100644
--- a/NumLib/ODESolver/ConvergenceCriterion.cpp
+++ b/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());
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
new file mode 100644
index 00000000000..dc94fb9b6e9
--- /dev/null
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.cpp
@@ -0,0 +1,121 @@
+/**
+ * \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
diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h
new file mode 100644
index 00000000000..d04e4f6e06f
--- /dev/null
+++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentResidual.h
@@ -0,0 +1,61 @@
+/**
+ * \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
--
GitLab