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)

	~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

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 )

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

Definition at line 92 of file NonlinearSolver.h.

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

◆ ~NonlinearSolver()

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

Definition at line 600 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 302 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 125 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 107 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 330 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;
        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);
            }
            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 160 of file NonlinearSolver.h.

◆ _convergence_criterion

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

private

Definition at line 135 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 144 of file NonlinearSolver.h.

◆ _equation_system

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

private

Definition at line 132 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 148 of file NonlinearSolver.h.

◆ _linear_solver

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

private

Definition at line 131 of file NonlinearSolver.h.

◆ _maxiter

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

private

maximum number of iterations

Definition at line 136 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 149 of file NonlinearSolver.h.

◆ _r_neq

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

private

non-equilibrium initial residuum.

Definition at line 146 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 152 of file NonlinearSolver.h.

◆ _recompute_jacobian

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

private

Definition at line 138 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 147 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 150 of file NonlinearSolver.h.

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

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

Detailed Description

Public Types

Public Member Functions

Private Attributes

Member Typedef Documentation

◆ System

Constructor & Destructor Documentation

◆ NonlinearSolver()

◆ ~NonlinearSolver()

Member Function Documentation

◆ calculateNonEquilibriumInitialResiduum()

◆ compensateNonEquilibriumInitialResiduum()

◆ setEquationSystem()

◆ solve()

Member Data Documentation

◆ _compensate_non_equilibrium_initial_residuum

◆ _convergence_criterion

◆ _damping

◆ _equation_system

◆ _J_id

◆ _linear_solver

◆ _maxiter

◆ _minus_delta_x_id

◆ _r_neq

◆ _r_neq_id

◆ _recompute_jacobian

◆ _res_id

◆ _x_new_id