NumLib::NonlinearSolver< NonlinearSolverTag::Newton > Class Reference

Detailed Description

Find a solution to a nonlinear equation using the Newton-Raphson method.

Definition at line 76 of file NonlinearSolver.h.

#include <NonlinearSolver.h>

Inheritance diagram for NumLib::NonlinearSolver< NonlinearSolverTag::Newton >:

Collaboration diagram for NumLib::NonlinearSolver< NonlinearSolverTag::Newton >:

Public Types
using	System = NonlinearSystem<NonlinearSolverTag::Newton>
	Type of the nonlinear equation system to be solved.

Public Member Functions
	NonlinearSolver (GlobalLinearSolver &linear_solver, int const maxiter, int const recompute_jacobian=1, double const damping=1.0, std::optional< double > const damping_reduction=std::nullopt)
	~NonlinearSolver ()
void	setEquationSystem (System &eq, ConvergenceCriterion &conv_crit)
void	calculateNonEquilibriumInitialResiduum (std::vector< GlobalVector * > const &x, std::vector< GlobalVector * > const &x_prev, int const process_id) override
NonlinearSolverStatus	solve (std::vector< GlobalVector * > &x, std::vector< GlobalVector * > const &x_prev, std::function< void(int, bool, std::vector< GlobalVector * > const &)> const &postIterationCallback, int const process_id) override
void	compensateNonEquilibriumInitialResiduum (bool const value)
Public Member Functions inherited from NumLib::NonlinearSolverBase
virtual	~NonlinearSolverBase ()=default

Private Attributes
GlobalLinearSolver &	_linear_solver
System *	_equation_system = nullptr
ConvergenceCriterion *	_convergence_criterion = nullptr
int const	_maxiter
	maximum number of iterations
int const	_recompute_jacobian = 1
double const	_damping
std::optional< double > const	_damping_reduction
GlobalVector *	_r_neq = nullptr
	non-equilibrium initial residuum.
std::size_t	_res_id = 0u
	ID of the residual vector.
std::size_t	_J_id = 0u
	ID of the Jacobian matrix.
std::size_t	_minus_delta_x_id = 0u
	ID of the $-\Delta x$ vector.
std::size_t	_x_new_id
	ID of the vector storing $x - (-\Delta x)$.
std::size_t	_r_neq_id = 0u
bool	_compensate_non_equilibrium_initial_residuum = false

Member Typedef Documentation

System

using NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::System = NonlinearSystem<NonlinearSolverTag::Newton>

Type of the nonlinear equation system to be solved.

Definition at line 81 of file NonlinearSolver.h.

Constructor & Destructor Documentation

NonlinearSolver()

NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::NonlinearSolver ( GlobalLinearSolver & linear_solver, int const maxiter, int const recompute_jacobian = 1, double const damping = 1.0, std::optional< double > const damping_reduction = std::nullopt )

inlineexplicit

Constructs a new instance.

Parameters

linear_solver	the linear solver used by this nonlinear solver.
maxiter	the maximum number of iterations used to solve the equation.
recompute_jacobian	recompute jacobian for the follow-up steps
damping	A positive damping factor.

See also

_damping

Parameters

damping_reduction

A parameter that reduces damping

See also

_damping_reduction

Definition at line 94 of file NonlinearSolver.h.

        : _linear_solver(linear_solver),
          _maxiter(maxiter),
          _recompute_jacobian(recompute_jacobian),
          _damping(damping),
          _damping_reduction(damping_reduction)
    {
    }

~NonlinearSolver()

NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::~NonlinearSolver

(

)

Definition at line 613 of file NonlinearSolver.cpp.

{
    if (_r_neq != nullptr)
    {
        NumLib::GlobalVectorProvider::provider.releaseVector(*_r_neq);
    }
}

References NumLib::GlobalVectorProvider::provider, and NumLib::VectorProvider::releaseVector().

Member Function Documentation

calculateNonEquilibriumInitialResiduum()

void NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::calculateNonEquilibriumInitialResiduum ( std::vector< GlobalVector * > const & x, std::vector< GlobalVector * > const & x_prev, int const process_id )

overridevirtual

Implements NumLib::NonlinearSolverBase.

Definition at line 304 of file NonlinearSolver.cpp.

{
    if (!_compensate_non_equilibrium_initial_residuum)
    {
        return;
    }
 
    INFO("Calculate non-equilibrium initial residuum.");
 
    _equation_system->assemble(x, x_prev, process_id);
    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector(_r_neq_id);
    _equation_system->getResidual(*x[process_id], *x_prev[process_id], *_r_neq);
 
    // Set the values of the selected entries of _r_neq, which are associated
    // with the equations that do not need initial residual compensation, to
    // zero.
    const auto selected_global_indices =
        _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
    std::vector<double> zero_entries(selected_global_indices.size(), 0.0);
 
    _r_neq->set(selected_global_indices, zero_entries);
 
    MathLib::LinAlg::finalizeAssembly(*_r_neq);
}

References MathLib::LinAlg::finalizeAssembly(), NumLib::VectorProvider::getVector(), INFO(), and NumLib::GlobalVectorProvider::provider.

compensateNonEquilibriumInitialResiduum()

void NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::compensateNonEquilibriumInitialResiduum ( bool const value )

inline

Definition at line 130 of file NonlinearSolver.h.

    {
        _compensate_non_equilibrium_initial_residuum = value;
    }

setEquationSystem()

void NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::setEquationSystem ( System & eq, ConvergenceCriterion & conv_crit )

inline

Set the nonlinear equation system that will be solved. TODO doc

Definition at line 112 of file NonlinearSolver.h.

    {
        _equation_system = &eq;
        _convergence_criterion = &conv_crit;
    }

solve()

NonlinearSolverStatus NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::solve ( std::vector< GlobalVector * > & x, std::vector< GlobalVector * > const & x_prev, std::function< void(int, bool, std::vector< GlobalVector * > const &)> const & postIterationCallback, int const process_id )

overridevirtual

Assemble and solve the equation system.

Parameters

x	in: the initial guess, out: the solution.
x_prev	previous time step solution.
postIterationCallback	called after each iteration if set.
process_id	usually used in staggered schemes.

Return values

true	if the equation system could be solved
false	otherwise

Implements NumLib::NonlinearSolverBase.

Definition at line 332 of file NonlinearSolver.cpp.

{
    namespace LinAlg = MathLib::LinAlg;
    auto& sys = *_equation_system;
 
    auto& res = NumLib::GlobalVectorProvider::provider.getVector(_res_id);
    auto& minus_delta_x =
        NumLib::GlobalVectorProvider::provider.getVector(_minus_delta_x_id);
    auto& J = NumLib::GlobalMatrixProvider::provider.getMatrix(_J_id);
 
    bool error_norms_met = false;
 
    // TODO be more efficient
    // init minus_delta_x to the right size
    LinAlg::copy(*x[process_id], minus_delta_x);
 
    _convergence_criterion->preFirstIteration();
 
    int iteration = 1;
#if !defined(USE_PETSC) && !defined(USE_LIS)
    int next_iteration_inv_jacobian_recompute = 1;
#endif
    for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
    {
        BaseLib::RunTime timer_dirichlet;
        double time_dirichlet = 0.0;
 
        BaseLib::RunTime time_iteration;
        INFO("Iteration #{:d} started.", iteration);
        time_iteration.start();
 
        timer_dirichlet.start();
        sys.computeKnownSolutions(*x[process_id], process_id);
        time_dirichlet += timer_dirichlet.elapsed();
 
        sys.preIteration(iteration, *x[process_id]);
 
        BaseLib::RunTime time_assembly;
        time_assembly.start();
        bool mpi_rank_assembly_ok = true;
        try
        {
            sys.assemble(x, x_prev, process_id);
        }
        catch (AssemblyException const& e)
        {
            ERR("Abort nonlinear iteration. Repeating timestep. Reason: {:s}",
                e.what());
            error_norms_met = false;
            iteration = _maxiter;
            mpi_rank_assembly_ok = false;
        }
        if (BaseLib::MPI::anyOf(!mpi_rank_assembly_ok))
        {
            break;
        }
        sys.getResidual(*x[process_id], *x_prev[process_id], res);
        sys.getJacobian(J);
        INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
 
        // Subtract non-equilibrium initial residuum if set
        if (_r_neq != nullptr)
            LinAlg::axpy(res, -1, *_r_neq);
 
        minus_delta_x.setZero();
 
        timer_dirichlet.start();
        sys.applyKnownSolutionsNewton(J, res, *x[process_id], minus_delta_x);
        time_dirichlet += timer_dirichlet.elapsed();
        INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
 
        if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
        {
            _convergence_criterion->checkResidual(res);
        }
 
        BaseLib::RunTime time_linear_solver;
        time_linear_solver.start();
#if !defined(USE_PETSC) && !defined(USE_LIS)
        auto linear_solver_behaviour = MathLib::LinearSolverBehaviour::REUSE;
        if (iteration == next_iteration_inv_jacobian_recompute)
        {
            linear_solver_behaviour =
                MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE;
            next_iteration_inv_jacobian_recompute =
                next_iteration_inv_jacobian_recompute + _recompute_jacobian;
        }
        if (!_linear_solver.compute(J, linear_solver_behaviour))
        {
            ERR("Newton: The linear solver failed in the compute() step.");
        }
        bool iteration_succeeded = _linear_solver.solve(res, minus_delta_x);
#else
        bool iteration_succeeded = _linear_solver.solve(J, res, minus_delta_x);
#endif
        INFO("[time] Linear solver took {:g} s.", time_linear_solver.elapsed());
 
        if (!iteration_succeeded)
        {
            ERR("Newton: The linear solver failed.");
        }
        else
        {
            // TODO could be solved in a better way
            // cf.
            // https://petsc.org/release/manualpages/Vec/VecWAXPY
 
            // Copy pointers, replace the one for the given process id.
            std::vector<GlobalVector*> x_new{x};
            x_new[process_id] =
                &NumLib::GlobalVectorProvider::provider.getVector(
                    *x[process_id], _x_new_id);
            double damping = _damping;
            if (_convergence_criterion->hasNonNegativeDamping())
            {
                damping = _convergence_criterion->getDampingFactor(
                    minus_delta_x, *x[process_id], _damping);
            }
            if (_damping_reduction)
            {
                damping =
                    damping +
                    (1 - damping) *
                        std::clamp(iteration / *_damping_reduction, 0.0, 1.0);
            }
            LinAlg::axpy(*x_new[process_id], -damping, minus_delta_x);
 
            if (postIterationCallback)
            {
                postIterationCallback(iteration, error_norms_met, x_new);
            }
 
            switch (sys.postIteration(*x_new[process_id]))
            {
                case IterationResult::SUCCESS:
                    break;
                case IterationResult::FAILURE:
                    ERR("Newton: The postIteration() hook reported a "
                        "non-recoverable error.");
                    iteration_succeeded = false;
                    break;
                case IterationResult::REPEAT_ITERATION:
                    INFO(
                        "Newton: The postIteration() hook decided that this "
                        "iteration has to be repeated.");
                    // TODO introduce some onDestroy hook.
                    NumLib::GlobalVectorProvider::provider.releaseVector(
                        *x_new[process_id]);
                    continue;  // That throws the iteration result away.
            }
 
            LinAlg::copy(*x_new[process_id],
                         *x[process_id]);  // copy new solution to x
            NumLib::GlobalVectorProvider::provider.releaseVector(
                *x_new[process_id]);
        }
 
        if (!iteration_succeeded)
        {
            // Don't compute further error norms, but break here.
            error_norms_met = false;
            break;
        }
 
        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[process_id]);
            }
 
            error_norms_met = _convergence_criterion->isSatisfied();
        }
 
        INFO("[time] Iteration #{:d} took {:g} s.", iteration,
             time_iteration.elapsed());
 
        if (error_norms_met)
        {
            break;
        }
 
        // Avoid increment of the 'iteration' if the error norms are not met,
        // but maximum number of iterations is reached.
        if (iteration >= _maxiter)
        {
            break;
        }
    }
 
    if (iteration > _maxiter)
    {
        ERR("Newton: Could not solve the given nonlinear system within {:d} "
            "iterations",
            _maxiter);
    }
 
    NumLib::GlobalMatrixProvider::provider.releaseMatrix(J);
    NumLib::GlobalVectorProvider::provider.releaseVector(res);
    NumLib::GlobalVectorProvider::provider.releaseVector(minus_delta_x);
 
    return {error_norms_met, iteration};
}

References BaseLib::MPI::anyOf(), MathLib::LinAlg::axpy(), MathLib::LinAlg::copy(), BaseLib::RunTime::elapsed(), ERR(), NumLib::FAILURE, NumLib::MatrixProvider::getMatrix(), NumLib::VectorProvider::getVector(), INFO(), NumLib::GlobalMatrixProvider::provider, NumLib::GlobalVectorProvider::provider, MathLib::RECOMPUTE_AND_STORE, NumLib::MatrixProvider::releaseMatrix(), NumLib::VectorProvider::releaseVector(), NumLib::REPEAT_ITERATION, MathLib::REUSE, BaseLib::RunTime::start(), and NumLib::SUCCESS.

Member Data Documentation

_compensate_non_equilibrium_initial_residuum

bool NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_compensate_non_equilibrium_initial_residuum = false

private

Enables computation of the non-equilibrium initial residuum $r_{\rm neq}$ before the first time step. The forces are zero if the external forces are in equilibrium with the initial state/initial conditions. During the simulation the new residuum reads $\tilde r = r - r_{\rm neq}$.

Definition at line 170 of file NonlinearSolver.h.

_convergence_criterion

ConvergenceCriterion* NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_convergence_criterion = nullptr

private

Definition at line 140 of file NonlinearSolver.h.

_damping

double const NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_damping

private

A positive damping factor. The default value 1.0 gives a non-damped Newton method. Common values are in the range 0.5 to 0.7 for somewhat conservative method and seldom become smaller than 0.2 for very conservative approach.

Definition at line 149 of file NonlinearSolver.h.

_damping_reduction

std::optional<double> const NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_damping_reduction

private

A parameter that reduces the damping by iteration / _damping_reduction If the expression is smaller than zero, then damping is the same for all iterations. If the expression is bigger than one it is clamped by one.

Definition at line 154 of file NonlinearSolver.h.

_equation_system

System* NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_equation_system = nullptr

private

Definition at line 137 of file NonlinearSolver.h.

_J_id

std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_J_id = 0u

private

ID of the Jacobian matrix.

Definition at line 158 of file NonlinearSolver.h.

_linear_solver

GlobalLinearSolver& NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_linear_solver

private

Definition at line 136 of file NonlinearSolver.h.

_maxiter

int const NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_maxiter

private

maximum number of iterations

Definition at line 141 of file NonlinearSolver.h.

_minus_delta_x_id

std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_minus_delta_x_id = 0u

private

ID of the $-\Delta x$ vector.

Definition at line 159 of file NonlinearSolver.h.

_r_neq

GlobalVector* NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_r_neq = nullptr

private

non-equilibrium initial residuum.

Definition at line 156 of file NonlinearSolver.h.

_r_neq_id

std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_r_neq_id = 0u

private

ID of the non-equilibrium initial residuum vector.

Definition at line 162 of file NonlinearSolver.h.

_recompute_jacobian

int const NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_recompute_jacobian = 1

private

Definition at line 143 of file NonlinearSolver.h.

_res_id

std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_res_id = 0u

private

ID of the residual vector.

Definition at line 157 of file NonlinearSolver.h.

_x_new_id

std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::_x_new_id

private

Initial value:

=
        0u

ID of the vector storing $x - (-\Delta x)$.

Definition at line 160 of file NonlinearSolver.h.

---

The documentation for this class was generated from the following files:

• NumLib/ODESolver/NonlinearSolver.h
• NumLib/ODESolver/NonlinearSolver.cpp