diff --git a/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp new file mode 100644 index 0000000000000000000000000000000000000000..5ab4ab6edbd8a525e97f9bd33f6f68c7da97f490 --- /dev/null +++ b/NumLib/ODESolver/ConvergenceCriterionDeltaX.cpp @@ -0,0 +1,73 @@ +/** + * \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 diff --git a/NumLib/ODESolver/ConvergenceCriterionDeltaX.h b/NumLib/ODESolver/ConvergenceCriterionDeltaX.h new file mode 100644 index 0000000000000000000000000000000000000000..6dae228950b2d423417110beaf9a01702e3359b8 --- /dev/null +++ b/NumLib/ODESolver/ConvergenceCriterionDeltaX.h @@ -0,0 +1,48 @@ +/** + * \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 diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponent.h b/NumLib/ODESolver/ConvergenceCriterionPerComponent.h new file mode 100644 index 0000000000000000000000000000000000000000..e1ff2230f208985eae9c47a11f7eb1049cd23457 --- /dev/null +++ b/NumLib/ODESolver/ConvergenceCriterionPerComponent.h @@ -0,0 +1,33 @@ +/** + * \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 diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp new file mode 100644 index 0000000000000000000000000000000000000000..8c9c9e925abb92b9c97a89421e6ad2f3997aee51 --- /dev/null +++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.cpp @@ -0,0 +1,108 @@ +/** + * \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 diff --git a/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.h b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.h new file mode 100644 index 0000000000000000000000000000000000000000..51601ee5817b01f28e1452417fd900e67329b17e --- /dev/null +++ b/NumLib/ODESolver/ConvergenceCriterionPerComponentDeltaX.h @@ -0,0 +1,56 @@ +/** + * \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