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 "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, bool, 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,
        int const recompute_jacobian = 1,
        double const damping = 1.0,
        std::optional<double> const damping_reduction = std::nullopt)
        : _linear_solver(linear_solver),
          _maxiter(maxiter),
          _recompute_jacobian(recompute_jacobian),
          _damping(damping),
          _damping_reduction(damping_reduction)
    {
    }
 
    ~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, bool, 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;
    int const _maxiter;  
 
    int const _recompute_jacobian = 1;

    double const _damping;

    std::optional<double> const _damping_reduction;
 
    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;
};

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

    explicit NonlinearSolver(GlobalLinearSolver& linear_solver,
                             const int maxiter)
        : _linear_solver(linear_solver), _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, bool, 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;
    const int _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
};

std::pair<std::unique_ptr<NonlinearSolverBase>, NonlinearSolverTag>
createNonlinearSolver(GlobalLinearSolver& linear_solver,
                      BaseLib::ConfigTree const& config);

 
}  // namespace NumLib