MatrixTranslator.h

Go to the documentation of this file.

#pragma once

#include <memory>

#include "TimeDiscretization.h"
#include "Types.h"

#include "MathLib/LinAlg/LinAlg.h"

namespace NumLib
{

template <ODESystemTag ODETag>
class MatrixTranslator;

template <>
class MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
{
public:
    virtual void computeA(GlobalMatrix const& M, GlobalMatrix const& K,
                          GlobalMatrix& A) const = 0;

    virtual void computeRhs(const GlobalMatrix& M, const GlobalMatrix& K,
                            const GlobalVector& b, GlobalVector& rhs) const = 0;

    virtual void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                                 GlobalVector const& b,
                                 GlobalVector const& x_new_timestep,
                                 GlobalVector const& xdot,
                                 GlobalVector& res) const = 0;

    virtual void computeJacobian(GlobalMatrix const& Jac_in,
                                 GlobalMatrix& Jac_out) const = 0;

    virtual void pushMatrices(GlobalMatrix const& /*M*/,
                              GlobalMatrix const& /*K*/,
                              GlobalVector const& /*b*/)
    {
    }

    virtual ~MatrixTranslator() = default;
};

template <ODESystemTag ODETag>
class MatrixTranslatorGeneral;

template <>
class MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>
    : public MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
{
public:
    MatrixTranslatorGeneral(TimeDiscretization const& timeDisc)
        : _time_disc(timeDisc)
    {
    }

    void computeA(GlobalMatrix const& M, GlobalMatrix const& K,
                  GlobalMatrix& A) const override;

    void computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
                    const GlobalVector& b, GlobalVector& rhs) const override;

    void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                         GlobalVector const& b,
                         GlobalVector const& x_new_timestep,
                         GlobalVector const& xdot,
                         GlobalVector& res) const override;

    void computeJacobian(GlobalMatrix const& Jac_in,
                         GlobalMatrix& Jac_out) const override;

private:
    TimeDiscretization const&
        _time_disc;  

    mutable std::size_t _tmp_id = 0u;
};

template <ODESystemTag ODETag>
class MatrixTranslatorForwardEuler;

template <>
class MatrixTranslatorForwardEuler<ODESystemTag::FirstOrderImplicitQuasilinear>
    : public MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
{
public:
    MatrixTranslatorForwardEuler(ForwardEuler const& timeDisc)
        : _fwd_euler(timeDisc)
    {
    }

    void computeA(GlobalMatrix const& M, GlobalMatrix const& /*K*/,
                  GlobalMatrix& A) const override;

    void computeRhs(const GlobalMatrix& M, const GlobalMatrix& K,
                    const GlobalVector& b, GlobalVector& rhs) const override;

    void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                         GlobalVector const& b,
                         GlobalVector const& x_new_timestep,
                         GlobalVector const& xdot,
                         GlobalVector& res) const override;

    void computeJacobian(GlobalMatrix const& Jac_in,
                         GlobalMatrix& Jac_out) const override;

private:
    ForwardEuler const& _fwd_euler;  

    mutable std::size_t _tmp_id = 0u;
};

template <ODESystemTag ODETag>
class MatrixTranslatorCrankNicolson;

template <>
class MatrixTranslatorCrankNicolson<ODESystemTag::FirstOrderImplicitQuasilinear>
    : public MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
{
public:
    MatrixTranslatorCrankNicolson(CrankNicolson const& timeDisc)
        : _crank_nicolson(timeDisc),
          _M_bar(NumLib::GlobalMatrixProvider::provider
                     .getMatrix()),
          _b_bar(
              NumLib::GlobalVectorProvider::provider.getVector())
    {
    }

    ~MatrixTranslatorCrankNicolson() override
    {
        NumLib::GlobalMatrixProvider::provider.releaseMatrix(
            _M_bar);
        NumLib::GlobalVectorProvider::provider.releaseVector(
            _b_bar);
    }

    void computeA(GlobalMatrix const& M, GlobalMatrix const& K,
                  GlobalMatrix& A) const override;

    void computeRhs(const GlobalMatrix& M, const GlobalMatrix& /*K*/,
                    const GlobalVector& b, GlobalVector& rhs) const override;

    void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                         GlobalVector const& b,
                         GlobalVector const& x_new_timestep,
                         GlobalVector const& xdot,
                         GlobalVector& res) const override;

    void computeJacobian(GlobalMatrix const& Jac_in,
                         GlobalMatrix& Jac_out) const override;

    void pushMatrices(GlobalMatrix const& M, GlobalMatrix const& K,
                      GlobalVector const& b) override;

private:
    CrankNicolson const& _crank_nicolson;

    GlobalMatrix&
        _M_bar;  
    GlobalVector& _b_bar;  

    mutable std::size_t _tmp_id = 0u;
};

template <ODESystemTag ODETag>
std::unique_ptr<MatrixTranslator<ODETag>> createMatrixTranslator(
    TimeDiscretization const& timeDisc)
{
    if (auto* fwd_euler = dynamic_cast<ForwardEuler const*>(&timeDisc))
    {
        return std::unique_ptr<MatrixTranslator<ODETag>>(
            new MatrixTranslatorForwardEuler<ODETag>(*fwd_euler));
    }
    if (auto* crank = dynamic_cast<CrankNicolson const*>(&timeDisc))
    {
        return std::unique_ptr<MatrixTranslator<ODETag>>(
            new MatrixTranslatorCrankNicolson<ODETag>(*crank));
    }

    return std::unique_ptr<MatrixTranslator<ODETag>>(
        new MatrixTranslatorGeneral<ODETag>(timeDisc));
}

}