Detailed Description

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

Definition at line 70 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, std::unique_ptr< NewtonStepStrategy > newton_strategy, int const recompute_jacobian=1)
	~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
int const	_maxiter
	maximum number of iterations
std::unique_ptr< NewtonStepStrategy >	_step_strategy
	Globalization / step-acceptance strategy (e.g. fixed damping).
ConvergenceCriterion *	_convergence_criterion = nullptr
	Convergence criterion used to terminate the Newton iteration.
int const	_recompute_jacobian
	Recompute Jacobian every this many steps.
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 75 of file NonlinearSolver.h.

Constructor & Destructor Documentation

◆ NonlinearSolver()

NumLib::NonlinearSolver< NonlinearSolverTag::Newton >::NonlinearSolver	(	GlobalLinearSolver &	linear_solver,
		int const	maxiter,
		std::unique_ptr< NewtonStepStrategy >	newton_strategy,
		int const	recompute_jacobian = 1 )

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.
newton_strategy	globalization / step-acceptance strategy (e.g. fixed damping, line search). Ownership is transferred to the solver.
recompute_jacobian	recompute the Jacobian every this many steps.

Definition at line 87 of file NonlinearSolver.h.

        : _linear_solver(linear_solver),
          _maxiter(maxiter),
          _step_strategy(std::move(newton_strategy)),
          _recompute_jacobian(recompute_jacobian)
    {
    }

References _linear_solver, _maxiter, _recompute_jacobian, and _step_strategy.

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

◆ ~NonlinearSolver()

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

Definition at line 605 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 327 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.
    auto selected_global_indices =
        _equation_system->getIndicesOfResiduumWithoutInitialCompensation();
 
#ifdef USE_PETSC
    // Ghost entry with global index 0 is encoded as -global_size
    // After abs(), it appears as global_size and must be converted back to 0
    auto const global_size = _r_neq->size();
    for (auto& idx : selected_global_indices)
    {
        if (idx == global_size)
        {
            idx = 0;
        }
    }
#endif
 
    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 125 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 and the convergence criterion used to terminate the iteration. Also forwards the criterion to the step strategy so that strategies that depend on residual or delta-x information (e.g. non-negative damping) can use it.

Definition at line 105 of file NonlinearSolver.h.

    {
        _equation_system = &eq;
        _convergence_criterion = &conv_crit;
        _step_strategy->setDampingPolicy(
            _convergence_criterion->dampingPolicy());
    }

References _convergence_criterion, _equation_system, and _step_strategy.

◆ 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 369 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();
 
    NewtonStepContext step_ctx{sys, x_prev, process_id};
 
    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);
            auto const step_result = _step_strategy->applyStep(
                *x[process_id], minus_delta_x, res, J, *x_new[process_id],
                step_ctx, iteration);
 
            if (step_result.step_length != 1.0)
            {
                INFO("Step length: {:g}", step_result.step_length);
            }
 
            if (!step_result.success)
            {
                ERR("Newton: step strategy failed.");
                iteration_succeeded = false;
            }
            else if (!step_result.x_new_is_set)
            {
                LinAlg::axpy(*x_new[process_id], -1.0, 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, _equation_system, _J_id, _linear_solver, _maxiter, _minus_delta_x_id, _r_neq, _recompute_jacobian, _res_id, _step_strategy, _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 159 of file NonlinearSolver.h.

Referenced by calculateNonEquilibriumInitialResiduum(), and compensateNonEquilibriumInitialResiduum().