FourierFlows / FourierFlows.jl

Tools for building fast, hackable, pseudospectral partial differential equation solvers on periodic domains
https://bit.ly/FourierFlows
MIT License
207 stars 29 forks source link

Fixes bug in `parsevalsum`/`parsevalsum2` #350

Closed navidcy closed 1 year ago

navidcy commented 1 year ago

This PR fixes a bug in parsevalsum/parsevalsum2. When rfft transforms are provided, then the mode that corresponds to k = nx/2 should only be summed once (and not twice) since its conjugate is not part of rfft output.

The previous tests did not catch this because the functions we were giving as input for the parsevalsum tests were smooth and thus had no Fourier amplitude at k = nx/2.

glwagner commented 1 year ago

Maybe another reference?

navidcy commented 1 year ago

What do you mean "another reference"?

navidcy commented 1 year ago
julia> using FourierFlows
[ Info: FourierFlows will use 6 threads

julia> grid = OneDGrid(nx=6, Lx=2π)
OneDimensionalGrid
  ├─────────── Device: CPU
  ├──────── FloatType: Float64
  ├────────── size Lx: 6.283185307179586
  ├──── resolution nx: 6
  ├── grid spacing dx: 1.0471975511965976
  ├─────────── domain: x ∈ [-3.141592653589793, 2.094395102393195]
  └─ aliased fraction: 0.3333333333333333

julia> grid.kr
4-element Vector{Float64}:
 0.0
 1.0
 2.0
 3.0

julia> grid.k
6-element Vector{Float64}:
  0.0
  1.0
  2.0
 -3.0
 -2.0
 -1.0

The mode k = -nx/2 = -3 (in the example above) only appears once in the full wavenumber grid grid.k. So when we only have the grid.kr wavenumber, we should not sum over kr = nx/2 twice.

navidcy commented 1 year ago

I discovered the bug when I was testing out something and was doing

fh = randn(Complex{eltype(grid)}, (grid.nkr, grid.nl))
f = irfft(deepcopy(f), grid.nx)

and then tried to confirm that parsevalsum2(fh, grid) was the same as sum(f.^2) * grid.dx * grid.dy...