Skip to content

Commit

Permalink
Merge pull request #356 from FourierFlows/ncc/docs-visit
Browse files Browse the repository at this point in the history
Adds comparison of filter for different `order` parameter
  • Loading branch information
navidcy authored Sep 7, 2023
2 parents 270c41d + 2bb7ffd commit a26ed3a
Show file tree
Hide file tree
Showing 3 changed files with 43 additions and 28 deletions.
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

0 comments on commit a26ed3a

Please sign in to comment.