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)
    {
    }

References _damping, _damping_reduction, _linear_solver, _maxiter, and _recompute_jacobian.

Referenced by ~NonlinearSolver(), calculateNonEquilibriumInitialResiduum(), and solve().

◆ ~NonlinearSolver()

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

Definition at line 637 of file NonlinearSolver.cpp.

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

References NonlinearSolver(), _r_neq, and NumLib::GlobalVectorProvider::provider.

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 322 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);
    _equation_system->setReleaseNodalForces(_r_neq, process_id);
 
    MathLib::LinAlg::finalizeAssembly(*_r_neq);
}

References NonlinearSolver(), _compensate_non_equilibrium_initial_residuum, _equation_system, _r_neq, _r_neq_id, MathLib::LinAlg::finalizeAssembly(), 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;
    }

References _compensate_non_equilibrium_initial_residuum.

◆ 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;
    }

References _convergence_criterion, and _equation_system.

◆ 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 351 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;
        }
 
        bool iteration_succeeded = false;
        if (!_linear_solver.compute(J, linear_solver_behaviour))
        {
            ERR("Newton: The linear solver failed in the compute() step.");
        }
        else
        {
            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 NonlinearSolver(), _convergence_criterion, _damping, _damping_reduction, _equation_system, _J_id, _linear_solver, _maxiter, _minus_delta_x_id, _r_neq, _recompute_jacobian, _res_id, _x_new_id, BaseLib::MPI::anyOf(), MathLib::LinAlg::axpy(), MathLib::LinAlg::copy(), BaseLib::RunTime::elapsed(), ERR(), NumLib::FAILURE, INFO(), NumLib::GlobalMatrixProvider::provider, NumLib::GlobalVectorProvider::provider, MathLib::RECOMPUTE_AND_STORE, 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.

Referenced by calculateNonEquilibriumInitialResiduum(), and compensateNonEquilibriumInitialResiduum().

◆ _convergence_criterion

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

private

Definition at line 140 of file NonlinearSolver.h.

Referenced by setEquationSystem(), and solve().

◆ _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.

Referenced by NonlinearSolver(), and solve().

◆ _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.

Referenced by NonlinearSolver(), and solve().