Documentation

approximations

BoussinesqApproximation(Ra, kappa=1, g=1, rho=1, alpha=1)

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
kappa	Function | Number	thermal diffusivity	1
g	Function | Number	gravitational acceleration	1
rho	Function | Number	reference density	1
alpha	Function | Number	coefficient of thermal expansion	1

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(u, p, T, bottom_id, top_id, degree=4)

Typical simulation diagnostics used in geodynamical simulations.

Parameters:

Name	Type	Description	Default
u	Function	Firedrake function for the velocity	required
p	Function	Firedrake function for the pressure	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
domain_volume	The numerical domain's volume
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

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

def __init__(
    self,
    u: Function,
    p: Function,
    T: Function,
    bottom_id: int,
    top_id: int,
    degree: int = 4
):
    self.mesh = extract_unique_domain(u)
    self.u = u
    self.p = p
    self.T = T

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

    self.n = FacetNormal(self.mesh)

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

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:
        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

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()

Advances solver.

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

def solve(self):
    """Advances solver."""
    if not self._solver_setup:
        self.setup_solver()
    t = 0  # not used atm
    self.ts.advance(t)

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=None, 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	Optional[dict[int, dict[str, int | float]]]	Dictionary of identifier-value pairs specifying boundary conditions	None
mu	Function | int | float	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 | float]]	Dictionary of solver parameters provided to PETSc	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

def __init__(
    self,
    z: fd.Function,
    T: fd.Function,
    approximation: BaseApproximation,
    bcs: Optional[dict[int, dict[str, int | float]]] = None,
    mu: fd.Function | int | float = 1,
    quad_degree: int = 6,
    cartesian: bool = True,
    solver_parameters: Optional[dict[str, str | float]] = 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.approximation = approximation
    self.mu = ensure_constant(mu)
    self.solver_parameters = solver_parameters
    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 self.solver_parameters is None:
        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 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 is setup only last minute
    # so people can overwrite parameters we've setup here
    self._solver_setup = False

setup_solver()

Sets up the solver.

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

def setup_solver(self):
    """Sets up the solver."""
    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": self.mu},
                                                 **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": self.mu},
                                                    **self.solver_kwargs)
    self._solver_setup = True

solve()

Solves the system.

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

def solve(self):
    """Solves the system."""
    if not self._solver_setup:
        self.setup_solver()
    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

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

@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, 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

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 = 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

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 = 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

@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, 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

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

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, 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

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

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()

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

@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

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()

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

Crack-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()