Integrate using a while loop. More...

Inheritance diagram for dg::AdaptiveTimeloop< ContainerType >:

Public Types
using	value_type = dg::get_value_type< ContainerType >

using	container_type = ContainerType

Public Types inherited from dg::aTimeloop< ContainerType >
using	value_type = dg::get_value_type< ContainerType >

using	container_type = ContainerType

Public Member Functions
	AdaptiveTimeloop ()=default
	no allocation More...

	AdaptiveTimeloop (std::function< void(value_type, const ContainerType &, value_type &, ContainerType &, value_type &)> step)
	Construct using a `std::function`. More...

template<class Adaptive , class ODE , class ErrorNorm = value_type( const container_type&), class ControlFunction = value_type(std::array<value_type,3>, std::array<value_type,3>, unsigned, unsigned)>
	AdaptiveTimeloop (Adaptive &&adapt, ODE &&ode, ControlFunction control, ErrorNorm norm, value_type rtol, value_type atol, value_type reject_limit=2)
	Bind the step function of a `dg::Adaptive` object. More...

template<class ... Params>
void	construct (Params &&...ps)
	Perfect forward parameters to one of the constructors. More...

void	set_dt (value_type dt)
	Set initial time-guess in integrate function. More...

template<class Domain >
void	integrate_in_domain (value_type t0, const ContainerType &u0, value_type &t1, ContainerType &u1, value_type dt, Domain &&domain, value_type eps_root)
	Integrate a differential equation within an integration domain. More...

virtual AdaptiveTimeloop *	clone () const
	Abstract copy method that returns a copy of *this on the heap. More...

Public Member Functions inherited from dg::aTimeloop< ContainerType >
void	integrate (value_type t0, const ContainerType &u0, value_type t1, ContainerType &u1)
	Integrate a differential equation between given bounds. More...

void	integrate (value_type &t0, const ContainerType &u0, value_type t1, ContainerType &u1, enum to mode)
	Build your own timeloop. More...

value_type	get_dt () const
	The current timestep. More...

virtual aTimeloop *	clone () const =0
	Abstract copy method that returns a copy of *this on the heap. More...

virtual	~aTimeloop ()

Additional Inherited Members
Protected Member Functions inherited from dg::aTimeloop< ContainerType >
	aTimeloop ()
	empty More...

	aTimeloop (const aTimeloop &)
	empty More...

aTimeloop &	operator= (const aTimeloop &)
	return *this More...

Detailed Description

template<class ContainerType>
struct dg::AdaptiveTimeloop< ContainerType >

Integrate using a while loop.

well suited for dg::Adaptive. The timeloop (for positive dt) corresponds to

t = t0; u1 = u0;
double dt = 1e-6; // some initial guess
while( t < t1)
{
    if( t + dt > t1)
        dt = t1 - t;
    step( t, u1, t, u1, dt);
}

In the integrate_at_least function the if statement is removed.

Note: the current timestep dt is saved by the class and re-used in the next call to integrate unless overwritten by set_dt.

Attention: The integrator may throw if it detects too small timesteps, too many failures, NaN, Inf, or other non-sanitary behaviour

See also: SinglestepTimeloop, MultistepTimeloop

Template Parameters

ContainerType

Any class for which a specialization of TensorTraits exists and which fulfills the requirements of the there defined data and execution policies derived from AnyVectorTag and AnyPolicyTag. Among others

dg::HVec (serial), dg::DVec (cuda / omp), dg::MHVec (mpi + serial) or dg::MDVec (mpi + cuda / omp)
std::vector<dg::DVec> (vector of shared device vectors), std::array<double, 4> (array of 4 doubles) or std::map < std::string, dg::DVec> ( a map of named vectors)
double (scalar) and other primitive types ...

If there are several ContainerTypes in the argument list, then TensorTraits must exist for all of them

See also: See The dg dispatch system for a detailed explanation of our type dispatch system

Member Typedef Documentation

◆ container_type

template<class ContainerType >

using dg::AdaptiveTimeloop< ContainerType >::container_type = ContainerType

◆ value_type

template<class ContainerType >

using dg::AdaptiveTimeloop< ContainerType >::value_type = dg::get_value_type<ContainerType>

Constructor & Destructor Documentation

◆ AdaptiveTimeloop() [1/3]

template<class ContainerType >

dg::AdaptiveTimeloop< ContainerType >::AdaptiveTimeloop ( )

default

no allocation

◆ AdaptiveTimeloop() [2/3]

template<class ContainerType >

dg::AdaptiveTimeloop< ContainerType >::AdaptiveTimeloop ( std::function< void(value_type, const ContainerType &, value_type &, ContainerType &, value_type &)> step )

inline

Construct using a std::function.

Parameters

step	Called in the timeloop as `step( t0, u1, t0, u1, dt)` . Has to advance the ode in-place by `dt` and suggest a new `dt` for the next step.

◆ AdaptiveTimeloop() [3/3]

template<class ContainerType >

template<class Adaptive , class ODE , class ErrorNorm = value_type( const container_type&), class ControlFunction = value_type(std::array<value_type,3>, std::array<value_type,3>, unsigned, unsigned)>

dg::AdaptiveTimeloop< ContainerType >::AdaptiveTimeloop	(	Adaptive &&	adapt,
		ODE &&	ode,
		ControlFunction	control,
		ErrorNorm	norm,
		value_type	rtol,
		value_type	atol,
		value_type	reject_limit = `2`
	)

inline

Bind the step function of a dg::Adaptive object.

Construct a lambda function that calls the step function of adapt with given parameters and stores it internally in a std::function

Template Parameters

Adaptive a type with a step function as in dg::Adaptive

Parameters

adapt If constructed in-place (rvalue), will be copied into the lambda. If an lvalue, then the lambda stores a reference

Attention: If adapt is an lvalue then you need to make sure that adapt remains valid to avoid a "dangling reference"

Template Parameters

ODE	The ExplicitRHS or tuple type that corresponds to what is inserted into the step member of the Stepper

Parameters

ode	rhs or tuple

Template Parameters

ControlFunction function or Functor called as dt' = control( {dt0, dt1, dt2}, {eps0, eps1, eps2}, embedded_order, order), where all parameters are of type value_type except the last two, which are unsigned.

Parameters

control The control function. For explicit and imex methods, dg::pid_control is a good choice with dg::ex_control or dg::imex_control as an alternative if too many steps fail. For implicit methods use the dg::im_control. The task of the control function is to compute a new timestep size based on the old timestep size, the order of the method and the past error(s). The behaviour of the controller is also changed by the set_reject_limit function

Template Parameters

ErrorNorm function or Functor of type value_type( const ContainerType&)

Parameters

norm	The error norm. Usually `dg::l2norm` is a good choice, but for very small vector sizes the time for the binary reproducible dot product might become a performance bottleneck. Then `dg::fast_l2norm` is a better choice.
rtol	the desired relative accuracy. Usually 1e-5 is a good choice.
atol	the desired absolute accuracy. Usually 1e-7 is a good choice.
reject_limit	the default value is 2. Sometimes even 2 is not enough, for example when the stepsize fluctuates very much and the stepper fails often then it may help to increase the reject_limit further (to 10 say).

Note: The error tolerance is computed such that on average every point in the solution vector fulfills \( |\delta_i| \approx r|u_i| + a\); Try not to mess with dt. The controller is best left alone and it does a very good job choosing timesteps. But how do I output my solution at certain (equidistant) timesteps? First, think about if you really, really need that. Why is it so bad to have output at non-equidistant timesteps? If you are still firm, then consider using an interpolation scheme (cf. dg::Extrapolation). Let choosing the timestep yourself be the very last option if the others are not viable; From Kennedy and Carpenter, Appl. num. Math., (2003): "Step-size control is a means by which accuracy, iteration, and to a lesser extent stability are controlled. The choice of (t) may be chosen from many criteria, among those are the (t) from the accuracy based step controller, the (t)_inviscid and (t)_viscous associated with the inviscid and viscous stability limits of the ERK, and the (t)_iter associated with iteration convergence."

Our control method is an error based control method. However, due to the CFL (or other stability) condition there might be a sharp barrier in the range of possible stepsizes, i.e the stability limited timestep where the error sharply increses leading to a large number of rejected steps. The pid-controller usually does a good job keeping the timestep "just right" but in future work we may want to consider stability based control methods as well. Furthermore, for implicit or semi-implicit methods the chosen timestep influences how fast the iterative solver converges. If the timestep is too large the solver might take too long to converge and a smaller timestep might be preferable. Even for explicit methods where in each timestep an elliptic equation has to be solved the timestep may influence the convergence rate. Currently, there is no communication implemented between these solvers and the controller, but we may consider it in future releases.

Attention: When you use the ARKStep in combination with the Adaptive time step algorithm pay attention to solve the implicit part with sufficient accuracy. Else, the error propagates into the time controller, which will then choose the timestep as if the implicit part was explicit i.e. far too small. This might have to do with stiffness-leakage [Kennedy and Carpenter, Appl. num. Math., (2003)]: "An essential requirement for the viability of stiff/nonstiff IMEX schemes is that the stiffness remains truely separable. If this were not the case then stiffness would leak out of the stiff terms and stiffen the nonstiff terms. It would manifest itself as a loss in stability or a forced reduction in stepsize of the nonstiff terms. A more expensive fully implicit approach might then be required, and hence, methods that leak substantial stiffness might best be avoided".

Note: The exact values of rtol and atol might not be too important. However, don't make rtol too large, 1e-1 say, since then the controller might get too close to the CFL barrier. The timestepper is still able to crash, mind, even though the chances of that happening are somewhat lower than in a fixed stepsize method.

Member Function Documentation

◆ clone()

template<class ContainerType >

virtual AdaptiveTimeloop * dg::AdaptiveTimeloop< ContainerType >::clone ( ) const

inlinevirtual

Abstract copy method that returns a copy of *this on the heap.

Returns: a copy of *this on the heap

See also: dg::ClonePtr

Implements dg::aTimeloop< ContainerType >.

◆ construct()

template<class ContainerType >

template<class ... Params>

void dg::AdaptiveTimeloop< ContainerType >::construct ( Params &&... ps )

inline

Perfect forward parameters to one of the constructors.

Template Parameters

Params deduced by the compiler

Parameters

ps	parameters forwarded to constructors

◆ integrate_in_domain()

template<class ContainerType >

template<class Domain >

void dg::AdaptiveTimeloop< ContainerType >::integrate_in_domain	(	value_type	t0,
		const ContainerType &	u0,
		value_type &	t1,
		ContainerType &	u1,
		value_type	dt,
		Domain &&	domain,
		value_type	eps_root
	)

Integrate a differential equation within an integration domain.

If the integration does not leave the domain then it is equivalent to a call to integrate, else a bisection algorithm is used to find the exact point of exit

Parameters

t0	initial time
u0	initial value at `t0`
t1	(read / write) end time; if the solution leaves the domain contains the last time where the solution still lies within the domain on output
u1	(write only) contains the result corresponding to t1 on output Use only if you know a good step-size.
dt	The initial timestep guess (if 0 the function chooses something for you). The exact value is not really important, the stepper does not even have to succeed. Usually the control function will very(!) quickly adapt the stepsize in just one or two steps (even if it's several orders of magnitude off in the beginning).
domain	a restriction of the solution space. The integrator checks after every step if the solution is still within the given domain `domain.contains(u1)`. If not, the integrator will bisect the exact domain boundary (up to the given tolerances) and return (t1, u1) that lies closest (but within) the domain boundary.
eps_root	Relative error of root finding algorithm `dt` < eps( t1 + 1)

Template Parameters

Domain Must have the contains(const ContainerType&) const member function returning true if the given solution is part of the domain, false else (can for example be dg::aRealTopology2d)

◆ set_dt()

template<class ContainerType >

void dg::AdaptiveTimeloop< ContainerType >::set_dt ( value_type dt )

inline

Set initial time-guess in integrate function.

Use only if you know a good step-size.

Parameters

dt The initial timestep guess (if 0 the function chooses something for you). The exact value is not really important, the stepper does not even have to succeed. Usually the control function will very(!) quickly adapt the stepsize in just one or two steps (even if it's several orders of magnitude off in the beginning).

The documentation for this struct was generated from the following file:

adaptive.h

Public Types

Public Member Functions

Additional Inherited Members

Detailed Description

Member Typedef Documentation

◆ container_type

◆ value_type

Constructor & Destructor Documentation

◆ AdaptiveTimeloop() [1/3]

◆ AdaptiveTimeloop() [2/3]

◆ AdaptiveTimeloop() [3/3]

Member Function Documentation

◆ clone()

◆ construct()

◆ integrate_in_domain()

◆ set_dt()