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)