Merge pull request #350 from FourierFlows/ncc/clarify

Fixes bug in `parsevalsum`/`parsevalsum2`

Author: navidcy

Project.toml

@@ -6,7 +6,7 @@ authors = ["Gregory L. Wagner <wagner.greg@gmail.com>", "Navid C. Constantinou <
 description = "Tools for building fast, hackable, pseudospectral partial differential equation solvers on periodic domains."
 documentation = "https://fourierflows.github.io/FourierFlowsDocumentation/stable/"
 repository = "https://github.com/FourierFlows/FourierFlows.jl"
-version = "0.10.2"
+version = "0.10.3"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -24,7 +24,7 @@ CUDA = "1, 2.4.2, 3.0.0 - 3.6.4, 3.7.1, 4"
 DocStringExtensions = "0.8, 0.9"
 FFTW = "1"
 Interpolations = "0.12, 0.13, 0.14"
-JLD2 = "^0.1, 0.2, 0.3, 0.4"
+JLD2 = "0.1, 0.2, 0.3, 0.4"
 Reexport = "0.2, 1"
 julia = "1.6"
 

src/utils.jl

@@ -97,18 +97,21 @@ end
 """
     parsevalsum2(uh, grid)
 
-Return `Σ |uh|²` on the `grid`, which is equal to the domain integral of `u`. More specifically, 
-it returns
+Return `Σ |uh|²` on the `grid`, which is equal to the domain integral `u` squared.
+For example on a 2D grid, `parsevalsum2` returns
+
 ```math
-\\sum_{𝐤} |û_{𝐤}|² L_x L_y = \\int u(𝐱)² \\, 𝖽x 𝖽y \\,,
+\\sum_{𝐤} |û_{𝐤}|² L_x L_y = \\iint u² \\, 𝖽x 𝖽y \\,,
 ```
+
 where ``û_{𝐤} =`` `uh` `` / (n_x e^{i 𝐤 ⋅ 𝐱₀})``. The elements of the vector ``𝐱₀`` are the
 left-most position in each direction, e.g., for a 2D grid `(grid.x[1], grid.y[1])`.
 """
 function parsevalsum2(uh, grid::TwoDGrid)
-  if size(uh, 1) == grid.nkr # uh is in conjugate symmetric form
-    U = sum(abs2, uh[1, :])           # k=0 modes
-    U += 2*sum(abs2, uh[2:end, :])    # sum k>0 modes twice
+  if size(uh, 1) == grid.nkr  # uh is in conjugate symmetric form
+    U = sum(abs2, uh[1, :])                  # k=0 modes
+    U += sum(abs2, uh[grid.nkr, :])          # k=nx/2 modes
+    U += 2 * sum(abs2, uh[2:grid.nkr-1, :])  # sum twice for 0 < k < nx/2 modes
   else # count every mode once
     U = sum(abs2, uh)
   end
@@ -119,9 +122,10 @@ function parsevalsum2(uh, grid::TwoDGrid)
 end
 
 function parsevalsum2(uh, grid::OneDGrid)
-  if size(uh, 1) == grid.nkr                 # uh is conjugate symmetric
-    U = sum(abs2, CUDA.@allowscalar uh[1])   # k=0 modes
-    U += @views 2 * sum(abs2, uh[2:end])     # sum k>0 modes twice
+  if size(uh, 1) == grid.nkr  # uh is conjugate symmetric
+    U = sum(abs2, CUDA.@allowscalar uh[1])          # k=0 mode
+    U += sum(abs2, CUDA.@allowscalar uh[grid.nkr])  # k=nx/2 mode
+    U += @views 2 * sum(abs2, uh[2:grid.nkr-1])     # sum twice for 0 < k < nx/2 modes
   else # count every mode once
     U = sum(abs2, uh)
   end
@@ -134,17 +138,20 @@ end
 """
     parsevalsum(uh, grid)
 
-Return `real(Σ uh)` on the `grid`, i.e.
+Return `real(Σ uh)` on the `grid`, i.e.,
+
 ```math
 ℜ [ \\sum_{𝐤} û_{𝐤} L_x L_y ] \\,,
+
 ```
 where ``û_{𝐤} =`` `uh` `` / (n_x e^{i 𝐤 ⋅ 𝐱₀})``. The elements of the vector ``𝐱₀`` are the
 left-most position in each direction, e.g., for a 2D grid `(grid.x[1], grid.y[1])`.
 """
 function parsevalsum(uh, grid::TwoDGrid)
-  if size(uh, 1) == grid.nkr    # uh is conjugate symmetric
-    U = sum(uh[1, :])           # k=0 modes
-    U += 2*sum(uh[2:end, :])    # sum k>0 modes twice
+  if size(uh, 1) == grid.nkr  # uh is conjugate symmetric
+    U = sum(uh[1, :])                  # k = 0 modes
+    U += sum(uh[grid.nkr, :])          # k = nx/2 modes
+    U += 2 * sum(uh[2:grid.nkr-1, :])  # sum twice for 0 < k < nx/2 modes
   else # count every mode once
     U = sum(uh)
   end
@@ -155,9 +162,10 @@ function parsevalsum(uh, grid::TwoDGrid)
 end
 
 function parsevalsum(uh, grid::OneDGrid)
-  if size(uh, 1) == grid.nkr    # uh is conjugate symmetric
-    U = uh[1]                   # k=0 mode
-    U += 2*sum(uh[2:end])       # sum k>0 modes twice
+  if size(uh, 1) == grid.nkr        # uh is conjugate symmetric
+    U = uh[1]                       # k=0 mode
+    U += uh[grid.nkr]               # k=nx/2 mode
+    U += 2 * sum(uh[2:grid.nkr-1])  # sum twice for 0 < k < nx/2 modes
   else # count every mode once
     U = sum(uh)
   end

test/runtests.jl

@@ -267,6 +267,13 @@ for dev in devices
     @test test_parsevalsums(f1, g; realvalued=false)      # Real valued f with fft
     @test test_parsevalsums(f2, g; realvalued=false)      # Complex valued f with fft
 
+    f3 = randn(eltype(g), (g.nx, g.ny))
+    @test test_parsevalsums(f3, g; realvalued=true)        # Real valued f with rfft
+    @test test_parsevalsums(f3, g; realvalued=false)       # Real valued f with fft
+
+    f4 = randn(Complex{eltype(g)}, (g.nx, g.ny))
+    @test test_parsevalsums(f4, g; realvalued=false)       # Complex valued f with fft
+
     @test test_jacobian(sinkl1, sinkl1, 0*sinkl1, g)      # Test J(a, a) = 0
     @test test_jacobian(sinkl1, sinkl2, Jsinkl1sinkl2, g) # Test J(sin1, sin2) = Jsin1sin2
     @test test_jacobian(expkl1, expkl2, Jexpkl1expkl2, g) # Test J(exp1, exp2) = Jexp1exps2
@@ -279,7 +286,13 @@ for dev in devices
     @test test_parsevalsums(f1, g; realvalued=true)        # Real valued f with rfft
     @test test_parsevalsums(f1, g; realvalued=false)       # Real valued f with fft
 
+    f2 = randn(eltype(g), (g.nx,))
+    @test test_parsevalsums(f2, g; realvalued=true)        # Real valued f with rfft
+    @test test_parsevalsums(f2, g; realvalued=false)       # Real valued f with fft
 
+    f3 = randn(Complex{eltype(g)}, (g.nx,))
+    @test test_parsevalsums(f3, g; realvalued=false)       # Complex valued f with fft
+
     # Radial spectrum tests. Note that ahρ = ∫ ah ρ dθ.
     n = 128; δ = n/10                                       # Parameters
     ahkl_isotropic(k, l) = exp(-(k^2 + l^2) / 2δ^2)                     # a  = exp(-ρ²/2δ²)