Detailed Description

Find a solution to a nonlinear equation using the Picard fixpoint iteration method.

Notation used throughout this class and its documentation: \( x_k \) is the current iterate and \( g \) is the Picard fixpoint map, i.e. \( g(x_k) \) is the solution of the linearized system \( A(x_k)\,x = b(x_k) \) (the result of one linear solve). A converged iteration is a fixpoint \( x = g(x) \), and one iteration reads \( x_{k+1} = g(x_k) \).

With an under-relaxation (damping) coefficient \( \beta \in (0, 1] \) the iteration becomes \( x_{k+1} = (1-\beta)\,x_k + \beta\,g(x_k) \); the value \( \beta = 1 \) recovers the plain Picard update above.

Definition at line 186 of file NonlinearSolver.h.

#include <NonlinearSolver.h>

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

[legend]

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

[legend]

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

Public Member Functions
	NonlinearSolver (GlobalLinearSolver &linear_solver, const int maxiter, const double damping)
	~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
double const	_damping
const int	_maxiter
	maximum number of iterations
GlobalVector *	_r_neq = nullptr
	non-equilibrium initial residuum.
std::size_t	_A_id = 0u
	ID of the \( A \) matrix.
std::size_t	_rhs_id = 0u
	ID of the right-hand side vector.
std::size_t	_x_new_id = 0u
std::size_t	_r_neq_id = 0u
bool	_compensate_non_equilibrium_initial_residuum = false

Member Typedef Documentation

◆ System

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

Type of the nonlinear equation system to be solved.

Definition at line 191 of file NonlinearSolver.h.

Constructor & Destructor Documentation

◆ NonlinearSolver()

NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::NonlinearSolver	(	GlobalLinearSolver &	linear_solver,
		const int	maxiter,
		const double	damping )

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.
damping	under-relaxation coefficient \( \beta \in (0, 1] \) applied to the Picard update (see class description); \( \beta = 1 \) disables damping.

Definition at line 202 of file NonlinearSolver.h.

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

References _damping, _linear_solver, and _maxiter.

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

◆ ~NonlinearSolver()

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

Definition at line 632 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::Picard >::calculateNonEquilibriumInitialResiduum	(	std::vector< GlobalVector * > const &	x,
		std::vector< GlobalVector * > const &	x_prev,
		int const	process_id )

overridevirtual

Implements NumLib::NonlinearSolverBase.

Definition at line 81 of file NonlinearSolver.cpp.

{
    if (!_compensate_non_equilibrium_initial_residuum)
    {
        return;
    }
 
    INFO("Calculate non-equilibrium initial residuum.");
 
    auto& A = NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
    auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(_rhs_id);
    _equation_system->assemble(x, x_prev, process_id);
    _equation_system->getA(A);
    _equation_system->getRhs(*x_prev[process_id], rhs);
 
    // r_neq = A * x - rhs
    _r_neq = &NumLib::GlobalVectorProvider::provider.getVector(_r_neq_id);
    MathLib::LinAlg::matMult(A, *x[process_id], *_r_neq);
    MathLib::LinAlg::axpy(*_r_neq, -1.0, rhs);  // res -= rhs
 
    // 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);
 
    NumLib::GlobalMatrixProvider::provider.releaseMatrix(A);
    NumLib::GlobalVectorProvider::provider.releaseVector(rhs);
}

References NonlinearSolver(), _A_id, _compensate_non_equilibrium_initial_residuum, _equation_system, _r_neq, _r_neq_id, _rhs_id, MathLib::LinAlg::axpy(), MathLib::LinAlg::finalizeAssembly(), INFO(), MathLib::LinAlg::matMult(), NumLib::GlobalMatrixProvider::provider, and NumLib::GlobalVectorProvider::provider.

◆ compensateNonEquilibriumInitialResiduum()

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

inline

Definition at line 230 of file NonlinearSolver.h.

    {
        _compensate_non_equilibrium_initial_residuum = value;
    }

References _compensate_non_equilibrium_initial_residuum.

◆ setEquationSystem()

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

inline

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

Definition at line 212 of file NonlinearSolver.h.

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

References _convergence_criterion, and _equation_system.

◆ solve()

NonlinearSolverStatus NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::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 133 of file NonlinearSolver.cpp.

{
    namespace LinAlg = MathLib::LinAlg;
    auto& sys = *_equation_system;
 
    if (_damping != 1.0 && sys.isLinear())
    {
        OGS_FATAL(
            "Damping (under-relaxation) is not compatible with a linear "
            "equation system: a single Picard step already yields the exact "
            "solution, so the damped iterate would be accepted as converged "
            "but wrong. Remove the 'damping' parameter for linear problems.");
    }
 
    auto& A = NumLib::GlobalMatrixProvider::provider.getMatrix(_A_id);
    auto& rhs = NumLib::GlobalVectorProvider::provider.getVector(_rhs_id);
 
    std::vector<GlobalVector*> x_new{x};
    x_new[process_id] =
        &NumLib::GlobalVectorProvider::provider.getVector(_x_new_id);
    LinAlg::copy(*x[process_id], *x_new[process_id]);  // set initial guess
 
    bool error_norms_met = false;
 
    _convergence_criterion->preFirstIteration();
 
    int iteration = 1;
    for (; iteration <= _maxiter; ++iteration, _convergence_criterion->reset())
    {
        BaseLib::RunTime timer_dirichlet;
        double time_dirichlet = 0.0;
 
        BaseLib::RunTime time_iteration;
        time_iteration.start();
 
        INFO("Iteration #{:d} started.", iteration);
        timer_dirichlet.start();
        auto& x_new_process = *x_new[process_id];
        LinAlg::setLocalAccessibleVector(x_new_process);
        sys.computeKnownSolutions(x_new_process, process_id);
        sys.applyKnownSolutions(x_new_process);
        time_dirichlet += timer_dirichlet.elapsed();
 
        sys.preIteration(iteration, x_new_process);
 
        BaseLib::RunTime time_assembly;
        time_assembly.start();
        sys.assemble(x_new, x_prev, process_id);
        sys.getA(A);
        sys.getRhs(*x_prev[process_id], rhs);
 
        // Normalize the linear equation system, if required
        if (sys.requiresNormalization() &&
            !_linear_solver.canSolveRectangular())
        {
            sys.getAandRhsNormalized(A, rhs);
            WARN(
                "The equation system is rectangular, but the current linear "
                "solver only supports square systems. "
                "The system will be normalized, which lead to a squared "
                "condition number and potential numerical issues. "
                "It is recommended to use a solver that supports rectangular "
                "equation systems for better numerical stability.");
        }
 
        INFO("[time] Assembly took {:g} s.", time_assembly.elapsed());
 
        // Subtract non-equilibrium initial residuum if set
        if (_r_neq != nullptr)
        {
            LinAlg::axpy(rhs, -1, *_r_neq);
        }
 
        auto const solver_needs_to_compute = sys.linearSolverNeedsToCompute();
        bool const solver_will_compute =
            _linear_solver.willCompute(solver_needs_to_compute);
 
        timer_dirichlet.start();
        sys.applyKnownSolutionsPicard(
            A, rhs, x_new_process,
            solver_will_compute
                ? MathLib::DirichletBCApplicationMode::COMPLETE_MATRIX_UPDATE
                : MathLib::DirichletBCApplicationMode::
                      FAST_INCOMPLETE_MATRIX_UPDATE);
        time_dirichlet += timer_dirichlet.elapsed();
        INFO("[time] Applying Dirichlet BCs took {:g} s.", time_dirichlet);
 
        if (!sys.isLinear() && _convergence_criterion->hasResidualCheck())
        {
            if (!solver_will_compute)
            {
                // !solver_will_compute means that the Dirichlet BC application
                // is incomplete (i.e., A not properly modified) and the
                // computed residual is wrong.
                OGS_FATAL(
                    "Logic error. The solver skips the compute step for a "
                    "non-linear equation system.");
            }
            GlobalVector res;
            LinAlg::matMult(A, x_new_process, res);  // res = A * x_new
            LinAlg::axpy(res, -1.0, rhs);            // res -= rhs
            _convergence_criterion->checkResidual(res);
        }
 
        bool iteration_succeeded = detail::solvePicard(
            _linear_solver, A, rhs, x_new_process, solver_needs_to_compute);
 
        if (iteration_succeeded)
        {
            // Notation (see NonlinearSolver<Picard> class doc):
            //   x_k      = x[process_id]  (iterate entering this step)
            //   g(x_k)   = x_new_process  (raw Picard output, this linear
            //   solve)
            // Under-relaxation with damping beta (active when beta != 1):
            //   x_{k+1} = x_k + beta*(g(x_k) - x_k)
            //           = (1-beta)*x_k + beta*g(x_k)
            if (_damping != 1.0)
            {
                LinAlg::scale(x_new_process, _damping);
                LinAlg::axpy(x_new_process, 1.0 - _damping, *x[process_id]);
            }
 
            if (postIterationCallback)
            {
                postIterationCallback(iteration, error_norms_met, x_new);
            }
 
            switch (sys.postIteration(x_new_process))
            {
                case IterationResult::SUCCESS:
                    // Don't copy here. The old x might still be used further
                    // below. Although currently it is not.
                    break;
                case IterationResult::FAILURE:
                    ERR("Picard: The postIteration() hook reported a "
                        "non-recoverable error.");
                    iteration_succeeded = false;
                    // Copy new solution to x.
                    // Thereby the failed solution can be used by the caller for
                    // debugging purposes.
                    LinAlg::copy(x_new_process, *x[process_id]);
                    break;
                case IterationResult::REPEAT_ITERATION:
                    INFO(
                        "Picard: The postIteration() hook decided that this "
                        "iteration has to be repeated.");
                    LinAlg::copy(
                        *x[process_id],
                        x_new_process);  // throw the iteration result away
                    continue;
            }
        }
 
        if (!iteration_succeeded)
        {
            // Don't compute error norms, break here.
            error_norms_met = false;
            break;
        }
 
        if (sys.isLinear())
        {
            error_norms_met = true;
        }
        else
        {
            if (_convergence_criterion->hasDeltaXCheck())
            {
                GlobalVector minus_delta_x(*x[process_id]);
                LinAlg::axpy(minus_delta_x, -1.0,
                             x_new_process);  // minus_delta_x = x - x_new
                _convergence_criterion->checkDeltaX(minus_delta_x,
                                                    x_new_process);
            }
 
            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_process, *x[process_id]);
 
        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("Picard: Could not solve the given nonlinear system within {:d} "
            "iterations",
            _maxiter);
    }
 
    NumLib::GlobalMatrixProvider::provider.releaseMatrix(A);
    NumLib::GlobalVectorProvider::provider.releaseVector(rhs);
    NumLib::GlobalVectorProvider::provider.releaseVector(*x_new[process_id]);
 
    return {error_norms_met, iteration};
}

References NonlinearSolver(), _A_id, _convergence_criterion, _damping, _equation_system, _linear_solver, _maxiter, _r_neq, _rhs_id, _x_new_id, MathLib::LinAlg::axpy(), MathLib::COMPLETE_MATRIX_UPDATE, MathLib::LinAlg::copy(), BaseLib::RunTime::elapsed(), ERR(), NumLib::FAILURE, INFO(), MathLib::LinAlg::matMult(), OGS_FATAL, NumLib::GlobalMatrixProvider::provider, NumLib::GlobalVectorProvider::provider, NumLib::REPEAT_ITERATION, MathLib::LinAlg::scale(), MathLib::LinAlg::setLocalAccessibleVector(), NumLib::detail::solvePicard(), BaseLib::RunTime::start(), NumLib::SUCCESS, and WARN().

Member Data Documentation

◆ _A_id

std::size_t NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::_A_id = 0u

private

ID of the \( A \) matrix.

Definition at line 250 of file NonlinearSolver.h.

Referenced by calculateNonEquilibriumInitialResiduum(), and solve().

◆ _compensate_non_equilibrium_initial_residuum

bool NumLib::NonlinearSolver< NonlinearSolverTag::Picard >::_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 259 of file NonlinearSolver.h.

Referenced by calculateNonEquilibriumInitialResiduum(), and compensateNonEquilibriumInitialResiduum().