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[^Hotta2016]; 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. 2. an Euler forward step with , then

3. one leapfrog time step with , then

4. 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	time step at
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[^Robert66][^Asselin72] 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[^Williams2009] then proposed an additional filter step to regain accuracy that is otherwise lost with a strong Robert-Asselin filter[^Amezcua2011][^Williams2011]. 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):

prognostic_spectral_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):

prognostic_steps(::NCycleLorenz) = 1
tendency_grid_steps(::NCycleLorenz) = 1
tendency_spectral_steps(::NCycleLorenz) = 2

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:

# 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:

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:

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

get_step(var, i) is the low-level accessor used by all of the above; for a variable with a step dimension it returns a view of step i.

SpeedyWeather.get_step Function

get_step(var) -> Any

Select step dimension from variable, when no step as 2nd argument provided select las index as this typically presents the "current" step (and not any previous ones). But this depends on the time stepping a variable with step dimension was created for.

source

get_step(
    var::LowerTriangularArray{T, 2, ArrayType} where ArrayType<:AbstractArray{T, 2},
    step::Integer
) -> LowerTriangularArray

Get the i-th step of a LowerTriangularArray as a view (wrapped into a LowerTriangularArray). "step" refers to the last dimension, for prognostic variables e.g. used for the leapfrog time step. This method is for a 2D spectral variable (horizontal only) with steps in the 3rd dimension.

source

get_step(
    var::LowerTriangularArray{T, 3, ArrayType} where ArrayType<:AbstractArray{T, 3},
    step::Integer
) -> LowerTriangularArray

Get the i-th step of a LowerTriangularArray as a view (wrapped into a LowerTriangularArray). "step" refers to the last dimension, for prognostic variables e.g. used for the leapfrog time step. This method is for a 3D spectral variable (horizontal + vertical) with steps in the 4rd dimension.

source

get_step(var::AbstractField{T, 2}, step::Integer) -> Any

Get the i-th step of a 3D field as a view (wrapped into the same type as the input variable). "step" refers to the last dimension, for prognostic variables e.g. used for the leapfrog time step. This method is for a 2D field (horizontal only) with steps in the 3rd dimension.

source

get_step(var::AbstractField{T, 3}, step::Integer) -> Any

Get the i-th step of a 4D field as a view (wrapped into the same type as the input variable). "step" refers to the last dimension, for prognostic variables e.g. used for the leapfrog time step. This method is for a 3D field (horizontal + vertical) with steps in the 4rd dimension.

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[^Hotta2016]).

Lorenz N-cycle

The Lorenz N-cycle NCycleLorenz is a semi-implicit time integration following Hotta et al. (2016)[^Hotta2016]. 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    ( )	Steps per period	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

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

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

Algorithm (per substep of a cycle) 2. G = w_F_E(x) + (1-w)_G (weighted tendency accumulation)

3. dx = (I - α_Δt_L_I)^(-1) * (G + L_I*x) (implicit solve)

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

source

SpeedyWeather.NCycleLorenzA Type

Version A: weights w_k = N/(N-k) for k=1,...,N-1; w_0 = 1

source

SpeedyWeather.NCycleLorenzB Type

Version B: weights w_k = N/k for k=1,...,N-1; w_0 = 1

source

SpeedyWeather.NCycleLorenzAB Type

Version AB: alternates A and B every N steps

source

SpeedyWeather.NCycleLorenzABBA Type

Version ABBA: uses A-B-B-A sequence (only for N=4, provides 4th order accuracy)

source

References

[^Robert66]: Robert, André. "The Integration of a Low Order Spectral Form of the Primitive Meteorological Equations." Journal of the Meteorological Society of Japan 44 (1966): 237-245.

[^Asselin72]: ASSELIN, R., 1972: Frequency Filter for Time Integrations. Mon. Wea. Rev., 100, 487-490, doi:10.1175/1520-0493(1972)100<0487:FFFTI>2.3.CO;2

[^Williams2009]: Williams, P. D., 2009: A Proposed Modification to the Robert-Asselin Time Filter. Mon. Wea. Rev., 137, 2538-2546, 10.1175/2009MWR2724.1.

[^Amezcua2011]: Amezcua, J., E. Kalnay, and P. D. Williams, 2011: The Effects of the RAW Filter on the Climatology and Forecast Skill of the SPEEDY Model. Mon. Wea. Rev., 139, 608-619, doi:10.1175/2010MWR3530.1.

[^Williams2011]: Williams, P. D., 2011: The RAW Filter: An Improvement to the Robert-Asselin Filter in Semi-Implicit Integrations. Mon. Wea. Rev., 139, 1996-2007, doi:10.1175/2010MWR3601.1.

[^Hotta2016]: Hotta, D., E. Kalnay, and P. Ullrich, 2016: A Semi-Implicit Modification to the Lorenz N-Cycle Scheme and Its Application for Integration of Meteorological Equations. Mon. Wea. Rev., 144, 2215-2233, doi:10.1175/MWR-D-15-0330.1.