diff --git a/dumux/experimental/timestepping/multistagemethods.hh b/dumux/experimental/timestepping/multistagemethods.hh index e8b1f39359a438ac0f1f31da26b6f34b705e3b5f..3533d7f7b2ba1e5935c1caa2bf0e12f2b185101f 100644 --- a/dumux/experimental/timestepping/multistagemethods.hh +++ b/dumux/experimental/timestepping/multistagemethods.hh @@ -3,6 +3,13 @@ // // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder // SPDX-License-Identifier: GPL-3.0-or-later +// SPDX-FileCopyrightInfo: 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. // /*! * \file @@ -188,6 +195,65 @@ private: std::array<Scalar, 5> paramD_; }; +/*! + * \brief Third order DIRK scheme + * \note see Alexander (2003) https://doi.org/10.1016/S0168-9274(03)00012-6 + */ +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_; +}; + } // end namespace MultiStage } // end namespace Dumux::Experimental diff --git a/test/experimental/timestepping/test_timestepmethods.cc b/test/experimental/timestepping/test_timestepmethods.cc index 41add0a105ad07c7f258232d08a20f35bb1bef51..b70da9e692c006c328388f4b9f40e71c27106811 100644 --- a/test/experimental/timestepping/test_timestepmethods.cc +++ b/test/experimental/timestepping/test_timestepmethods.cc @@ -188,6 +188,7 @@ int main(int argc, char* argv[]) testIntegration(std::make_shared<ExplicitEuler<Scalar>>(), 4.9917e-03); testIntegration(std::make_shared<ImplicitEuler<Scalar>>(), 5.0083e-03); testIntegration(std::make_shared<Theta<Scalar>>(0.5), 8.3333e-06); + testIntegration(std::make_shared<DIRKThirdOrderAlexander<Scalar>>(), 7.9214e-09); testIntegration(std::make_shared<RungeKuttaExplicitFourthOrder<Scalar>>(), 3.4829e-12); return 0;