# gusto package¶

## Submodules¶

class gusto.advection.NoAdvection(state, field, equation=None, *, solver_parameters=None, limiter=None)[source]

Bases: gusto.advection.Advection

An non-advection scheme that does nothing.

apply(x_in, x_out)[source]

Function takes x as input, computes L(x) as defined by the equation, and returns x_out as output.

Parameters: x – Function object, the input Function. x_out – Function object, the output Function.
lhs()[source]
rhs()[source]
update_ubar(xn, xnp1, alpha)[source]
class gusto.advection.ForwardEuler(state, field, equation=None, *, subcycles=None, solver_parameters=None, limiter=None)[source]

Bases: gusto.advection.ExplicitAdvection

Class to implement the forward Euler timestepping scheme: y_(n+1) = y_n + dt*L(y_n) where L is the advection operator

apply_cycle(x_in, x_out)[source]

Function takes x as input, computes L(x) as defined by the equation, and returns x_out as output.

Parameters: x – Function object, the input Function. x_out – Function object, the output Function.
lhs[source]
rhs[source]
class gusto.advection.SSPRK3(state, field, equation=None, *, subcycles=None, solver_parameters=None, limiter=None)[source]

Bases: gusto.advection.ExplicitAdvection

Class to implement the Strongly Structure Preserving Runge Kutta 3-stage timestepping method: y^1 = y_n + L(y_n) y^2 = (3/4)y_n + (1/4)(y^1 + L(y^1)) y_(n+1) = (1/3)y_n + (2/3)(y^2 + L(y^2)) where subscripts indicate the timelevel, superscripts indicate the stage number and L is the advection operator.

apply_cycle(x_in, x_out)[source]

Function takes x as input, computes L(x) as defined by the equation, and returns x_out as output.

Parameters: x – Function object, the input Function. x_out – Function object, the output Function.
lhs[source]
rhs[source]
solve_stage(x_in, stage)[source]
class gusto.advection.ThetaMethod(state, field, equation, theta=0.5, solver_parameters=None)[source]

Bases: gusto.advection.Advection

Class to implement the theta timestepping method: y_(n+1) = y_n + dt*(theta*L(y_n) + (1-theta)*L(y_(n+1))) where L is the advection operator.

apply(x_in, x_out)[source]

Function takes x as input, computes L(x) as defined by the equation, and returns x_out as output.

Parameters: x – Function object, the input Function. x_out – Function object, the output Function.
lhs[source]
rhs[source]

## gusto.configuration module¶

Some simple tools for making model configuration nicer.

class gusto.configuration.TimesteppingParameters(**kwargs)[source]

Bases: gusto.configuration.Configuration

Timestepping parameters for Gusto

alpha = 0.5
dt = None
maxi = 1
maxk = 4
class gusto.configuration.OutputParameters(**kwargs)[source]

Bases: gusto.configuration.Configuration

Output parameters for Gusto

checkpoint = True
dirname = None
dump_diagnostics = True
dumpfreq = 1
dumplist = None
dumplist_latlon = []
log_level = 30

log_level for logger, can be DEBUG, INFO or WARNING. Takes default value “warning”

perturbation_fields = []

List of fields for computing perturbations

point_data = []
project_fields = False

Should the output fields be interpolated or projected to a linear space? Default is interpolation.

steady_state_error_fields = []

List of fields to dump error fields for steady state simulation

class gusto.configuration.CompressibleParameters(**kwargs)[source]

Bases: gusto.configuration.Configuration

Physical parameters for Compressible Euler

L_v0 = 2500000.0
N = 0.01
R_d = 287.0
R_v = 461.0
T_0 = 273.15
c_pl = 4186.0
c_pv = 1885.0
c_vv = 1424.0
cp = 1004.5
cv = 717.0
g = 9.810616
kappa = 0.2857142857142857
p_0 = 100000.0
w_sat1 = 380.3
w_sat2 = -17.27
w_sat3 = 35.86
w_sat4 = 610.9
class gusto.configuration.ShallowWaterParameters(**kwargs)[source]

Bases: gusto.configuration.Configuration

Physical parameters for 3d Compressible Euler

H = None
Omega = 7.292e-05
g = 9.80616
class gusto.configuration.EadyParameters(**kwargs)[source]

Bases: gusto.configuration.Configuration

H = None
L = None
Nsq = 2.5e-05
dbdy = -1e-07
deltax = None
deltaz = None
f = None
fourthorder = False
class gusto.configuration.CompressibleEadyParameters(**kwargs)[source]

N = 0.005
Pi0 = 0.0
dthetady = -3e-06
g = 10.0
theta_surf = 300.0
class gusto.configuration.EmbeddedDGOptions(**kwargs)[source]

Bases: gusto.configuration.AdvectionOptions

embedding_space = None
name = 'embedded_dg'
class gusto.configuration.RecoveredOptions(**kwargs)[source]

Bases: gusto.configuration.AdvectionOptions

boundary_method = None
broken_space = None
embedding_space = None
name = 'recovered'
recovered_space = None

## gusto.diagnostics module¶

class gusto.diagnostics.Diagnostics(*fields)[source]

Bases: object

available_diagnostics = ['min', 'max', 'rms', 'l2', 'total']
static l2(f)[source]
static max(f)[source]
static min(f)[source]
register(*fields)[source]
static rms(f)[source]
static total(f)[source]
class gusto.diagnostics.CourantNumber(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'CourantNumber'
setup(state)[source]
class gusto.diagnostics.VelocityX(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'VelocityX'
setup(state)[source]
class gusto.diagnostics.VelocityZ(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'VelocityZ'
setup(state)[source]
class gusto.diagnostics.VelocityY(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'VelocityY'
setup(state)[source]
class gusto.diagnostics.Gradient(name)[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name

The name of this diagnostic field

setup(state)[source]
class gusto.diagnostics.RichardsonNumber(density_field, factor=1.0)[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'RichardsonNumber'
setup(state)[source]
class gusto.diagnostics.Energy(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

kinetic(u, factor=None)[source]

Computes 0.5*dot(u, u) with an option to multiply rho

class gusto.diagnostics.KineticEnergy(required_fields=())[source]
compute(state)[source]

Compute the diagnostic field from the current state

name = 'KineticEnergy'
class gusto.diagnostics.CompressibleKineticEnergy(required_fields=())[source]
compute(state)[source]

Compute the diagnostic field from the current state

name = 'CompressibleKineticEnergy'
class gusto.diagnostics.ExnerPi(reference=False)[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name

The name of this diagnostic field

setup(state)[source]
class gusto.diagnostics.Sum(field1, field2)[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name

The name of this diagnostic field

setup(state)[source]
class gusto.diagnostics.Difference(field1, field2)[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name

The name of this diagnostic field

setup(state)[source]
class gusto.diagnostics.SteadyStateError(state, name)[source]
name

The name of this diagnostic field

class gusto.diagnostics.Perturbation(name)[source]
name

The name of this diagnostic field

class gusto.diagnostics.PotentialVorticity(required_fields=())[source]

Bases: gusto.diagnostics.Vorticity

Diagnostic field for potential vorticity.

name = 'PotentialVorticity'
setup(state)[source]

Solver for potential vorticity. Solves a weighted mass system to generate the potential vorticity from known velocity and depth fields.

Parameters: state – The state containing model.
class gusto.diagnostics.Theta_e(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'Theta_e'
setup(state)[source]
class gusto.diagnostics.InternalEnergy(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'InternalEnergy'
setup(state)[source]
class gusto.diagnostics.Dewpoint(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'Dewpoint'
setup(state)[source]
class gusto.diagnostics.Temperature(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'Temperature'
setup(state)[source]
class gusto.diagnostics.Theta_d(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'Theta_d'
setup(state)[source]
class gusto.diagnostics.RelativeHumidity(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'RelativeHumidity'
setup(state)[source]
class gusto.diagnostics.HydrostaticImbalance(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'HydrostaticImbalance'
setup(state)[source]
class gusto.diagnostics.RelativeVorticity(required_fields=())[source]

Bases: gusto.diagnostics.Vorticity

Diagnostic field for relative vorticity.

name = 'RelativeVorticity'
setup(state)[source]

Solver for relative vorticity.

Parameters: state – The state containing model.
class gusto.diagnostics.AbsoluteVorticity(required_fields=())[source]

Bases: gusto.diagnostics.Vorticity

Diagnostic field for absolute vorticity.

name = 'AbsoluteVorticity'
setup(state)[source]

Solver for absolute vorticity.

Parameters: state – The state containing model.
class gusto.diagnostics.ShallowWaterKineticEnergy(required_fields=())[source]
compute(state)[source]

Compute the diagnostic field from the current state

name = 'ShallowWaterKineticEnergy'
class gusto.diagnostics.ShallowWaterPotentialEnergy(required_fields=())[source]
compute(state)[source]

Compute the diagnostic field from the current state

name = 'ShallowWaterPotentialEnergy'
class gusto.diagnostics.ShallowWaterPotentialEnstrophy(base_field_name='PotentialVorticity')[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name

The name of this diagnostic field

class gusto.diagnostics.Precipitation(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'Precipitation'
setup(state)[source]

## gusto.diffusion module¶

class gusto.diffusion.InteriorPenalty(state, V, kappa, mu, bcs=None)[source]

Bases: gusto.diffusion.Diffusion

Interior penalty diffusion method

Parameters: state – State object. V – Function space of diffused field direction – list containing directions in which function space mu: the penalty weighting function, which is

:recommended to be proportional to 1/dx :arg: kappa: strength of diffusion :arg: bcs: (optional) a list of boundary conditions to apply

apply(x_in, x_out)[source]

Function takes x as input, computes F(x) and returns x_out as output.

Parameters: x – Function object, the input Function. x_out – Function object, the output Function.

class gusto.eady_diagnostics.KineticEnergyY(required_fields=())[source]
compute(state)[source]

Computes the kinetic energy of the y component

name = 'KineticEnergyY'
class gusto.eady_diagnostics.CompressibleKineticEnergyY(required_fields=())[source]
compute(state)[source]

Computes the kinetic energy of the y component

name = 'CompressibleKineticEnergyY'
class gusto.eady_diagnostics.EadyPotentialEnergy(required_fields=())[source]
compute(state)[source]

Compute the diagnostic field from the current state

name = 'EadyPotentialEnergy'
class gusto.eady_diagnostics.CompressibleEadyPotentialEnergy(required_fields=())[source]
compute(state)[source]

Compute the diagnostic field from the current state

name = 'CompressibleEadyPotentialEnergy'
class gusto.eady_diagnostics.GeostrophicImbalance(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'GeostrophicImbalance'
setup(state)[source]
class gusto.eady_diagnostics.TrueResidualV(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'TrueResidualV'
setup(state)[source]
class gusto.eady_diagnostics.SawyerEliassenU(required_fields=())[source]

Bases: gusto.diagnostics.DiagnosticField

compute(state)[source]

Compute the diagnostic field from the current state

name = 'SawyerEliassenU'
setup(state)[source]

## gusto.forcing module¶

class gusto.forcing.CompressibleForcing(state, euler_poincare=True, linear=False, extra_terms=None, moisture=None)[source]

Bases: gusto.forcing.Forcing

Forcing class for compressible Euler equations.

apply(scaling, x_in, x_nl, x_out, **kwargs)[source]

Function takes x as input, computes F(x_nl) and returns x_out = x + scale*F(x_nl) as output.

Parameters: x_in – Function object x_nl – Function object x_out – Function object implicit – forcing stage for sponge and hydrostatic terms, if present
gravity_term()[source]
pressure_gradient_term()[source]
theta_forcing()[source]
class gusto.forcing.IncompressibleForcing(state, euler_poincare=True, linear=False, extra_terms=None, moisture=None)[source]

Bases: gusto.forcing.Forcing

Forcing class for incompressible Euler Boussinesq equations.

apply(scaling, x_in, x_nl, x_out, **kwargs)[source]

Function takes x as input, computes F(x_nl) and returns x_out = x + scale*F(x_nl) as output.

Parameters: x_in – Function object x_nl – Function object x_out – Function object implicit – forcing stage for sponge and hydrostatic terms, if present
gravity_term()[source]
pressure_gradient_term()[source]
class gusto.forcing.EadyForcing(state, euler_poincare=True, linear=False, extra_terms=None, moisture=None)[source]

Forcing class for Eady Boussinesq equations.

apply(scaling, x_in, x_nl, x_out, **kwargs)[source]

Function takes x as input, computes F(x_nl) and returns x_out = x + scale*F(x_nl) as output.

Parameters: x_in – Function object x_nl – Function object x_out – Function object implicit – forcing stage for sponge and hydrostatic terms, if present
forcing_term()[source]
class gusto.forcing.CompressibleEadyForcing(state, euler_poincare=True, linear=False, extra_terms=None, moisture=None)[source]

Forcing class for compressible Eady equations.

apply(scaling, x_in, x_nl, x_out, **kwargs)[source]

Function takes x as input, computes F(x_nl) and returns x_out = x + scale*F(x_nl) as output.

Parameters: x_in – Function object x_nl – Function object x_out – Function object implicit – forcing stage for sponge and hydrostatic terms, if present
forcing_term()[source]
class gusto.forcing.ShallowWaterForcing(state, euler_poincare=True, linear=False, extra_terms=None, moisture=None)[source]

Bases: gusto.forcing.Forcing

coriolis_term()[source]
pressure_gradient_term()[source]
topography_term()[source]

## gusto.initialisation_tools module¶

A module containing some tools for computing initial conditions, such as balanced initial conditions.

gusto.initialisation_tools.latlon_coords(mesh)[source]
gusto.initialisation_tools.sphere_to_cartesian(mesh, u_zonal, u_merid)[source]
gusto.initialisation_tools.incompressible_hydrostatic_balance(state, b0, p0, top=False, params=None)[source]
gusto.initialisation_tools.compressible_hydrostatic_balance(state, theta0, rho0, pi0=None, top=False, pi_boundary=Constant(FiniteElement('Real', None, 0), 0), water_t=None, solve_for_rho=False, params=None)[source]

Compute a hydrostatically balanced density given a potential temperature profile. By default, this uses a vertically-oriented hybridization procedure for solving the resulting discrete systems.

Parameters: state – The State object. theta0 – :class:.Functioncontaining the potential temperature. rho0 – Function to write the initial density into. top – If True, set a boundary condition at the top. Otherwise, set

it at the bottom. :arg pi_boundary: a field or expression to use as boundary data for pi on the top or bottom as specified. :arg water_t: the initial total water mixing ratio field.

gusto.initialisation_tools.remove_initial_w(u, Vv)[source]
gusto.initialisation_tools.eady_initial_v(state, p0, v)[source]
gusto.initialisation_tools.compressible_eady_initial_v(state, theta0, rho0, v)[source]
gusto.initialisation_tools.calculate_Pi0(state, theta0, rho0)[source]
gusto.initialisation_tools.saturated_hydrostatic_balance(state, theta_e, water_t, pi0=None, top=False, pi_boundary=Constant(FiniteElement('Real', None, 0), 1), max_outer_solve_count=40, max_theta_solve_count=5, max_inner_solve_count=3)[source]

Given a wet equivalent potential temperature, theta_e, and the total moisture content, water_t, compute a hydrostatically balance virtual potential temperature, dry density and water vapour profile.

The general strategy is to split up the solving into two steps: 1) finding rho to balance the theta profile 2) finding theta_v and r_v to get back theta_e and saturation We iteratively solve these steps until we (hopefully) converge to a solution.

Parameters: state – The State object. theta_e – The initial wet equivalent potential temperature profile. water_t – The total water pseudo-mixing ratio profile. pi0 – Optional function to put exner pressure into. top – If True, set a boundary condition at the top, otherwise it will be at the bottom. pi_boundary – The value of pi on the specified boundary. max_outer_solve_count – Max number of outer iterations for balance solver. max_theta_solve_count – Max number of iterations for theta solver (middle part of solve). max_inner_solve_count – Max number of iterations on the inner most loop for the water vapour solver.
gusto.initialisation_tools.unsaturated_hydrostatic_balance(state, theta_d, H, pi0=None, top=False, pi_boundary=Constant(FiniteElement('Real', None, 0), 2), max_outer_solve_count=40, max_inner_solve_count=20)[source]

Given vertical profiles for dry potential temperature and relative humidity compute hydrostatically balanced virtual potential temperature, dry density and water vapour profiles.

The general strategy is to split up the solving into two steps: 1) finding rho to balance the theta profile 2) finding theta_v and r_v to get back theta_d and H We iteratively solve these steps until we (hopefully) converge to a solution.

Parameters: state – The State object. theta_d – The initial dry potential temperature profile. H – The relative humidity profile. pi0 – Optional function to put exner pressure into. top – If True, set a boundary condition at the top, otherwise it will be at the bottom. pi_boundary – The value of pi on the specified boundary. max_outer_solve_count – Max number of iterations for outer loop of balance solver. max_inner_solve_count – Max number of iterations for inner loop of balanace solver.

## gusto.limiters module¶

class gusto.limiters.ThetaLimiter(space)[source]

Bases: object

A vertex based limiter for fields in the DG1xCG2 space, i.e. temperature variables. This acts like the vertex-based limiter implemented in Firedrake, but in addition corrects the central nodes to prevent new maxima or minima forming.

Initialise limiter

Parameters: space – the space in which theta lies.

It should be the DG1xCG2 space.

apply(field)[source]

The application of the limiter to the theta-space field.

check_midpoint_values(field)[source]

Checks the midpoint field values are less than the maximum and more than the minimum values. Amends them to the average if they are not.

copy_vertex_values(field)[source]

Copies the vertex values from temperature space to Q1DG space which only has vertices.

copy_vertex_values_back(field)[source]

Copies the vertex values back from the Q1DG space to the original temperature space.

remap_to_embedded_space(field)[source]

Remap from DG space to embedded DG space.

class gusto.limiters.NoLimiter[source]

Bases: object

A blank limiter that does nothing.

Initialise the blank limiter.

apply(field)[source]

The application of the blank limiter.

## gusto.linear_solvers module¶

class gusto.linear_solvers.CompressibleSolver(state, quadrature_degree=None, solver_parameters=None, overwrite_solver_parameters=False, moisture=None)[source]

Bases: gusto.linear_solvers.TimesteppingSolver

Timestepping linear solver object for the compressible equations in theta-pi formulation with prognostic variables u,rho,theta.

This solver follows the following strategy: (1) Analytically eliminate theta (introduces error near topography) (2) Solve resulting system for (u,rho) using a Schur preconditioner (3) Reconstruct theta

Parameters: state – a State object containing everything else. degree (quadrature) – tuple (q_h, q_v) where q_h is the required quadrature degree in the horizontal direction and q_v is that in the vertical direction (optional) (moisture) – solver parameters overwrite_solver_parameters – boolean, if True use only the solver_parameters that have been passed in, if False then update the default solver parameters with the solver_parameters passed in. (optional) – list of names of moisture fields.
solve()[source]

Apply the solver with rhs state.xrhs and result state.dy.

solver_parameters = {'fieldsplit_0': {'ksp_type': 'preonly', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'}, 'fieldsplit_1': {'ksp_type': 'preonly', 'mg_levels': {'ksp_chebyshev_esteig': True, 'ksp_max_it': 1, 'ksp_type': 'chebyshev', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'}, 'pc_type': 'gamg'}, 'ksp_gmres_restart': 50, 'ksp_max_it': 100, 'ksp_type': 'gmres', 'pc_fieldsplit_schur_fact_type': 'FULL', 'pc_fieldsplit_schur_precondition': 'selfp', 'pc_fieldsplit_type': 'schur', 'pc_type': 'fieldsplit'}
class gusto.linear_solvers.IncompressibleSolver(state, solver_parameters=None, overwrite_solver_parameters=False)[source]

Bases: gusto.linear_solvers.TimesteppingSolver

Timestepping linear solver object for the incompressible Boussinesq equations with prognostic variables u, p, b.

This solver follows the following strategy: (1) Analytically eliminate b (introduces error near topography) (2) Solve resulting system for (u,p) using a hybrid-mixed method (3) Reconstruct b

Parameters: state – a State object containing everything else. solver_parameters – (optional) Solver parameters. overwrite_solver_parameters – boolean, if True use only the solver_parameters that have been passed in, if False then update the default solver parameters with the solver_parameters passed in.
solve()[source]

Apply the solver with rhs state.xrhs and result state.dy.

solver_parameters = {'hybridization': {'ksp_rtol': 1e-08, 'ksp_type': 'cg', 'mg_levels': {'ksp_max_it': 2, 'ksp_type': 'chebyshev', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'}, 'pc_type': 'gamg'}, 'ksp_type': 'preonly', 'mat_type': 'matfree', 'pc_python_type': 'firedrake.HybridizationPC', 'pc_type': 'python'}
class gusto.linear_solvers.ShallowWaterSolver(state, solver_parameters=None, overwrite_solver_parameters=False)[source]

Bases: gusto.linear_solvers.TimesteppingSolver

Timestepping linear solver object for the nonlinear shallow water equations with prognostic variables u and D. The linearized system is solved using a hybridized-mixed method.

solve()[source]

Apply the solver with rhs state.xrhs and result state.dy.

solver_parameters = {'hybridization': {'ksp_rtol': 1e-08, 'ksp_type': 'cg', 'mg_levels': {'ksp_max_it': 2, 'ksp_type': 'chebyshev', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'}, 'pc_type': 'gamg'}, 'ksp_type': 'preonly', 'mat_type': 'matfree', 'pc_python_type': 'firedrake.HybridizationPC', 'pc_type': 'python'}
class gusto.linear_solvers.HybridizedCompressibleSolver(state, quadrature_degree=None, solver_parameters=None, overwrite_solver_parameters=False, moisture=None)[source]

Bases: gusto.linear_solvers.TimesteppingSolver

Timestepping linear solver object for the compressible equations in theta-pi formulation with prognostic variables u, rho, and theta.

This solver follows the following strategy:

1. Analytically eliminate theta (introduces error near topography)
(2a) Formulate the resulting mixed system for u and rho using a
hybridized mixed method. This breaks continuity in the linear perturbations of u, and introduces a new unknown on the mesh interfaces approximating the average of the Exner pressure perturbations. These trace unknowns also act as Lagrange multipliers enforcing normal continuity of the “broken” u variable.
(2b) Statically condense the block-sparse system into a single system
for the Lagrange multipliers. This is the only globally coupled system requiring a linear solver.
(2c) Using the computed trace variables, we locally recover the
broken velocity and density perturbations. This is accomplished in two stages: (i): Recover rho locally using the multipliers. (ii): Recover “broken” u locally using rho and the multipliers.
(2d) Project the “broken” velocity field into the HDiv-conforming
space using local averaging.
1. Reconstruct theta
Parameters: state – a State object containing everything else. degree (quadrature) – tuple (q_h, q_v) where q_h is the required quadrature degree in the horizontal direction and q_v is that in the vertical direction. (optional) (moisture) – solver parameters for the trace system. overwrite_solver_parameters – boolean, if True use only the solver_parameters that have been passed in, if False then update. the default solver parameters with the solver_parameters passed in. (optional) – list of names of moisture fields.
solve()[source]

Apply the solver with rhs state.xrhs and result state.dy.

solver_parameters = {'ksp_rtol': 1e-08, 'ksp_type': 'gmres', 'pc_type': 'bjacobi', 'sub_pc_type': 'ilu'}

## gusto.physics module¶

class gusto.physics.Condensation(state, iterations=1)[source]

Bases: gusto.physics.Physics

The process of condensation of water vapour into liquid water and evaporation of liquid water into water vapour, with the associated latent heat changes. The parametrization follows that used in Bryan and Fritsch (2002).

Parameters: state – State. object. iterations – number of iterations to do of condensation scheme per time step.
apply()[source]

Function computes the value of specific fields at the next time step.

class gusto.physics.Fallout(state, moments=<AdvectedMoments.M3: 1>, limit=True)[source]

Bases: gusto.physics.Physics

The fallout process of hydrometeors.

:arg state :class: .State. object. :arg moments: an AdvectedMoments Enum object, indicating which

rainfall scheme to use. Current valid values are: AdvectedMoments.M0 – advects all rain at constant speed; AdvectedMoments.M3 – advects mean mass of droplet distribution. The default value is AdvectedMoments.M3.
Parameters: limit – if True (the default value), applies a limiter to the rainfall advection.
apply()[source]

Function computes the value of specific fields at the next time step.

class gusto.physics.Coalescence(state, accretion=True, accumulation=True)[source]

Bases: gusto.physics.Physics

The process of the coalescence of cloud droplets to form rain droplets. These parametrizations come from Klemp and Wilhelmson (1978).

Parameters: state – State. object. accretion – Boolean which determines whether the accretion process is used. accumulation – Boolean which determines whether the accumulation process is used.
apply()[source]

Function computes the value of specific fields at the next time step.

class gusto.physics.Evaporation(state)[source]

Bases: gusto.physics.Physics

The process of evaporation of rain into water vapour with the associated latent heat change. This parametrization comes from Klemp and Wilhelmson (1978).

Parameters: state – State. object.
apply()[source]

Function computes the value of specific fields at the next time step.

class gusto.physics.AdvectedMoments[source]

Bases: enum.Enum

An Enum object storing moments of the rain drop size distribution. This is then used in deciding which to be advected in the Fallout step of the model.

M0 = 0
M3 = 1

## gusto.preconditioners module¶

A module containing some specialized preconditioners for Gusto applications.

class gusto.preconditioners.VerticalHybridizationPC[source]

Bases: firedrake.preconditioners.base.PCBase

A Slate-based python preconditioner for solving the hydrostatic pressure equation (after rewriting as a mixed vertical HDiv x L2 system). This preconditioner hybridizes a mixed system in the vertical direction. This means that velocities are rendered discontinuous in the vertical and Lagrange multipliers are introduced on the top/bottom facets to weakly enforce continuity through the top/bottom faces of each cell.

This PC assembles a statically condensed problem for the multipliers and inverts the resulting system using the provided solver options. The original unknowns are recovered element-wise by solving local linear systems.

All elimination and recovery kernels are generated using the Slate DSL in Firedrake.

Create a PC context suitable for PETSc.

Matrix free preconditioners should inherit from this class and implement:

apply(pc, x, y)[source]

We solve the forward eliminated problem for the approximate traces of the scalar solution (the multipliers) and reconstruct the “broken flux and scalar variable.”

Lastly, we project the broken solutions into the mimetic non-broken finite element space.

applyTranspose(pc, x, y)[source]

Apply the transpose of the preconditioner.

initialize(pc)[source]

Set up the problem context. Takes the original mixed problem and transforms it into the equivalent hybrid-mixed system.

A KSP object is created for the Lagrange multipliers on the top/bottom faces of the mesh cells.

update(pc)[source]

Update by assembling into the operator. No need to reconstruct symbolic objects.

view(pc, viewer=None)[source]

Viewer calls for the various configurable objects in this PC.

## gusto.recovery module¶

The recovery operators used for lowest-order advection schemes.

class gusto.recovery.Averager(v, v_out)[source]

Bases: object

An object that ‘recovers’ a low order field (e.g. in DG0) into a higher order field (e.g. in CG1). The code is essentially that of the Firedrake Projector object, using the “average” method, and could possibly be replaced by it if it comes into the master branch.

Parameters: v – the ufl.Expr or Function to project. v_out – Function to put the result in.
project()[source]

Apply the recovery.

class gusto.recovery.Boundary_Recoverer(v0, v1, v_out, method='density')[source]

Bases: object

An object that performs a recovery process at the domain boundaries that has second order accuracy. This is necessary because the Averager object does not recover a field with sufficient accuracy at the boundaries.

The strategy is to minimise the curvature of the function in the boundary cells, subject to the constraints of conserved mass and continuity on the interior facets. The quickest way to perform this is by using the analytic solution and a parloop.

Currently this is only implemented for the (DG0, DG1, CG1) set of spaces, and only on a PeriodicIntervalMesh or ‘PeriodicUnitIntervalMesh that has been extruded.

Parameters: v0 – the function providing the mass conservation constraints. Should be in DG0 and the initial function before the recovery process. v1 – the continuous function providing the continuity constraints. Should be in CG1 and is the field output by the initial recovery process. v_out – the function to be output. Should be in DG1. method – string giving the method used for the recovery. Valid options are ‘density’ and ‘physics’.
apply()[source]
class gusto.recovery.Recoverer(v_in, v_out, VDG=None, boundary_method=None)[source]

Bases: object

An object that ‘recovers’ a field from a low order space (e.g. DG0) into a higher order space (e.g. CG1). This encompasses the process of interpolating first to a the right space before using the Averager object, and also automates the boundary recovery process. If no boundary method is specified, this simply performs the action of the :class: Averager.

Parameters: v_in – the ufl.Expr or Function to project. (e.g. a VDG0 function) v_out – Function to put the result in. (e.g. a CG1 function) VDG – optional FunctionSpace. If not None, v_in is interpolated to this space first before recovery happens. boundary_method – a string defining which type of method needs to be used at the boundaries. Valid options are ‘density’, ‘velocity’ or ‘physics’.
extract_scalar()[source]
project()[source]

Perform the fully specified recovery.

restore_vector()[source]

## gusto.state module¶

class gusto.state.State(mesh, vertical_degree=None, horizontal_degree=1, family='RT', Coriolis=None, sponge_function=None, hydrostatic=None, timestepping=None, output=None, parameters=None, diagnostics=None, fieldlist=None, diagnostic_fields=None)[source]

Bases: object

Build a model state to keep the variables in, and specify parameters.

Parameters: mesh – The Mesh to use. vertical_degree – integer, required for vertically extruded meshes.

Specifies the degree for the pressure space in the vertical (the degrees for other spaces are inferred). Defaults to None. :arg horizontal_degree: integer, the degree for spaces in the horizontal (specifies the degree for the pressure space, other spaces are inferred) defaults to 1. :arg family: string, specifies the velocity space family to use. Options: “RT”: The Raviart-Thomas family (default, recommended for quads) “BDM”: The BDM family “BDFM”: The BDFM family :arg Coriolis: (optional) Coriolis function. :arg sponge_function: (optional) Function specifying a sponge layer. :arg timestepping: class containing timestepping parameters :arg output: class containing output parameters :arg parameters: class containing physical parameters :arg diagnostics: class containing diagnostic methods :arg fieldlist: list of prognostic field names :arg diagnostic_fields: list of diagnostic field classes

dump(t=0, pickup=False)[source]

Dump output :arg t: the current model time (default is zero). :arg pickup: recover state from the checkpointing file if true, otherwise dump and checkpoint to disk. (default is False).

initialise(initial_conditions)[source]

Initialise state variables

Parameters: initial_conditions – An iterable of pairs (field_name, pointwise_value)
set_reference_profiles(reference_profiles)[source]

Initialise reference profiles

Parameters: reference_profiles – An iterable of pairs (field_name, interpolatory_value)
setup_diagnostics()[source]

setup_dump(tmax, pickup=False)[source]

Setup dump files Check for existence of directory so as not to overwrite output files Setup checkpoint file

Parameters: tmax – model stop time pickup – recover state from the checkpointing file if true,

otherwise dump and checkpoint to disk. (default is False).

## gusto.thermodynamics module¶

Some thermodynamic expressions to help declutter the code.

gusto.thermodynamics.theta(parameters, T, p)[source]

Returns an expression for dry potential temperature theta in K.

Parameters: parameters – a CompressibleParameters object. T – temperature in K. p – pressure in Pa.
gusto.thermodynamics.pi(parameters, rho, theta_v)[source]

Returns an expression for the Exner pressure.

Parameters: parameters – a CompressibleParameters object. rho – the dry density of air in kg / m^3. theta – the potential temperature (or the virtual potential temperature for wet air), in K.
gusto.thermodynamics.pi_rho(parameters, rho, theta_v)[source]

Returns an expression for the derivative of Exner pressure with respect to density.

Parameters: parameters – a CompressibleParameters object. rho – the dry density of air in kg / m^3. theta – the potential temperature (or the virtual potential temperature for wet air), in K.
gusto.thermodynamics.pi_theta(parameters, rho, theta_v)[source]

Returns an expression for the deriavtive of Exner pressure with respect to potential temperature.

Parameters: parameters – a CompressibleParameters object. rho – the dry density of air in kg / m^3. theta – the potential temperature (or the virtual potential temperature for wet air), in K.
gusto.thermodynamics.p(parameters, pi)[source]

Returns an expression for the pressure in Pa from the Exner Pi.

Parameters: parameters – a CompressibleParameters object. pi – the Exner pressure.
gusto.thermodynamics.T(parameters, theta_v, pi, r_v=None)[source]

Returns an expression for temperature T in K.

Parameters: parameters – a CompressibleParameters object. theta_v – the virtual potential temperature in K. pi – the Exner pressure. r_v – the mixing ratio of water vapour.
gusto.thermodynamics.rho(parameters, theta_v, pi)[source]

Returns an expression for the dry density rho in kg / m^3 from the (virtual) potential temperature and Exner pressure.

Parameters: parameters – a CompressibleParameters object. theta_v – the virtual potential temperature in K. pi – the Exner pressure.
gusto.thermodynamics.r_sat(parameters, T, p)[source]

Returns an expression from Tetens’ formula for the saturation mixing ratio of water vapour.

Parameters: parameters – a CompressibleParameters object. T – the temperature in K. p – the pressure in Pa.
gusto.thermodynamics.Lv(parameters, T)[source]

Returns an expression for the latent heat of vaporisation of water.

Parameters: parameters – a CompressibleParameters object. T – the temperature in K.
gusto.thermodynamics.theta_e(parameters, T, p, r_v, r_t)[source]

Returns an expression for the wet equivalent potential temperature in K.

Parameters: parameters – a CompressibleParameters object. T – the temperature in K. p – the pressure in Pa. r_v – the mixing ratio of water vapour. r_t – the total mixing ratio of water.
gusto.thermodynamics.internal_energy(parameters, rho, T, r_v=0.0, r_l=0.0)[source]

Returns an expression for the (possibly wet) internal energy density in J.

Parameters: parameters – a CompressibleParameters object. rho – the dry density in kg / m^3. T – the temperature in K. r_v – the mixing ratio of water vapour. r_l – the mixing ratio of all forms of liquid water.
gusto.thermodynamics.RH(parameters, r_v, T, p)[source]

Returns an expression for the relative humidity.

Parameters: parameters – a CompressibleParameters object. r_v – the mixing ratio of water vapour. T – the temperature in K. p – the pressure in Pa.
gusto.thermodynamics.e_sat(parameters, T)[source]

Returns an expression for the saturated partial pressure of water vapour as a function of T, based on Tetens’ formula.

Parameters: parameters – a CompressibleParameters object. T – the temperature in K.
gusto.thermodynamics.r_v(parameters, H, T, p)[source]

Returns an expression for the mixing ratio of water vapour from the relative humidity, pressure and temperature.

Parameters: parameters – a CompressibleParameters object. H – the relative humidity (as a decimal). T – the temperature in K. p – the pressure in Pa.
gusto.thermodynamics.T_dew(parameters, p, r_v)[source]

Returns the dewpoint temperature as a function of temperature and the water vapour mixing ratio.

Parameters: parameters – a CompressibleParameters object. T – the temperature in K. r_v – the water vapour mixing ratio.

## gusto.timeloop module¶

class gusto.timeloop.CrankNicolson(state, advected_fields, linear_solver, forcing, diffused_fields=None, physics_list=None, prescribed_fields=None)[source]

Bases: gusto.timeloop.BaseTimestepper

This class implements a Crank-Nicolson discretisation, with Strang splitting and auxilliary semi-Lagrangian advection.

Parameters: state – a State object advected_fields – iterable of (field_name, scheme) pairs indicating the fields to advect, and the Advection to use. linear_solver – a TimesteppingSolver object forcing – a Forcing object diffused_fields – optional iterable of (field_name, scheme) pairs indictaing the fields to diffusion, and the Diffusion to use. physics_list – optional list of classes that implement physics schemes prescribed_fields – an order list of tuples, pairing a field name with a function that returns the field as a function of time.
passive_advection

Advected fields that are not part of the semi implicit step are passively advected

semi_implicit_step()[source]

Implement the semi implicit step for the timestepping scheme.

class gusto.timeloop.AdvectionDiffusion(state, advected_fields=None, diffused_fields=None, physics_list=None, prescribed_fields=None)[source]

Bases: gusto.timeloop.BaseTimestepper

This class implements a timestepper for the advection-diffusion equations. No semi implicit step is required.

passive_advection

semi_implicit_step()[source]

Implement the semi implicit step for the timestepping scheme.

## gusto.transport_equation module¶

class gusto.transport_equation.LinearAdvection(state, V, qbar, ibp=None, equation_form='advective', solver_params=None)[source]

Bases: gusto.transport_equation.TransportEquation

Class for linear transport equation.

Parameters: state – State object. V – :class:.FunctionSpace object. The function space that q lives in. qbar – The reference function that the equation has been linearised about. It is assumed that the reference velocity is zero and the ubar below is the nonlinear advecting velocity 0.5*(u’^(n+1) + u’(n))) ibp – string, stands for ‘integrate by parts’ and can take the value None, “once” or “twice”. Defaults to “once”. equation_form – (optional) string, can take the values ‘continuity’, which means the equation is in continuity form L(q) = div(u’*qbar), or ‘advective’, which means the equation is in advective form L(q) = u’ dot grad(qbar). Default is “advective” solver_params – (optional) dictionary of solver parameters to pass to the linear solver.
advection_term(q)[source]
class gusto.transport_equation.AdvectionEquation(state, V, *, ibp='once', equation_form='advective', vector_manifold=False, solver_params=None, outflow=False)[source]

Bases: gusto.transport_equation.TransportEquation

Class for discretisation of the transport equation.

Parameters: state – State object. V – :class:.FunctionSpace object. The function space that q lives in. ibp – string, stands for ‘integrate by parts’ and can take the value None, “once” or “twice”. Defaults to “once”. equation_form – (optional) string, can take the values ‘continuity’, which means the equation is in continuity form L(q) = div(u*q), or ‘advective’, which means the equation is in advective form L(q) = u dot grad(q). Default is “advective” vector_manifold – Boolean. If true adds extra terms that are needed for

advecting vector equations on manifolds. :arg solver_params: (optional) dictionary of solver parameters to pass to the

linear solver.
Parameters: outflow – Boolean specifying whether advected quantity can be advected out of domain.
advection_term(q)[source]
class gusto.transport_equation.EmbeddedDGAdvection(state, V, ibp='once', equation_form='advective', vector_manifold=False, solver_params=None, outflow=False, options=None)[source]

Class for the transport equation, using an embedded DG advection scheme.

Parameters: state – State object. V – :class:.FunctionSpace object. The function space that q lives in. ibp – (optional) string, stands for ‘integrate by parts’ and can take the value None, “once” or “twice”. Defaults to “once”. equation_form – (optional) string, can take the values ‘continuity’, which means the equation is in continuity form L(q) = div(u*q), or ‘advective’, which means the equation is in advective form L(q) = u dot grad(q). Default is “advective” vector_manifold – Boolean. If true adds extra terms that are needed for

advecting vector equations on manifolds. :arg solver_params: (optional) dictionary of solver parameters to pass to the

linear solver.
Parameters: outflow – Boolean specifying whether advected quantity can be advected out of domain. options – an instance of the AdvectionOptions class specifying which options to use with the embedded DG scheme.
class gusto.transport_equation.SUPGAdvection(state, V, ibp='twice', equation_form='advective', supg_params=None, solver_params=None, outflow=False)[source]

Class for the transport equation.

Parameters: state – State object. V – :class:.FunctionSpace object. The function space that q lives in. ibp – (optional) string, stands for ‘integrate by parts’ and can take the value None, “once” or “twice”. Defaults to “twice” since we commonly use this scheme for parially continuous spaces, in which case we don’t want to take a derivative of the test function. If using for a fully continuous space, we don’t integrate by parts at all (so you can set ibp=None). equation_form – (optional) string, can take the values ‘continuity’, which means the equation is in continuity form L(q) = div(u*q), or ‘advective’, which means the equation is in advective form L(q) = u dot grad(q). Default is “advective” supg_params – (optional) dictionary of parameters for the SUPG method. Can contain: ‘ax’, ‘ay’, ‘az’, which specify the coefficients in the x, y, z directions respectively ‘dg_direction’, which can be ‘horizontal’ or ‘vertical’, and specifies the direction in which the function space is discontinuous so that we can apply DG upwinding in this direction. Appropriate defaults are provided for these parameters, in particular, the space is assumed to be continuous. solver_params – (optional) dictionary of solver parameters to pass to the linear solver. outflow – Boolean specifying whether advected quantity can be advected out of domain.
class gusto.transport_equation.VectorInvariant(state, V, *, ibp='once', solver_params=None)[source]

Bases: gusto.transport_equation.TransportEquation

Class defining the vector invariant form of the vector advection equation.

Parameters: state – State object. V – Function space ibp – (optional) string, stands for ‘integrate by parts’ and can take the value None, “once” or “twice”. Defaults to “once”. solver_params – (optional) dictionary of solver parameters to pass to the linear solver.
advection_term(q)[source]
class gusto.transport_equation.EulerPoincare(state, V, *, ibp='once', solver_params=None)[source]

Class defining the Euler-Poincare form of the vector advection equation.

Parameters: state – State object. V – Function space ibp – string, stands for ‘integrate by parts’ and can take the value None, “once” or “twice”. Defaults to “once”. solver_params – (optional) dictionary of solver parameters to pass to the linear solver.
advection_term`(q)[source]