NonlinearSolver.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) OpenGeoSys Community (opengeosys.org)
// SPDX-License-Identifier: BSD-3-Clause
 
#pragma once
 
#include <memory>
#include <utility>
 
#include "ConvergenceCriterion.h"
#include "MathLib/LinAlg/GlobalLinearSolverType.h"
#include "NewtonStepStrategy.h"
#include "NonlinearSolverStatus.h"
#include "NonlinearSystem.h"
#include "Types.h"
 
namespace BaseLib
{
class ConfigTree;
}
 
// TODO Document in the ODE solver lib, which matrices and vectors that are
// passed around as method arguments are guaranteed to be of the right size
// (and zeroed out) and which are not.
 
namespace NumLib
{
class NonlinearSolverBase
{
public:
    virtual void calculateNonEquilibriumInitialResiduum(
        std::vector<GlobalVector*> const& x,
        std::vector<GlobalVector*> const& x_prev, int const process_id) = 0;

    virtual NonlinearSolverStatus solve(
        std::vector<GlobalVector*>& x,
        std::vector<GlobalVector*> const& x_prev,
        std::function<void(int, std::vector<GlobalVector*> const&)> const&
            postIterationCallback,
        int const process_id) = 0;
 
    virtual ~NonlinearSolverBase() = default;
};


template <NonlinearSolverTag NLTag>
class NonlinearSolver;

template <>
class NonlinearSolver<NonlinearSolverTag::Newton> final
    : public NonlinearSolverBase
{
public:
    using System = NonlinearSystem<NonlinearSolverTag::Newton>;

    explicit NonlinearSolver(GlobalLinearSolver& linear_solver,
                             int const maxiter,
                             std::unique_ptr<NewtonStepStrategy>
                                 newton_strategy,
                             int const recompute_jacobian = 1)
        : _linear_solver(linear_solver),
          _maxiter(maxiter),
          _step_strategy(std::move(newton_strategy)),
          _recompute_jacobian(recompute_jacobian)
    {
    }
 
    ~NonlinearSolver();

    void setEquationSystem(System& eq, ConvergenceCriterion& conv_crit)
    {
        _equation_system = &eq;
        _convergence_criterion = &conv_crit;
        _step_strategy->setDampingPolicy(
            _convergence_criterion->dampingPolicy());
    }
 
    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, std::vector<GlobalVector*> const&)> const&
            postIterationCallback,
        int const process_id) override;
 
    void compensateNonEquilibriumInitialResiduum(bool const value)
    {
        _compensate_non_equilibrium_initial_residuum = value;
    }
 
    void setTikhonovLambda(double const lambda, int const starting_iteration)
    {
        _tikhonov_lambda = lambda;
        _tikhonov_starting_iteration = starting_iteration;
    }
 
private:
    GlobalLinearSolver& _linear_solver;
    System* _equation_system = nullptr;
 
    int const _maxiter;  

    std::unique_ptr<NewtonStepStrategy> _step_strategy;

    ConvergenceCriterion* _convergence_criterion = nullptr;
 
    int const _recompute_jacobian =
        1;  
 
    GlobalVector* _r_neq = nullptr;      
    std::size_t _res_id = 0u;            
    std::size_t _J_id = 0u;              
    std::size_t _minus_delta_x_id = 0u;  
    std::size_t _x_new_id =
        0u;  
    std::size_t _r_neq_id = 0u;  

    bool _compensate_non_equilibrium_initial_residuum = false;
    double _tikhonov_lambda = 0.0;  
    int _tikhonov_starting_iteration =
        0;  
};

template <>
class NonlinearSolver<NonlinearSolverTag::Picard> final
    : public NonlinearSolverBase
{
public:
    using System = NonlinearSystem<NonlinearSolverTag::Picard>;

    explicit NonlinearSolver(GlobalLinearSolver& linear_solver,
                             int const maxiter, int const anderson_depth,
                             double const damping)
        : _linear_solver(linear_solver),
          _damping(damping),
          _anderson_depth(anderson_depth),
          _maxiter(maxiter)
    {
    }
 
    ~NonlinearSolver();

    void setEquationSystem(System& eq, ConvergenceCriterion& conv_crit)
    {
        _equation_system = &eq;
        _convergence_criterion = &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, std::vector<GlobalVector*> const&)> const&
            postIterationCallback,
        int const process_id) override;
 
    void compensateNonEquilibriumInitialResiduum(bool const value)
    {
        _compensate_non_equilibrium_initial_residuum = value;
    }
 
private:
    GlobalLinearSolver& _linear_solver;
    System* _equation_system = nullptr;
 
    // TODO doc
    ConvergenceCriterion* _convergence_criterion = nullptr;

    double const _damping;

    int const _anderson_depth;
 
    int const _maxiter;  
 
    GlobalVector* _r_neq = nullptr;  
    std::size_t _A_id = 0u;          
    std::size_t _rhs_id = 0u;        
    std::size_t _x_new_id = 0u;  
    std::size_t _r_neq_id = 0u;  
 
    // clang-format off
    bool _compensate_non_equilibrium_initial_residuum = false;
    // clang-format on
};

 
}  // namespace NumLib