Oceananigans.jl

Appearance

Time-stepping ​

Available time steppers ​

The TimeSteppers module provides three generic, explicit time-stepping schemes. Importantly, the module handles only the time integration of prognostic fields, the tendency computation is the responsibility of each model's implementation.

Time stepper availability in Oceananigans models

  • NonhydrostaticModel: QuasiAdamsBashforth2, RungeKutta3TimeStepper (default)

  • ShallowWaterModel: QuasiAdamsBashforth2TimeStepper, RungeKutta3 (default)

  • HydrostaticFreeSurfaceModel: QuasiAdamsBashforth2TimeStepper (default), SplitRungeKuttaTimeStepper

Quasi-Adams-Bashforth second order ​

The QuasiAdamsBashforth2TimeStepper approximates the time integral of tendencies via

where is a parameter. Ascher et al. (1995) suggests that   is optimal. With the additional parameter, the scheme is formally first-order accurate but offers improved stability properties. The default   provides a reasonable balance between accuracy and stability;   recovers the standard second-order Adams-Bashforth method.

The scheme requires storing tendencies from the previous time step, but only requires one tendency evaluation, making it computationally-efficient compared to multi-stage methods. The first time step automatically uses forward Euler (  ) since no previous tendencies exist.

Runge-Kutta third order ​

The RungeKutta3TimeStepper implements a low-storage, third-order Runge-Kutta scheme following Le and Moin (1991). The scheme advances the state through three substeps per time step:

with default coefficients  ,  ,  ,   , and    ( ).

This scheme requires three tendency evaluations per time step but provides higher accuracy and better stability for problems with oscillatory dynamics.

Split Runge-Kutta ​

The SplitRungeKuttaTimeStepper implements a Runge-Kutta scheme suitable for split-explicit computations that follows the implementation detailed by Wicker and Skamarock (2002). At the beginning of each time step the prognostic fields are cached, and subsequent substeps compute:

where is the cached initial state and are stage coefficients. The user can specify an arbitrary number of stages with custom coefficients. The default three-stage scheme uses  . This time stepper is used by HydrostaticFreeSurfaceModel for split-explicit treatment of the barotropic and baroclinic modes.

Usage ​

NonhydrostaticModel and ShallowWaterModel ​

The NonhydrostaticModel and ShallowWaterModel support two time steppers: :QuasiAdamsBashforth2 and :RungeKutta3 (default). The timestepper can be specified through a keyword argument in the model constructor using the symbols shown in the examples below or by building the timestepper manually.

julia
using Oceananigans

# Use QuasiAdamsBashforth2 time stepper
grid  = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; timestepper=:QuasiAdamsBashforth2)
model.timestepper

# output
QuasiAdamsBashforth2TimeStepper{Float64}
├── χ: 0.1
└── implicit_solver: nothing
julia
using Oceananigans

# Use RungeKutta3 time stepper (default)
grid  = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = NonhydrostaticModel(grid; timestepper=:RungeKutta3)
model.timestepper

# output
RungeKutta3TimeStepper{Float64}
├── γ: (0.5333333333333333, 0.4166666666666667, 0.75)
├── ζ: (-0.2833333333333333, -0.4166666666666667)
└── implicit_solver: nothing

HydrostaticFreeSurfaceModel ​

The HydrostaticFreeSurfaceModel supports :QuasiAdamsBashforth2 (default) and SplitRungeKuttaTimeStepper variants. The QuasiAdamsBashforth2 timestepper can be constructed as for the above or by manually building the timestepper and passing it to the model:

julia
using Oceananigans
using Oceananigans.TimeSteppers: QuasiAdamsBashforth2TimeStepper

# Use QuasiAdamsBashforth2 time stepper (default)
grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
timestepper = QuasiAdamsBashforth2TimeStepper(χ = 0.125)
model = HydrostaticFreeSurfaceModel(grid; timestepper)
model.timestepper

# output
QuasiAdamsBashforth2TimeStepper{Float64}
├── χ: 0.125
└── implicit_solver: nothing

The HydrostaticFreeSurfaceModel supports a SplitRungeKutta time stepper with 2 to 5 stages using convenience constructors timestepper = :SplitRungeKuttaN where N is the desired number of stages:

julia
using Oceananigans

# Construct models with convenience split Runge-Kutta constructor
grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
model = HydrostaticFreeSurfaceModel(grid; timestepper=:SplitRungeKutta2)
model = HydrostaticFreeSurfaceModel(grid; timestepper=:SplitRungeKutta3)
model = HydrostaticFreeSurfaceModel(grid; timestepper=:SplitRungeKutta4)
model = HydrostaticFreeSurfaceModel(grid; timestepper=:SplitRungeKutta5)
model.timestepper

# output
SplitRungeKuttaTimeStepper
├── stages: 5
├── β: (5, 4, 3, 2, 1)
└── implicit_solver: nothing

For custom coefficients, it is also possible to construct a SplitRungeKuttaTimeStepper directly:

julia
using Oceananigans
using Oceananigans.TimeSteppers: SplitRungeKuttaTimeStepper

grid = RectilinearGrid(size=(4, 4, 4), extent=(1, 1, 1))
timestepper = SplitRungeKuttaTimeStepper(coefficients=(4.3, 3.2, 2.1, 1.0))
model = HydrostaticFreeSurfaceModel(grid; timestepper)
model.timestepper

# output
SplitRungeKuttaTimeStepper
├── stages: 4
├── β: (4.3, 3.2, 2.1, 1.0)
└── implicit_solver: nothing

Extending time steppers for custom models ​

To use the existing time steppers with a new model type, implement the following methods. For QuasiAdamsBashforth2TimeStepper:

julia
ab2_step!(model::MyModel, Δt, callbacks)           # advance fields one AB2 step
cache_previous_tendencies!(model::MyModel)         # store Gⁿ → G⁻

For RungeKutta3TimeStepper:

julia
rk3_substep!(model::MyModel, Δt, γⁿ, ζⁿ, callbacks) # advance one RK3 substep
cache_previous_tendencies!(model::MyModel)          # store Gⁿ → G⁻

For SplitRungeKuttaTimeStepper:

julia
cache_current_fields!(model::MyModel)               # store U → Ψ⁻ at step start
rk_substep!(model::MyModel, Δτ, callbacks)          # advance one substep

All models must also implement update_state!(model, callbacks) to fill halo regions and compute any diagnostic quantities after each step or substep.

The fractional step method ​

With the pressure decomposition as discussed, the momentum evolves via:

where, e.g., for the non-hydrostatic model (ignoring background velocities and surface-wave effects)

collects all terms on the right side of the momentum equation , except the contribution of non-hydrostatic pressure .

The time-integral of the momentum equation from time step at   to time step   at is:

where the superscript and   imply evaluation at and , such that   . The crux of the fractional step method is to treat the pressure term implicitly using the approximation

while treating the rest of the terms on the right hand side of explicitly. The implicit treatment of pressure ensures that the velocity field obtained at time step   is divergence-free.

To effect such a fractional step method, we define an intermediate velocity field such that

The integral on the right of the equation for may be approximated by any of the time steppers described in Available time steppers.

Combining the equation for and the time integral of the non-hydrostatic pressure yields

Taking the divergence of fractional step equation and requiring that    yields a Poisson equation for the kinematic pressure at time-step  :

With in hand we can invert to get and then is computed via the fractional step equation .

Tracers are stepped forward explicitly via

where

and the same time-stepping scheme used for the momentum tendencies is applied to evaluate the integral of .