Time stepping

This page describes the mathematics and internal implementation of SpeedyWeather.jl's time integration schemes. For how to create and configure a time stepper and pass it to the model constructor, see Time integration in the Usage section.

SpeedyWeather.jl supports several time integration schemes, selected by passing a time_stepping component to the model constructor:

Leapfrog, a 2-step leapfrog scheme with a Robert-Asselin and Williams filter (the default for ShallowWaterModel, PrimitiveDryModel and PrimitiveWetModel), described below.
NCycleLorenz, a family of semi-implicit Lorenz N-cycle schemes ((Hotta et al., 2016); the default for BarotropicModel).

All schemes share a common framework (see Time steppers and variable steps) in which the time stepper decides, for every model component, which stored step of each variable to read or write. This decouples the dynamical core and parameterizations from the time-stepping bookkeeping: e.g. leapfrog stores two steps of the prognostic variables, while the Lorenz N-cycle stores only one but keeps a second tendency for its weighted accumulation.

Leapfrog

SpeedyWeather.jl's default time integration is the Leapfrog time integration, which, for relative vorticity , is in its simplest form

meaning we step from the previous time step , leapfrogging over the current time step to the next time step by evaluating the tendencies on the right-hand side at the current time step . The time stepping is done in spectral space (see LowerTriangularArrays for how spectral coefficients are stored). Once the right-hand side is evaluated, leapfrogging is a linear operation, meaning that its simply applied to every spectral coefficient as one would evaluate it on every grid point in grid-point models.

For the Leapfrog time integration two time steps of the prognostic variables have to be stored, and . Time step is used to evaluate the tendencies which are then added to in a step that also swaps the indices for the next time step and , so that no additional memory than two time steps have to be stored at the same time.

Leapfrog initialisation

The Leapfrog time integration has to be initialized with an Euler forward step in order to have a second time step available when starting from to actually leapfrog over. SpeedyWeather.jl therefore does two initial time steps that are different from the leapfrog time steps that follow and that have been described above.

an Euler forward step with , then
one leapfrog time step with , then
leapfrog with till the end

This is particularly done in a way that after 2. we have at and at available so that 3. can start the leapfrogging without any offset from the intuitive spacing . The following schematic can be useful

	time at step	time at step
Initial conditions
1: Euler	(T)
2: Leapfrog with		(T)
3 to : Leapfrog with		(T)

The time step that is used to evaluate the tendencies is denoted with (T). It is always the time step furthest in time that is available.

Before the time integration starts, the initial conditions – which only occupy step – are also copied into step . This is purely a bookkeeping convenience: it lets which_prognostic_step always return step (2) as the one to evaluate the right-hand side on, without a special case for the very first (Euler) step where step would otherwise not yet exist.

The initial Euler step is not filtered, see Robert-Asselin and Williams filter below. On the first two steps (Euler, and the first leapfrog step) Leapfrog's update_prognostic! therefore disables the filter weights (both are 0), and enables them for every step after that.

Robert-Asselin and Williams filter

The standard leapfrog time integration is often combined with a Robert-Asselin filter (Robert, 1966; Asselin, 1972) to dampen a computational mode. The idea is to start with a standard leapfrog step to obtain the next time step but then to correct the current time step by applying a filter which dampens the computational mode. The filter looks like a discrete Laplacian in time with a stencil, and so, maybe unsurprisingly, is efficient to filter out a "grid-scale oscillation" in time, aka the computational mode. Let be the unfiltered variable and be the filtered variable, the right-hand side tendency, then the standard leapfrog step is

Meaning we start with a filtered variable at the previous time step , evaluate the tendency based on the current time step to obtain an unfiltered next time step . We then filter the current time step (which will become on the next iteration)

by adding a discrete Laplacian with coefficient to it, evaluated from the available filtered and unfiltered time steps centred around : is not available anymore because it was overwritten by the filtering at the previous iteration, are not filtered yet when applying the Laplacian. The filter parameter is typically chosen between 0.01-0.2, with stronger filtering for higher values.

Williams (2009) then proposed an additional filter step to regain accuracy that is otherwise lost with a strong Robert-Asselin filter (Amezcua et al., 2011; Williams, 2011). Now let be unfiltered, be once filtered, and twice filtered, then

with the Williams filter parameter . For we're back with the Robert-Asselin filter (the first two lines).

The Laplacian in the parentheses is often called a displacement, meaning that the filtered value is displaced (or corrected) in the direction of the two surrounding time steps. The Williams filter now also applies the same displacement, but in the opposite direction to the next time step as a correction step (line 3 above) for a once-filtered value which will then be twice-filtered by the Robert-Asselin filter on the next iteration. For more details see the referenced publications.

The initial Euler step (see Leapfrog initialisation) is not filtered. Both the the Robert-Asselin and Williams filter are then switched on for all following leapfrog time steps.

Implementation-wise, update_prognostic! for Leapfrog folds both filters into two weights and (both zero on the unfiltered initial steps) and applies the leapfrog step, the Robert-Asselin and the Williams filter in a single kernel per spectral coefficient (and vertical layer).

Time steppers and variable steps

Different time integration schemes need to store a different number of past states of the prognostic variables and/or tendencies. Leapfrog needs the two spectral steps and ; the Lorenz N-cycle needs only one prognostic step but a second tendency to accumulate weighted tendencies. SpeedyWeather handles this generically: prognostic variables and tendencies carry an extra Step dimension (the last dimension of their underlying array), and the time stepper decides how many steps to allocate and which step each model component should read or write.

The number of steps is requested by the time stepper through prognostic_steps and tendency_steps (with prognostic_grid_steps, prognostic_spectral_steps, tendency_grid_steps, tendency_spectral_steps to distinguish grid/spectral and dispatch over the model). For example leapfrog requests two spectral prognostic steps but, in the primitive-equation models, also two grid steps (because the parameterizations are evaluated at the previous grid state):

julia

prognostic_steps(::AbstractLeapfrog) = 2
prognostic_grid_steps(::AbstractLeapfrog, ::PrimitiveEquation) = 2
tendency_steps(::AbstractLeapfrog) = 1

The Lorenz N-cycle, by contrast, only needs one prognostic step (the state is updated in place) but two spectral tendency steps – one for the explicit tendency , one for the weighted accumulation – and only one grid tendency step ( is only ever needed in spectral space):

julia

prognostic_steps(::NCycleLorenz) = 1
tendency_steps(::NCycleLorenz) = 2          # also for ocean and land variables, stepped directly on the grid
tendency_grid_steps(::NCycleLorenz) = 1     # atmospheric grid tendencies only hold F

Throughout the dynamical core and parameterizations a variable is then accessed with get_prognostic_step / get_tendency_step, which return a view of the appropriate step:

julia

# in a model component, e.g. the spectral→grid transform or a tendency term
vor  = get_prognostic_step(vars.prognostic.vorticity, time_stepping, component)
ζtend = get_tendency_step(vars.tendencies.vorticity,  time_stepping, component)

Which step is returned is decided by the time stepper via which_prognostic_step / which_tendency_step, dispatched on the variable, the time stepper, the component and (optionally) the model — so a scheme can choose a different step per process. The fallback is step 1, and leapfrog for instance overrides it to read the current (2nd) step for transforms and the nonlinear dynamical core, but the previous (1st) step for the linear terms and horizontal diffusion:

julia

which_prognostic_step(var, ::AbstractLeapfrog, ::AbstractSpectralTransform)       = 2  # current
which_prognostic_step(var, ::AbstractLeapfrog, ::LinearDynamicalCore)             = 1  # previous
which_prognostic_step(var, ::AbstractLeapfrog, ::AbstractHorizontalDiffusion)     = 1  # previous

For the Lorenz N-cycle, only the explicit tendency is reset to zero every substep for accumulation into ; the fallback step 1 covers this and every other component that only ever needs one tendency step:

julia

which_tendency_step(var, ::AbstractNCycleLorenz, ::ResetTendencies) = 1

get_step(var, i) is the low-level accessor used by all of the above. Which array dimension is the step dimension is taken from the variable's dimension tag (see Array dimensions): get_step selects the dimension tagged as time T, so for a variable tagged XYT, XYZT, LMT or LMZT it returns a view of step i. A variable without a time tag (XY, XYZ, LM, LMZ) has no step dimension to select from — the full variable is returned as a view and i is ignored. Model components can therefore call it unconditionally, whether or not the variable they are handed has steps. get_steps(var) returns all steps as a tuple of views (a 1-tuple for a variable without a step dimension), and get_steps(var, Val(N)) does the same with a compile-time length, which is what Enzyme needs.

Dispatching on whether a variable has a step dimension at all is done with the AbstractFieldWithTime / LowerTriangularArrayWithTime aliases, as nsteps does:

julia

nsteps(var::Union{AbstractFieldWithTime, LowerTriangularArrayWithTime}) = size(var, ndims(var))
nsteps(::Union{AbstractField, LowerTriangularArray}) = 1

SpeedyWeather.get_step Function

julia

get_step(var) -> Any

Select the step dimension from a variable, when no step is provided as 2nd argument select the last step as this typically represents the "current" step (and not any previous ones). But this depends on the time stepping a variable with step dimension was created for. Variables without a step dimension have a single step, so get_step(var) returns the full variable as a view.

source

julia

get_step(
    var::LowerTriangularArray{T, 1, A, S, SpeedyWeatherInternals.ArrayDimensions.LM},
    step::Integer
) -> Union{LowerTriangularArray{_A, _B, <:AbstractArray{_A, _B}, <:AbstractSpectrum, SpeedyWeatherInternals.ArrayDimensions.LM} where _B, LowerTriangularArray{_A, N, <:AbstractArray{_A, N}, <:AbstractSpectrum, SpeedyWeatherInternals.ArrayDimensions.LM} where N} where _A

Get the i-th step of a LowerTriangularArray as a view (wrapped into a LowerTriangularArray). "step" refers to the time dimension T in the variable's ArrayDimensions, for prognostic variables e.g. used for the leapfrog time step. These methods are for spectral variables WITHOUT a time dimension (LM, LMZ), which have no step dimension to select from: the full array is returned as a view and step is ignored (it should be 1).

source

julia

get_step(
    var::LowerTriangularArray{T, 2, A, S, SpeedyWeatherInternals.ArrayDimensions.LMT},
    step::Integer
) -> Union{LowerTriangularArray{_A, _B, <:AbstractArray{_A, _B}, <:AbstractSpectrum, SpeedyWeatherInternals.ArrayDimensions.LM} where _B, LowerTriangularArray{_A, N, <:AbstractArray{_A, N}, <:AbstractSpectrum, SpeedyWeatherInternals.ArrayDimensions.LM} where N} where _A

Get the i-th step of a LowerTriangularArray as a view (wrapped into a LowerTriangularArray). "step" refers to the time dimension T in the variable's ArrayDimensions, for prognostic variables e.g. used for the leapfrog time step. These methods are for spectral variables WITH a time dimension (LMT horizontal + time, LMZT horizontal + vertical + time), whose last dimension is the step dimension.

source

julia

get_step(
    var::AbstractField{T, 1, A, G, SpeedyWeatherInternals.ArrayDimensions.XY},
    step::Integer
) -> Any

Get the i-th step of a field as a view (wrapped into the same type as the input variable). "step" refers to the time dimension T in the variable's ArrayDimensions, for prognostic variables e.g. used for the leapfrog time step. These methods are for fields WITHOUT a time dimension (XY, XYZ), which have no step dimension to select from: the full array is returned as a view and step is ignored (it should be 1).

source

julia

get_step(
    var::AbstractField{T, 1, A, G, SpeedyWeatherInternals.ArrayDimensions.XYT},
    step::Integer
) -> Any

Get the i-th step of a field as a view (wrapped into the same type as the input variable). "step" refers to the time dimension T in the variable's ArrayDimensions, for prognostic variables e.g. used for the leapfrog time step. These methods are for fields WITH a time dimension (XYT horizontal + time, XYZT horizontal + vertical + time), whose last dimension is the step dimension.

source

SpeedyWeather.get_steps Function

Get all steps of a variable as a tuple of views (wrapped into the same type as the input variable) as defined by get_step. "Steps" refer to the step dimension, i.e. the time dimension T in the variable's ArrayDimensions (e.g. XYT, LMZT), for prognostic variables e.g. used for the leapfrog time step. Variables without a time dimension have no steps to splat, so a 1-tuple holding the full variable (as a view, consistent with get_step) is returned.

source

When writing a new time stepper you implement the *_steps methods (how many steps to store), the which_*_step methods (which step each component reads/writes) and an update_prognostic! method (how a tendency advances the state); the rest of the model is agnostic to the scheme.

Time steppers can also reorder how diffusion and the semi-implicit correction are applied via diffusion_and_implicit!: Leapfrog applies the implicit correction first, then horizontal diffusion, while the Lorenz N-cycle applies diffusion first, then the implicit correction (consistent with the ordering in Hotta et al. (2016)).

Ocean, sea ice and land

Dynamic ocean, sea ice and land components (e.g. SlabOcean, ThermodynamicSeaIce, LandBucketTemperature, LandBucketMoisture, SnowModel) write tendencies into vars.tendencies.ocean and vars.tendencies.land. The variables in these namespaces are time stepped with their own time stepper within the time stepping setup model.time_stepping, SpeedyWeather.time_stepper(model.time_stepping, :ocean) and (model.time_stepping, :land) (sea ice lives in the :ocean namespace). By default this is the time stepper of the atmosphere, but Leapfrog has ocean and land fields which default to EulerForward, which steps the variables with Euler forward over :

julia

time_stepping = Leapfrog(spectral_grid)                             # Euler forward for ocean and land
time_stepping = Leapfrog(spectral_grid, land = nothing)             # leapfrog land as the atmosphere
model = PrimitiveWetModel(spectral_grid; time_stepping)

NCycleLorenz has the same ocean and land fields, defaulting to nothing, i.e. the same N-cycle as the atmosphere. Any time stepper can be used for ocean and land, e.g. Leapfrog(spectral_grid, land = NCycleLorenz(spectral_grid)).

The variables are allocated with the number of steps of their time stepper (prognostic_steps and tendency_steps, see SpeedyWeather.get_nsteps(time_stepping, :land)), so with EulerForward ocean and land variables have only one step, with NCycleLorenz one prognostic step and two tendency steps. The time step of ocean and land is always set to the of the leapfrog (also with set!(model, Δt=...)), and on the first two leapfrog steps, which each only advance the clock by , also Euler steps with .

Why not leapfrog? The tendencies are computed from surface fluxes (sensible and latent heat, longwave radiation) that the parameterizations evaluate at the previous (1st) step, and they contain stiff relaxation terms, e.g. the sensible heat flux feedback on the thin top soil layer has timescales of about an hour in strong winds. Leapfrogging those over from is an Euler step of , which halves the stability limit compared to an Euler step of and made long integrations unstable. Terms that are meant to act "within one time step" (e.g. restoring sea surface temperatures to freezing when sea ice forms) use SpeedyWeather.time_step(model.time_stepping, :ocean, clock) for the same reason.

SpeedyWeather.EulerForward Type

Euler forward time stepping, used by default for the ocean and land namespaces with Leapfrog (see its ocean and land fields). Defined by the following fields

Δt_at_T32::Any: [OPTION] Time step for T32, scale linearly to spectral resolution truncation
adjust_with_output::Any: [OPTION] Adjust Δt_at_T32 with the output interval to output exactly after integer time steps
Δt_millisec::Any: [DERIVED] Time step Δt in milliseconds at specified resolution
Δt::Any: [DERIVED] Time step Δt [s] at specified resolution

source

Lorenz N-cycle

The Lorenz N-cycle NCycleLorenz is a semi-implicit time integration following Hotta et al. (2016). Over a cycle of substeps it advances the state , per substep, as

where are the explicit tendencies, the linear (gravity-wave) terms treated implicitly, and a substep-dependent weight coefficient. Only one prognostic step is stored (the state is updated in place), but two tendencies are kept: (the current tendency) and (the weighted accumulation that carries memory across substeps), hence tendency_spectral_steps(::NCycleLorenz) = 2. In SpeedyWeather's implementation the implicit solve is folded into the same implicit_correction! step used by Leapfrog (see Semi-implicit time stepping for the primitive equations and Semi-implicit time integration for the shallow water model – both implicit solves are shared between time steppers), so update_prognostic! for NCycleLorenz itself only performs the weighted accumulation and the explicit state update

for every spectral coefficient (and vertical layer), with the implicit correction applied to beforehand by diffusion_and_implicit!.

Four weight variants are available, following the naming in Hotta et al. (2016), selected via the variant option:

Variant	Weight for substep ()	Notes
`NCycleLorenzA` (default)
`NCycleLorenzB`
`NCycleLorenzAB`	alternates a full cycle of `A`, then a full cycle of `B`	2 cycles per period
`NCycleLorenzABBA`	sequence `A`, `B`, `B`, `A`	4th-order accurate for ; other run but without that accuracy guarantee

"Steps per period" is the number of substeps after which the weight sequence repeats. These substeps can be interpreted in two ways: as sub-stages of a single step, the way a Runge-Kutta method would use them, where the full time step is only completed once all substeps of a period have run; or as a multi-step method with changing weights that restarts every period, where every substep already advances the state by a full (as in the update above). SpeedyWeather.jl uses the latter, so the computational cost of the Lorenz N-cycle is independent of but, like leapfrog, inversely proportional to .

The cycle length is set with steps (3 or 4 recommended, 4 is more stable). The current substep is current_substep(L, clock) = mod(clock.step_counter, L.steps), and the weight for that substep is computed by weight_coefficient, dispatched on the variant.

SpeedyWeather.NCycleLorenz Type

julia

NCycleLorenz{NF, V, ...} <: AbstractTimeStepper

A semi-implicit Lorenz N-cycle time integration scheme following Hotta et al. (2016).

Algorithm (per substep of a cycle)

G = w_F_E(x) + (1-w)_G (weighted tendency accumulation)
dx = (I - α_Δt_L_I)^(-1) * (G + L_I*x) (implicit solve)
x = x + Δt*dx (state update)

steps::Any: [OPTION] Number of steps N in a cycle (3 or 4 recommended, 4 is more stable)
variant::Any: [OPTION] Variant: NCycleLorenzA() (default), B, AB, or ABBA
Δt_at_T32::Any: [OPTION] Time step for T32 resolution, scale linearly with resolution
adjust_with_output::Any: [OPTION] Adjust Δt_at_T32 with the interval to reach interval exactly
Δt_millisec::Any: [DERIVED] Time step Δt in milliseconds at specified resolution
Δt::Any: [DERIVED] Time step Δt [s] at specified resolution
ocean::Any: [OPTION] Time stepper for the ocean (incl sea ice) variables, nothing for the same N-cycle as the atmosphere
land::Any: [OPTION] Time stepper for the land variables, nothing for the same N-cycle as the atmosphere