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

Terrarium.initialize — Method

initialize(
    model::Terrarium.AbstractModel{NF, Grid} where Grid<:(Terrarium.AbstractLandGrid{NF});
    ...
) -> ModelIntegrator{_A, _B, _C, _D, var"#s179", StateVariables{NF, prognames, closurenames, auxnames, inputnames, nsnames, ProgFields, TendFields, AuxFields, InputFields, Namespaces, Cache, ClockType}, Clock{_A1, _B1, Float64, Int64, Int64}, @NamedTuple{}, InputSources{NF1, name, Sources}} where {_A, _B<:Oceananigans.Architectures.AbstractArchitecture, _C<:(Terrarium.AbstractLandGrid{_A}), _D<:Terrarium.AbstractTimeStepper{_A}, var"#s179"<:Terrarium.AbstractModel{_A, _C}, NF, prognames, closurenames, auxnames, inputnames, nsnames, ProgFields, TendFields, AuxFields, InputFields, Namespaces, Cache, ClockType, _A1, _B1, NF1, name, Sources<:Tuple{Vararg{InputSource{NF1}}}}
initialize(
    model::Terrarium.AbstractModel{NF, Grid} where Grid<:(Terrarium.AbstractLandGrid{NF}),
    params;
    clock,
    inputs,
    boundary_conditions,
    initializers,
    fields
) -> ModelIntegrator{_A, _B, _C, _D, var"#s179", StateVariables{NF, prognames, closurenames, auxnames, inputnames, nsnames, ProgFields, TendFields, AuxFields, InputFields, Namespaces, Cache, ClockType}, Clock{_A1, _B1, Float64, Int64, Int64}, @NamedTuple{}, InputSources{NF1, name, Sources}} where {_A, _B<:Oceananigans.Architectures.AbstractArchitecture, _C<:(Terrarium.AbstractLandGrid{_A}), _D<:Terrarium.AbstractTimeStepper{_A}, var"#s179"<:Terrarium.AbstractModel{_A, _C}, NF, prognames, closurenames, auxnames, inputnames, nsnames, ProgFields, TendFields, AuxFields, InputFields, Namespaces, Cache, ClockType, _A1, _B1, NF1, name, Sources<:Tuple{Vararg{InputSource{NF1}}}}

Creates and initializes a ModelIntegrator for the given model with input variables populated by the given inputs and optionally params . InputSources can be specified via the inputs keyword argument. This method allocates all necessary Fields for the state variables and subsequently calls initialize!(::ModelIntegrator).

Note that this method is not type stable and thus should not be called from Enzyme autodiff. To reinitialize the model for an existing state, use initialize!(state, model).

See the docstring for initialize(::AbstractModel) for further details.

source

This will return a ModelIntegrator:

Terrarium.ModelIntegrator — Type

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

Represents a "integrator" for a simulation of a given model. ModelIntegrator consists of a clock, a model, and an initialized StateVariables data structure, as well as any relevant inputs provided by a corresponding InputProvider. The ModelIntegrator implements the Oceananigans.AbstractModel interface and can thus be treated as a "model" in Oceananigans Simulations and output reading/writing utilities.

source

A single iteration of the ModelIntegrator roughly involves four steps:

Invoke update_inputs! to populate all input Fields from their respective InputSources based on the current clock time,
Invoke compute_auxiliary! on all model components to derive auxiliary state variables from the current prognostic state,
Invoke compute_tendencies! on all model components to calculate tendencies for all prognostic variables,
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))
model = SoilModel(grid, timestepper=ForwardEuler(Float32))
integrator = initialize(model)

Integrator of SoilModel{Float32, ColumnGrid{Float32, CPU, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}}, SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, Terrarium.ExplicitTwoPhaseHeatConduction, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, NoFlow, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKLinear{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 Oceananigans.Grids.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 OffsetArray(::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). Terrarium currently provides two explicit time steppers, ForwardEuler and Heun:

Terrarium.ForwardEuler — Type

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

Simple forward Euler time stepping scheme.

source

Terrarium.Heun — Type

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

Simple forward 2nd order Heun / improved Euler time stepping scheme.

source

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!.

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, ColumnGrid{Float32, CPU, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}}, SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, Terrarium.ExplicitTwoPhaseHeatConduction, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, NoFlow, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKLinear{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 (81.121 μs)
[ 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 (123.016 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, ColumnGrid{Float32, CPU, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}}, ForwardEuler{Float32}, SoilModel{Float32, ColumnGrid{Float32, CPU, Oceananigans.Grids.RectilinearGrid{Float32, Oceananigans.Grids.Periodic, Oceananigans.Grids.Flat, Oceananigans.Grids.Bounded, Oceananigans.Grids.StaticVerticalDiscretization{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}}, SoilEnergyWaterCarbon{Float32, SoilStratigraphy{Float32, 1, Tuple{ConstantSoilHorizon{Float32, :soil, ConstantSoilPorosity{Float32}}}}, SoilThermodynamics{Float32, Terrarium.ExplicitTwoPhaseHeatConduction, SoilEnergyTemperatureClosure, SoilThermalProperties{Float32, FreeWater, InverseQuadratic, SoilThermalConductivities{Float32}}}, SoilHydrology{Float32, NoFlow, SoilSaturationPressureClosure, SoilHydraulicsSURFEX{Float32, BrooksCorey{FreezeCurves.SoilWaterVolume{Float32, Float32, Float32}, Float32, Float32}, UnsatKLinear{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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, UnitRange{Int64}}, OffsetArrays.OffsetArray{Float32, 3, 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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 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, 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{OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}, OffsetArrays.OffsetVector{Float32, Vector{Float32}}}, Float32, Float32, OffsetArrays.OffsetVector{Float32, StepRangeLen{Float32, Float64, Float64, Int64}}, Nothing, CPU, Oceananigans.Grids.GridSize{1, 1, 10, 1, 0, 3}}, Tuple{Colon, Colon, Colon}, OffsetArrays.OffsetArray{Float32, 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{}}}
├── 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. Attach it to a Simulation like any other callback,

integrator = initialize(SoilModel(grid))
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))
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 is relevant for processes involving advection.

Richards equation is stiff near saturation

As the soil saturates, $\partial\theta/\partial\psi \to 0$ and the hydraulic timescale tends to zero. Explicit integration of the Richardson–Richards equation is therefore extremely stiff in the saturated limit, which is the regime that motivates implicit time stepping (planned). The TimeStepWizard's min_change bound prevents Δt from collapsing to zero in a single adjustment, but a saturated explicit run may still require very small steps.

Reactant

Adaptive stepping applies to host-driven Simulations only. Under a ReactantState grid the time step is compiled into the traced run! loop, so a host-side wizard callback does not take effect there.

Oceananigans.Diagnostics.cell_diffusion_timescale — Method

cell_diffusion_timescale(integrator::ModelIntegrator) -> Any

Return the minimum diffusive stability timescale $τ = Δz² / D$ (seconds) over all grid cells of the integrator, where D is the largest effective diffusivity among the model's diffusive processes. This is the diagnostic consumed by the Oceananigans TimeStepWizard to enforce a diffusive Courant–Friedrichs–Lewy (CFL) constraint of the form $Δt ≈ \mathrm{diffusive\_cfl} · τ$.

The generic fallback (for models without an implemented diffusive timescale) returns Inf, i.e. no diffusive restriction, mirroring the Oceananigans infinite_diffusion_timescale convention.

source

Oceananigans.Advection.cell_advection_timescale — Method

cell_advection_timescale(
    integrator::ModelIntegrator{NF, Arch, Grid, TimeStepper, Model} where {Arch<:Oceananigans.Architectures.AbstractArchitecture, Grid<:(Terrarium.AbstractLandGrid{NF}), TimeStepper<:Terrarium.AbstractTimeStepper{NF}, Model<:Terrarium.AbstractModel{NF, Grid}}
) -> Any

Return the advective stability timescale of the integrator. Land models currently transport no quantity advectively, so this is always Inf (no advective restriction); it exists only so that a plain TimeStepWizard (which always evaluates an advective timescale) can be applied to a Terrarium Simulation without providing a custom cell_advection_timescale.

source

Terrarium.compute_thermal_diffusion_timescale — Function

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

Kernel function returning the thermal diffusion timescale $Δz² C / κ$ at cell i, j, k, with the bulk thermal conductivity κ and heat capacity C computed from the local soil composition.

source

Terrarium.compute_hydraulic_diffusion_timescale — Function

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

Kernel function returning the hydraulic diffusion timescale $Δz² (∂θ/∂ψ) / K$ at cell i, j, k. The hydraulic conductivity K is evaluated at the cell centre from the local soil composition; the specific moisture capacity ∂θ/∂ψ is the analytic derivative of the soil-water retention curve (from FreezeCurves) evaluated at the matric potential ψₘ, reconstructed from the stored total pressure_head by removing its hydrostatic and elevation components. Cells with zero hydraulic conductivity (fully dry or frozen) impose no restriction and return Inf.

source

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_existing = true,
    including = [:grid], # save the grid with the output
    schedule = TimeInterval(2hours),
)

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

# Run the simulation again
run!(sim)

Writing output to /tmp/jl_uZqTGRGfSj.jld2
[ Info: Initializing simulation...
┌ Warning: error ArgumentError("a group or dataset named grid is already present within this group") thrown when trying to serialize the grid at serialized/grid
└ @ Oceananigans.OutputWriters ~/.julia/packages/Oceananigans/RPT5V/src/OutputWriters/jld2_writer.jl:251
[ Info:     ... simulation initialization complete (4.265 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (43.542 ms).
[ Info: Simulation is stopping after running for 4.341 seconds.
[ Info: Simulation time 1 day equals or exceeds stop time 1 day.

Note

Always call Terrarium.initialize!(integrator) before calling run!(sim) when re-running a simulation; otherwise, the simulation will not run due to the stopping condition already being satisfied.

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 Oceananigans.Grids.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 OffsetArray(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:

Schedule type	Description
`TimeInterval(Δt)`	Write every `Δt` seconds of simulation time
`IterationInterval(n)`	Write every `n` timesteps
`AveragedTimeInterval(Δt)`	Write time-averaged outputs over windows of `Δt` seconds

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