From c7399b02740fad7ce8b575aa07184444b54672e3 Mon Sep 17 00:00:00 2001
From: Christoph Lehmann <christoph.lehmann@ufz.de>
Date: Fri, 19 Aug 2016 10:51:43 +0200
Subject: [PATCH] [NL] nonlinear solver uses conv criteria
---
NumLib/ODESolver/NonlinearSolver.cpp | 100 +++++++++++++--------------
NumLib/ODESolver/NonlinearSolver.h | 33 ++++++---
NumLib/ODESolver/TimeLoopSingleODE.h | 16 +++--
3 files changed, 84 insertions(+), 65 deletions(-)
diff --git a/NumLib/ODESolver/NonlinearSolver.cpp b/NumLib/ODESolver/NonlinearSolver.cpp
index c8569a2881a..0ab6d7e5bf1 100644
--- a/NumLib/ODESolver/NonlinearSolver.cpp
+++ b/NumLib/ODESolver/NonlinearSolver.cpp
@@ -19,6 +19,7 @@
#include "BaseLib/RunTime.h"
#include "MathLib/LinAlg/LinAlg.h"
#include "NumLib/DOF/GlobalMatrixProviders.h"
+#include "ConvergenceCriterion.h"
namespace NumLib
{
@@ -32,7 +33,7 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
GlobalVector& x,
std::function<void(unsigned, GlobalVector const&)> const& postIterationCallback)
{
- namespace LinAlg= MathLib::LinAlg;
+ namespace LinAlg = MathLib::LinAlg;
auto& sys = *_equation_system;
auto& A =
@@ -51,6 +52,8 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
unsigned iteration = 1;
for (; iteration <= _maxiter; ++iteration)
{
+ _convergence_criterion->reset();
+
BaseLib::RunTime time_iteration;
time_iteration.start();
@@ -69,6 +72,13 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
sys.applyKnownSolutionsPicard(A, rhs, x_new);
INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
+ if (_convergence_criterion->hasResidualCheck()) {
+ GlobalVector res;
+ LinAlg::matMult(A, x_new, res); // res = A * x_new
+ LinAlg::axpy(res, -1.0, rhs); // res -= rhs
+ _convergence_criterion->checkResidual(res);
+ }
+
BaseLib::RunTime time_linear_solver;
time_linear_solver.start();
bool iteration_succeeded = _linear_solver.solve(A, rhs, x_new);
@@ -114,14 +124,18 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
break;
}
- auto const norm_x = LinAlg::norm2(x);
- // x is used as delta_x in order to compute the error.
- LinAlg::aypx(x, -1.0, x_new); // x = _x_new - x
- auto const error_dx = LinAlg::norm2(x);
- INFO(
- "Picard: Iteration #%u |dx|=%.4e, |x|=%.4e, |dx|/|x|=%.4e,"
- " tolerance(dx)=%.4e",
- iteration, error_dx, norm_x, error_dx / norm_x, _tol);
+ if (sys.isLinear()) {
+ error_norms_met = true;
+ } else {
+ if (_convergence_criterion->hasDeltaXCheck()) {
+ GlobalVector minus_delta_x(x);
+ LinAlg::axpy(minus_delta_x, -1.0,
+ x_new); // minus_delta_x = x - x_new
+ _convergence_criterion->checkDeltaX(minus_delta_x, x_new);
+ }
+
+ error_norms_met = _convergence_criterion->isSatisfied();
+ }
// Update x s.t. in the next iteration we will compute the right delta x
LinAlg::copy(x_new, x);
@@ -129,17 +143,8 @@ bool NonlinearSolver<NonlinearSolverTag::Picard>::solve(
INFO("[time] Iteration #%u took %g s.", iteration,
time_iteration.elapsed());
- if (error_dx < _tol)
- {
- error_norms_met = true;
- break;
- }
-
- if (sys.isLinear())
- {
- error_norms_met = true;
+ if (error_norms_met)
break;
- }
}
if (iteration > _maxiter)
@@ -190,6 +195,8 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
unsigned iteration = 1;
for (; iteration <= _maxiter; ++iteration)
{
+ _convergence_criterion->reset();
+
BaseLib::RunTime time_iteration;
time_iteration.start();
@@ -208,7 +215,8 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
sys.applyKnownSolutionsNewton(J, res, minus_delta_x);
INFO("[time] Applying Dirichlet BCs took %g s.", time_dirichlet.elapsed());
- auto const error_res = LinAlg::norm2(res);
+ if (_convergence_criterion->hasResidualCheck())
+ _convergence_criterion->checkResidual(res);
BaseLib::RunTime time_linear_solver;
time_linear_solver.start();
@@ -266,28 +274,22 @@ bool NonlinearSolver<NonlinearSolverTag::Newton>::solve(
break;
}
- auto const error_dx = LinAlg::norm2(minus_delta_x);
- auto const norm_x = LinAlg::norm2(x);
- INFO(
- "Newton: Iteration #%u |dx|=%.4e, |r|=%.4e, |x|=%.4e, "
- "|dx|/|x|=%.4e,"
- " tolerance(dx)=%.4e",
- iteration, error_dx, error_res, norm_x, error_dx / norm_x, _tol);
+ if (sys.isLinear()) {
+ error_norms_met = true;
+ } else {
+ if (_convergence_criterion->hasDeltaXCheck()) {
+ // Note: x contains the new solution!
+ _convergence_criterion->checkDeltaX(minus_delta_x, x);
+ }
+
+ error_norms_met = _convergence_criterion->isSatisfied();
+ }
INFO("[time] Iteration #%u took %g s.", iteration,
time_iteration.elapsed());
- if (error_dx < _tol)
- {
- error_norms_met = true;
- break;
- }
-
- if (sys.isLinear())
- {
- error_norms_met = true;
+ if (error_norms_met)
break;
- }
}
if (iteration > _maxiter)
@@ -313,27 +315,23 @@ createNonlinearSolver(GlobalLinearSolver& linear_solver,
//! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__type}
auto const type = config.getConfigParameter<std::string>("type");
- //! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__tol}
- auto const tol = config.getConfigParameter<double>("tol");
//! \ogs_file_param{prj__nonlinear_solvers__nonlinear_solver__max_iter}
auto const max_iter = config.getConfigParameter<unsigned>("max_iter");
- //! \ogs_file_param_special{prj__nonlinear_solvers__nonlinear_solver__Picard}
- if (type == "Picard")
- {
+ if (type == "Picard") {
auto const tag = NonlinearSolverTag::Picard;
using ConcreteNLS = NonlinearSolver<tag>;
- return std::make_pair(std::unique_ptr<AbstractNLS>(new ConcreteNLS{
- linear_solver, tol, max_iter}),
- tag);
- }
- else if (type == "Newton")
- {
+ return std::make_pair(
+ std::unique_ptr<AbstractNLS>(
+ new ConcreteNLS{linear_solver, max_iter}),
+ tag);
+ } else if (type == "Newton") {
auto const tag = NonlinearSolverTag::Newton;
using ConcreteNLS = NonlinearSolver<tag>;
- return std::make_pair(std::unique_ptr<AbstractNLS>(new ConcreteNLS{
- linear_solver, tol, max_iter}),
- tag);
+ return std::make_pair(
+ std::unique_ptr<AbstractNLS>(
+ new ConcreteNLS{linear_solver, max_iter}),
+ tag);
}
OGS_FATAL("Unsupported nonlinear solver type");
}
diff --git a/NumLib/ODESolver/NonlinearSolver.h b/NumLib/ODESolver/NonlinearSolver.h
index 4a7de970df3..1a16fc98fdf 100644
--- a/NumLib/ODESolver/NonlinearSolver.h
+++ b/NumLib/ODESolver/NonlinearSolver.h
@@ -14,6 +14,7 @@
#include <utility>
#include <logog/include/logog.hpp>
+#include "ConvergenceCriterion.h"
#include "NonlinearSystem.h"
#include "Types.h"
@@ -88,14 +89,21 @@ public:
* \param maxiter the maximum number of iterations used to solve the
* equation.
*/
- explicit NonlinearSolver(GlobalLinearSolver& linear_solver, double const tol,
- const unsigned maxiter)
- : _linear_solver(linear_solver), _tol(tol), _maxiter(maxiter)
+ explicit NonlinearSolver(
+ GlobalLinearSolver& linear_solver,
+ const unsigned maxiter)
+ : _linear_solver(linear_solver),
+ _maxiter(maxiter)
{
}
//! Set the nonlinear equation system that will be solved.
- void setEquationSystem(System& eq) { _equation_system = &eq; }
+ //! TODO doc
+ void setEquationSystem(System& eq, ConvergenceCriterion& conv_crit)
+ {
+ _equation_system = &eq;
+ _convergence_criterion = &conv_crit;
+ }
void assemble(GlobalVector const& x) const override;
bool solve(GlobalVector& x,
@@ -106,7 +114,8 @@ private:
GlobalLinearSolver& _linear_solver;
System* _equation_system = nullptr;
- const double _tol; //!< tolerance of the solver
+ // TODO doc
+ ConvergenceCriterion* _convergence_criterion = nullptr;
const unsigned _maxiter; //!< maximum number of iterations
double const _alpha =
@@ -139,14 +148,19 @@ public:
* \param maxiter the maximum number of iterations used to solve the
* equation.
*/
- explicit NonlinearSolver(GlobalLinearSolver& linear_solver, double const tol,
+ explicit NonlinearSolver(GlobalLinearSolver& linear_solver,
const unsigned maxiter)
- : _linear_solver(linear_solver), _tol(tol), _maxiter(maxiter)
+ : _linear_solver(linear_solver), _maxiter(maxiter)
{
}
//! Set the nonlinear equation system that will be solved.
- void setEquationSystem(System& eq) { _equation_system = &eq; }
+ //! TODO doc
+ void setEquationSystem(System& eq, ConvergenceCriterion& conv_crit)
+ {
+ _equation_system = &eq;
+ _convergence_criterion = &conv_crit;
+ }
void assemble(GlobalVector const& x) const override;
bool solve(GlobalVector& x,
@@ -157,7 +171,8 @@ private:
GlobalLinearSolver& _linear_solver;
System* _equation_system = nullptr;
- const double _tol; //!< tolerance of the solver
+ // TODO doc
+ ConvergenceCriterion* _convergence_criterion = nullptr;
const unsigned _maxiter; //!< maximum number of iterations
std::size_t _A_id = 0u; //!< ID of the \f$ A \f$ matrix.
diff --git a/NumLib/ODESolver/TimeLoopSingleODE.h b/NumLib/ODESolver/TimeLoopSingleODE.h
index 228a5153786..29df9d6e26f 100644
--- a/NumLib/ODESolver/TimeLoopSingleODE.h
+++ b/NumLib/ODESolver/TimeLoopSingleODE.h
@@ -36,13 +36,18 @@ public:
* \param linear_solver the linear solver used to solve the linearized ODE
* system.
* \param nonlinear_solver The solver to be used to resolve nonlinearities.
+ * \param convergence_criterion The convergence criterion used by the
+ * nonlinear solver.
*/
- TimeLoopSingleODE(TDiscODESys& ode_sys,
- std::unique_ptr<GlobalLinearSolver>&& linear_solver,
- std::unique_ptr<NLSolver>&& nonlinear_solver)
+ TimeLoopSingleODE(
+ TDiscODESys& ode_sys,
+ std::unique_ptr<GlobalLinearSolver>&& linear_solver,
+ std::unique_ptr<NLSolver>&& nonlinear_solver,
+ std::unique_ptr<ConvergenceCriterion>&& convergence_criterion)
: _ode_sys(ode_sys),
_linear_solver(std::move(linear_solver)),
- _nonlinear_solver(std::move(nonlinear_solver))
+ _nonlinear_solver(std::move(nonlinear_solver)),
+ _convergence_criterion(std::move(convergence_criterion))
{
}
@@ -69,6 +74,7 @@ private:
TDiscODESys& _ode_sys;
std::unique_ptr<GlobalLinearSolver> _linear_solver;
std::unique_ptr<NLSolver> _nonlinear_solver;
+ std::unique_ptr<ConvergenceCriterion> _convergence_criterion;
};
//! @}
@@ -87,7 +93,7 @@ bool TimeLoopSingleODE<NLTag>::loop(const double t0, GlobalVector const& x0,
time_disc.setInitialState(t0, x0); // push IC
- _nonlinear_solver->setEquationSystem(_ode_sys);
+ _nonlinear_solver->setEquationSystem(_ode_sys, *_convergence_criterion);
if (time_disc.needsPreload())
{
--
GitLab