TimeDiscretizedODESystem.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 "MatrixTranslator.h"
#include "NonlinearSystem.h"
#include "ODESystem.h"
#include "TimeDiscretization.h"
 
namespace NumLib
{

template <NonlinearSolverTag NLTag>
class TimeDiscretizedODESystemBase : public NonlinearSystem<NLTag>
{
public:
    virtual TimeDiscretization& getTimeDiscretization() = 0;
};

template <ODESystemTag ODETag, NonlinearSolverTag NLTag>
class TimeDiscretizedODESystem;

template <>
class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                               NonlinearSolverTag::Newton>
    final : public TimeDiscretizedODESystemBase<NonlinearSolverTag::Newton>
{
public:
    static const ODESystemTag ODETag =
        ODESystemTag::FirstOrderImplicitQuasilinear;

    using ODE = ODESystem<ODETag, NonlinearSolverTag::Newton>;
    using MatTrans = MatrixTranslator<ODETag>;
    using TimeDisc = TimeDiscretization;

    explicit TimeDiscretizedODESystem(const int process_id, ODE& ode,
                                      TimeDisc& time_discretization);
 
    ~TimeDiscretizedODESystem() override;
 
    void assemble(std::vector<GlobalVector*> const& x_new_timestep,
                  std::vector<GlobalVector*> const& x_prev,
                  int const process_id) override;

    std::vector<GlobalIndexType>
    getIndicesOfResiduumWithoutInitialCompensation() const override
    {
        return _ode.getIndicesOfResiduumWithoutInitialCompensation();
    }
 
    void setReleaseNodalForces(GlobalVector const* r_neq,
                               int const process_id) override
    {
        _ode.setReleaseNodalForces(r_neq, process_id);
    }
 
    void getResidual(GlobalVector const& x_new_timestep,
                     GlobalVector const& x_prev,
                     GlobalVector& res) const override;
 
    void getJacobian(GlobalMatrix& Jac) const override;
 
    void computeKnownSolutions(GlobalVector const& x,
                               int const process_id) override;
 
    void applyKnownSolutions(GlobalVector& x) const override;
 
    void applyKnownSolutionsNewton(GlobalMatrix& Jac, GlobalVector& res,
                                   GlobalVector const& x,
                                   GlobalVector& minus_delta_x) const override;
 
    void applyKnownSolutionsPETScSNES(GlobalMatrix& Jac, GlobalVector& res,
                                      GlobalVector& x) const override;
 
    void updateConstraints(GlobalVector& lower, GlobalVector& upper,
                           int const process_id) override
    {
        _ode.updateConstraints(lower, upper, process_id);
    }
 
    bool isLinear() const override { return _ode.isLinear(); }
 
    bool requiresNormalization() const override
    {
        return _ode.requiresNormalization();
    }
 
    void preIteration(const unsigned iter, GlobalVector const& x) override
    {
        _ode.preIteration(iter, x);
    }
 
    IterationResult postIteration(GlobalVector const& x) override
    {
        return _ode.postIteration(x);
    }
 
    TimeDisc& getTimeDiscretization() override { return _time_disc; }
    MathLib::MatrixSpecifications getMatrixSpecifications(
        const int process_id) const override
    {
        return _ode.getMatrixSpecifications(process_id);
    }
 
private:
    ODE& _ode;             
    TimeDisc& _time_disc;  

    std::unique_ptr<MatTrans> _mat_trans;
 
    using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
    std::vector<NumLib::IndexValueVector<Index>> const* _known_solutions =
        nullptr;  
 
    GlobalMatrix* _Jac;  
    GlobalVector* _b;    
 
    std::size_t _Jac_id = 0u;  
    std::size_t _b_id = 0u;    
};

template <>
class TimeDiscretizedODESystem<ODESystemTag::FirstOrderImplicitQuasilinear,
                               NonlinearSolverTag::Picard>
    final : public TimeDiscretizedODESystemBase<NonlinearSolverTag::Picard>
{
public:
    static const ODESystemTag ODETag =
        ODESystemTag::FirstOrderImplicitQuasilinear;

    using ODE = ODESystem<ODETag, NonlinearSolverTag::Picard>;
    using MatTrans = MatrixTranslator<ODETag>;
    using TimeDisc = TimeDiscretization;

    explicit TimeDiscretizedODESystem(const int process_id, ODE& ode,
                                      TimeDisc& time_discretization);
 
    ~TimeDiscretizedODESystem() override;
 
    void assemble(std::vector<GlobalVector*> const& x_new_timestep,
                  std::vector<GlobalVector*> const& x_prev,
                  int const process_id) override;

    std::vector<GlobalIndexType>
    getIndicesOfResiduumWithoutInitialCompensation() const override
    {
        return _ode.getIndicesOfResiduumWithoutInitialCompensation();
    }
 
    void setReleaseNodalForces(GlobalVector const* r_neq,
                               int const process_id) override
    {
        _ode.setReleaseNodalForces(r_neq, process_id);
    }
 
    void getA(GlobalMatrix& A) const override
    {
        _mat_trans->computeA(*_M, *_K, A);
    }
 
    void getRhs(GlobalVector const& x_prev, GlobalVector& rhs) const override
    {
        _mat_trans->computeRhs(*_M, *_K, *_b, x_prev, rhs);
    }
 
    void getAandRhsNormalized(GlobalMatrix& A, GlobalVector& rhs) const override
    {
        _mat_trans->normalizeAandRhs(A, rhs);
    }
 
    void computeKnownSolutions(GlobalVector const& x,
                               int const process_id) override;
 
    void applyKnownSolutions(GlobalVector& x) const override;
 
    void applyKnownSolutionsPicard(
        GlobalMatrix& A, GlobalVector& rhs, GlobalVector& x,
        MathLib::DirichletBCApplicationMode const mode) const override;
 
    bool isLinear() const override { return _ode.isLinear(); }
 
    bool requiresNormalization() const override
    {
        return _ode.requiresNormalization();
    }
 
    void preIteration(const unsigned iter, GlobalVector const& x) override
    {
        _ode.preIteration(iter, x);
    }
 
    IterationResult postIteration(GlobalVector const& x) override
    {
        return _ode.postIteration(x);
    }
 
    TimeDisc& getTimeDiscretization() override { return _time_disc; }
    MathLib::MatrixSpecifications getMatrixSpecifications(
        const int process_id) const override
    {
        return _ode.getMatrixSpecifications(process_id);
    }
 
    MathLib::LinearSolverBehaviour linearSolverNeedsToCompute() const override
    {
        if (_ode.shouldLinearSolverComputeOnlyUponTimestepChange() &&
            _time_disc.getCurrentTimeIncrement() !=
                _time_disc.getPreviousTimeIncrement())
        {
            return MathLib::LinearSolverBehaviour::RECOMPUTE_AND_STORE;
        }
        if (_ode.shouldLinearSolverComputeOnlyUponTimestepChange() &&
            _time_disc.getCurrentTimeIncrement() ==
                _time_disc.getPreviousTimeIncrement())
        {
            return MathLib::LinearSolverBehaviour::REUSE;
        }
        return MathLib::LinearSolverBehaviour::RECOMPUTE;
    }
 
private:
    ODE& _ode;             
    TimeDisc& _time_disc;  

    std::unique_ptr<MatTrans> _mat_trans;
 
    using Index = MathLib::MatrixVectorTraits<GlobalMatrix>::Index;
    std::vector<NumLib::IndexValueVector<Index>> const* _known_solutions =
        nullptr;  
 
    GlobalMatrix* _M;  
    GlobalMatrix* _K;  
    GlobalVector* _b;  
 
    std::size_t _M_id = 0u;  
    std::size_t _K_id = 0u;  
    std::size_t _b_id = 0u;  
};

}  // namespace NumLib