Time stepper

time_stepper

This module provides several classes to perform integration of time-dependent equations. Users choose if they require an explicit or implicit time integrator, and they instantiate one of the implemented algorithm class, for example, ERKEuler, by providing relevant parameters defined in the parent class (i.e. ERKGeneric or DIRKGeneric). Then, they call the advance method to request a solver update.

TimeIntegratorBase

Bases: ABC

Defines the API for all time integrators.

advance(t, update_forcings=None) abstractmethod

Advances equations for one time step.

Parameters:

Name	Type	Description	Default
t	float	Current simulation time	required
update_forcings	Optional[Function]	Firedrake function used to update any time-dependent boundary conditions	None

Source code in g-adopt/gadopt/time_stepper.py

@abstractmethod
def advance(self, t: float, update_forcings: Optional[firedrake.Function] = None):
    """Advances equations for one time step.

    Arguments:
      t: Current simulation time
      update_forcings: Firedrake function used to update any time-dependent boundary conditions

    """
    pass

initialize(init_solution) abstractmethod

Initialises the time integrator.

Parameters:

Name	Type	Description	Default
init_solution		Firedrake function representing the initial solution.	required

Source code in g-adopt/gadopt/time_stepper.py

@abstractmethod
def initialize(self, init_solution):
    """Initialises the time integrator.

    Arguments:
      init_solution: Firedrake function representing the initial solution.

    """
    pass

TimeIntegrator(equation, solution, dt, solution_old=None, solver_parameters=None, strong_bcs=None)

Bases: TimeIntegratorBase

Time integrator object that marches a single equation.

Parameters:

Name	Type	Description	Default
equation	Equation	G-ADOPT equation to integrate	required
solution	Function	Firedrake function representing the equation's solution	required
dt	float	Integration time step	required
solution_old	Optional[Function]	Firedrake function representing the equation's solution at the previous timestep	None
solver_parameters	Optional[dict[str, Any]]	Dictionary of solver parameters provided to PETSc	None
strong_bcs	Optional[list[DirichletBC]]	List of Firedrake Dirichlet boundary conditions	None

Source code in g-adopt/gadopt/time_stepper.py

def __init__(
    self,
    equation: Equation,
    solution: firedrake.Function,
    dt: float,
    solution_old: Optional[firedrake.Function] = None,
    solver_parameters: Optional[dict[str, Any]] = None,
    strong_bcs: Optional[list[firedrake.DirichletBC]] = None,
):
    super(TimeIntegrator, self).__init__()

    self.equation = equation
    self.test = firedrake.TestFunction(solution.function_space())
    self.solution = solution
    self.dt = float(dt)
    self.dt_const = ensure_constant(dt)
    self.solution_old = solution_old or firedrake.Function(solution, name='Old'+solution.name())

    # unique identifier used in solver
    self.name = '-'.join([self.__class__.__name__,
                          self.equation.__class__.__name__])

    self.solver_parameters = {}
    if solver_parameters:
        self.solver_parameters.update(solver_parameters)

    self.strong_bcs = strong_bcs or []
    self.hom_bcs = [bci.__class__(bci.function_space(), 0, bci.sub_domain) for bci in self.strong_bcs]

RungeKuttaTimeIntegrator(equation, solution, dt, solution_old=None, solver_parameters=None, strong_bcs=None)

Bases: TimeIntegrator

Abstract base class for all Runge-Kutta time integrators

Source code in g-adopt/gadopt/time_stepper.py

def __init__(
    self,
    equation: Equation,
    solution: firedrake.Function,
    dt: float,
    solution_old: Optional[firedrake.Function] = None,
    solver_parameters: Optional[dict[str, Any]] = None,
    strong_bcs: Optional[list[firedrake.DirichletBC]] = None,
):
    super(TimeIntegrator, self).__init__()

    self.equation = equation
    self.test = firedrake.TestFunction(solution.function_space())
    self.solution = solution
    self.dt = float(dt)
    self.dt_const = ensure_constant(dt)
    self.solution_old = solution_old or firedrake.Function(solution, name='Old'+solution.name())

    # unique identifier used in solver
    self.name = '-'.join([self.__class__.__name__,
                          self.equation.__class__.__name__])

    self.solver_parameters = {}
    if solver_parameters:
        self.solver_parameters.update(solver_parameters)

    self.strong_bcs = strong_bcs or []
    self.hom_bcs = [bci.__class__(bci.function_space(), 0, bci.sub_domain) for bci in self.strong_bcs]

get_final_solution() abstractmethod

Evaluates the final solution

Source code in g-adopt/gadopt/time_stepper.py

@abstractmethod
def get_final_solution(self):
    """Evaluates the final solution"""
    pass

solve_stage(i_stage, t, update_forcings=None) abstractmethod

Solves a single stage of step from t to t+dt. All functions that the equation depends on must be at right state corresponding to each sub-step.

Source code in g-adopt/gadopt/time_stepper.py

@abstractmethod
def solve_stage(self, i_stage, t, update_forcings=None):
    """Solves a single stage of step from t to t+dt.
    All functions that the equation depends on must be at right state
    corresponding to each sub-step.

    """
    pass

advance(t, update_forcings=None)

Advances equations for one time step.

Source code in g-adopt/gadopt/time_stepper.py

def advance(self, t, update_forcings=None):
    """Advances equations for one time step."""
    if not self._initialized:
        self.initialize(self.solution)
    for i in range(self.n_stages):
        self.solve_stage(i, t, update_forcings)
    self.get_final_solution()

ERKGeneric(equation, solution, dt, solution_old=None, bnd_conditions=None, solver_parameters={}, strong_bcs=None)

Bases: RungeKuttaTimeIntegrator

Generic explicit Runge-Kutta time integrator.

Implements the Butcher form. All terms in the equation are treated explicitly.

Parameters:

Name	Type	Description	Default
equation	Equation	G-ADOPT equation to solve	required
solution	Function	Firedrake function reperesenting the equation's solution	required
dt	float	Integration time step	required
solution_old	Optional[Function]	Firedrake function representing the equation's solution at the previous timestep	None
bnd_conditions	Optional[dict[int, dict[str, Number]]]	Dictionary of boundary conditions passed to the equation	None
solver_parameters	Optional[dict[str, Any]]	Dictionary of solver parameters provided to PETSc	{}
strong_bcs	Optional[list[DirichletBC]]	List of Firedrake Dirichlet boundary conditions	None

Source code in g-adopt/gadopt/time_stepper.py

def __init__(
    self,
    equation: Equation,
    solution: firedrake.Function,
    dt: float,
    solution_old: Optional[firedrake.Function] = None,
    bnd_conditions: Optional[dict[int, dict[str, Number]]] = None,
    solver_parameters: Optional[dict[str, Any]] = {},
    strong_bcs: Optional[list[firedrake.DirichletBC]] = None,
):
    super(ERKGeneric, self).__init__(
        equation, solution, dt, solution_old, solver_parameters, strong_bcs
    )
    self._initialized = False
    V = solution.function_space()
    assert V == equation.trial_space

    self.tendency = []
    for i in range(self.n_stages):
        k = firedrake.Function(V, name='tendency{:}'.format(i))
        self.tendency.append(k)

    # fully explicit evaluation
    trial = firedrake.TrialFunction(V)
    self.a_rk = self.equation.mass(trial)
    self.l_rk = self.dt_const * self.equation.residual(self.solution)

    self._nontrivial = self.l_rk != 0

    # construct expressions for stage solutions
    if self._nontrivial:
        self.sol_expressions = []
        for i_stage in range(self.n_stages):
            sol_expr = sum(map(operator.mul, self.tendency[:i_stage], self.a[i_stage][:i_stage]))
            self.sol_expressions.append(sol_expr)
        self.final_sol_expr = sum(map(operator.mul, self.tendency, self.b))

    self.update_solver()

update_solver()

Create solver objects

Source code in g-adopt/gadopt/time_stepper.py

def update_solver(self):
    """Create solver objects"""
    if self._nontrivial:
        self.solver = []
        for i in range(self.n_stages):
            prob = firedrake.LinearVariationalProblem(self.a_rk, self.l_rk, self.tendency[i], bcs=self.hom_bcs)
            solver = firedrake.LinearVariationalSolver(prob, options_prefix=self.name + '_k{:}'.format(i),
                                                       solver_parameters=self.solver_parameters)
            self.solver.append(solver)

update_solution(i_stage)

Computes the solution of the i-th stage

Tendencies must have been evaluated first.

Source code in g-adopt/gadopt/time_stepper.py

def update_solution(self, i_stage):
    """Computes the solution of the i-th stage

    Tendencies must have been evaluated first.

    """
    self.solution.assign(self.solution_old)
    if self._nontrivial and i_stage > 0:
        self.solution += self.sol_expressions[i_stage]

solve_tendency(i_stage, t, update_forcings=None)

Evaluates the tendency of i-th stage

Source code in g-adopt/gadopt/time_stepper.py

def solve_tendency(self, i_stage, t, update_forcings=None):
    """Evaluates the tendency of i-th stage"""
    if self._nontrivial:
        if update_forcings is not None:
            update_forcings(t + self.c[i_stage]*self.dt)
        self.solver[i_stage].solve()

DIRKGeneric(equation, solution, dt, solution_old=None, bnd_conditions=None, solver_parameters={}, strong_bcs=None, terms_to_add='all')

Bases: RungeKuttaTimeIntegrator

Generic implementation of Diagonally Implicit Runge Kutta schemes.

All derived classes must define the Butcher tableau coefficients :attr:a, :attr:b, :attr:c.

Parameters:

Name	Type	Description	Default
equation	Equation	G-ADOPT equation to solve	required
solution	Function	Firedrake function reperesenting the equation's solution	required
dt	float	Integration time step	required
solution_old	Optional[Function]	Firedrake function representing the equation's solution at the previous timestep	None
bnd_conditions	Optional[dict[int, dict[str, Number]]]	Dictionary of boundary conditions passed to the equation	None
solver_parameters	Optional[dict[str, Any]]	Dictionary of solver parameters provided to PETSc	{}
strong_bcs	Optional[list[DirichletBC]]	List of Firedrake Dirichlet boundary conditions	None
terms_to_add	Optional[str | list[str]]	Defines which terms of the equation are to be added to this solver. Default 'all' implies ['implicit', 'explicit', 'source'].	'all'

Source code in g-adopt/gadopt/time_stepper.py

def __init__(
    self,
    equation: Equation,
    solution: firedrake.Function,
    dt: float,
    solution_old: Optional[firedrake.Function] = None,
    bnd_conditions: Optional[dict[int, dict[str, Number]]] = None,
    solver_parameters: Optional[dict[str, Any]] = {},
    strong_bcs: Optional[list[firedrake.DirichletBC]] = None,
    terms_to_add: Optional[str | list[str]] = "all",
):
    super(DIRKGeneric, self).__init__(
        equation, solution, dt, solution_old, solver_parameters, strong_bcs
    )
    self.solver_parameters.setdefault('snes_type', 'newtonls')
    self._initialized = False

    fs = solution.function_space()
    assert fs == equation.trial_space

    mixed_space = len(fs) > 1

    # Allocate tendency fields
    self.k = []
    for i in range(self.n_stages):
        fname = '{:}_k{:}'.format(self.name, i)
        self.k.append(firedrake.Function(fs, name=fname))

    # construct variational problems
    self.F = []
    if not mixed_space:
        for i in range(self.n_stages):
            for j in range(i+1):
                if j == 0:
                    u = self.solution_old + self.a[i][j]*self.dt_const*self.k[j]
                else:
                    u += self.a[i][j]*self.dt_const*self.k[j]
            self.F.append(self.equation.mass(self.k[i]) - self.equation.residual(u))
    else:
        # solution must be split before computing sum
        # pass components to equation in a list
        for i in range(self.n_stages):
            for j in range(i+1):
                if j == 0:
                    u = []  # list of components in the mixed space
                    for s, k in zip(firedrake.split(self.solution_old), firedrake.split(self.k[j])):
                        u.append(s + self.a[i][j]*self.dt_const*k)
                else:
                    for l, k in enumerate(firedrake.split(self.k[j])):
                        u[l] += self.a[i][j]*self.dt_const*k
            self.F.append(self.equation.mass(self.k[i]) - self.equation.residual(u))
    self.update_solver()

    # construct expressions for stage solutions
    self.sol_expressions = []
    for i_stage in range(self.n_stages):
        sol_expr = sum(map(operator.mul, self.k[:i_stage+1], self.dt_const*self.a[i_stage][:i_stage+1]))
        self.sol_expressions.append(sol_expr)
    self.final_sol_expr = self.solution_old + sum(map(operator.mul, self.k, self.dt_const*self.b))

update_solver()

Create solver objects

Source code in g-adopt/gadopt/time_stepper.py

def update_solver(self):
    """Create solver objects"""
    self.solver = []
    for i in range(self.n_stages):
        p = firedrake.NonlinearVariationalProblem(self.F[i], self.k[i], bcs=self.hom_bcs)
        sname = '{:}_stage{:}_'.format(self.name, i)
        self.solver.append(
            firedrake.NonlinearVariationalSolver(
                p, solver_parameters=self.solver_parameters,
                options_prefix=sname))

update_solution(i_stage)

Updates solution to i_stage sub-stage.

Tendencies must have been evaluated first.

Source code in g-adopt/gadopt/time_stepper.py

def update_solution(self, i_stage):
    """Updates solution to i_stage sub-stage.

    Tendencies must have been evaluated first.

    """
    self.solution.assign(self.solution_old + self.sol_expressions[i_stage])

solve_tendency(i_stage, t, update_forcings=None)

Evaluates the tendency of i-th stage

Source code in g-adopt/gadopt/time_stepper.py

def solve_tendency(self, i_stage, t, update_forcings=None):
    """Evaluates the tendency of i-th stage"""
    if i_stage == 0:
        # NOTE solution may have changed in coupled system
        for bci in self.strong_bcs:
            bci.apply(self.solution)
        self.solution_old.assign(self.solution)
    if not self._initialized:
        raise ValueError('Time integrator {:} is not initialized'.format(self.name))
    if update_forcings is not None:
        update_forcings(t + self.c[i_stage]*self.dt)
    self.solver[i_stage].solve()

AbstractRKScheme()

Bases: ABC

Abstract class for defining Runge-Kutta schemes.

Derived classes must define the Butcher tableau (arrays :attr:a, :attr:b, :attr:c) and the CFL number (:attr:cfl_coeff).

Currently only explicit or diagonally implicit schemes are supported.

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

cfl_coeff abstractmethod property

CFL number of the scheme

Value 1.0 corresponds to Forward Euler time step.

a = np.array(self.a) abstractmethod instance-attribute property

Runge-Kutta matrix :math:a_{i,j} of the Butcher tableau

b = np.array(self.b) abstractmethod instance-attribute property

weights :math:b_{i} of the Butcher tableau

c = np.array(self.c) abstractmethod instance-attribute property

nodes :math:c_{i} of the Butcher tableau

ForwardEulerAbstract()

Bases: AbstractRKScheme

Forward Euler method

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

ERKLSPUM2Abstract()

Bases: AbstractRKScheme

ERKLSPUM2, 3-stage, 2nd order Explicit Runge Kutta method

From IMEX RK scheme (17) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

ERKLPUM2Abstract()

Bases: AbstractRKScheme

ERKLPUM2, 3-stage, 2nd order Explicit Runge Kutta method

From IMEX RK scheme (20) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

SSPRK33Abstract()

Bases: AbstractRKScheme

3rd order Strong Stability Preserving Runge-Kutta scheme, SSP(3,3).

This scheme has Butcher tableau

.. math:: \begin{array}{c|ccc} 0 & \ 1 & 1 \ 1/2 & 1/4 & 1/4 & \ \hline & 1/6 & 1/6 & 2/3 \end{array}

CFL coefficient is 1.0

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs3p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs4p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs5p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs6p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs7p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs8p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs9p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

eSSPRKs10p3Abstract()

Bases: AbstractRKScheme

Explicit SSP Runge-Kutta method with nondecreasing abscissas. See Isherwood, Grant, and Gottlieb (2018).

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

BackwardEulerAbstract()

Bases: AbstractRKScheme

Backward Euler method

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

ImplicitMidpointAbstract()

Bases: AbstractRKScheme

Implicit midpoint method, second order.

This method has the Butcher tableau

.. math:: \begin{array}{c|c} 0.5 & 0.5 \ \hline & 1.0 \end{array}

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

CrankNicolsonAbstract()

Bases: AbstractRKScheme

Crank-Nicolson scheme

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

DIRK22Abstract()

Bases: AbstractRKScheme

2-stage, 2nd order, L-stable Diagonally Implicit Runge Kutta method

This method has the Butcher tableau

.. math:: \begin{array}{c|cc} \gamma & \gamma & 0 \ 1 & 1-\gamma & \gamma \ \hline & 1/2 & 1/2 \end{array}

with :math:\gamma = (2 + \sqrt{2})/2.

From DIRK(2,3,2) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

DIRK23Abstract()

Bases: AbstractRKScheme

2-stage, 3rd order Diagonally Implicit Runge Kutta method

This method has the Butcher tableau

.. math:: \begin{array}{c|cc} \gamma & \gamma & 0 \ 1-\gamma & 1-2\gamma & \gamma \ \hline & 1/2 & 1/2 \end{array}

with :math:\gamma = (3 + \sqrt{3})/6.

From DIRK(2,3,3) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

DIRK33Abstract()

Bases: AbstractRKScheme

3-stage, 3rd order, L-stable Diagonally Implicit Runge Kutta method

From DIRK(3,4,3) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

DIRK43Abstract()

Bases: AbstractRKScheme

4-stage, 3rd order, L-stable Diagonally Implicit Runge Kutta method

From DIRK(4,4,3) IMEX scheme in Ascher et al. (1997)

Ascher et al. (1997). Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, 25:151-167. http://dx.doi.org/10.1137/0732037

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

DIRKLSPUM2Abstract()

Bases: AbstractRKScheme

DIRKLSPUM2, 3-stage, 2nd order, L-stable Diagonally Implicit Runge Kutta method

From IMEX RK scheme (17) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

DIRKLPUM2Abstract()

Bases: AbstractRKScheme

DIRKLPUM2, 3-stage, 2nd order, L-stable Diagonally Implicit Runge Kutta method

From IMEX RK scheme (20) in Higureras et al. (2014).

Higueras et al (2014). Optimized strong stability preserving IMEX Runge-Kutta methods. Journal of Computational and Applied Mathematics 272(2014) 116-140. http://dx.doi.org/10.1016/j.cam.2014.05.011

Source code in g-adopt/gadopt/time_stepper.py

def __init__(self):
    super(AbstractRKScheme, self).__init__()
    self.a = np.array(self.a)
    self.b = np.array(self.b)
    self.c = np.array(self.c)

    assert not np.triu(self.a, 1).any(), 'Butcher tableau must be lower diagonal'
    assert np.allclose(np.sum(self.a, axis=1), self.c), 'Inconsistent Butcher tableau: Row sum of a is not c'

    self.n_stages = len(self.b)
    self.butcher = np.vstack((self.a, self.b))

    self.is_implicit = np.diag(self.a).any()
    self.is_dirk = np.diag(self.a).all()

shu_osher_butcher(α_or_λ, β_or_μ)

Generate arrays composing the Butcher tableau of a Runge-Kutta method from the coefficient arrays of the equivalent, original or modified, Shu-Osher form. Code adapted from RK-Opt written in MATLAB by David Ketcheson. See also Ketcheson, Macdonald, and Gottlieb (2009).

Function arguments: α_or_λ : array_like, shape (n + 1, n) β_or_μ : array_like, shape (n + 1, n)

Source code in g-adopt/gadopt/time_stepper.py

def shu_osher_butcher(α_or_λ, β_or_μ):
    """
    Generate arrays composing the Butcher tableau of a Runge-Kutta method from the
    coefficient arrays of the equivalent, original or modified, Shu-Osher form.
    Code adapted from RK-Opt written in MATLAB by David Ketcheson.
    See also Ketcheson, Macdonald, and Gottlieb (2009).

    Function arguments:
    α_or_λ : array_like, shape (n + 1, n)
    β_or_μ : array_like, shape (n + 1, n)
    """

    X = np.identity(α_or_λ.shape[1]) - α_or_λ[:-1]
    A = np.linalg.solve(X, β_or_μ[:-1])
    b = np.transpose(β_or_μ[-1] + np.dot(α_or_λ[-1], A))
    c = np.sum(A, axis=1)
    return A, b, c