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.

Quasi-Adams-Bashforth second order

The QuasiAdamsBashforth2TimeStepper approximates the time integral of tendencies via

\[ \begin{equation} \label{eq:adams-bashforth} \int_{t_n}^{t_{n+1}} G \, \mathrm{d} t \approx \Delta t \left [ \left ( \tfrac{3}{2} + \chi \right ) G^n - \left ( \tfrac{1}{2} + \chi \right ) G^{n-1} \right ] \, , \end{equation}\]

where $\chi$ is a parameter. Ascher et al. (1995) suggests that $\chi = \tfrac{1}{8}$ is optimal. With the additional $\chi$ parameter, the scheme is formally first-order accurate but offers improved stability properties. The default $\chi = 0.1$ provides a reasonable balance between accuracy and stability; $\chi = 0$ 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 ($\chi = -1/2$) 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:

\[U^{m+1} = U^m + \Delta t \left( \gamma^m G^m + \zeta^m G^{m-1} \right)\]

with default coefficients $\gamma^1 = 8/15$, $\gamma^2 = 5/12$, $\gamma^3 = 3/4$, $\zeta^2 = -17/60$, and $\zeta^3 = -5/12$ ($\zeta^1 = 0$).

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:

\[U^{m+1} = U^0 + \frac{\Delta t}{\beta^m} G^m\]

where $U^0$ is the cached initial state and $\beta$ are stage coefficients. The user can specify an arbitrary number of stages with custom coefficients. The default three-stage scheme uses $\beta = (3, 2, 1)$. 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.

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

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:

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:

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:

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:

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

For RungeKutta3TimeStepper:

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

For SplitRungeKuttaTimeStepper:

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:

\[ \begin{equation} \label{eq:momentum-time-derivative} \partial_t \boldsymbol{v} = \boldsymbol{G}_{\boldsymbol{v}} - \boldsymbol{\nabla} p_{\rm{non}} \, , \end{equation}\]

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

\[\boldsymbol{G}_{\boldsymbol{v}} \equiv - \boldsymbol{\nabla}_h p_{\rm{hyd}} - \left ( \boldsymbol{v} \boldsymbol{\cdot} \boldsymbol{\nabla} \right ) \boldsymbol{v} - \boldsymbol{f} \times \boldsymbol{v} + \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{\tau} + \boldsymbol{F}_{\boldsymbol{v}}\]

collects all terms on the right side of the momentum equation \eqref{eq:momentum-time-derivative}, except the contribution of non-hydrostatic pressure $\boldsymbol{\nabla} p_{\rm{non}}$.

The time-integral of the momentum equation \eqref{eq:momentum-time-derivative} from time step $n$ at $t = t_n$ to time step $n+1$ at $t_{n+1}$ is:

\[ \begin{equation} \label{eq:momentum-time-integral} \boldsymbol{v}^{n+1} - \boldsymbol{v}^n = \int_{t_n}^{t_{n+1}} \Big [ - \boldsymbol{\nabla} p_{\rm{non}} + \boldsymbol{G}_{\boldsymbol{v}} \Big ] \, \mathrm{d} t \, , \end{equation}\]

where the superscript $n$ and $n+1$ imply evaluation at $t_n$ and $t_{n+1}$, such that $\boldsymbol{v}^n \equiv \boldsymbol{v}(t=t_n)$. The crux of the fractional step method is to treat the pressure term $\boldsymbol{\nabla} p_{\rm{non}}$ implicitly using the approximation

\[ \begin{align} \label{eq:pnon_implicit} \int_{t_n}^{t_{n+1}} \boldsymbol{\nabla} p_{\rm{non}} \, \mathrm{d} t \approx \Delta t \boldsymbol{\nabla} p_{\rm{non}}^{n+1} \, , \end{align}\]

while treating the rest of the terms on the right hand side of \eqref{eq:momentum-time-integral} explicitly. The implicit treatment of pressure ensures that the velocity field obtained at time step $n+1$ is divergence-free.

To effect such a fractional step method, we define an intermediate velocity field $\boldsymbol{v}^\star$ such that

\[ \begin{equation} \label{eq:intermediate-velocity-field} \boldsymbol{v}^\star - \boldsymbol{v}^n = \int_{t_n}^{t_{n+1}} \boldsymbol{G}_{\boldsymbol{v}} \, \mathrm{d} t \, , \end{equation}\]

The integral on the right of the equation for $\boldsymbol{v}^\star$ may be approximated by any of the time steppers described in Available time steppers.

Combining the equation \eqref{eq:intermediate-velocity-field} for $\boldsymbol{v}^\star$ and the time integral of the non-hydrostatic pressure \eqref{eq:pnon_implicit} yields

\[ \begin{equation} \label{eq:fractional-step} \boldsymbol{v}^{n+1} - \boldsymbol{v}^\star = - \Delta t \boldsymbol{\nabla} p_{\rm{non}}^{n+1} \, . \end{equation}\]

Taking the divergence of fractional step equation \eqref{eq:fractional-step} and requiring that $\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}^{n+1} = 0$ yields a Poisson equation for the kinematic pressure $p_{\rm{non}}$ at time-step $n+1$:

\[ \begin{equation} \label{eq:pressure-poisson} \nabla^2 p_{\rm{non}}^{n+1} = \frac{\boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{v}^{\star}}{\Delta t} \, . \end{equation}\]

With $\boldsymbol{v}^\star$ in hand we can invert \eqref{eq:pressure-poisson} to get $p_{\rm{non}}^{n+1}$ and then $\boldsymbol{v}^{n+1}$ is computed via the fractional step equation \eqref{eq:fractional-step}.

Tracers are stepped forward explicitly via

\[ \begin{equation} \label{eq:tracer-timestep} c^{n+1} - c^n = \int_{t_n}^{t_{n+1}} G_c \, \mathrm{d} t \, , \end{equation}\]

where

\[ \begin{equation} G_c \equiv - \boldsymbol{\nabla} \boldsymbol{\cdot} \left ( \boldsymbol{v} c \right ) - \boldsymbol{\nabla} \boldsymbol{\cdot} \boldsymbol{q}_c + F_c \, , \end{equation}\]

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