diff --git a/NumLib/ODESolver/ConvergenceCriterion.cpp b/NumLib/ODESolver/ConvergenceCriterion.cpp index 74b5cf521ff17c203914611bebf2657c316a4aaa..7e9fba73c6da33f5464360ee76b6a6b119d557cb 100644 --- a/NumLib/ODESolver/ConvergenceCriterion.cpp +++ b/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); } diff --git a/NumLib/ODESolver/ConvergenceCriterionResidual.cpp b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp new file mode 100644 index 0000000000000000000000000000000000000000..ff0d719c078b8aff1554489b287bfa7991f67fe9 --- /dev/null +++ b/NumLib/ODESolver/ConvergenceCriterionResidual.cpp @@ -0,0 +1,76 @@ +/** + * \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 diff --git a/NumLib/ODESolver/ConvergenceCriterionResidual.h b/NumLib/ODESolver/ConvergenceCriterionResidual.h new file mode 100644 index 0000000000000000000000000000000000000000..105fd7bb1a031840faa82cead084fa322b5b5c5d --- /dev/null +++ b/NumLib/ODESolver/ConvergenceCriterionResidual.h @@ -0,0 +1,51 @@ +/** + * \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