We propose to implement a low pass filter updater to gkeyll.
Motivation
We identified aliasing rising from the twist-and-shift boundary condition and we want to be able to filter the waves that cannot be resolved on the grid. Note that this aliasing appears in the x-direction, which is not periodic.
Method
Since we cannot expect periodicity in the direction, we propose here to implement a discrete sinc filter combined with a Blackman window. The discrete impulse response for a sample index $k \in [-M, M]$ (where $M$ is the half-width of the filter) is:
$$h_{\text{sinc}}[k] = \text{sinc}(2 f_c k) = \frac{\sin(2 \pi f_c k)}{2 \pi f_c k}$$
To prevent spectral leakage (the Gibbs phenomenon) caused by truncating the sinc function, we apply a Blackman window function $w[k]$ over the same domain $k \in [-M, M]$,
$$w[k] = 0.42 + 0.5 \cos\left(\frac{\pi k}{M}\right) + 0.08 \cos\left(\frac{2 \pi k}{M}\right),$$
where $0.42$, $0.5$, and $0.08$ are standard coefficients.
The practical filter coefficients are obtained by multiplying the ideal sinc response by the window function element-wise:
$$h[k] = h_{\text{sinc}}[k] \cdot w[k]$$
To guarantee that the filtering process preserves the overall mean amplitude of the wave, we normalize the discrete weights so they sum exactly to $1$,
$$\hat{h}[k] = \frac{h[k]}{\sum_{j=-M}^{M} h[j]}.$$
The filtered field writes
$$\hat{\phi}(x) = \sum_{k=-M}^{M} \hat{h}[k] \phi(x - k\Delta x),$$
Implementation
We plan to create a new updater in the same fashion as dg_interpolate in core/zero. This will make it accessible to the entire code base since it is pretty generic. The updater will have the usual _new, _advance, and _release methods. Here is a sketch of it's interface:
/**
* @param dir Direction along which to filter.
* @param half_width Stencil half-width M in cells (stencil spans 2M+1 cells).
* @param cutoff_wavelength Cutoff wavelength (physical units) of the filter.
* @param basis DG basis of the filtered field.
* @param grid Grid the filtered field is defined on.
* @param range Range to filter in (a sub-range of the range the fields are
* defined on). Donor cells are taken from this range only: the
* stencil is truncated (and renormalized) at its boundaries
* along dir.
* @param use_gpu bool to determine if on GPU.
* @return New filter updater.
*/
struct gkyl_dg_lowpass_filter*
gkyl_dg_lowpass_filter_new(int dir, int half_width, double cutoff_wavelength,
const struct gkyl_basis *basis, const struct gkyl_rect_grid *grid,
const struct gkyl_range *range, bool use_gpu);
Possible caveat
The behavior of the filter is not consistent at the boundary, and this is a tradeoff we need to accept. It will mainly not filter well at M/2 number of cell close to the filtering range boundaries. We will see what consequences this have in some examples. We can hope that the boundary will be dominated by boundary conditions like Dirichlet but nothing is certain.
In the case of clopen simulations, one could think of making the filter going inside the closed field line region so that we do not have boundary effects at the separatrix at least.
We propose to implement a low pass filter updater to gkeyll.
Motivation
We identified aliasing rising from the twist-and-shift boundary condition and we want to be able to filter the waves that cannot be resolved on the grid. Note that this aliasing appears in the x-direction, which is not periodic.
Method
Since we cannot expect periodicity in the direction, we propose here to implement a discrete sinc filter combined with a Blackman window. The discrete impulse response for a sample index$k \in [-M, M]$ (where $M$ is the half-width of the filter) is:
To prevent spectral leakage (the Gibbs phenomenon) caused by truncating the sinc function, we apply a Blackman window function$w[k]$ over the same domain $k \in [-M, M]$ ,
where$0.42$ , $0.5$ , and $0.08$ are standard coefficients.
The practical filter coefficients are obtained by multiplying the ideal sinc response by the window function element-wise:
To guarantee that the filtering process preserves the overall mean amplitude of the wave, we normalize the discrete weights so they sum exactly to$1$ ,
The filtered field writes
Implementation
We plan to create a new updater in the same fashion as
dg_interpolateincore/zero. This will make it accessible to the entire code base since it is pretty generic. The updater will have the usual_new,_advance, and_releasemethods. Here is a sketch of it's interface:Possible caveat
The behavior of the filter is not consistent at the boundary, and this is a tradeoff we need to accept. It will mainly not filter well at M/2 number of cell close to the filtering range boundaries. We will see what consequences this have in some examples. We can hope that the boundary will be dominated by boundary conditions like Dirichlet but nothing is certain.
In the case of clopen simulations, one could think of making the filter going inside the closed field line region so that we do not have boundary effects at the separatrix at least.