butcher_tableau.h

Go to the documentation of this file.

//    |  /           |
//    ' /   __| _` | __|  _ \   __|
//    . \  |   (   | |   (   |\__ `
//   _|\_\_|  \__,_|\__|\___/ ____/
//                   Multi-Physics
//
//  License:         BSD License
//                   Kratos default license: kratos/license.txt
//
//  Main authors:    Eduard Gomez
//
//
 
#pragma once
 
// System includes
 
// External includes
 
// Project includes
#include "containers/array_1d.h"
 
namespace Kratos
{
 
 
 
 
 
 
 
 
 
 
template<typename Derived, unsigned int TOrder, unsigned int TSubstepCount>
class ButcherTableau
{
 
public:
 
    using VectorType = std::vector<double>;
 
    /* Using the following constructs allows us to multiply parts of vectors with parts of matrices
     * while avoiding BOOST's size checks. This is useful to skip multiplications by zero, since
     * for all explicit runge-kutta methods a_ij = 0 for i>j
     */
 
    using RowType = std::vector<double>;
    using MatrixType = std::vector<RowType>;
 
    static constexpr unsigned int Order() {return TOrder;}
    static constexpr unsigned int SubstepCount() {return TSubstepCount; }
 
 
    virtual ~ButcherTableau() = default;
 
 
 
 
    std::tuple<RowType::const_iterator, RowType::const_iterator> GetMatrixRow(const unsigned int SubStepIndex) const
    {
        return {GetMatrixRowBegin(SubStepIndex), GetMatrixRowEnd(SubStepIndex)};
    }
 
    RowType::const_iterator GetMatrixRowBegin(const unsigned int SubStepIndex) const
    {
        KRATOS_DEBUG_ERROR_IF(SubStepIndex == 0) << "Provided substep is 0. This must be greater than or equal to 1." << std::endl;
        return mA[SubStepIndex - 1].begin();
    }
 
    RowType::const_iterator GetMatrixRowEnd(const unsigned int SubStepIndex) const
    {
        KRATOS_DEBUG_ERROR_IF(SubStepIndex == 0) << "Provided substep is 0. This must be greater than or equal to 1." << std::endl;
        return mA[SubStepIndex - 1].begin() + SubStepIndex;
    }
 
    constexpr const VectorType& GetWeights() const
    {
        return mB;
    }
 
    constexpr double GetIntegrationTheta(const unsigned int SubStepIndex) const
    {
        KRATOS_DEBUG_ERROR_IF(SubStepIndex == 0) << "Provided substep is 0. This must be greater than or equal to 1." << std::endl;
        return mC[SubStepIndex - 1];
    }
 
 
    static std::string Name()
    {
        return Derived::Name();
    }
 
    virtual std::string Info() const = 0;
 
protected:
 
 
 
    const MatrixType mA = Derived::GenerateRKMatrix();     // Runge-Kutta matrix
    const VectorType mB = Derived::GenerateWeights();      // Weights vector
    const VectorType mC = Derived::GenerateThetasVector(); // Nodes vector
 
 
 
 
 
 
 
 
 
private:
 
 
 
 
 
 
 
 
 
 
 
 
 
};
 
 
 
 
 
class ButcherTableauForwardEuler : public ButcherTableau<ButcherTableauForwardEuler, 1, 1>
{
public:
    typedef ButcherTableau<ButcherTableauForwardEuler, 1, 1> BaseType;
 
    static const BaseType::MatrixType GenerateRKMatrix()
    {
        return {};
    }
 
    static const BaseType::VectorType GenerateWeights()
    {
        return VectorType {1.0};
    }
 
    static const BaseType::VectorType GenerateThetasVector()
    {
        return VectorType {0.0};
    }
 
    static std::string Name()
    {
        return "butcher_tableau_forward_euler";
    }
 
    std::string Info() const override
    {
        return "ButcherTableauForwardEuler";
    }
};
 
 
class ButcherTableauMidPointMethod : public ButcherTableau<ButcherTableauMidPointMethod, 2, 2>
{
public:
    typedef ButcherTableau<ButcherTableauMidPointMethod, 2, 2> BaseType;
 
    static const BaseType::MatrixType GenerateRKMatrix()
    {
        MatrixType A(1, RowType(1));
        A[0][0] = 0.5;
        return A;
    }
 
    static const BaseType::VectorType GenerateWeights()
    {
        return VectorType {0.0, 1.0};
    }
 
    static const BaseType::VectorType GenerateThetasVector()
    {
        return VectorType {0.0, 0.5};
    }
 
    static std::string Name()
    {
        return "butcher_tableau_midpoint_method";
    }
 
    std::string Info() const override
    {
        return "ButcherTableauMidPointMethod";
    }
};
 
class ButcherTableauRK3TVD : public ButcherTableau<ButcherTableauRK3TVD, 3, 3>
{
public:
    typedef ButcherTableau<ButcherTableauRK3TVD, 3, 3> BaseType;
 
    static const BaseType::MatrixType GenerateRKMatrix()
    {
        MatrixType A(2, RowType(2));
        A[0][0] = 1;
        A[1][0] = 0.25;
        A[0][1] = 0.0;
        A[1][1] = 0.25;
        return A;
    }
 
    static const BaseType::VectorType GenerateWeights()
    {
        return VectorType {1.0 / 6.0,
                           1.0 / 6.0,
                           2.0 / 3.0};
    }
 
    static const BaseType::VectorType GenerateThetasVector()
    {
        return VectorType {0.0,
                           1.0,
                           0.5};
    }
 
    static std::string Name()
    {
        return "butcher_tableau_RK3TVD";
    }
 
    std::string Info() const override
    {
        return "ButcherTableauRK3TVD";
    }
};
 
 
class ButcherTableauRK4 : public ButcherTableau<ButcherTableauRK4, 4, 4>
{
public:
    typedef ButcherTableau<ButcherTableauRK4, 4, 4> BaseType;
    static const BaseType::MatrixType GenerateRKMatrix()
    {
        MatrixType A(3, RowType(3,0.0));
        A[0][0] = 0.5;
        A[1][1] = 0.5;
        A[2][2] = 1.0;
        return A;
    }
 
    static const BaseType::VectorType GenerateWeights()
    {
        return VectorType {1.0 / 6.0,
                           1.0 / 3.0,
                           1.0 / 3.0,
                           1.0 / 6.0};
    }
 
    static const BaseType::VectorType GenerateThetasVector()
    {
        return VectorType {0.0,
                           0.5,
                           0.5,
                           1.0};
    }
 
    static std::string Name()
    {
        return "butcher_tableau_RK4";
    }
 
    std::string Info() const override
    {
        return "ButcherTableauRK4";
    }
};
 
 
} /* namespace Kratos.*/