diff --git a/gusto/time_discretisation.py b/gusto/time_discretisation.py
index ec6428228..4a0dcd4d6 100644
--- a/gusto/time_discretisation.py
+++ b/gusto/time_discretisation.py
@@ -7,7 +7,7 @@
from abc import ABCMeta, abstractmethod, abstractproperty
from firedrake import (Function, TestFunction, NonlinearVariationalProblem,
- NonlinearVariationalSolver, DirichletBC, Constant)
+ NonlinearVariationalSolver, DirichletBC, split, Constant)
from firedrake.formmanipulation import split_form
from firedrake.utils import cached_property
@@ -22,9 +22,9 @@
import numpy as np
-__all__ = ["ForwardEuler", "BackwardEuler", "ExplicitMultistage", "SSPRK3", "RK4",
- "Heun", "ThetaMethod", "TrapeziumRule", "BDF2", "TR_BDF2", "Leapfrog",
- "AdamsMoulton", "AdamsBashforth"]
+__all__ = ["ForwardEuler", "BackwardEuler", "ExplicitMultistage", "ImplicitMultistage",
+ "SSPRK3", "RK4", "Heun", "ThetaMethod", "TrapeziumRule", "BDF2", "TR_BDF2",
+ "Leapfrog", "AdamsMoulton", "AdamsBashforth", "ImplicitMidpoint", "QinZhang"]
def wrapper_apply(original_apply):
@@ -243,6 +243,144 @@ def apply(self, x_out, x_in):
pass
+class ImplicitMultistage(TimeDiscretisation):
+ """
+ A class for implementing general diagonally implicit multistage (Runge-Kutta)
+ methods based on its Butcher tableau.
+
+ Unlike the explicit method, all upper diagonal a_ij elements are non-zero for implicit methods.
+
+ There are three steps to move from the current solution, y^n, to the new one, y^{n+1}
+
+ For each i = 1, s in an s stage method
+ we have the intermediate solutions: \n
+ y_i = y^n + dt*(a_i1*k_1 + a_i2*k_2 + ... + a_ii*k_i) \n
+ We compute the gradient at the intermediate location, k_i = F(y_i) \n
+
+ At the last stage, compute the new solution by: \n
+ y^{n+1} = y^n + dt*(b_1*k_1 + b_2*k_2 + .... + b_s*k_s) \n
+
+ """
+ # ---------------------------------------------------------------------------
+ # Butcher tableau for a s-th order
+ # diagonally implicit scheme:
+ # c_0 | a_00 0 . 0
+ # c_1 | a_10 a_11 . 0
+ # . | . . . .
+ # . | . . . .
+ # c_s | a_s0 a_s1 . a_ss
+ # -------------------------
+ # | b_1 b_2 ... b_s
+ #
+ #
+ # The corresponding square 'butcher_matrix' is:
+ #
+ # [a_00 0 . 0 ]
+ # [a_10 a_11 . 0 ]
+ # [ . . . . ]
+ # [ b_0 b_1 . b_s]
+ # ---------------------------------------------------------------------------
+
+ def __init__(self, domain, butcher_matrix, field_name=None,
+ solver_parameters=None, limiter=None, options=None,):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ butcher_matrix (numpy array): A matrix containing the coefficients of
+ a butcher tableau defining a given Runge Kutta time discretisation.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
+ super().__init__(domain, field_name=field_name,
+ solver_parameters=solver_parameters,
+ limiter=limiter, options=options)
+ self.butcher_matrix = butcher_matrix
+ self.nStages = int(np.shape(self.butcher_matrix)[1])
+
+ def setup(self, equation, apply_bcs=True, *active_labels):
+ """
+ Set up the time discretisation based on the equation.
+
+ Args:
+ equation (:class:`PrognosticEquation`): the model's equation.
+ *active_labels (:class:`Label`): labels indicating which terms of
+ the equation to include.
+ """
+
+ super().setup(equation, apply_bcs, *active_labels)
+
+ self.k = [Function(self.fs) for i in range(self.nStages)]
+
+ def lhs(self):
+ return super().lhs
+
+ def rhs(self):
+ return super().rhs
+
+ def solver(self, stage):
+ residual = self.residual.label_map(
+ lambda t: t.has_label(time_derivative),
+ map_if_true=drop,
+ map_if_false=replace_subject(self.xnph, self.idx),
+ )
+ mass_form = self.residual.label_map(
+ lambda t: t.has_label(time_derivative),
+ map_if_false=drop)
+ residual += mass_form.label_map(all_terms,
+ replace_subject(self.x_out, self.idx))
+
+ problem = NonlinearVariationalProblem(residual.form, self.x_out, bcs=self.bcs)
+
+ solver_name = self.field_name+self.__class__.__name__ + "%s" % (stage)
+ return NonlinearVariationalSolver(problem, solver_parameters=self.solver_parameters,
+ options_prefix=solver_name)
+
+ @cached_property
+ def solvers(self):
+ solvers = []
+ for stage in range(self.nStages):
+ solvers.append(self.solver(stage))
+ return solvers
+
+ def solve_stage(self, x0, stage):
+ self.x1.assign(x0)
+ for i in range(stage):
+ self.x1.assign(self.x1 + self.butcher_matrix[stage, i]*self.dt*self.k[i])
+
+ if self.limiter is not None:
+ self.limiter.apply(self.x1)
+
+ if self.idx is None and len(self.fs) > 1:
+ self.xnph = tuple([self.dt*self.butcher_matrix[stage, stage]*a + b for a, b in zip(split(self.x_out), split(self.x1))])
+ else:
+ self.xnph = self.x1 + self.butcher_matrix[stage, stage]*self.dt*self.x_out
+ solver = self.solvers[stage]
+ solver.solve()
+
+ self.k[stage].assign(self.x_out)
+
+ def apply(self, x_out, x_in):
+
+ for i in range(self.nStages):
+ self.solve_stage(x_in, i)
+
+ x_out.assign(x_in)
+ for i in range(self.nStages):
+ x_out.assign(x_out + self.butcher_matrix[self.nStages, i]*self.dt*self.k[i])
+
+ if self.limiter is not None:
+ self.limiter.apply(x_out)
+
+
class ExplicitTimeDiscretisation(TimeDiscretisation):
"""Base class for explicit time discretisations."""
@@ -365,50 +503,78 @@ class ExplicitMultistage(ExplicitTimeDiscretisation):
A class for implementing general explicit multistage (Runge-Kutta)
methods based on its Butcher tableau.
- A Butcher tableau is formed in the following way for a s-th order explicit scheme:
-
- c_0 | a_00
- c_1 | a_10 a_11
- . | a_20 a_21 a_22
- . | . . . .
- -------------------------
- c_s | b_1 b_2 ... b_s
+ A Butcher tableau is formed in the following way for a s-th order explicit scheme: \n
- The gradient function has no time-dependence, so c elements are not needed.
- A square 'butcher_matrix' is defined in each class that uses
- the ExplicitMultiStage structure,
-
- [a_00 0 . 0 ]
- [a_10 a_11 . 0 ]
- [ . . . . ]
- [ b_0 b_1 . b_{s-1}]
-
- All upper diagonal a_ij elements are zero for explicit methods.
+ All upper diagonal a_ij elements are zero for explicit methods. We exclude the first
+ row of the butcher tableau from our butcher matrix as the row is always zeros.
There are three steps to move from the current solution, y^n, to the new one, y^{n+1}
- For an s stage method,
- At iteration i (from 0 to s-1)
- An intermediate location is computed as y_i = y^n + sum{over j less than i} (dt*a_ij*k_i)
- Compute the gradient at the intermediate location, k_i = F(y_i)
+ For each i = 1, s in an s stage method
+ we have the intermediate solutions: \n
+ y_i = y^n + dt*(a_i1*k_1 + a_i2*k_2 + ... + a_i{i-1}*k_{i-1}) \n
+ We compute the gradient at the intermediate location, k_i = F(y_i) \n
At the last stage, compute the new solution by:
- y^{n+1} = y^n + sum_{j from 0 to s-1} (dt*b_i*k_i)
+ y^{n+1} = y^n + dt*(b_1*k_1 + b_2*k_2 + .... + b_s*k_s) \n
"""
-
- def __init__(self, domain, field_name=None, fixed_subcycles=None,
+ # ---------------------------------------------------------------------------
+ # Butcher tableau for a s-th order
+ # explicit scheme:
+ # c_0 | 0 0 . 0
+ # c_1 | a_10 0 . 0
+ # . | . . . .
+ # . | . . . .
+ # c_s | a_s0 a_s1 . 0
+ # -------------------------
+ # | b_1 b_2 ... b_s
+ #
+ #
+ # The corresponding square 'butcher_matrix' is:
+ #
+ # [a_10 0 . 0 ]
+ # [a_20 a_21 . 0 ]
+ # [ . . . . ]
+ # [ b_0 b_1 . b_s]
+ # ---------------------------------------------------------------------------
+
+ def __init__(self, domain, butcher_matrix, field_name=None, fixed_subcycles=None,
subcycle_by_courant=None, solver_parameters=None,
- limiter=None, options=None, butcher_matrix=None):
-
+ limiter=None, options=None):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ butcher_matrix (numpy array): A matrix containing the coefficients of
+ a butcher tableau defining a given Runge Kutta time discretisation.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ fixed_subcycles (int, optional): the fixed number of sub-steps to
+ perform. This option cannot be specified with the
+ `subcycle_by_courant` argument. Defaults to None.
+ subcycle_by_courant (float, optional): specifying this option will
+ make the scheme perform adaptive sub-cycling based on the
+ Courant number. The specified argument is the maximum Courant
+ for one sub-cycle. Defaults to None, in which case adaptive
+ sub-cycling is not used. This option cannot be specified with the
+ `fixed_subcycles` argument.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
super().__init__(domain, field_name=field_name,
fixed_subcycles=fixed_subcycles,
subcycle_by_courant=subcycle_by_courant,
solver_parameters=solver_parameters,
limiter=limiter, options=options)
- if butcher_matrix is not None:
- self.butcher_matrix = butcher_matrix
- self.nbutcher = int(np.shape(self.butcher_matrix)[0])
+ self.butcher_matrix = butcher_matrix
+ self.nbutcher = int(np.shape(self.butcher_matrix)[0])
@property
def nStages(self):
@@ -500,44 +666,88 @@ class ForwardEuler(ExplicitMultistage):
"""
Implements the forward Euler timestepping scheme.
- The forward Euler method for operator F is the most simple explicit scheme:
- k0 = F[y^n]
- y^(n+1) = y^n + dt*k0
+ The forward Euler method for operator F is the most simple explicit scheme: \n
+ k0 = F[y^n] \n
+ y^(n+1) = y^n + dt*k0 \n
"""
def __init__(self, domain, field_name=None, fixed_subcycles=None,
subcycle_by_courant=None, solver_parameters=None,
- limiter=None, options=None, butcher_matrix=None):
- super().__init__(domain, field_name=field_name,
+ limiter=None, options=None):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ fixed_subcycles (int, optional): the fixed number of sub-steps to
+ perform. This option cannot be specified with the
+ `subcycle_by_courant` argument. Defaults to None.
+ subcycle_by_courant (float, optional): specifying this option will
+ make the scheme perform adaptive sub-cycling based on the
+ Courant number. The specified argument is the maximum Courant
+ for one sub-cycle. Defaults to None, in which case adaptive
+ sub-cycling is not used. This option cannot be specified with the
+ `fixed_subcycles` argument.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
+ butcher_matrix = np.array([1.]).reshape(1, 1)
+ super().__init__(domain, butcher_matrix, field_name=field_name,
fixed_subcycles=fixed_subcycles,
subcycle_by_courant=subcycle_by_courant,
solver_parameters=solver_parameters,
- limiter=limiter, options=options,
- butcher_matrix=butcher_matrix)
- self.butcher_matrix = np.array([1.]).reshape(1, 1)
- self.nbutcher = int(np.shape(self.butcher_matrix)[0])
+ limiter=limiter, options=options)
class SSPRK3(ExplicitMultistage):
u"""
Implements the 3-stage Strong-Stability-Preserving Runge-Kutta method
- for solving ∂y/∂t = F(y). It can be written as:
+ for solving ∂y/∂t = F(y). It can be written as: \n
- k0 = F[y^n]
- k1 = F[y^n + dt*k1]
- k2 = F[y^n + (1/4)*dt*(k0+k1)]
- y^(n+1) = y^n + (1/6)*dt*(k0 + k1 + 4*k2)
+ k0 = F[y^n] \n
+ k1 = F[y^n + dt*k1] \n
+ k2 = F[y^n + (1/4)*dt*(k0+k1)] \n
+ y^(n+1) = y^n + (1/6)*dt*(k0 + k1 + 4*k2) \n
"""
def __init__(self, domain, field_name=None, fixed_subcycles=None,
subcycle_by_courant=None, solver_parameters=None,
- limiter=None, options=None, butcher_matrix=None):
- super().__init__(domain, field_name=field_name,
+ limiter=None, options=None):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ fixed_subcycles (int, optional): the fixed number of sub-steps to
+ perform. This option cannot be specified with the
+ `subcycle_by_courant` argument. Defaults to None.
+ subcycle_by_courant (float, optional): specifying this option will
+ make the scheme perform adaptive sub-cycling based on the
+ Courant number. The specified argument is the maximum Courant
+ for one sub-cycle. Defaults to None, in which case adaptive
+ sub-cycling is not used. This option cannot be specified with the
+ `fixed_subcycles` argument.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
+ butcher_matrix = np.array([[1., 0., 0.], [1./4., 1./4., 0.], [1./6., 1./6., 2./3.]])
+ super().__init__(domain, butcher_matrix, field_name=field_name,
fixed_subcycles=fixed_subcycles,
subcycle_by_courant=subcycle_by_courant,
solver_parameters=solver_parameters,
- limiter=limiter, options=options,
- butcher_matrix=butcher_matrix)
- self.butcher_matrix = np.array([[1., 0., 0.], [1./4., 1./4., 0.], [1./6., 1./6., 2./3.]])
- self.nbutcher = int(np.shape(self.butcher_matrix)[0])
+ limiter=limiter, options=options)
class RK4(ExplicitMultistage):
@@ -545,27 +755,49 @@ class RK4(ExplicitMultistage):
Implements the classic 4-stage Runge-Kutta method.
The classic 4-stage Runge-Kutta method for solving ∂y/∂t = F(y). It can be
- written as:
+ written as: \n
- k0 = F[y^n]
- k1 = F[y^n + 1/2*dt*k1]
- k2 = F[y^n + 1/2*dt*k2]
- k3 = F[y^n + dt*k3]
- y^(n+1) = y^n + (1/6) * dt * (k0 + 2*k1 + 2*k2 + k3)
+ k0 = F[y^n] \n
+ k1 = F[y^n + 1/2*dt*k1] \n
+ k2 = F[y^n + 1/2*dt*k2] \n
+ k3 = F[y^n + dt*k3] \n
+ y^(n+1) = y^n + (1/6) * dt * (k0 + 2*k1 + 2*k2 + k3) \n
- where superscripts indicate the time-level.
+ where superscripts indicate the time-level. \n
"""
def __init__(self, domain, field_name=None, fixed_subcycles=None,
subcycle_by_courant=None, solver_parameters=None,
- limiter=None, options=None, butcher_matrix=None):
- super().__init__(domain, field_name=field_name,
+ limiter=None, options=None):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ fixed_subcycles (int, optional): the fixed number of sub-steps to
+ perform. This option cannot be specified with the
+ `subcycle_by_courant` argument. Defaults to None.
+ subcycle_by_courant (float, optional): specifying this option will
+ make the scheme perform adaptive sub-cycling based on the
+ Courant number. The specified argument is the maximum Courant
+ for one sub-cycle. Defaults to None, in which case adaptive
+ sub-cycling is not used. This option cannot be specified with the
+ `fixed_subcycles` argument.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
+ butcher_matrix = np.array([[0.5, 0., 0., 0.], [0., 0.5, 0., 0.], [0., 0., 1., 0.], [1./6., 1./3., 1./3., 1./6.]])
+ super().__init__(domain, butcher_matrix, field_name=field_name,
fixed_subcycles=fixed_subcycles,
subcycle_by_courant=subcycle_by_courant,
solver_parameters=solver_parameters,
- limiter=limiter, options=options,
- butcher_matrix=butcher_matrix)
- self.butcher_matrix = np.array([[0.5, 0., 0., 0.], [0., 0.5, 0., 0.], [0., 0., 1., 0.], [1./6., 1./3., 1./3., 1./6.]])
- self.nbutcher = int(np.shape(self.butcher_matrix)[0])
+ limiter=limiter, options=options)
class Heun(ExplicitMultistage):
@@ -573,31 +805,55 @@ class Heun(ExplicitMultistage):
Implements Heun's method.
The 2-stage Runge-Kutta scheme known as Heun's method,for solving
- ∂y/∂t = F(y). It can be written as:
+ ∂y/∂t = F(y). It can be written as: \n
- y_1 = F[y^n]
- y^(n+1) = (1/2)y^n + (1/2)F[y_1]
+ y_1 = F[y^n] \n
+ y^(n+1) = (1/2)y^n + (1/2)F[y_1] \n
where superscripts indicate the time-level and subscripts indicate the stage
number.
"""
def __init__(self, domain, field_name=None, fixed_subcycles=None,
subcycle_by_courant=None, solver_parameters=None,
- limiter=None, options=None, butcher_matrix=None):
- super().__init__(domain, field_name, fixed_subcycles=fixed_subcycles,
+ limiter=None, options=None):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ fixed_subcycles (int, optional): the fixed number of sub-steps to
+ perform. This option cannot be specified with the
+ `subcycle_by_courant` argument. Defaults to None.
+ subcycle_by_courant (float, optional): specifying this option will
+ make the scheme perform adaptive sub-cycling based on the
+ Courant number. The specified argument is the maximum Courant
+ for one sub-cycle. Defaults to None, in which case adaptive
+ sub-cycling is not used. This option cannot be specified with the
+ `fixed_subcycles` argument.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
+ butcher_matrix = np.array([[1., 0.], [0.5, 0.5]])
+ super().__init__(domain, butcher_matrix, field_name=field_name,
+ fixed_subcycles=fixed_subcycles,
subcycle_by_courant=subcycle_by_courant,
solver_parameters=solver_parameters,
limiter=limiter, options=options)
- self.butcher_matrix = np.array([[1., 0.], [0.5, 0.5]])
- self.nbutcher = np.shape(self.butcher_matrix)[0]
class BackwardEuler(TimeDiscretisation):
"""
Implements the backward Euler timestepping scheme.
- The backward Euler method for operator F is the most simple implicit scheme:
- y^(n+1) = y^n + dt*F[y^(n+1)].
+ The backward Euler method for operator F is the most simple implicit scheme: \n
+ y^(n+1) = y^n + dt*F[y^(n+1)]. \n
"""
def __init__(self, domain, field_name=None, solver_parameters=None,
limiter=None, options=None):
@@ -672,9 +928,9 @@ class ThetaMethod(TimeDiscretisation):
Implements the theta implicit-explicit timestepping method, which can
be thought as a generalised trapezium rule.
- The theta implicit-explicit timestepping method for operator F is written as
- y^(n+1) = y^n + dt*(1-theta)*F[y^n] + dt*theta*F[y^(n+1)]
- for off-centring parameter theta.
+ The theta implicit-explicit timestepping method for operator F is written as: \n
+ y^(n+1) = y^n + dt*(1-theta)*F[y^n] + dt*theta*F[y^(n+1)] \n
+ for off-centring parameter theta. \n
"""
def __init__(self, domain, theta, field_name=None,
@@ -754,9 +1010,9 @@ class TrapeziumRule(ThetaMethod):
Implements the trapezium rule timestepping method, also commonly known as
Crank Nicholson.
- The trapezium rule timestepping method for operator F is written as
- y^(n+1) = y^n + dt/2*F[y^n] + dt/2*F[y^(n+1)].
- It is equivalent to the "theta" method with theta = 1/2.
+ The trapezium rule timestepping method for operator F is written as: \n
+ y^(n+1) = y^n + dt/2*F[y^n] + dt/2*F[y^(n+1)]. \n
+ It is equivalent to the "theta" method with theta = 1/2. \n
"""
def __init__(self, domain, field_name=None, solver_parameters=None,
@@ -818,10 +1074,10 @@ def setup(self, equation, apply_bcs=True, *active_labels):
class BDF2(MultilevelTimeDiscretisation):
"""
- Implements the implicit multistep BDF2 timestepping method
+ Implements the implicit multistep BDF2 timestepping method.
- The BDF2 timestepping method for operator F is written as
- y^(n+1) = (4/3)*y^n - (1/3)*y^(n-1) + (2/3)*dt*F[y^(n+1)]
+ The BDF2 timestepping method for operator F is written as: \n
+ y^(n+1) = (4/3)*y^n - (1/3)*y^(n-1) + (2/3)*dt*F[y^(n+1)] \n
"""
@property
@@ -919,10 +1175,10 @@ class TR_BDF2(TimeDiscretisation):
Implements the two stage implicit TR-BDF2 time stepping method, with a
trapezoidal stage (TR) followed by a second order backwards difference stage (BDF2).
- The TR-BDF2 time stepping method for operator F is written as
- y^(n+g) = y^n + dt*g/2*F[y^n] + dt*g/2*F[y^(n+g)] (TR stage)
- y^(n+1) = 1/(g(2-g))*y^(n+g) - (1-g)**2/(g(2-g))*y^(n) + (1-g)/(2-g)*dt*F[y^(n+1)] (BDF2 stage)
- for an off-centring parameter g (gamma).
+ The TR-BDF2 time stepping method for operator F is written as: \n
+ y^(n+g) = y^n + dt*g/2*F[y^n] + dt*g/2*F[y^(n+g)] (TR stage) \n
+ y^(n+1) = 1/(g(2-g))*y^(n+g) - (1-g)**2/(g(2-g))*y^(n) + (1-g)/(2-g)*dt*F[y^(n+1)] (BDF2 stage) \n
+ for an off-centring parameter g (gamma). \n
"""
def __init__(self, domain, gamma, field_name=None,
solver_parameters=None, options=None):
@@ -1047,8 +1303,8 @@ class Leapfrog(MultilevelTimeDiscretisation):
"""
Implements the multistep Leapfrog timestepping method.
- The Leapfrog timestepping method for operator F is written as
- y^(n+1) = y^(n-1) + 2*dt*F[y^n]
+ The Leapfrog timestepping method for operator F is written as: \n
+ y^(n+1) = y^(n-1) + 2*dt*F[y^n] \n
"""
@property
def nlevels(self):
@@ -1122,10 +1378,11 @@ def apply(self, x_out, *x_in):
class AdamsBashforth(MultilevelTimeDiscretisation):
"""
- Implements the explicit multistep Adams-Bashforth timestepping method of general order up to 5
+ Implements the explicit multistep Adams-Bashforth timestepping
+ method of general order up to 5.
- The general AB timestepping method for operator F is written as
- y^(n+1) = y^n + dt*(b_0*F[y^(n)] + b_1*F[y^(n-1)] + b_2*F[y^(n-2)] + b_3*F[y^(n-3)] + b_4*F[y^(n-4)])
+ The general AB timestepping method for operator F is written as: \n
+ y^(n+1) = y^n + dt*(b_0*F[y^(n)] + b_1*F[y^(n-1)] + b_2*F[y^(n-2)] + b_3*F[y^(n-3)] + b_4*F[y^(n-4)]) \n
"""
def __init__(self, domain, order, field_name=None,
solver_parameters=None, options=None):
@@ -1252,10 +1509,11 @@ def apply(self, x_out, *x_in):
class AdamsMoulton(MultilevelTimeDiscretisation):
"""
- Implements the implicit multistep Adams-Moulton timestepping method of general order up to 5
+ Implements the implicit multistep Adams-Moulton
+ timestepping method of general order up to 5
- The general AM timestepping method for operator F is written as
- y^(n+1) = y^n + dt*(b_0*F[y^(n+1)] + b_1*F[y^(n)] + b_2*F[y^(n-1)] + b_3*F[y^(n-2)])
+ The general AM timestepping method for operator F is written as \n
+ y^(n+1) = y^n + dt*(b_0*F[y^(n+1)] + b_1*F[y^(n)] + b_2*F[y^(n-1)] + b_3*F[y^(n-2)]) \n
"""
def __init__(self, domain, order, field_name=None,
solver_parameters=None, options=None):
@@ -1397,3 +1655,70 @@ def apply(self, x_out, *x_in):
self.x[n].assign(x_in[n])
solver.solve()
x_out.assign(self.x_out)
+
+
+class ImplicitMidpoint(ImplicitMultistage):
+ u"""
+ Implements the Implicit Midpoint method as a 1-stage Runge Kutta method.
+
+ The method, for solving
+ ∂y/∂t = F(y), can be written as: \n
+
+ k0 = F[y^n + 0.5*dt*k0] \n
+ y^(n+1) = y^n + dt*k0 \n
+ """
+ def __init__(self, domain, field_name=None, solver_parameters=None,
+ limiter=None, options=None):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
+ butcher_matrix = np.array([[0.5], [1.]])
+ super().__init__(domain, butcher_matrix, field_name,
+ solver_parameters=solver_parameters,
+ limiter=limiter, options=options)
+
+
+class QinZhang(ImplicitMultistage):
+ u"""
+ Implements Qin and Zhang's two-stage, 2nd order, implicit Runge–Kutta method.
+
+ The method, for solving
+ ∂y/∂t = F(y), can be written as: \n
+
+ k0 = F[y^n + 0.25*dt*k0] \n
+ k1 = F[y^n + 0.5*dt*k0 + 0.25*dt*k1] \n
+ y^(n+1) = y^n + 0.5*dt*(k0 + k1) \n
+ """
+ def __init__(self, domain, field_name=None, solver_parameters=None,
+ limiter=None, options=None):
+ """
+ Args:
+ domain (:class:`Domain`): the model's domain object, containing the
+ mesh and the compatible function spaces.
+ field_name (str, optional): name of the field to be evolved.
+ Defaults to None.
+ solver_parameters (dict, optional): dictionary of parameters to
+ pass to the underlying solver. Defaults to None.
+ limiter (:class:`Limiter` object, optional): a limiter to apply to
+ the evolving field to enforce monotonicity. Defaults to None.
+ options (:class:`AdvectionOptions`, optional): an object containing
+ options to either be passed to the spatial discretisation, or
+ to control the "wrapper" methods, such as Embedded DG or a
+ recovery method. Defaults to None.
+ """
+ butcher_matrix = np.array([[0.25, 0], [0.5, 0.25], [0.5, 0.5]])
+ super().__init__(domain, butcher_matrix, field_name,
+ solver_parameters=solver_parameters,
+ limiter=limiter, options=options)
diff --git a/integration-tests/model/test_time_discretisation.py b/integration-tests/model/test_time_discretisation.py
index ed882bcd9..6f114fe90 100644
--- a/integration-tests/model/test_time_discretisation.py
+++ b/integration-tests/model/test_time_discretisation.py
@@ -8,7 +8,7 @@ def run(timestepper, tmax, f_end):
return norm(timestepper.fields("f") - f_end) / norm(f_end)
-@pytest.mark.parametrize("scheme", ["ssprk", "implicit_midpoint",
+@pytest.mark.parametrize("scheme", ["ssprk", "TrapeziumRule", "ImplicitMidpoint", "QinZhang",
"RK4", "Heun", "BDF2", "TR_BDF2", "AdamsBashforth", "Leapfrog", "AdamsMoulton"])
def test_time_discretisation(tmpdir, scheme, tracer_setup):
if (scheme == "AdamsBashforth"):
@@ -27,8 +27,12 @@ def test_time_discretisation(tmpdir, scheme, tracer_setup):
if scheme == "ssprk":
transport_scheme = SSPRK3(domain)
- elif scheme == "implicit_midpoint":
+ elif scheme == "TrapeziumRule":
transport_scheme = TrapeziumRule(domain)
+ elif scheme == "ImplicitMidpoint":
+ transport_scheme = ImplicitMidpoint(domain)
+ elif scheme == "QinZhang":
+ transport_scheme = QinZhang(domain)
elif scheme == "RK4":
transport_scheme = RK4(domain)
elif scheme == "Heun":