Skip to content

GitLab

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

Merge branch 'feature/timestep-methods' into 'master' 

Feature/timestep methods

See merge request !2296
parents d3208f29 4afa8261
Hide whitespace changes
Inline Side-by-side
dumux/common/pdesolver.hh
View file @ c7192d51
......@@ -64,16 +64,19 @@ using AssemblerVariables = typename VariablesChooser<Assembler>::Type;
* and has a method solve that linearizes (if not already linear), assembles, solves and updates
* given an initial solution producing a new solution.
*
* \tparam Assembler A PDE linearized system assembler
* \tparam LinearSolver A linear system solver
* \tparam A Assembler for linearized system of the PDE
* \tparam LS Linear system solver
*/
template<class Assembler, class LinearSolver>
template<class A, class LS>
class PDESolver
{
using Scalar = typename Assembler::Scalar;
using Scalar = typename A::Scalar;
using TimeLoop = TimeLoopBase<Scalar>;
public:
//! export the assembler and linear solver types
using Assembler = A;
using LinearSolver = LS;
//! export the type of variables that represent a numerical solution
using Variables = Detail::AssemblerVariables<Assembler>;
......@@ -112,8 +115,6 @@ public:
solve(vars);
}
protected:
/*!
* \brief Access the assembler
*/
......@@ -132,6 +133,8 @@ protected:
const LinearSolver& linearSolver() const
{ return *linearSolver_; }
protected:
/*!
* \brief Access the linear solver
*/
......
dumux/timestepping/CMakeLists.txt
View file @ c7192d51
install(FILES
timelevel.hh
multistagemethods.hh
multistagetimestepper.hh
DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}/dumux/timestepping)
dumux/timestepping/multistagemethods.hh 0 → 100644
View file @ c7192d51
// -*- 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
dumux/timestepping/multistagetimestepper.hh 0 → 100644
View file @ c7192d51
// -*- 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
test/CMakeLists.txt
View file @ c7192d51
......@@ -11,6 +11,7 @@ add_subdirectory(nonlinear)
add_subdirectory(porenetwork)
add_subdirectory(porousmediumflow)
add_subdirectory(discretization)
add_subdirectory(timestepping)
# if Python bindings are enabled, include Python binding tests
if(DUNE_ENABLE_PYTHONBINDINGS)
......
test/timestepping/CMakeLists.txt 0 → 100644
View file @ c7192d51
dumux_add_test(SOURCES test_timestepmethods.cc LABELS unit timestepping)
test/timestepping/test_timestepmethods.cc 0 → 100644
View file @ c7192d51
#include <config.h>
#include <cmath>
#include <cassert>
#include <iostream>
#include <algorithm>
#include <utility>
#include <memory>
#include <array>
#include <dune/common/float_cmp.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/parallel/mpihelper.hh>
#include <dumux/io/format.hh>
#include <dumux/common/variables.hh>
#include <dumux/nonlinear/newtonsolver.hh>
#include <dumux/timestepping/timelevel.hh>
#include <dumux/timestepping/multistagemethods.hh>
#include <dumux/timestepping/multistagetimestepper.hh>
/*
This tests the time integration methods by solving the
linear ODE du/dt = exp(t), where u is the unknown
and t is the time. We use the initial condition u_0 = 1,
and thus, the exact solution is u_e = exp(t) - 1.
*/
namespace Dumux {
class ScalarAssembler
{
public:
using Scalar = double;
using SolutionVector = Scalar;
using ResidualType = Scalar;
using JacobianMatrix = Scalar;
using Variables = Experimental::Variables<SolutionVector>;
using StageParams = Experimental::MultiStageParams<Scalar>;
void setLinearSystem() {}
JacobianMatrix& jacobian() { return jac_; }
ResidualType& residual() { return res_; }
void assembleResidual(const Variables& vars)
{
res_ = 0.0;
const auto storage = [] (const auto& stageVars)
{ return stageVars.dofs(); };
const auto source = [] (const auto& stageVars)
{ using std::exp; return exp(stageVars.timeLevel().current()); };
for (std::size_t k = 0; k < stageParams_->size(); ++k)
{
const bool isCurrent = k == stageParams_->size() - 1;
const auto& curVars = isCurrent ? vars : prevStageVariables_[k];
if (!stageParams_->skipTemporal(k))
res_ += stageParams_->temporalWeight(k)*storage(curVars);
if (!stageParams_->skipSpatial(k))
res_ -= stageParams_->spatialWeight(k)*source(curVars);
}
}
void assembleJacobianAndResidual(const Variables& vars)
{
assembleResidual(vars);
jac_ = 1.0;
}
void prepareStage(Variables& variables,
std::shared_ptr<const StageParams> params)
{
stageParams_ = params;
const auto curStage = params->size() - 1;
// we keep track of previous stages, they are needed for residual assembly
prevStageVariables_.push_back(variables);
// Now we update the time level of the given grid variables
const auto t = params->timeAtStage(curStage);
const auto prevT = params->timeAtStage(0);
const auto dtFraction = params->timeStepFraction(curStage);
variables.updateTime(Experimental::TimeLevel{t, prevT, dtFraction});
}
void clearStages()
{
prevStageVariables_.clear();
stageParams_.reset();
}
private:
ResidualType res_;
JacobianMatrix jac_;
std::vector<Variables> prevStageVariables_;
std::shared_ptr<const StageParams> stageParams_;
};
class ScalarLinearSolver
{
public:
void setResidualReduction(double residualReduction) {}
template<class Vector>
bool solve(const double& A, Vector& x, const Vector& b) const
{ x[0] = b[0]/A; return true; }
double norm(const double residual) const
{
using std::abs;
return abs(residual);
}
};
} // end namespace Dumux
int main(int argc, char* argv[])
{
using namespace Dumux;
// maybe initialize MPI
Dune::MPIHelper::instance(argc, argv);
// initialize parameters
// TODO Try to remove this once the Newton does not depend on global default parameters anymore (#1003)
Parameters::init(argc, argv);
using Assembler = ScalarAssembler;
using LinearSolver = ScalarLinearSolver;
using NewtonSolver = NewtonSolver<Assembler, LinearSolver, DefaultPartialReassembler>;
auto assembler = std::make_shared<Assembler>();
auto linearSolver = std::make_shared<LinearSolver>();
auto newtonSolver = std::make_shared<NewtonSolver>(assembler, linearSolver);
// initial solution and variables
using Scalar = typename Assembler::Scalar;
using Variables = typename Assembler::Variables;
using SolutionVector = typename Variables::SolutionVector;
const auto exact = [] (const Scalar t)
{
using std::exp;
return exp(t) - 1.0;
};
const auto computeError = [&] (const Variables& vars)
{
using std::abs;
const auto time = vars.timeLevel().current();
const auto exactSol = exact(time);
const auto absErr = abs(vars.dofs()-exactSol);
const auto relErr = absErr/exactSol;
return std::make_pair(absErr, relErr);
};
const auto testIntegration = [&] (auto method, Scalar expectedRelError)
{
std::cout << "\n-- Integration with " << method->name() << ":\n\n";
SolutionVector x = 0.0;
Variables vars(x);
using TimeStepper = Experimental::MultiStageTimeStepper<NewtonSolver>;
TimeStepper timeStepper(newtonSolver, method);
const Scalar dt = 0.01;
timeStepper.step(vars, /*time*/0.0, dt);
const auto [abs, rel] = computeError(vars);
std::cout << "\n"