Commit 8bf983ed authored by Dennis Gläser's avatar Dennis Gläser
Browse files

[timeintegration] Add a multistage time stepping scheme

Based on initial draft by Timo Koch.
parent d3208f29
install(FILES
timelevel.hh
multistagemethods.hh
multistagetimestepper.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/timestepping)
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief Parameters for different multistage time stepping methods
* \note See e.g. https://en.wikipedia.org/wiki/List_of_Runge%E2%80%93Kutta_methods
*/
#ifndef DUMUX_TIMESTEPPING_MULTISTAGE_METHODS_HH
#define DUMUX_TIMESTEPPING_MULTISTAGE_METHODS_HH
#include <cmath>
#include <string>
#include <array>
namespace Dumux::Experimental {
/*!
* \brief Abstract interface for one-step multi-stage method parameters in Shu/Osher form.
*
* This implementation is based on the Shu/Osher form from:
* Chi W. Shu and Stanley Osher. Efficient implementation of essentially
* non- oscillatory shock-capturing schemes. J. Comput. Phys., 77:439–471, 1988.
* https://doi.org/10.1016/0021-9991(88)90177-5.
* To this end Eq. (2.12) is extended for implicit schemes.
*
* We consider the general PDE form
*
* \f[
* \begin{equation}
* \frac{\partial M(x)}{\partial t} - R(x, t) = 0,
* \end{equation}
* \f]
*
* where \f$ M(x)\f$ is the temporal operator/residual in terms of the solution \f$ x \f$,
* and \f$ R(x)\f$ is the discrete spatial operator/residual.
* \f$ M(x)\f$ usually corresponds to the conserved quanity (e.g. mass), and \f$ R(x)\f$
* contains the rest of the residual. We can then construct \f$ m \f$-stage time discretization methods.
* Integrating from time \f$ t^n\f$ to \f$ t^{n+1}\f$ with time step size \f$ \Delta t^n\f$, we solve:
*
* \f[
* \begin{aligned}
* x^{(0)} &= u^n\\
* \sum_{k=0}^i \left[ \alpha_{ik} M\left(x^{(k)}, t^n + d_k\Delta t^n\right)
* + \beta_{ik}\Delta t^n R \left(x^{(k)}, t^n+d_k\Delta t^n \right)\right] &= 0 & \forall i \in \{1,\ldots,m\} \\
* x^{n+1} &= x^{(m)}
* \end{aligned}
* \f]
* where \f$ x^{(k)} \f$ denotes the intermediate solution at stage \f$ k \f$.
* Dependent on the number of stages \f$ m \f$, and the coefficients \f$ \alpha, \beta, d\f$,
* schemes with different properties and order of accuracy can be constructed.
*/
template<class Scalar>
class MultiStageMethod
{
public:
virtual bool implicit () const = 0;
virtual std::size_t numStages () const = 0;
//! weights of the temporal operator residual (\f$ \alpha_{ik} \f$)
virtual Scalar temporalWeight (std::size_t i, std::size_t k) const = 0;
//! weights of the spatial operator residual (\f$ \beta_{ik} \f$)
virtual Scalar spatialWeight (std::size_t i, std::size_t k) const = 0;
//! time step weights for each stage (\f$ d_k \f$)
virtual Scalar timeStepWeight (std::size_t k) const = 0;
virtual std::string name () const = 0;
virtual ~MultiStageMethod() = default;
};
//! Multi-stage time stepping scheme implementations
namespace MultiStage {
/*!
* \brief A theta time stepping scheme
* theta=1.0 is an implicit Euler scheme,
* theta=0.0 an explicit Euler scheme,
* theta=0.5 is a Cranck-Nicholson scheme
*/
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_;
};
/*!
* \brief An explicit Euler time stepping scheme
*/
template<class Scalar>
class ExplicitEuler final : public Theta<Scalar>
{
public:
ExplicitEuler() : Theta<Scalar>(0.0) {}
std::string name () const final
{ return "explicit Euler"; }
};
/*!
* \brief An implicit Euler time stepping scheme
*/
template<class Scalar>
class ImplicitEuler final : public Theta<Scalar>
{
public:
ImplicitEuler() : Theta<Scalar>(1.0) {}
std::string name () const final
{ return "implicit Euler"; }
};
/*!
* \brief Classical explicit fourth order Runge-Kutta scheme
*/
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_;
};
} // end namespace MultiStage
} // end namespace Dumux::Experimental
#endif
// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
// vi: set et ts=4 sw=4 sts=4:
/*****************************************************************************
* See the file COPYING for full copying permissions. *
* *
* This program is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* This program is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with this program. If not, see <http://www.gnu.org/licenses/>. *
*****************************************************************************/
/*!
* \file
* \brief A time stepper performing a single time step of a transient simulation
*/
#ifndef DUMUX_TIMESTEPPING_MULTISTAGE_TIMESTEPPER_HH
#define DUMUX_TIMESTEPPING_MULTISTAGE_TIMESTEPPER_HH
#include <memory>
#include <vector>
#include <cmath>
namespace Dumux::Experimental {
//! forward declaration
template<class Scalar>
class MultiStageMethod;
//! Data object for the parameters of a given stage
template<class Scalar>
class MultiStageParams
{
struct Params {
Scalar alpha, betaDt, timeAtStage, dtFraction;
bool skipTemporal, skipSpatial;
};
public:
//! Extract params for stage i from method m
MultiStageParams(const MultiStageMethod<Scalar>& m, std::size_t i, const Scalar t, const Scalar dt)
: size_(i+1)
{
params_.resize(size_);
for (std::size_t k = 0; k < size_; ++k)
{
auto& p = params_[k];
p.alpha = m.temporalWeight(i, k);
p.betaDt = m.spatialWeight(i, k)*dt;
p.timeAtStage = t + m.timeStepWeight(k)*dt;
p.dtFraction = m.timeStepWeight(k);
using std::abs;
p.skipTemporal = (abs(p.alpha) < 1e-6);
p.skipSpatial = (abs(p.betaDt) < 1e-6);
}
}
std::size_t size () const
{ return size_; }
//! weights of the temporal operator residual (\f$ \alpha_{ik} \f$)
Scalar temporalWeight (std::size_t k) const
{ return params_[k].alpha; }
//! weights of the spatial operator residual (\f$ \beta_{ik} \f$)
Scalar spatialWeight (std::size_t k) const
{ return params_[k].betaDt; }
//! the time at which we have to evaluate the operators
Scalar timeAtStage (std::size_t k) const
{ return params_[k].timeAtStage; }
//! the fraction of a time step corresponding to the k-th stage
Scalar timeStepFraction (std::size_t k) const
{ return params_[k].dtFraction; }
//! If \f$ \alpha_{ik} = 0\f$
Scalar skipTemporal (std::size_t k) const
{ return params_[k].skipTemporal; }
//! If \f$ \beta_{ik} = 0\f$
Scalar skipSpatial (std::size_t k) const
{ return params_[k].skipSpatial; }
private:
std::size_t size_;
std::vector<Params> params_;
};
/*!
* \brief Time stepping with a multi-stage method
* \note We limit ourselves to "diagonally" implicit multi-stage methods where solving
* a stage can only depend on the values of the same stage and stages before
* but not future stages (which would require solving larger linear systems)
*/
template<class PDESolver>
class MultiStageTimeStepper
{
using Variables = typename PDESolver::Variables;
using Scalar = typename Variables::Scalar;
using StageParams = MultiStageParams<Scalar>;
public:
/*!
* \brief The constructor
* \param pdeSolver Solver class for solving a PDE in each stage
* \param msMethod The multi-stage method which is to be used for time integration
* \todo TODO: Add time step control if the pde solver doesn't converge
*/
MultiStageTimeStepper(std::shared_ptr<PDESolver> pdeSolver,
std::shared_ptr<const MultiStageMethod<Scalar>> msMethod)
: pdeSolver_(pdeSolver)
, msMethod_(msMethod)
{}
/*!
* \brief Advance one time step of the given time loop
* \params vars The variables object at the current time level.
* \param t The current time level
* \param dt The time step size to be used
* \note We expect the time level in vars to correspond to the given time `t`
* \todo: TODO: Add time step control if the pde solver doesn't converge
*/
void step(Variables& vars, const Scalar t, const Scalar dt)
{
// make sure there are no traces of previous stages
pdeSolver_->assembler().clearStages();
for (auto stageIdx = 1UL; stageIdx <= msMethod_->numStages(); ++stageIdx)
{
// extract parameters for this stage from the time stepping method
auto stageParams = std::make_shared<StageParams>(*msMethod_, stageIdx, t, dt);
// prepare the assembler for this stage
pdeSolver_->assembler().prepareStage(vars, stageParams);
// assemble & solve
pdeSolver_->solve(vars);
}
// clear traces of previously registered stages
pdeSolver_->assembler().clearStages();
}
/*!
* \brief Set/change the time step method
*/
void setMethod(std::shared_ptr<const MultiStageMethod<Scalar>> msMethod)
{ msMethod_ = msMethod; }
private:
std::shared_ptr<PDESolver> pdeSolver_;
std::shared_ptr<const MultiStageMethod<Scalar>> msMethod_;
};
} // end namespace Dumux::Experimental
#endif
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment