Skip to content

Documentation

approximations

BoussinesqApproximation(Ra, *, rho=1, alpha=1, T0=0, g=1, RaB=0, delta_rho=1, kappa=1, H=0)

Bases: BaseApproximation

Expressions for the Boussinesq approximation.

Density variations are considered small and only affect the buoyancy term. Reference parameters are typically constant. Viscous dissipation is neglected (Di << 1).

Parameters:

Name Type Description Default
Ra Function | Number

Rayleigh number

required
rho Function | Number

reference density

1
alpha Function | Number

coefficient of thermal expansion

1
T0 Function | Number

reference temperature

0
g Function | Number

gravitational acceleration

1
RaB Function | Number

compositional Rayleigh number; product of the Rayleigh and buoyancy numbers

0
delta_rho Function | Number

compositional density difference from the reference density

1
kappa Function | Number

thermal diffusivity

1
H Function | Number

internal heating rate

0
Note

The thermal diffusivity, gravitational acceleration, reference density, and coefficient of thermal expansion are normally kept at 1 when non-dimensionalised.

ExtendedBoussinesqApproximation(Ra, Di, mu=1, H=None, cartesian=True, **kwargs)

Bases: BoussinesqApproximation

Expressions for the extended Boussinesq approximation.

Extends the Boussinesq approximation by including viscous dissipation and work against gravity (both scaled with Di).

Parameters:

Name Type Description Default
Ra Number

Rayleigh number

required
Di Number

Dissipation number

required
mu Number

dynamic viscosity

1
H Optional[Number]

volumetric heat production

None
cartesian bool
  • True: gravity is assumed to point in the negative z-direction
  • False: gravity is assumed to point radially inward
True

Other Parameters:

Name Type Description
kappa Number

thermal diffusivity

g Number

gravitational acceleration

rho Number

reference density

alpha Number

coefficient of thermal expansion

Note

The thermal diffusivity, gravitational acceleration, reference density, and coefficient of thermal expansion are normally kept at 1 when non-dimensionalised.

TruncatedAnelasticLiquidApproximation(Ra, Di, Tbar=0, chi=1, cp=1, gamma0=1, cp0=1, cv0=1, **kwargs)

Bases: ExtendedBoussinesqApproximation

Truncated Anelastic Liquid Approximation

Compressible approximation. Excludes linear dependence of density on pressure (chi)

Parameters:

Name Type Description Default
Ra Number

Rayleigh number

required
Di Number

Dissipation number

required
Tbar Function | Number

reference temperature. In the diffusion term we use Tbar + T (i.e. T is the pertubartion) - default 0

0
chi Function | Number

reference isothermal compressibility

1
cp Function | Number

reference specific heat at constant pressure

1
gamma0 Function | Number

Gruneisen number (in pressure-dependent buoyancy term)

1
cp0 Function | Number

specific heat at constant pressure, reference for entire Mantle (in pressure-dependent buoyancy term)

1
cv0 Function | Number

specific heat at constant volume, reference for entire Mantle (in pressure-dependent buoyancy term)

1

Other Parameters:

Name Type Description
rho Number

reference density

alpha Number

reference thermal expansion coefficient

mu Number

viscosity used in viscous dissipation

H Number

volumetric heat production - default 0

cartesian bool
  • True: gravity points in negative z-direction
  • False: gravity points radially inward
kappa Number

diffusivity

g Number

gravitational acceleration

Note

The keyword arguments may be depth-dependent, but default to 1 if not supplied.

AnelasticLiquidApproximation(Ra, Di, Tbar=0, chi=1, cp=1, gamma0=1, cp0=1, cv0=1, **kwargs)

Bases: TruncatedAnelasticLiquidApproximation

Anelastic Liquid Approximation

Compressible approximation. Includes linear dependence of density on pressure (chi)

diagnostics

GeodynamicalDiagnostics(z, T, bottom_id, top_id, degree=4)

Typical simulation diagnostics used in geodynamical simulations.

Parameters:

Name Type Description Default
z Function

Firedrake function for the mixed Stokes function space

required
T Function

Firedrake function for the temperature

required
bottom_id int

bottom boundary identifier

required
top_id int

top boundary identifier

required
degree int

degree of the polynomial approximation

4
Note

All the diagnostics are returned as a float value.

Methods:

Name Description
u_rms

Root-mean squared velocity

u_rms_top

Root-mean squared velocity along the top boundary

Nu_top

Nusselt number at the top boundary

Nu_bottom

Nusselt number at the bottom boundary

T_avg

Average temperature in the domain

ux_max

Maximum velocity (optionally over a given boundary)

Source code in g-adopt/gadopt/diagnostics.py
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
def __init__(
    self,
    z: Function,
    T: Function,
    bottom_id: int,
    top_id: int,
    degree: int = 4,
):
    mesh = extract_unique_domain(z)

    self.u, self.p, *_ = z.subfunctions
    self.T = T

    self.dx = dx(domain=mesh, degree=degree)
    self.ds = (
        CombinedSurfaceMeasure(mesh, degree)
        if T.function_space().extruded
        else ds(mesh)
    )
    self.ds_t = self.ds(top_id)
    self.ds_b = self.ds(bottom_id)

    self.n = FacetNormal(mesh)

    self.domain_volume = assemble(Constant(1) * self.dx)
    self.top_surface = assemble(Constant(1) * self.ds_t)
    self.bottom_surface = assemble(Constant(1) * self.ds_b)

energy_solver

iterative_energy_solver_parameters: dict[str, Any] = {'mat_type': 'aij', 'snes_type': 'ksponly', 'ksp_type': 'gmres', 'ksp_rtol': 1e-05, 'pc_type': 'sor'} module-attribute

Default solver parameters for iterative solvers

direct_energy_solver_parameters: dict[str, Any] = {'mat_type': 'aij', 'snes_type': 'ksponly', 'ksp_type': 'preonly', 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps'} module-attribute

Default solver parameters for direct solvers

EnergySolver(T, u, approximation, delta_t, timestepper, bcs=None, solver_parameters=None, su_advection=False)

Timestepper solver for the energy equation.

Parameters:

Name Type Description Default
T Function

Firedrake function for the temperature

required
u Function

Firedrake function for the velocity

required
approximation BaseApproximation

G-ADOPT base approximation describing the system of equations

required
delta_t Constant

the simulation time step

required
bcs Optional[dict[int, dict[str, Number]]]

dictionary of identifier-value pairs specifying boundary conditions

None
solver_parameters Optional[dict[str, Any]]

additional solver parameters provided to PETSc

None
su_advection bool

whether to use of the streamline-upwind scheme

False
Source code in g-adopt/gadopt/energy_solver.py
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
def __init__(
    self,
    T: Function,
    u: Function,
    approximation: BaseApproximation,
    delta_t: Constant,
    timestepper: RungeKuttaTimeIntegrator,
    bcs: Optional[dict[int, dict[str, Number]]] = None,
    solver_parameters: Optional[dict[str, Any]] = None,
    su_advection: bool = False,
):
    self.T = T
    self.Q = T.function_space()
    self.mesh = self.Q.mesh()
    self.delta_t = delta_t
    rhocp = approximation.rhocp()
    self.eq = EnergyEquation(self.Q, self.Q, rhocp=rhocp)
    self.fields = {
        'diffusivity': ensure_constant(approximation.kappa()),
        'reference_for_diffusion': approximation.Tbar,
        'source': approximation.energy_source(u),
        'velocity': u,
        'advective_velocity_scaling': rhocp
    }
    sink = approximation.linearized_energy_sink(u)
    if sink:
        self.fields['absorption_coefficient'] = sink

    if su_advection:
        if not is_continuous(self.Q):
            raise TypeError("SU advection requires a continuous function space.")

        log("Using SU advection")
        # SU(PG) ala Donea & Huerta:
        # Columns of Jacobian J are the vectors that span the quad/hex
        # which can be seen as unit-vectors scaled with the dx/dy/dz in that direction (assuming physical coordinates x,y,z aligned with local coordinates)
        # thus u^T J is (dx * u , dy * v)
        # and following (2.44c) Pe = u^T J / 2 kappa
        # beta(Pe) is the xibar vector in (2.44a)
        # then we get artifical viscosity nubar from (2.49)

        J = Function(TensorFunctionSpace(self.mesh, 'DQ', 1), name='Jacobian').interpolate(Jacobian(self.mesh))
        kappa = self.fields['diffusivity'] + 1e-12  # Set lower bound for diffusivity in case zero diffusivity specified for pure advection.
        vel = self.fields['velocity']
        Pe = absv(dot(vel, J)) / (2*kappa)  # Calculate grid peclet number
        nubar = su_nubar(vel, J, Pe)  # Calculate SU artifical diffusion

        self.fields['su_nubar'] = nubar

    if solver_parameters is None:
        if self.mesh.topological_dimension() == 2:
            self.solver_parameters = direct_energy_solver_parameters.copy()
            if INFO >= log_level:
                # not really "informative", but at least we get a 1-line message saying we've passed the energy solve
                self.solver_parameters['ksp_converged_reason'] = None
        else:
            self.solver_parameters = iterative_energy_solver_parameters.copy()
            if DEBUG >= log_level:
                self.solver_parameters['ksp_monitor'] = None
            elif INFO >= log_level:
                self.solver_parameters['ksp_converged_reason'] = None
    else:
        self.solver_parameters = solver_parameters
    apply_strongly = is_continuous(T)
    self.strong_bcs = []
    self.weak_bcs = {}
    bcs = bcs or {}
    for id, bc in bcs.items():
        weak_bc = {}
        for type, value in bc.items():
            if type == 'T':
                if apply_strongly:
                    self.strong_bcs.append(DirichletBC(self.Q, value, id))
                else:
                    weak_bc['q'] = value
            else:
                weak_bc[type] = value
        self.weak_bcs[id] = weak_bc

    self.timestepper = timestepper
    self.T_old = Function(self.Q)
    # solver is setup only last minute
    # so people can overwrite parameters we've setup here
    self._solver_setup = False

setup_solver()

Sets up timestepper and associated solver.

Source code in g-adopt/gadopt/energy_solver.py
137
138
139
140
141
142
143
def setup_solver(self):
    """Sets up timestepper and associated solver."""
    self.ts = self.timestepper(self.eq, self.T, self.fields, self.delta_t,
                               bnd_conditions=self.weak_bcs, solution_old=self.T_old,
                               strong_bcs=self.strong_bcs,
                               solver_parameters=self.solver_parameters)
    self._solver_setup = True

solve(t=0, update_forcings=None)

Advances solver.

Source code in g-adopt/gadopt/energy_solver.py
145
146
147
148
149
def solve(self, t=0, update_forcings=None):
    """Advances solver."""
    if not self._solver_setup:
        self.setup_solver()
    self.ts.advance(t, update_forcings)

stokes_integrators

iterative_stokes_solver_parameters = {'mat_type': 'matfree', 'ksp_type': 'preonly', 'pc_type': 'fieldsplit', 'pc_fieldsplit_type': 'schur', 'pc_fieldsplit_schur_type': 'full', 'fieldsplit_0': {'ksp_type': 'cg', 'ksp_rtol': 1e-05, 'ksp_max_it': 1000, 'pc_type': 'python', 'pc_python_type': 'gadopt.SPDAssembledPC', 'assembled_pc_type': 'gamg', 'assembled_mg_levels_pc_type': 'sor', 'assembled_pc_gamg_threshold': 0.01, 'assembled_pc_gamg_square_graph': 100, 'assembled_pc_gamg_coarse_eq_limit': 1000, 'assembled_pc_gamg_mis_k_minimum_degree_ordering': True}, 'fieldsplit_1': {'ksp_type': 'fgmres', 'ksp_rtol': 0.0001, 'ksp_max_it': 200, 'pc_type': 'python', 'pc_python_type': 'firedrake.MassInvPC', 'Mp_pc_type': 'ksp', 'Mp_ksp_ksp_rtol': 1e-05, 'Mp_ksp_ksp_type': 'cg', 'Mp_ksp_pc_type': 'sor'}} module-attribute

Default solver parameters for iterative solvers

direct_stokes_solver_parameters = {'mat_type': 'aij', 'ksp_type': 'preonly', 'pc_type': 'lu', 'pc_factor_mat_solver_type': 'mumps'} module-attribute

Default solver parameters for direct solvers

newton_stokes_solver_parameters = {'snes_type': 'newtonls', 'snes_linesearch_type': 'l2', 'snes_max_it': 100, 'snes_atol': 1e-10, 'snes_rtol': 1e-05} module-attribute

Default solver parameters for non-linear systems

StokesSolver(z, T, approximation, bcs={}, mu=1, quad_degree=6, cartesian=True, solver_parameters=None, J=None, constant_jacobian=False, **kwargs)

Solves the Stokes system.

Parameters:

Name Type Description Default
z Function

Firedrake function representing the mixed Stokes system

required
T Function

Firedrake function representing the temperature

required
approximation BaseApproximation

Approximation describing the system of equations

required
bcs dict[int, dict[str, Number]]

Dictionary of identifier-value pairs specifying boundary conditions

{}
mu Function | Number

Firedrake function representing dynamic viscosity

1
quad_degree int

Quadrature degree. Default value is 2p + 1, where p is the polynomial degree of the trial space

6
cartesian bool

Whether to use Cartesian coordinates

True
solver_parameters Optional[dict[str, str | Number] | str]

Either a dictionary of PETSc solver parameters or a string specifying a default set of parameters defined in G-ADOPT

None
J Optional[Function]

Firedrake function representing the Jacobian of the system

None
constant_jacobian bool

Whether the Jacobian of the system is constant

False
Source code in g-adopt/gadopt/stokes_integrators.py
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
def __init__(
    self,
    z: fd.Function,
    T: fd.Function,
    approximation: BaseApproximation,
    bcs: dict[int, dict[str, Number]] = {},
    mu: fd.Function | Number = 1,
    quad_degree: int = 6,
    cartesian: bool = True,
    solver_parameters: Optional[dict[str, str | Number] | str] = None,
    J: Optional[fd.Function] = None,
    constant_jacobian: bool = False,
    **kwargs,
):
    self.Z = z.function_space()
    self.mesh = self.Z.mesh()
    self.test = fd.TestFunctions(self.Z)
    self.equations = StokesEquations(self.Z, self.Z, quad_degree=quad_degree,
                                     compressible=approximation.compressible)
    self.solution = z
    self.solution_old = None
    self.approximation = approximation
    self.mu = ensure_constant(mu)
    self.J = J
    self.constant_jacobian = constant_jacobian
    self.linear = not depends_on(self.mu, self.solution)

    self.solver_kwargs = kwargs
    u, p = fd.split(self.solution)
    self.k = upward_normal(self.Z.mesh(), cartesian)
    self.fields = {
        'velocity': u,
        'pressure': p,
        'viscosity': self.mu,
        'interior_penalty': fd.Constant(2.0),  # allows for some wiggle room in imposition of weak BCs
                                               # 6.25 matches C_ip=100. in "old" code for Q2Q1 in 2d.
        'source': self.approximation.buoyancy(p, T) * self.k,
        'rho_continuity': self.approximation.rho_continuity(),
    }

    self.weak_bcs = {}
    self.strong_bcs = []
    for id, bc in bcs.items():
        weak_bc = {}
        for bc_type, value in bc.items():
            if bc_type == 'u':
                self.strong_bcs.append(fd.DirichletBC(self.Z.sub(0), value, id))
            elif bc_type == 'ux':
                self.strong_bcs.append(fd.DirichletBC(self.Z.sub(0).sub(0), value, id))
            elif bc_type == 'uy':
                self.strong_bcs.append(fd.DirichletBC(self.Z.sub(0).sub(1), value, id))
            elif bc_type == 'uz':
                self.strong_bcs.append(fd.DirichletBC(self.Z.sub(0).sub(2), value, id))
            else:
                weak_bc[bc_type] = value
        self.weak_bcs[id] = weak_bc

    self.F = 0
    for test, eq, u in zip(self.test, self.equations, fd.split(self.solution)):
        self.F -= eq.residual(test, u, u, self.fields, bcs=self.weak_bcs)

    if isinstance(solver_parameters, dict):
        self.solver_parameters = solver_parameters
    else:
        if self.linear:
            self.solver_parameters = {"snes_type": "ksponly"}
        else:
            self.solver_parameters = newton_stokes_solver_parameters.copy()

        if INFO >= log_level:
            self.solver_parameters["snes_monitor"] = None

        if isinstance(solver_parameters, str):
            match solver_parameters:
                case "direct":
                    self.solver_parameters.update(direct_stokes_solver_parameters)
                case "iterative":
                    self.solver_parameters.update(
                        iterative_stokes_solver_parameters
                    )
                case _:
                    raise ValueError(
                        f"Solver type '{solver_parameters}' not implemented."
                    )
        elif self.mesh.topological_dimension() == 2 and cartesian:
            self.solver_parameters.update(direct_stokes_solver_parameters)
        else:
            self.solver_parameters.update(iterative_stokes_solver_parameters)

            if DEBUG >= log_level:
                self.solver_parameters['fieldsplit_0']['ksp_converged_reason'] = None
                self.solver_parameters['fieldsplit_1']['ksp_monitor'] = None
            elif INFO >= log_level:
                self.solver_parameters['fieldsplit_1']['ksp_converged_reason'] = None

    # solver object is set up later to permit editing default solver parameters
    self._solver_setup = False

setup_solver()

Sets up the solver.

Source code in g-adopt/gadopt/stokes_integrators.py
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
def setup_solver(self):
    """Sets up the solver."""
    # mu used in MassInvPC:
    mu_over_rho = self.mu / self.approximation.rho_continuity()
    if self.constant_jacobian:
        z_tri = fd.TrialFunction(self.Z)
        F_stokes_lin = fd.replace(self.F, {self.solution: z_tri})
        a, L = fd.lhs(F_stokes_lin), fd.rhs(F_stokes_lin)
        self.problem = fd.LinearVariationalProblem(a, L, self.solution,
                                                   bcs=self.strong_bcs,
                                                   constant_jacobian=True)
        self.solver = fd.LinearVariationalSolver(self.problem,
                                                 solver_parameters=self.solver_parameters,
                                                 options_prefix=self.name,
                                                 appctx={"mu": mu_over_rho},
                                                 **self.solver_kwargs)
    else:
        self.problem = fd.NonlinearVariationalProblem(self.F, self.solution,
                                                      bcs=self.strong_bcs, J=self.J)
        self.solver = fd.NonlinearVariationalSolver(self.problem,
                                                    solver_parameters=self.solver_parameters,
                                                    options_prefix=self.name,
                                                    appctx={"mu": mu_over_rho},
                                                    **self.solver_kwargs)
    self._solver_setup = True

solve()

Solves the system.

Source code in g-adopt/gadopt/stokes_integrators.py
262
263
264
265
266
267
def solve(self):
    """Solves the system."""
    if not self._solver_setup:
        self.setup_solver()
    self.solution_old = self.solution.copy(deepcopy=True)
    self.solver.solve()

create_stokes_nullspace(Z, closed=True, rotational=False, translations=None)

Create a null space for the mixed Stokes system.

Parameters:

Name Type Description Default
Z WithGeometry

Firedrake mixed function space associated with the Stokes system

required
closed bool

Whether to include a constant pressure null space

True
rotational bool

Whether to include all rotational modes

False
translations Optional[list[int]]

List of translations to include

None

Returns:

Type Description
MixedVectorSpaceBasis

A Firedrake mixed vector space basis incorporating the null space components

Source code in g-adopt/gadopt/stokes_integrators.py
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
def create_stokes_nullspace(
    Z: fd.functionspaceimpl.WithGeometry,
    closed: bool = True,
    rotational: bool = False,
    translations: Optional[list[int]] = None,
) -> fd.nullspace.MixedVectorSpaceBasis:
    """Create a null space for the mixed Stokes system.

    Arguments:
      Z: Firedrake mixed function space associated with the Stokes system
      closed: Whether to include a constant pressure null space
      rotational: Whether to include all rotational modes
      translations: List of translations to include

    Returns:
      A Firedrake mixed vector space basis incorporating the null space components

    """
    X = fd.SpatialCoordinate(Z.mesh())
    dim = len(X)
    V, W = Z.subfunctions

    if rotational:
        if dim == 2:
            rotV = fd.Function(V).interpolate(fd.as_vector((-X[1], X[0])))
            basis = [rotV]
        elif dim == 3:
            x_rotV = fd.Function(V).interpolate(fd.as_vector((0, -X[2], X[1])))
            y_rotV = fd.Function(V).interpolate(fd.as_vector((X[2], 0, -X[0])))
            z_rotV = fd.Function(V).interpolate(fd.as_vector((-X[1], X[0], 0)))
            basis = [x_rotV, y_rotV, z_rotV]
        else:
            raise ValueError("Unknown dimension")
    else:
        basis = []

    if translations:
        for tdim in translations:
            vec = [0] * dim
            vec[tdim] = 1
            basis.append(fd.Function(V).interpolate(fd.as_vector(vec)))

    if basis:
        V_nullspace = fd.VectorSpaceBasis(basis, comm=Z.mesh().comm)
        V_nullspace.orthonormalize()
    else:
        V_nullspace = V

    if closed:
        p_nullspace = fd.VectorSpaceBasis(constant=True, comm=Z.mesh().comm)
    else:
        p_nullspace = W

    return fd.MixedVectorSpaceBasis(Z, [V_nullspace, p_nullspace])

time_stepper

Timestepper implementation, mostly copied from Thetis.

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
17
18
19
20
21
22
23
24
25
26
@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
28
29
30
31
32
33
34
35
36
@abstractmethod
def initialize(self, init_solution):
    """Initialises the time integrator.

    Arguments:
      init_solution: Firedrake function representing the initial solution.

    """
    pass

TimeIntegrator(equation, solution, fields, 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 BaseEquation

G-ADOPT equation to integrate

required
solution Function

Firedrake function representing the equation's solution

required
fields dict[str, Function | Constant]

Dictionary of Firedrake fields passed to the equation

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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def __init__(
    self,
    equation: BaseEquation,
    solution: firedrake.Function,
    fields: dict[str, firedrake.Function | firedrake.Constant],
    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.fields = fields
    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 = [firedrake.DirichletBC(bci.function_space(), 0, bci.sub_domain) for bci in self.strong_bcs]

RungeKuttaTimeIntegrator(equation, solution, fields, 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
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
def __init__(
    self,
    equation: BaseEquation,
    solution: firedrake.Function,
    fields: dict[str, firedrake.Function | firedrake.Constant],
    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.fields = fields
    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 = [firedrake.DirichletBC(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
89
90
91
92
@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
 94
 95
 96
 97
 98
 99
100
101
@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
103
104
105
106
107
108
109
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, fields, 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 BaseEquation

G-ADOPT equation to solve

required
solution Function

Firedrake function reperesenting the equation's solution

required
fields dict[str, Function | Constant]

Dictionary of Firedrake fields passed to the equation

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
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
def __init__(
        self,
        equation: BaseEquation,
        solution: firedrake.Function,
        fields: dict[str, firedrake.Function | firedrake.Constant],
        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, fields, 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_term(self.test, trial)
    self.l_rk = self.dt_const*self.equation.residual(self.test, self.solution, self.solution, self.fields, bnd_conditions)

    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
168
169
170
171
172
173
174
175
176
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
182
183
184
185
186
187
188
189
190
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
192
193
194
195
196
197
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, fields, 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 BaseEquation

G-ADOPT equation to solve

required
solution Function

Firedrake function reperesenting the equation's solution

required
fields dict[str, Function | Constant]

Dictionary of Firedrake fields passed to the equation

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
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
def __init__(
        self,
        equation: BaseEquation,
        solution: firedrake.Function,
        fields: dict[str, firedrake.Function | firedrake.Constant],
        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, fields, 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_term(self.test, self.k[i]) -
                          self.equation.residual(self.test, u, self.solution_old, fields, bnd_conditions))
    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_term(self.test, self.k[i]) -
                          self.equation.residual(self.test, u, self.solution_old, fields, bnd_conditions))
    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
293
294
295
296
297
298
299
300
301
302
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
308
309
310
311
312
313
314
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
316
317
318
319
320
321
322
323
324
325
326
327
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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()

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

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

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

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

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

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

cfl_coeff()

CFL number of the scheme

Value 1.0 corresponds to Forward Euler time step.

Source code in g-adopt/gadopt/time_stepper.py
364
365
366
367
368
369
370
371
@abstractproperty
def cfl_coeff(self):
    """CFL number of the scheme

    Value 1.0 corresponds to Forward Euler time step.

    """
    pass

ForwardEulerAbstract()

Bases: AbstractRKScheme

Forward Euler method

Source code in g-adopt/gadopt/time_stepper.py
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
373
374
375
376
377
378
379
380
381
382
383
384
385
386
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
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
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