MatrixTranslator.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 "MathLib/LinAlg/LinAlg.h"
#include "TimeDiscretization.h"
#include "Types.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, const GlobalVector& x_prev,
                            GlobalVector& rhs) const = 0;

    virtual void normalizeAandRhs(GlobalMatrix& A, GlobalVector& b) const = 0;

    virtual void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                                 GlobalVector const& b, double dt,
                                 GlobalVector const& x_curr,
                                 GlobalVector const& x_prev,
                                 GlobalVector& res) const = 0;

    virtual void computeJacobian(GlobalMatrix const& Jac_in,
                                 GlobalMatrix& Jac_out) const = 0;
 
    virtual ~MatrixTranslator() = default;
};

template <ODESystemTag ODETag>
class MatrixTranslatorGeneral;

template <>
class MatrixTranslatorGeneral<ODESystemTag::FirstOrderImplicitQuasilinear>
    : public MatrixTranslator<ODESystemTag::FirstOrderImplicitQuasilinear>
{
public:
    explicit 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, const GlobalVector& x_prev,
                    GlobalVector& rhs) const override;

    void normalizeAandRhs(GlobalMatrix& A, GlobalVector& b) const override;

    void computeResidual(GlobalMatrix const& M, GlobalMatrix const& K,
                         GlobalVector const& b, double dt,
                         GlobalVector const& x_curr, GlobalVector const& x_prev,
                         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>
std::unique_ptr<MatrixTranslator<ODETag>> createMatrixTranslator(
    TimeDiscretization const& timeDisc)
{
    return std::unique_ptr<MatrixTranslator<ODETag>>(
        new MatrixTranslatorGeneral<ODETag>(timeDisc));
}

}  // namespace NumLib