Add a filter and upsampler to gk twist-shift bc

[Antoinehoff]

This PR is related to DR #1083 and gkylcas PR gkeyllorg/gkylcas#106

Modal filter application

For simplicity we apply the filter on each DG modal coefficient separately. This holds as long as the filter width is a multiple of the cell width. Since the basis elements follow the property

$$\psi_i^{(p)}(x)=\psi^{(p)}(x-x_j),$$

i.e. they are all the same except for a shift that center them in the middle of their respective cell $i$, the introduction of a shift[+] $k\Delta x$ yields[/+][-] which consist in a multiple of the cell width $k\Delta x$ yields[/-]

$$\psi_i^{(p)}(x - k \Delta x) = \psi^{(p)}(x-x_i - k\Delta x) = \psi^{(p)}(x-x_{k+i}) = \psi_{i+k}^{(p)}(x).$$

Consequently we can write

$$\phi(x - k\Delta x) = \sum_{i=1}^N \sum_{p=0}^{P} \phi_i^{(p)} \psi_{i+k}^{(p)}(x),$$

where we do not care about the $i+k>N$ limit for now (this can be resolved by a reflection at the boundary which preserves M0).
Now we see that the coefficient of order $q$ of the filtered DG rerpresentation in a cell $j$ is

$$\hat{\phi}_{j}^{(q)} = \int dx \hat{\phi}(x) \psi_j^{(q)}(x) = \sum_{k=-M}^{M} \hat{h}[k] \sum_{i=1}^N \sum_{p=0}^{P} \phi_i^{(p)} \int dx \psi_{i+k}^{(p)}(x)\psi_j^{(q)} = \sum_{k=-M}^{M} \hat{h}[k] \phi_{j-k}^{(q)}.$$

Caveat

The filter is done on the DG modal coefficients directly and cannot prevent oscillations smaller than the cell width. This is not a problem for our main application, twist-shift anti-aliasing, because we filter a field that is previously upscaled to a resolution expected to be large enough so that the twist-shift does not introduce grid scale mode.

Single mode example

We show here the example of a single two period $k_y$ mode as a "donor", the data that need to be twisted. Below, you see from left to right: i) the donor on the original resolution, ii) the shifted version of the donor on a 4x resolution (the truth), iii) the shifted version using the original Gkeyll TSBC, iv) the new Gkeyll TSBC that includes supersampling and filter.

In the original TSBC we see a reversed kx>0 wave appearing in the 2.4<x<3.0 region. This is the aliasing we want to get rid of. A finite perturbation remains at the border of the domain. This is part of the limit of a discrete filter capability. At the boundary, the filter cannot be applied entirely as it is applied to neighbors cells. Here we reflect the data to apply the filter but this causes a reduction of the damping capability.

The Fourier spectral analysis shows that the new TSBC is indeed removing that branch and mimics better the spectrum of the truth. The higher modes, rising from the shift, that cannot be resolved are simply dissipated similarly to some local GK codes (GS2, GENE, stella, Gyacomo).

Two mode example

We consider here a higher resolution and a two-mode input. The first mode is sheared outside of the resolution that can be resolved by the original grid. We see that the TSBC without anti-aliasing filter is presenting a checkerboard pattern while the anti-aliased version presents a much smoother wave, though boundary effects remain.