Time-stepping

FourierFlows.jl includes several time-stepping algorithms.

Most of the time-stepping algorithms are fully explicit schemes: ForwardEulerTimeStepper, AB3TimeStepper, RK4TimeStepper, and LSRK54TimeStepper but also implemented is the ETDRK4TimeStepper.

High-wavenumber filtering

Most of the time steppers also come with their Filtered equivalents: FilteredForwardEulerTimeStepper, FilteredAB3TimeStepper, FilteredRK4TimeStepper, FilteredLSRK54TimeStepper, and FilteredETDRK4TimeStepper.

The filtered time steppers apply a high-wavenumber filter to the solution at the end of each step. The motivation behind filtering is to remove enstrophy accumulating at high wavenumbers and creating noise at grid-scale level.

The high-wavenumber filter used in the filtered timesteppers is:

\[\mathrm{filter}(\boldsymbol{k}) = \begin{cases} 1 & \quad |\boldsymbol{k}| ≤ k_{\textrm{cutoff}} \, ,\\ \exp{ \left [- \alpha (|\boldsymbol{k}| - k_{\textrm{cutoff}})^p \right]} & \quad |\boldsymbol{k}| > k_{\textrm{cutoff}} \, . \end{cases}\]

For fluid equations with quadratic non-linearities it makes sense to choose a cutoff wavenumber at 2/3 of the highest wavenumber resolved in our domain, $k_{\textrm{cutoff}} = \tfrac{2}{3} k_{\textrm{max}}$ (see discussion in Aliasing section).

Given the order $p$, we calculate the coefficient $\alpha$ so that the the filter value that corresponds to the highest allowed wavenumber in our domain is a small value, $\delta$, taken to be close to machine precision.

That is:

\[\alpha = \frac{- \log\delta}{(k_{\textrm{max}} - k_{\textrm{cutoff}})^p} \ .\]

The above filter originates from Canuto et al. (1988). In geophysical turbulence applications it was used by LaCasce (1996) and later by Arbic & Flierl (2003).

Using the default parameters provided by the filtered time steppers (see FourierFlows.makefilter), the filter has the following form:

using FourierFlows, CairoMakie

K = 0:0.01:1 # non-dimensional wavenumber k * dx / π

filter = FourierFlows.makefilter(K)

fig = Figure()
ax = Axis(fig[1, 1], xlabel = "k dx / π", ylabel = "filter", aspect=2.5, xticks=0:0.2:1)

lines!(ax, K, filter)