Time stepping

Overview

Terrarium explicitly separates process computations in compute_auxiliary! and compute_tendencies! from the choice of time stepping scheme. As a general rule, only models can be configured for timestepping. A model can be initialized for timestepping via

initialize(model::Terrarium.AbstractModel{NF}; ...) where NF
initialize(
    model::Terrarium.AbstractModel{NF},
    params;
    clock::Clock,
    inputs::InputSource,
    boundary_conditions,
    initializers,
    fields
) where NF

This will return a ModelIntegrator:

struct ModelIntegrator{NF, Arch<:Oceananigans.Architectures.AbstractArchitecture, Grid<:Oceananigans.Grids.AbstractGrid, TimeStepper<:Terrarium.AbstractTimeStepper, Model<:Terrarium.AbstractModel, StateVars<:Terrarium.AbstractStateVariables, ClockType<:Clock, Inits<:NamedTuple, Inputs<:InputSources, InitialTime} <: Oceananigans.AbstractModel{TimeStepper, Arch}

A single iteration of the ModelIntegrator roughly involves four steps:

1. Invoke update_inputs! to populate all input Fields from their respective InputSources based on the current clock time,
2. Invoke compute_auxiliary! on all model components to derive auxiliary state variables from the current prognostic state,
3. Invoke compute_tendencies! on all model components to calculate tendencies for all prognostic variables,
4. Apply the tendencies computed in step (3) to update the prognostic state based on the selected timestepping scheme. Higher order explicit or implicit timesteppers may repeat steps 1-3 multiple times within a single time step.

As an example, let's consider the construction and initialization of a SoilModel:

arch = CPU()
grid = ColumnGrid(arch, Float32, ExponentialSpacing(N=10))
soil = SoilEnergyWaterCarbon(Float32; hydrology = SoilHydrology(Float32, NoFlow())) # immobile soil water
model = SoilModel(grid; soil, timestepper = ForwardEuler(Float32))
integrator = initialize(model)

Integrator of SoilModel{Float32, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, NoFlow, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKBrooksCorey{Float32}}, Nothing}, ConstantSoilCarbonDensity{Float32}}, DefaultInitializer{Float32}, ForwardEuler{Float32}} with timestepper ForwardEuler{Float32}(300.0f0)
├── Current time: 0.0
├── StateVariables{Float32}(clock = Clock{Float32, Float64}(time=0 seconds, iteration=0, last_Δt=Inf days), prognostic = (:internal_energy,), auxiliary = (:temperature, :liquid_water_fraction, :ground_temperature, :saturation_water_ice, :water_table, :hydraulic_conductivity), inputs = (), namespaces = (:soil,), timestepper_cache = EmptyCache)

Here integrator corresponds to a ModelIntegrator configured for a ForwardEuler time stepping scheme that was set in the SoilModel. State variable Fields can be accessed via integrator.state:

integrator.state.temperature   # current temperature Field

1×1×10 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×10 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CPU with 1×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 3×1×16 StaticStridesOffsetArray(::Array{Float32, 3}, 0:2, 1:1, -2:13) with eltype Float32 with indices 0:2×1:1×-2:13
    └── max=0.0, min=0.0, mean=0.0

Time stepping schemes

Each model carries a single timestepper (set via the timestepper keyword on the model constructor). LandModel, and SoilModel with Richards-equation water flow, default to the implicit-explicit IMEX(ForwardEuler(NF; Δt), BackwardEuler(NF)), with a time step of 15 minutes for LandModel and 1 hour for SoilModel (see Implicit-explicit time stepping); all other models, including SoilModel without water flow, default to ForwardEuler(NF), whose default time step is 300 s. Terrarium currently provides two explicit time steppers, ForwardEuler and Heun:

struct ForwardEuler{NF} <: Terrarium.AbstractTimeStepper{NF}

struct Heun{NF} <: Terrarium.AbstractTimeStepper{NF}

Both are constructed with a default timestep Δt in seconds. Δt can be overridden at run time by specifying it when calling run! or timestep!.

Implicit-explicit time stepping

Explicit schemes are cheap and simple, but for diffusive processes the stable time step shrinks quadratically with the grid spacing (see Adaptive time stepping). A thin surface layer or a fine vertical resolution can therefore force steps of seconds to minutes, even when the quantity of interest evolves over hours. Terrarium addresses this with an implicit-explicit (IMEX) time stepper, IMEX, which integrates selected prognostic variables with an implicit scheme and all remaining variables with an explicit one.

For an implicitly integrated variable $u$, the tendency is split into the terms treated explicitly, the explicit tendency $G_\text{exp}$, and the terms treated implicitly, the implicit tendency $G_\text{imp}$ (e.g. the vertical diffusion operator):

\[\frac{\partial u}{\partial t} = G_\text{exp}(u) + G_\text{imp}(u)\,.\]

The explicit tendency holds everything else, such as sources, sinks, and fluxes through the boundaries of the column. For explicitly integrated variables, $G_\text{imp}(u) = 0$ and thus the whole tendency is explicit.

struct IMEX{NF, E<:Terrarium.AbstractTimeStepper{NF}, I<:Terrarium.AbstractTimeStepper{NF}} <: AbstractIMEX{NF}

Structure of an IMEX step

A single IMEX step of size $\Delta t$ consists of the following stages:

1. update_cache! on the implicit sub-stepper's cache stores what it needs from $u^n$ (for BackwardEuler, a copy of $u^n$ of each implicitly routed variable) before the explicit step overwrites it.
2. The explicit sub-stepper's step (explicit_step!) is applied to all prognostic variables. It evaluates the inputs, auxiliary variables, and tendencies itself; for implicitly routed variables, the tendency holds only the explicitly treated terms $G_\text{exp}$, and the implicit operator is omitted (see Routing variables below). The result is $u^+ = u^n + \Delta t\, G_\text{exp}$, where $G_\text{exp}$ is $G_\text{exp}(u^n)$ at $t^n$ for ForwardEuler, and the average of the tendencies at $u^n$ and at the predictor for Heun. For explicitly routed variables, this is the final result. For implicitly routed variables, it is the right-hand side of the implicit system.
3. For the implicitly routed variables, implicit_step! of the implicit sub-stepper solves

   \[u^{n+1} = u^+ + \Delta t\, G_\text{imp}(u^{n+1})\]

   in place.
4. closure!, which also enforces any constraints of the model on its prognostic state, is applied once, and the clock ticks once.

This is sequential (Lie) operator splitting: the IMEX composes the steps of its two partners without depending on how either is done.

IMEX(ForwardEuler(NF; Δt), BackwardEuler(NF)) is the default timestepper of LandModel (Δt = 900) and of SoilModel with Richards-equation water flow (Δt = 3600). Under IMEX, the snow meltwater drainage of a LandModel is implicit as well (see Implicit meltwater drainage), while the snow heat conduction, vegetation, surface hydrology, and surface energy balance remain explicit and are integrated with an explicit integrator (by default ForwardEuler) at the same step; the shorter default time step of LandModel is due to the explicit snow heat conduction (which will later also be treated implicitly). Value boundary conditions on temperature or pressure_head (e.g. PrescribedSurfaceTemperature) currently require an explicit timestepper instead (see Limitations).

With ForwardEuler as the explicit partner, an IMEX step costs one evaluation of the explicit right-hand side plus the implicit solve. Heun evaluates it twice and reduces the error of the explicitly integrated processes, but the splitting keeps the scheme first order overall. The sub-steppers are given in the order IMEX(explicit, implicit), which is checked when the integrator is initialized. The Δt of the explicit partner is the default time step of the IMEX. The implicit sub-stepper is an AbstractTimeStepper whose behavior is further specialized by an AbstractImplicitSolver:

abstract type AbstractImplicitSolver

update_cache!(cache::AbstractTimeStepperCache, integrator::ModelIntegrator, Δt, names::Tuple)

update_cache!(
    cache::Terrarium.BackwardEulerCache,
    integrator::ModelIntegrator,
    Δt,
    names::Tuple
)

explicit_step!(integrator::ModelIntegrator, timestepper::AbstractTimeStepper, Δt, names::Tuple)

implicit_step!(integrator::ModelIntegrator, timestepper::AbstractTimeStepper, Δt, names::Tuple)

implicit_step!(solver::AbstractImplicitSolver, operator::AbstractLinearOperator, workspaces::NamedTuple, state, model::AbstractModel, Δt)
implicit_step!(solver::AbstractImplicitSolver, operator::AbstractLinearOperator, workspace::ImplicitSolverWorkspace, state, model::AbstractModel, Δt)

implicit_step!(
    solver::Picard,
    operator::Terrarium.SplitOperator,
    workspaces::NamedTuple,
    state,
    model::Terrarium.AbstractModel,
    Δt
)

implicit_step!(
    solver::Picard,
    operator::AbstractLinearOperator,
    workspace::Terrarium.ImplicitSolverWorkspace,
    state,
    model::Terrarium.AbstractModel,
    Δt
)

implicit_step!(
    integrator::ModelIntegrator,
    timestepper::BackwardEuler,
    Δt,
    names::Tuple
)

Routing variables

Which variables are integrated implicitly is determined by the timestepping trait, evaluated once when the integrator is initialized as timestepping(var, model, imex). The default for every variable is Explicit. SoilModel and LandModel both specialize it so that internal_energy is Implicit under any AbstractIMEX timestepper, via timestepping(::AbstractVariable{:internal_energy}, ::SoilModel, ::AbstractIMEX) = Implicit(). The soil saturation saturation_water_ice is Implicit as well whenever the soil hydrology follows the Richards equation (SoilHydrology with RichardsEq), and Explicit otherwise, as given by default_timestepping(hydrology). LandModel and SnowModel also route snow_water_equivalent implicitly, for the meltwater drainage (see Implicit meltwater drainage); the snow energy, vegetation, and surface hydrology variables always remain explicit.

To change the routing, define a more specific method for your own model type or timestepper type. For example, to keep soil energy explicit under a custom IMEX subtype MyIMEX, define timestepping(::AbstractVariable{:internal_energy}, ::SoilModel, ::MyIMEX) = Explicit(). The resolved class of each prognostic variable is stored in the type of the timestepper cache, so that routing is resolved at compile time and kernels remain type stable.

timestepping(timestepper::AbstractTimeStepper)::Timestepping

timestepping(var::AbstractVariable, model::AbstractModel, timestepper::AbstractTimeStepper)::Timestepping

timestepping(state::StateVariables, ::Val{name})
timestepping(state::StateVariables, var::AbstractVariable{name})

timestepping(timestepper::Optional{AbstractTimeStepper}, state::StateVariables, ::Val{name})

timestepping(cache::IMEXCache)

Processes see the resolved class as a trailing Timestepping argument of compute_tendencies!. Under Implicit(), the soil thermodynamics omit the diffusive and advective heat fluxes through interior cell faces from the tendency (see Implicit vertical heat conduction), and the Richards-equation hydrology omits the Darcy fluxes through interior cell faces (see Implicit Richards equation).

Implicit Euler with a Picard solver

The only implicit scheme currently provided is first-order backward Euler, BackwardEuler, with a Picard nonlinear solver. When the integrator is initialized, BackwardEuler queries the model for the linearized operator of all variables routed to it with implicit_operator(model, state, names) and stores this single operator in its cache. For SoilModel and LandModel, it is a SplitOperator that holds a batched vertical TridiagonalOperator for each routed soil variable and solves saturation_water_ice before internal_energy (operator splitting; see Solvers for implicit time stepping and Coupling of water and energy).

BackwardEuler also resolves, for each variable, the field its linear system is solved for: the variable of its closure relation (as given by varname of the closure), i.e. the temperature for internal_energy and the pressure_head for saturation_water_ice, and for a variable without closure its implicit_variable, e.g. the snow_liquid_fraction for snow_water_equivalent. In each step, BackwardEuler stores the explicit increment $u^+ - u^n = \Delta t\, G_\text{exp}$ of every implicit variable in the ImplicitSolverWorkspace of the variable and resets all of their prognostic fields to $u^n$, and then calls the solver-level implicit_step!(solver, operator, workspaces, state, model, Δt) of its solver once; when a variable is solved, the solver forms its explicit stage from its current field and the increment. The solver iterates over linearizations of the operator, which it obtains from the process that owns each variable through the process hooks:

struct BackwardEuler{NF, S<:Terrarium.AbstractImplicitSolver} <: Terrarium.AbstractTimeStepper{NF}

uⁿ⁺¹ = u⁺ + Δt G_imp(uⁿ⁺¹),    u⁺ = uⁿ + Δt G_exp,

struct Picard <: Terrarium.AbstractImplicitSolver

The implicit steps of the soil variables are described with their processes: vertical heat conduction in Implicit vertical heat conduction, the Richards equation in Implicit Richards equation, and their coupling in Coupling of water and energy. Both follow the same pattern. The backward Euler step is linearized in the variable of the closure relation (the temperature or the total hydraulic head) with coefficients lagged at the latest iterate, which gives a tridiagonal system per column (see Tridiagonal operator). The conserved variable is then updated from the implicit fluxes through the interior faces rather than from the linearization, so the column-integrated energy and water change only by the boundary fluxes and sources in $G_\text{exp}$, independent of the number of iterations.

With iterations = 1, this is the linearly implicit scheme used by many land surface models [2]. The default is iterations = 2. With iterations = M > 1, the lagged coefficients are re-evaluated at the new iterate (after re-applying the closure) and the solve is repeated, which is a Picard iteration on the nonlinear system [3]. A single step is still only first order accurate in $\Delta t$, but additional iterations reduce the error from the lagged capacity near freezing fronts.

Process hooks

Each iteration of the Picard solver calls four hooks for each implicit variable, which the model forwards to the process that owns the variable, so that processes depend on neither the solver nor the timestepper. The hooks receive the explicit stage $u^+$ and dispatch on an Iteration{name, target, final}, which names the variable, the field target that its linear system is solved for (see Implicit Euler with a Picard solver), and whether the iteration is the last one:

• linearize!(operator, rhs, state, model, explicit_stage, Δt, ::Iteration{name}) evaluates the lagged coefficients of operator and the right-hand side rhs at the current iterate, where explicit_stage holds $u^+$;
• solve!(target, operator, rhs, Δt) solves the linear system for the field target; it dispatches on the operator type and needs no process-specific method;
• implicit_update!(state, model, operator, explicit_stage, Δt, ::Iteration{name}) updates the prognostic variable conservatively, from the fluxes implied by the solution;
• implicit_closure!(state, model, ::Iteration{name, target, final}) re-evaluates the closure, where final is true only in the last iteration.

After the last iteration, the prognostic field must hold $u^{n+1}$, with the closure variables derived from it consistent. A hook may read another implicitly integrated variable from its prognostic field: since all implicit variables are reset to $u^n$ before the first solve, it holds $u^n$ if that variable is solved later and $u^{n+1}$ if it was solved earlier, never the explicit stage.

Iteration{name, target, final}(number)

julia> iteration = Terrarium.Iteration(:internal_energy; target = :temperature)
Terrarium.Iteration{:internal_energy, :temperature, false}(1)

linearize!(operator::AbstractLinearOperator, rhs, state, model::AbstractModel, explicit_stage, Δt, ::Iteration{name})

linearize!(
    operator::TridiagonalOperator,
    rhs,
    state,
    grid,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    explicit_saturation,
    Δt
) where NF

linearize!(
    operator::TridiagonalOperator,
    rhs,
    state,
    grid,
    energy::SoilThermodynamics,
    soil::Terrarium.AbstractSoil,
    constants::PhysicalConstants,
    explicit_energy,
    Δt
)

linearize!(
    operator::DiagonalOperator,
    rhs,
    state,
    grid,
    snow::SingleLayerSnow,
    constants::PhysicalConstants,
    explicit_water,
    Δt,
    ::Terrarium.Iteration{:snow_water_equivalent}
)

implicit_update!(state, model::AbstractModel, operator::AbstractLinearOperator, explicit_stage, Δt, ::Iteration{name})

implicit_closure!(state, model::AbstractModel, ::Iteration{name, target, final})

Example

We configure a SoilModel with an IMEX(ForwardEuler(Float32), BackwardEuler(Float32)) timestepper, a thawed initial column, and a constant cooling flux at the surface via SoilHeatFlux (positive upward, i.e. out of the soil). The default time step of the IMEX is that of its explicit partner. The time step of one hour is more than ten times the explicit stability limit of this column (roughly six minutes):

solver = Picard(Float32; iterations = 2)
timestepper = IMEX(ForwardEuler(Float32; Δt = 3600), BackwardEuler(Float32; solver))
imex_grid = ColumnGrid(CPU(), Float32, ExponentialSpacing(Δz_min = 0.02, Δz_max = 0.5, N = 20))
imex_model = SoilModel(imex_grid; soil, timestepper)
imex_integrator = initialize(
    imex_model;
    boundary_conditions = SoilHeatFlux(80.0f0),
    initializers = (temperature = 0.5f0, saturation_water_ice = 0.8f0),
)
for _ in 1:24
    timestep!(imex_integrator, 3600.0)
end
imex_integrator.state.temperature   # temperature after 24 steps of one hour

1×1×20 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×20 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CPU with 1×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 3×1×26 StaticStridesOffsetArray(::Array{Float32, 3}, 0:2, 1:1, -2:23) with eltype Float32 with indices 0:2×1:1×-2:23
    └── max=0.5, min=-0.999971, mean=0.21328

The surface cools below 0 °C and a freezing front enters the column. The same configuration under ForwardEuler would require a time step of a few minutes.

Choosing Δt and iterations

The implicit step is unconditionally stable for the diffusion operator (the linear system is an M-matrix, see Tridiagonal operator, so the temperature obeys a discrete maximum principle), but it is only first order accurate. The choice of Δt is therefore governed by accuracy and by the time scales of the explicitly treated terms (e.g. boundary forcing), not by the diffusive CFL. A useful reference is the explicit limit cell_diffusion_timescale; steps of tens to hundreds of times that limit are stable, and the error of the solution decreases linearly when Δt is halved.

Phase change is the main source of error. With iterations = 1, the lagged apparent heat capacity overestimates the cooling of cells that begin to freeze within a step. Using iterations = 2 reduces the maximum temperature error of a freezing front by roughly a factor of ten in the tests, at the cost of one additional coefficient evaluation, tridiagonal solve, and closure per step. More iterations are not guaranteed to reduce the error further, since the regularized capacity of a cell switches on and off as the iterates enter and leave the phase change. For thawed columns (no phase change), a single iteration is sufficient.

Limitations

• Only the vertical heat conduction, the advection of heat by the soil water flux, and the Richards-equation Darcy flux in the soil can be integrated implicitly. The advection is first-order upwind and therefore numerically diffusive at large advective Courant numbers. Infiltrating water enters at the temperature of the top soil cell rather than at the skin or rain temperature, water removed by evapotranspiration carries no heat, and the surface excess pool stores no heat, i.e. the heat of the excess water leaves the system with it.
• Value (Dirichlet) boundary conditions on temperature, such as PrescribedSurfaceTemperature, are rejected when internal_energy is implicit, with an informative error at initialization. Use flux boundary conditions on internal_energy, such as SoilHeatFlux and GeothermalHeatFlux, or an explicit timestepper. Gradient (Neumann) boundary conditions on temperature are supported; the flux $-\kappa\, \partial T / \partial z$ they prescribe is treated explicitly, at $t^n$, like a flux boundary condition.
• Boundary fluxes are evaluated at $t^n$ and held fixed over the step; only the diffusive fluxes through interior faces are implicit.
• With active soil hydrology (e.g. RichardsEq), the energy and water equations are coupled by operator splitting (see Coupling of water and energy): the Richards step evaluates the hydraulic conductivity at the liquid water fraction at $t^n$, so the effect of freezing and thawing on the water flow is lagged by one step. This affects only accuracy; energy and water conservation are unaffected.
• The redistribution of water by adjust_saturation_profile! moves only the sensible heat of liquid water, also when the water it moves is frozen.
• Boundary Darcy fluxes, e.g. gravity drainage, are evaluated at $t^n$ and are not limited by the water available in the boundary cell. At very large steps they can overdraw a column so strongly that it alternates between empty and saturated; the water solution is then inaccurate, and since the hydrology closure moves the excess water of a cell (and its heat) before the heat step advects the heat of the water that caused it, the temperature is then no longer guaranteed to stay within its initial range. A TimeStepWizard keeps the step below this bound through the explicit water outflow limit of cell_advection_timescale (see Adaptive time stepping with IMEX).
• A prescribed infiltration flux on the saturation, such as InfiltrationFlux, forces water into the top cell whatever its state, also when it is frozen and conducts no water (the surface runoff of a LandModel instead limits the infiltration to what the top cell can take). Water then accumulates in frozen surface cells, and when they thaw within a step of several hours, the default two Picard iterations can fail to converge: the iterates alternately drain and overfill the cells near the surface, water is lost, and the solution diverges. In a test with 5 mm/day into freezing soil, this happened at steps of eight hours. Keep implicit steps to a few hours in this situation, or limit the prescribed infiltration to the capacity of the soil.
• The meltwater drained by the implicit snow step reaches the soil in the following step (see Implicit meltwater drainage).
• Value (Dirichlet) boundary conditions on pressure_head are rejected when saturation_water_ice is implicit; use flux conditions on saturation_water_ice (e.g. InfiltrationFlux) or gradient conditions on pressure_head (e.g. FreeDrainage).
• Only the FreeWater freeze curve is supported, since the regularized apparent heat capacity is only implemented for it.
• The scheme is first order accurate in time.

Reactant and Enzyme

IMEX time stepping with the Picard solver and the TridiagonalOperator is tested both under Reactant (ReactantState grids, compiled with XLA) and for reverse-mode differentiation with Enzyme. The validity checks of the implicit operators run on the host when the integrator is initialized (in implicit_operator), so that the kernels contain no error paths. As with other timesteppers, adaptive stepping with a host-side TimeStepWizard does not apply to compiled Reactant runs.

Transient simulations

ModelIntegrator defines dispatches for timestep! and run! that allow for transient, open-ended simulations starting from the initialized state. We can use timestep! to take a single time step,

timestep!(integrator)

and run! to advance over a fixed number steps or a specified period,

# Advance by a fixed number of steps
run!(integrator; steps = 100)

# Advance for a given wall-clock period
run!(integrator; period = Day(30), Δt = 3600.0)

Integrator of SoilModel{Float32, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, NoFlow, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKBrooksCorey{Float32}}, Nothing}, ConstantSoilCarbonDensity{Float32}}, DefaultInitializer{Float32}, ForwardEuler{Float32}} with timestepper ForwardEuler{Float32}(300.0f0)
├── Current time: 2.6223e6
├── StateVariables{Float32}(clock = Clock{Float32, Float64}(time=30.351 days, iteration=821, last_Δt=1 hour), prognostic = (:internal_energy,), auxiliary = (:temperature, :liquid_water_fraction, :ground_temperature, :saturation_water_ice, :water_table, :hydraulic_conductivity), inputs = (), namespaces = (:soil,), timestepper_cache = EmptyCache)

Setting up a Simulation

The timestep! and run! methods allow us to control the integrate the model forward in time, but we only have access to the transient model state at each time step. In order for the model to be useful, we also need to be able to define a finite time period over which to run a simulation and save outputs during that time period.

Fortunately, Terrarium's ModelIntegrator implements the Oceananigans model interface. As a result, we can take advantage of the built-in Oceananigans infrastructure for configuring and running Simulations:

sim = Simulation(integrator; stop_time = 24*3600.0, Δt = 300.0)
run!(sim)

[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (31.304 ms)
[ Info: Executing initial time step...
[ Info: Simulation is stopping after running for 0 seconds.
[ Info: Simulation time 30.354 days equals or exceeds stop time 1 day.
[ Info:     ... initial time step complete (91.093 ms).

Note that stop_time is by default expressed in the same units as the clock, which is here seconds. However, Oceananigans.Units provides convenient conversion factors for other time units (e.g. days, hours):

using Oceananigans.Units: days

sim = Simulation(integrator; stop_time = 1days, Δt = 300.0)

Simulation of ModelIntegrator{Float32, CPU, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, ForwardEuler{Float32}, SoilModel{Float32, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, NoFlow, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKBrooksCorey{Float32}}, Nothing}, ConstantSoilCarbonDensity{Float32}}, DefaultInitializer{Float32}, ForwardEuler{Float32}}, StateVariables{Float32, (:internal_energy,), (), (:temperature, :liquid_water_fraction, :ground_temperature, :saturation_water_ice, :water_table, :hydraulic_conductivity), (), (:soil,), Tuple{Field{Center, Center, Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 16), (-1, 0, -3), Array{Float32, 3}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Nothing, @NamedTuple{south_and_north::Nothing, bottom_and_top::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 16)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:16)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}}, Tuple{Field{Center, Center, Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 16), (-1, 0, -3), Array{Float32, 3}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Nothing, @NamedTuple{south_and_north::Nothing, bottom_and_top::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 16)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:16)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}}, Tuple{Field{Center, Center, Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 16), (-1, 0, -3), Array{Float32, 3}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Nothing, @NamedTuple{south_and_north::Nothing, bottom_and_top::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 16)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:16)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}, Field{Center, Center, Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 16), (-1, 0, -3), Array{Float32, 3}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Nothing, @NamedTuple{south_and_north::Nothing, bottom_and_top::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 16)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:16)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}, Field{Center, Center, Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, UnitRange{Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 1), (-1, 0, 9), SubArray{Float32, 3, Array{Float32, 3}, Tuple{Base.Slice{Base.OneTo{Int64}}, Base.Slice{Base.OneTo{Int64}}, UnitRange{Int64}}, true}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, Nothing, @NamedTuple{bottom_and_top::Nothing, south_and_north::Nothing, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{bottom_and_top::Tuple{Nothing, Nothing}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}, Field{Center, Center, Center, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 16), (-1, 0, -3), Array{Float32, 3}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, Nothing, @NamedTuple{south_and_north::Nothing, bottom_and_top::KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_bottom_and_top_halo!)}, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 16)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:16)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{south_and_north::Tuple{Nothing, Nothing}, bottom_and_top::Tuple{NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}, NoFluxBoundaryCondition{Flux{Oceananigans.Utils.ExplicitTimeDiscretization}}}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}, Field{Center, Center, Nothing, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 1), (-1, 0, 0), Array{Float32, 3}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, Nothing, @NamedTuple{bottom_and_top::Nothing, south_and_north::Nothing, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 1)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:1)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{bottom_and_top::Tuple{Nothing, Nothing}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}, Field{Center, Center, Face, Nothing, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (17,), (-4,), Vector{Float32}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (16,), (-3,), Vector{Float32}}}, Float32, Float32, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 1, (3,), (-1,), StepRangeLen{Float32, Float64, Float64, Int64}}, Float32, Float32, Nothing, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, Oceananigans.Utils.StaticStridesOffsetArray{Float32, 3, (3, 1, 17), (-1, 0, -3), Array{Float32, 3}}, Float32, Oceananigans.BoundaryConditions.FieldBoundaryConditions{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Nothing, Nothing, Nothing, Nothing, Nothing, @NamedTuple{bottom_and_top::Nothing, south_and_north::Nothing, west_and_east::Oceananigans.BoundaryConditions.PeriodicFillHalo{KernelAbstractions.Kernel{KernelAbstractions.CPU, KernelAbstractions.NDIteration.StaticSize{(1, 17)}, Oceananigans.Utils.OffsetStaticSize{(1:1, 1:17)}, typeof(Oceananigans.BoundaryConditions.cpu__fill_periodic_west_and_east_halo!)}, 1, 1}}, @NamedTuple{bottom_and_top::Tuple{Nothing, Nothing}, south_and_north::Tuple{Nothing, Nothing}, west_and_east::Tuple{Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}, Oceananigans.BoundaryConditions.BoundaryCondition{Oceananigans.BoundaryConditions.Periodic, Nothing}}}}, Nothing, Nothing}}, Tuple{}, Tuple{StateVariables{Float32, (), (), (), (), (), Tuple{}, Tuple{}, Tuple{}, Tuple{}, Tuple{}, Terrarium.EmptyCache{Float32}, Clock{Float32, Float32, Float64, Int64, Int64}}}, Terrarium.EmptyCache{Float32}, Clock{Float32, Float32, Float64, Int64, Int64}}, Clock{Float32, Float32, Float64, Int64, Int64}, @NamedTuple{}, InputSources{Float32, nothing, Tuple{}}, Float32}
├── Next time step: 5 minutes
├── run_wall_time: 0 seconds
├── run_wall_time / iteration: 0 seconds
├── stop_time: 1 day
├── stop_iteration: Inf
├── wall_time_limit: Inf
├── minimum_relative_step: 0.0
├── callbacks: OrderedDict with 3 entries:
│   ├── stop_time_exceeded => Callback of stop_time_exceeded on IterationInterval(1)
│   ├── stop_iteration_exceeded => Callback of stop_iteration_exceeded on IterationInterval(1)
│   └── wall_time_limit_exceeded => Callback of wall_time_limit_exceeded on IterationInterval(1)
└── output_writers: OrderedDict with no entries

Callbacks

Terrarium simulations inherit the full Oceananigans Callback machinery. A callback is a function that is called by Simulation at a given schedule during run!:

using Oceananigans: Callback, IterationInterval

function print_progress(sim)
    t = sim.model.clock.time
    iter = sim.model.clock.iteration
    @info "Iteration $iter, time $t s"
end

sim.callbacks[:progress] = Callback(print_progress, IterationInterval(100))

Callback of print_progress on IterationInterval(100)

The callback function receives the Simulation object, giving it access to the full integrator via sim.model, and to the current state via sim.model.state.

Adaptive time stepping

Some land components, such as the nonlinear advection-diffusion equations for soil energy, water, and carbon transport, are subject to stability constraints when using explicit time-stepping schemes like ForwardEuler or Heun. A suitable constraint is the diffusive Courant–Friedrichs–Lewy (CFL) condition. The standard formulation of the CFL for advective transport places a limit on the wave propagation speed based on the grid cell length $\Delta x$,

\[\frac{ \lvert v \rvert \Delta t}{\Delta x} \leq C\,.\]

While diffusive systems do not have an independent velocity $v$, they do have a characteristic length scale based on the diffusion coefficient $D$,

\[\ell \sim \sqrt{D \Delta t}\]

which we can understand as the distance over which information gets dispersed in a single time step $\Delta t$. Plugging that in as the numerator in CFL yields,

\[\frac{\sqrt{D \Delta t}}{\Delta x} \leq C\]

and solving for $\Delta t$ we get,

\[\Delta t \leq \frac{C^2 \Delta x^2}{D}\,.\]

Since $C$ is a constant, we can simply redefine $C_{\text{diff}} = C^2$. The choice of $C$ can be determined by von Neumann stability analysis for any given time stepping scheme; in the case of ForwardEuler, this leads to $C = \frac{1}{2}$.

Terrarium defines the diffusive CFL via the Oceananigans diagnostic cell_diffusion_timescale which is computed as the minimum over all grid cells of the diffusive timescale $\tau = \Delta z^2 / D$. For the soil model it is the minimum of the thermal timescale $\Delta z^2 C / \kappa$ (heat conduction, with thermal conductivity $\kappa$ and volumetric heat capacity $C$) and, for a Richards-equation hydrology, the hydraulic timescale $\Delta z^2 (\partial\theta/\partial\psi) / K$ (with hydraulic conductivity $K$ and specific moisture capacity $\partial\theta/\partial\psi$).

Because a ModelIntegrator implements the Oceananigans model interface, this diagnostic can be consumed directly by the Oceananigans TimeStepWizard callback, which rescales simulation.Δt each time it is invoked to target a chosen diffusive CFL number. The following examples use ForwardEuler explicitly, since an IMEX timestepper imposes no diffusive restriction from implicitly integrated terms (see Adaptive time stepping with IMEX). Attach it to a Simulation like any other callback,

integrator = initialize(SoilModel(grid; soil, timestepper = ForwardEuler(Float32)))
sim = Simulation(integrator; stop_iteration = 50, Δt = 60.0)

wizard = TimeStepWizard(diffusive_cfl = 0.5, max_change = 1.1)
sim.callbacks[:wizard] = Callback(wizard, IterationInterval(5))

run!(sim)
sim.Δt   # adapted time step (s)

171.187f0

or more succinctly via conjure_time_step_wizard!:

integrator = initialize(SoilModel(grid; soil, timestepper = ForwardEuler(Float32)))
sim = Simulation(integrator; stop_iteration = 50, Δt = 60.0)
conjure_time_step_wizard!(sim, show_progress = false)
run!(sim)
sim.Δt   # adapted time step (s)

572.59076f0

Oceananigans also defines a companion advective diagnostic cell_advection_timescale, which the wizard multiplies by its cfl (0.2 by default). For a soil with Richards-equation hydrology, it is the minimum of two time scales of explicit transport terms:

• the stability limit $C \Delta z / (C_\text{w} |j_\text{w}|)$ of the explicit upwind advection of heat by the Darcy flux $j_\text{w}$ (with bulk and water volumetric heat capacities $C$ and $C_\text{w}$), Inf without water flow;
• the explicit water outflow limit water_outflow_timescale, $\phi (\xi - \xi_\text{res}) \Delta z / q_\text{out}$, the time in which the explicit terms of the saturation tendency would drain a cell of porosity $\phi$ and saturation $\xi$ to its residual saturation $\xi_\text{res}$. The outflow $q_\text{out}$ (m/s) is the net loss of water of the cell per unit area through these terms: the Darcy fluxes that are integrated explicitly, the flux boundary conditions on the saturation (e.g. InfiltrationFlux), evapotranspiration, and the user-defined forcing. Cells without net explicit outflow impose no limit.

The wizard calls these diagnostics after each step, when the auxiliary variables have been updated but the tendencies are those of the beginning of the step. The outflow limit is therefore evaluated directly from the current state by kernel functions (compute_water_outflow_timescale) rather than read from the stored tendency, so that it measures the explicit outflow of the next step.

Adaptive time stepping with IMEX

Implicitly integrated operators impose no stability restriction, so the wizard is not throttled by them: under an IMEX timestepper, cell_diffusion_timescale is Inf for implicit heat conduction and implicit Richards flow, and the advective limit of the heat advection is Inf when internal_energy is implicit. The step is then bounded by the terms that remain explicit and by max_Δt. Without water flow, only max_Δt remains:

wizard_integrator = initialize(SoilModel(imex_grid; soil, timestepper = IMEX(ForwardEuler(Float32), BackwardEuler(Float32))))
sim = Simulation(wizard_integrator; stop_iteration = 50, Δt = 60.0)
wizard = TimeStepWizard(diffusive_cfl = 0.5, max_change = 1.5, max_Δt = 6 * 3600.0)
sim.callbacks[:wizard] = Callback(wizard, IterationInterval(1))
run!(sim)
sim.Δt   # reaches `max_Δt`, far beyond the explicit diffusive limit of a few minutes

21600.0f0

With a Richards-equation hydrology, the Darcy fluxes through the boundary faces of the column (e.g. gravity drainage through a FreeDrainage or unit-gradient condition on pressure_head), the flux boundary conditions, and the sinks stay explicit at $t^n$ (see Implicit Richards equation). A step longer than the outflow limit would remove more water from a cell than it holds, and the hydrology closure would then move water from other cells to fill the deficit, so the column alternates between drained and saturated cells. The outflow limit keeps the wizard below this bound. In a draining sandy column, it can be as restrictive as the explicit diffusive limit of the Richards equation, since the boundary drainage carries the full flux through the column; the cfl of the wizard sets the fraction of the available water that a step may remove. Under explicit time stepping of the saturation, the outflow limit includes the interior Darcy fluxes, but the diffusive limit is then usually the more restrictive one.

If internal_energy is routed explicitly while saturation_water_ice is implicit (see the warning in Coupling of water and energy), the advective limit of the explicit heat advection and the diffusive limit of explicit heat conduction both apply.

cell_diffusion_timescale(integrator::ModelIntegrator)

cell_advection_timescale(integrator::ModelIntegrator)

compute_thermal_diffusion_timescale(
    i,
    j,
    k,
    grid,
    fields,
    energy::SoilThermodynamics,
    hydrology::Terrarium.AbstractSoilHydrology,
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
)

compute_hydraulic_diffusion_timescale(
    i,
    j,
    k,
    grid,
    fields,
    hydrology::SoilHydrology{NF, RichardsEq, SaturationClosure, SoilHydraulics} where {SaturationClosure<:Terrarium.AbstractSoilWaterClosure, SoilHydraulics<:(Terrarium.AbstractSoilHydraulics{NF})},
    strat::Terrarium.AbstractStratigraphy,
    bgc::Terrarium.AbstractSoilBiogeochemistry
) where NF

Output writers

We can make use of Oceananigans.OutputWriters to flexibly save output to disk in a variety of formats.

The simplest choice output writer is JLD2Writer which saves selected Fields to a .jld2 file at a specified schedule:

using Oceananigans: JLD2Writer, TimeInterval
using Oceananigans.Units: hours

output_file = "$(tempname()).jld2"
println("Writing output to $(output_file)")

sim = Simulation(integrator; stop_time = 24*3600.0, Δt = 300.0)
sim.output_writers[:soil] = JLD2Writer(
    integrator,
    (temperature = integrator.state.temperature,
     saturation  = integrator.state.saturation_water_ice);
    filename = output_file,
    overwrite_files = true,
    schedule = TimeInterval(2hours),
)

# Re-initialize integrator
Terrarium.initialize!(integrator)

# Run the simulation again
run!(sim)

Writing output to /tmp/jl_H0bmJibO3u.jld2
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (1.102 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (20.258 ms).
[ Info: Simulation is stopping after running for 1.148 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.

The saved file can be read back as a FieldTimeSeries for post-processing:

using Oceananigans: FieldTimeSeries

temperature_series = FieldTimeSeries(output_file, "temperature")
temperature_series[end]  # Extract temperature Field at the last saved time

1×1×10 Field{Center, Center, Center} on RectilinearGrid on CPU
├── grid: 1×1×10 RectilinearGrid{Float32, Periodic, Flat, Bounded} on CPU with 1×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 3×1×16 StaticStridesOffsetArray(view(::Array{Float32, 4}, :, :, :, 13), 0:2, 1:1, -2:13) with eltype Float32 with indices 0:2×1:1×-2:13
    └── max=0.0, min=0.0, mean=0.0

Output writers can also accept subtypes AbstractSchedule via the schedule keyword. Schedules determine the frequency and aggregation of the output Fields:

Multiple output writers can be added to the same simulation, e.g. to save different variables at different frequencies:

sim.output_writers[:fast]   = JLD2Writer(integrator, (temperature = ...,); schedule = TimeInterval(1hours),  ...)
sim.output_writers[:daily]  = JLD2Writer(integrator, (pressure_head = ...,); schedule = AveragedTimeInterval(1days), ...)

References