Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Adds comparison of filter for different order parameter #356

Merged
merged 8 commits into from
Sep 7, 2023
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion docs/make.jl
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ using
##### Generate examples
#####

@show bib_filepath = joinpath(@__DIR__, "src/references.bib")
bib_filepath = joinpath(@__DIR__, "src/references.bib")
bib = CitationBibliography(bib_filepath)

const EXAMPLES_DIR = joinpath(@__DIR__, "..", "examples")
Expand Down
62 changes: 38 additions & 24 deletions docs/src/timestepping.md
Original file line number Diff line number Diff line change
Expand Up @@ -16,37 +16,46 @@ For example, `"RK4"` for the Runge-Kutta 4th-order time stepper.
Most of the time steppers also come with their `Filtered` equivalents: [`FilteredForwardEulerTimeStepper`](@ref), [`FilteredAB3TimeStepper`](@ref), [`FilteredRK4TimeStepper`](@ref), [`FilteredLSRK54TimeStepper`](@ref), and [`FilteredETDRK4TimeStepper`](@ref).

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 motivation behind filtering is to preclude enstrophy from accumulating at high wavenumbers and
creating noise at grid-scale level.

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

```math
\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}
\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](@ref aliasing)).
For fluid equations with quadratic non-linearities it makes sense to choose a cut-off wavenumber
at 2/3 of the highest wavenumber that is resolved in our domain,
``k_{\textrm{cutoff}} = \tfrac{2}{3} k_{\textrm{max}}`` (see discussion in [Aliasing section](@ref aliasing)).

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:
Given the order ``p``, we choose coefficient ``\alpha`` so that the filter value that corresponds
to the highest allowed wavenumber in our domain is a small number ``\delta``, usually taken to be
close to machine precision. That is:

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

The above filter originates from the book by [Canuto-etal-1987](@cite). In geophysical turbulence
applications it was used by [LaCasce-1996](@cite) and later by [Arbic-Flierl-2004](@cite).
The above-mentioned filter form originates from the book by [Canuto-etal-1987](@cite).
In geophysical turbulence applications it was used by [LaCasce-1996](@cite) and later
by [Arbic-Flierl-2004](@cite).

!!! warning "Not too steep, not too shallow"
Care should be taken if one decides to fiddle with the filter parameters. Changing
the order ``p`` affects how steeply the filter falls off. Lower order values imply
that the filter might fall off too quickly and may lead to Gibbs artifacts; higher
order value implies that the filter might fall off too slow and won't suffice to
remove enstrophy accumulating at the grid scale.

Using the default parameters provided by the filtered time steppers (see
[`FourierFlows.makefilter`](@ref)), the filter has the following form:
The filter using the default parameters provided by the filtered time steppers (see
[`FourierFlows.makefilter`](@ref)) is depicted below. The same plot also compares how
the filter changes when we vary the order parameter ``p`` while keeping everything
else the same.

```@setup 1
using CairoMakie
Expand All @@ -55,16 +64,21 @@ set_theme!(Theme(linewidth = 3, fontsize = 20))
```

```@example 1
using FourierFlows, CairoMakie

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

filter = FourierFlows.makefilter(K)
K = 0:0.001:1 # non-dimensional wavenumber k * dx / π

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

vlines!(ax, 2/3, color = (:gray, 0.4), linewidth = 6, label = "cutoff wavenumber |𝐤| dx / π = 2/3 (default)")

lines!(ax, K, makefilter(K), linewidth = 4, label = "order 4 (default)")
lines!(ax, K, makefilter(K, order = 1), linestyle = :dash, label = "order 1")
lines!(ax, K, makefilter(K, order = 10), linestyle = :dot, label = "order 10")

lines!(ax, K, filter)
axislegend(position = :lb)

current_figure() # hide
```
7 changes: 4 additions & 3 deletions src/domains.jl
Original file line number Diff line number Diff line change
Expand Up @@ -485,7 +485,7 @@ one-dimensional grid, the non-dimensional wavenumber `K` is
K = k * dx / π
```

and thus it takes values `K` ``[-1, 1]``.
and thus take values in ``[-1, 1]``.

For `K ≤ innerK` the filter is inactive, i.e., equal to 1. For `K > innerK`,
the filter decays exponentially to remove high-wavenumber content from
Expand All @@ -495,8 +495,9 @@ the solution, i.e.,
filter(K) = exp(- decay * (K - innerK)^order)
```

For a given `order` and , the `decay` rate is determined so that the filter value at the
outer wavenumber `outerK` is `tol`, where `tol` is a small number, close to machine precision.
For a given `order`, the `decay` rate is determined so that the filter value at the
outer wavenumber `outerK` is `tol`, where `tol` is a small number, close to machine
precision.

```julia
decay = - log(tol) / (outerK - innerK)^order
Expand Down