multistagemethods.hh

Go to the documentation of this file.

// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
//
// SPDX-FileCopyrightText: Copyright © DuMux Project contributors, see AUTHORS.md in root folder
// SPDX-License-Identifier: GPL-3.0-or-later
// SPDX-FileCopyrightText: Copyright © dune-pdelab developers, see LICENSE.md (permalink below)
// SPDX-License-Identifier: GPL-3.0-or-later
// The code is based on the implementation of time stepping method parameters
// in dune-pdelab (https://archive.softwareheritage.org/swh:1:cnt:9c3d72412f8e4d48d84a0090b2bb4362b5a0d843)
// licensed under GPL-3.0-or-later, see their LICENSE.md for a full list of copyright holders at
// https://archive.softwareheritage.org/swh:1:cnt:b11b484e74eefe20c74f7043309b1b02853df2eb.
// Modifications (different interface naming, comments, types) are licensed under GPL-3.0-or-later.
//
#ifndef DUMUX_TIMESTEPPING_MULTISTAGE_METHODS_HH
#define DUMUX_TIMESTEPPING_MULTISTAGE_METHODS_HH
 
#include <cmath>
#include <string>
#include <array>
 
namespace Dumux::Experimental {

template<class Scalar>
class MultiStageMethod
{
public:
    virtual bool implicit () const = 0;
 
    virtual std::size_t numStages () const = 0;

    virtual Scalar temporalWeight (std::size_t i, std::size_t k) const = 0;

    virtual Scalar spatialWeight (std::size_t i, std::size_t k) const = 0;

    virtual Scalar timeStepWeight (std::size_t k) const = 0;
 
    virtual std::string name () const = 0;
 
    virtual ~MultiStageMethod() = default;
};

namespace MultiStage {

template<class Scalar>
class Theta : public MultiStageMethod<Scalar>
{
public:
    explicit Theta(const Scalar theta)
    : paramAlpha_{{-1.0, 1.0}}
    , paramBeta_{{1.0-theta, theta}}
    , paramD_{{0.0, 1.0}}
    {}
 
    bool implicit () const final
    { return paramBeta_[1] > 0.0; }
 
    std::size_t numStages () const final
    { return 1; }
 
    Scalar temporalWeight (std::size_t, std::size_t k) const final
    { return paramAlpha_[k]; }
 
    Scalar spatialWeight (std::size_t, std::size_t k) const final
    { return paramBeta_[k]; }
 
    Scalar timeStepWeight (std::size_t k) const final
    { return paramD_[k]; }
 
    std::string name () const override
    { return "theta scheme"; }
 
private:
    std::array<Scalar, 2> paramAlpha_;
    std::array<Scalar, 2> paramBeta_;
    std::array<Scalar, 2> paramD_;
};

template<class Scalar>
class ExplicitEuler final : public Theta<Scalar>
{
public:
    ExplicitEuler() : Theta<Scalar>(0.0) {}
 
    std::string name () const final
    { return "explicit Euler"; }
};

template<class Scalar>
class ImplicitEuler final : public Theta<Scalar>
{
public:
    ImplicitEuler() : Theta<Scalar>(1.0) {}
 
    std::string name () const final
    { return "implicit Euler"; }
};

template<class Scalar>
class RungeKuttaExplicitSecondOrderHeun final : public MultiStageMethod<Scalar>
{
public:
    RungeKuttaExplicitSecondOrderHeun()
    : paramAlpha_{{{-1.0, 1.0, 0.0},
                   {-1.0, 0.0, 1.0}}}
    , paramBeta_{{{1.0, 0.0, 0.0},
                  {0.5, 0.5, 0.0}}}
    , paramD_{{0.0, 1.0, 1.0}}
    {}
 
    bool implicit () const final
    { return false; }
 
    std::size_t numStages () const final
    { return 2; }
 
    Scalar temporalWeight (std::size_t i, std::size_t k) const final
    { return paramAlpha_[i-1][k]; }
 
    Scalar spatialWeight (std::size_t i, std::size_t k) const final
    { return paramBeta_[i-1][k]; }
 
    Scalar timeStepWeight (std::size_t k) const final
    { return paramD_[k]; }
 
    std::string name () const final
    { return "explicit Runge-Kutta 2nd order (Heun)"; }
 
private:
    std::array<std::array<Scalar, 3>, 2> paramAlpha_;
    std::array<std::array<Scalar, 3>, 2> paramBeta_;
    std::array<Scalar, 3> paramD_;
};

template<class Scalar>
class RungeKuttaExplicitThirdOrder final : public MultiStageMethod<Scalar>
{
public:
    RungeKuttaExplicitThirdOrder()
    : paramAlpha_{{{-1.0, 1.0, 0.0, 0.0},
                   {-1.0, 0.0, 1.0, 0.0},
                   {-1.0, 0.0, 0.0, 1.0}}}
    , paramBeta_{{{0.5, 0.0, 0.0, 0.0},
                  {-1.0, 2.0, 0.0, 0.0},
                  {1.0/6.0, 2.0/3.0, 1.0/6.0, 0.0}}}
    , paramD_{{0.0, 0.5, 1.0, 1.0}}
    {}
 
    bool implicit () const final
    { return false; }
 
    std::size_t numStages () const final
    { return 3; }
 
    Scalar temporalWeight (std::size_t i, std::size_t k) const final
    { return paramAlpha_[i-1][k]; }
 
    Scalar spatialWeight (std::size_t i, std::size_t k) const final
    { return paramBeta_[i-1][k]; }
 
    Scalar timeStepWeight (std::size_t k) const final
    { return paramD_[k]; }
 
    std::string name () const final
    { return "explicit Runge-Kutta 3rd order"; }
 
private:
    std::array<std::array<Scalar, 4>, 3> paramAlpha_;
    std::array<std::array<Scalar, 4>, 3> paramBeta_;
    std::array<Scalar, 4> paramD_;
};

template<class Scalar>
class RungeKuttaExplicitFourthOrder final : public MultiStageMethod<Scalar>
{
public:
    RungeKuttaExplicitFourthOrder()
    : paramAlpha_{{{-1.0, 1.0, 0.0, 0.0, 0.0},
                   {-1.0, 0.0, 1.0, 0.0, 0.0},
                   {-1.0, 0.0, 0.0, 1.0, 0.0},
                   {-1.0, 0.0, 0.0, 0.0, 1.0}}}
    , paramBeta_{{{0.5, 0.0, 0.0, 0.0, 0.0},
                  {0.0, 0.5, 0.0, 0.0, 0.0},
                  {0.0, 0.0, 1.0, 0.0, 0.0},
                  {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0, 0.0}}}
    , paramD_{{0.0, 0.5, 0.5, 1.0, 1.0}}
    {}
 
    bool implicit () const final
    { return false; }
 
    std::size_t numStages () const final
    { return 4; }
 
    Scalar temporalWeight (std::size_t i, std::size_t k) const final
    { return paramAlpha_[i-1][k]; }
 
    Scalar spatialWeight (std::size_t i, std::size_t k) const final
    { return paramBeta_[i-1][k]; }
 
    Scalar timeStepWeight (std::size_t k) const final
    { return paramD_[k]; }
 
    std::string name () const final
    { return "explicit Runge-Kutta 4th order"; }
 
private:
    std::array<std::array<Scalar, 5>, 4> paramAlpha_;
    std::array<std::array<Scalar, 5>, 4> paramBeta_;
    std::array<Scalar, 5> paramD_;
};
 

template<class Scalar>
class DIRKSecondOrderAlexander final : public MultiStageMethod<Scalar>
{
public:
    DIRKSecondOrderAlexander()
    {
        using std::sqrt;
        const Scalar gamma = 1.0 - 1.0/sqrt(2.0);
 
        paramD_ = {{0.0, gamma, 1.0}};
        paramAlpha_ = {{
            {-1.0, 1.0, 0.0},
            {-1.0, 0.0, 1.0}
        }};
        paramBeta_ = {{
            {0.0, gamma, 0.0},
            {0.0, 1.0-gamma, gamma}
        }};
    }
 
    bool implicit () const final
    { return true; }
 
    std::size_t numStages () const final
    { return 2; }
 
    Scalar temporalWeight (std::size_t i, std::size_t k) const final
    { return paramAlpha_[i-1][k]; }
 
    Scalar spatialWeight (std::size_t i, std::size_t k) const final
    { return paramBeta_[i-1][k]; }
 
    Scalar timeStepWeight (std::size_t k) const final
    { return paramD_[k]; }
 
    std::string name () const final
    { return "diagonally implicit Runge-Kutta 2nd order (Alexander)"; }
 
private:
    std::array<std::array<Scalar, 3>, 2> paramAlpha_;
    std::array<std::array<Scalar, 3>, 2> paramBeta_;
    std::array<Scalar, 3> paramD_;
};

template<class Scalar>
class DIRKThirdOrderAlexander final : public MultiStageMethod<Scalar>
{
public:
    DIRKThirdOrderAlexander()
    {
        constexpr Scalar alpha = []{
            // Newton iteration for alpha
            Scalar alpha = 0.4358665215; // initial guess
            for (int i = 0; i < 10; ++i)
                alpha = alpha - (alpha*(alpha*alpha-3.0*(alpha-0.5))-1.0/6.0)/(3.0*alpha*(alpha-2.0)+1.5);
            return alpha;
        }();
 
        constexpr Scalar tau2 = (1.0+alpha)*0.5;
        constexpr Scalar b1 = -(6.0*alpha*alpha -16.0*alpha + 1.0)*0.25;
        constexpr Scalar b2 = (6*alpha*alpha - 20.0*alpha + 5.0)*0.25;
 
        paramD_ = {{0.0, alpha, tau2, 1.0}};
        paramAlpha_ = {{
            {-1.0, 1.0, 0.0, 0.0},
            {-1.0, 0.0, 1.0, 0.0},
            {-1.0, 0.0, 0.0, 1.0}
        }};
        paramBeta_ = {{
            {0.0, alpha, 0.0, 0.0},
            {0.0, tau2-alpha, alpha, 0.0},
            {0.0, b1, b2, alpha}
        }};
    }
 
    bool implicit () const final
    { return true; }
 
    std::size_t numStages () const final
    { return 3; }
 
    Scalar temporalWeight (std::size_t i, std::size_t k) const final
    { return paramAlpha_[i-1][k]; }
 
    Scalar spatialWeight (std::size_t i, std::size_t k) const final
    { return paramBeta_[i-1][k]; }
 
    Scalar timeStepWeight (std::size_t k) const final
    { return paramD_[k]; }
 
    std::string name () const final
    { return "diagonally implicit Runge-Kutta 3rd order (Alexander)"; }
 
private:
    std::array<std::array<Scalar, 4>, 3> paramAlpha_;
    std::array<std::array<Scalar, 4>, 3> paramBeta_;
    std::array<Scalar, 4> paramD_;
};

template<class Scalar>
class QinZhangSymplecticDIRK final : public MultiStageMethod<Scalar>
{
public:
    QinZhangSymplecticDIRK()
    : paramAlpha_{{
        {-1.0, 1.0, 0.0, 0.0},
        {-1.0, 0.0, 1.0, 0.0},
        {-1.0, 0.0, 0.0, 1.0}
    }}
    , paramBeta_{{
        {0.0, 0.25, 0.0,  0.0},   // stage 1: a_11 = 1/4
        {0.0, 0.5,  0.25, 0.0},   // stage 2: a_21 = 1/2, a_22 = 1/4
        {0.0, 0.5,  0.5,  0.0}    // update:  b_1  = 1/2, b_2  = 1/2
    }}
    , paramD_{{0.0, 0.25, 0.75, 1.0}}
    {}
 
    bool implicit () const final
    { return true; }
 
    std::size_t numStages () const final
    { return 3; }
 
    Scalar temporalWeight (std::size_t i, std::size_t k) const final
    { return paramAlpha_[i-1][k]; }
 
    Scalar spatialWeight (std::size_t i, std::size_t k) const final
    { return paramBeta_[i-1][k]; }
 
    Scalar timeStepWeight (std::size_t k) const final
    { return paramD_[k]; }
 
    std::string name () const final
    { return "symplectic diagonally implicit Runge-Kutta 2nd order (Qin-Zhang)"; }
 
private:
    std::array<std::array<Scalar, 4>, 3> paramAlpha_;
    std::array<std::array<Scalar, 4>, 3> paramBeta_;
    std::array<Scalar, 4> paramD_;
};
 
} // end namespace MultiStage
} // end namespace Dumux::Experimental
 
#endif