Merge pull request #323 from FourierFlows/ncc/add-palinstrophy

Add palinstrophy diagnostic in `TwoDNavierStokes`

Project.toml

@@ -2,7 +2,7 @@ name = "GeophysicalFlows"
 uuid = "44ee3b1c-bc02-53fa-8355-8e347616e15e"
 license = "MIT"
 authors = ["Navid C. Constantinou <[email protected]>", "Gregory L. Wagner <[email protected]>", "and co-contributors"]
-version = "0.15.1"
+version = "0.15.2"
 
 [deps]
 CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba"
@@ -18,15 +18,15 @@ StaticArrays = "90137ffa-7385-5640-81b9-e52037218182"
 Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
 
 [compat]
-CUDA = "^1, ^2.4.2, 3.0.0 - 3.6.4, ^3.7.1, 4"
-DocStringExtensions = "^0.8, 0.9"
-FFTW = "^1"
-FourierFlows = "^0.10.0"
-JLD2 = "^0.1, ^0.2, ^0.3, ^0.4"
-Reexport = "^0.2, ^1"
-SpecialFunctions = "^0.10, ^1, 2"
-StaticArrays = "^0.12, ^1"
-julia = "^1.6"
+CUDA = "1, 2.4.2, 3.0.0 - 3.6.4, 3.7.1, 4"
+DocStringExtensions = "0.8, 0.9"
+FFTW = "1"
+FourierFlows = "0.10.3"
+JLD2 = "0.1, 0.2, 0.3, 0.4"
+Reexport = "0.2, 1"
+SpecialFunctions = "0.10, 1, 2"
+StaticArrays = "0.12, 1"
+julia = "1.6"
 
 [extras]
 BenchmarkTools = "6e4b80f9-dd63-53aa-95a3-0cdb28fa8baf"

src/twodnavierstokes.jl

@@ -382,7 +382,7 @@ end
 """
     enstrophy(prob)
 
-Returns the problem's (`prob`) domain-averaged enstrophy,
+Return the problem's (`prob`) domain-averaged enstrophy,
 
 ```math
 \\int \\frac1{2} ζ² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |ζ̂|² ,
@@ -392,6 +392,25 @@ where ``ζ`` is the relative vorticity.
 """
 @inline enstrophy(prob) = 1 / (2 * prob.grid.Lx * prob.grid.Ly) * parsevalsum(abs2.(prob.sol), prob.grid)
 
+"""
+    palinstrophy(prob)
+
+Return the problem's (`prob`) domain-averaged palinstrophy,
+
+```math
+\\int \\frac1{2} |{\\bf ∇} ζ|² \\frac{𝖽x 𝖽y}{L_x L_y} = \\sum_{𝐤} \\frac1{2} |𝐤|² |ζ̂|² ,
+```
+
+where ``ζ`` is the relative vorticity.
+"""
+@inline function palinstrophy(prob)
+  sol, vars, grid = prob.sol, prob.vars, prob.grid
+  palinstrophyh = vars.uh # use vars.uh as scratch variable
+
+  @. palinstrophyh = 1 / 2 * grid.Krsq * abs2(sol)
+  return 1 / (grid.Lx * grid.Ly) * parsevalsum(palinstrophyh, grid)
+end
+
 """
     energy_dissipation(prob, ξ, νξ)
 

test/runtests.jl

@@ -46,7 +46,7 @@ for dev in devices
     @test test_twodnavierstokes_stochasticforcing_energybudget(dev)
     @test test_twodnavierstokes_deterministicforcing_enstrophybudget(dev)
     @test test_twodnavierstokes_stochasticforcing_enstrophybudget(dev)
-    @test test_twodnavierstokes_energyenstrophy(dev)
+    @test test_twodnavierstokes_energyenstrophypalinstrophy(dev)
     @test test_twodnavierstokes_problemtype(dev, Float32)
     @test TwoDNavierStokes.nothingfunction() == nothing
   end

test/test_twodnavierstokes.jl

@@ -257,19 +257,46 @@ function test_twodnavierstokes_advection(dt, stepper, dev::Device=CPU(); n=128,
   isapprox(prob.vars.ζ, ζf, rtol=rtol_twodnavierstokes)
 end
 
-function test_twodnavierstokes_energyenstrophy(dev::Device=CPU())
+function test_twodnavierstokes_energyenstrophypalinstrophy(dev::Device=CPU())
   nx, Lx  = 128, 2π
   ny, Ly  = 126, 3π
 
   grid = TwoDGrid(dev; nx, Lx, ny, Ly)
   x, y = gridpoints(grid)
 
   k₀, l₀ = 2π/grid.Lx, 2π/grid.Ly # fundamental wavenumbers
-  ψ₀ = @. sin(2k₀*x)*cos(2l₀*y) + 2sin(k₀*x)*cos(3l₀*y)
-  ζ₀ = @. -((2k₀)^2+(2l₀)^2)*sin(2k₀*x)*cos(2l₀*y) - (k₀^2+(3l₀)^2)*2sin(k₀*x)*cos(3l₀*y)
+  ψ₀ = @.                        sin(2k₀ * x) * cos(2l₀ * y) +  2sin(k₀ * x) * cos(3l₀ * y)
+  ζ₀ = @. -((2k₀)^2 + (2l₀)^2) * sin(2k₀ * x) * cos(2l₀ * y) - (k₀^2 + (3l₀)^2) * 2sin(k₀ * x) * cos(3l₀ * y)
 
-  energy_calc = 29/9
-  enstrophy_calc = 2701/162
+  #=
+
+  We can analytically compute the energy, enstrophy, and palinstrophy
+  that corresponds to the above initial condition using Mathematica
+
+    k0 = 2 Pi/Lx; l0 = 2 Pi/Ly;
+
+    psi[x_, y_] = Sin[2 k0 x] Cos[2 l0 y] + 2 Sin[k0 x] Cos[3 l0 y];
+    u[x, y]     = -D[psi[x, y], y];
+    v[x, y]     =  D[psi[x, y], x];
+    zeta[x, y]  =  D[v[x, y], x] - D[u[x, y], y];
+
+    E0 = Integrate[
+        1/2 (u[x, y]^2 + v[x, y]^2), {x, 0, Lx}, {y, 0, Ly}] / (Lx Ly);
+    E0 /. {Lx -> 2 Pi, Ly -> 3 Pi}
+
+    Z0 = Integrate[1/2 zeta[x, y]^2, {x, 0, Lx}, {y, 0, Ly}] / (Lx Ly);
+    Z0 /. {Lx -> 2 Pi, Ly -> 3 Pi}
+
+    P0 = Integrate[
+        1/2 (D[zeta[x, y], x]^2 + D[zeta[x, y], y]^2), {x, 0, Lx}, {y, 0, 
+        Ly}] / (Lx Ly);
+    P0 /. {Lx -> 2 Pi, Ly -> 3 Pi}
+
+  =#
+
+  energy_analytic       = 29/9
+  enstrophy_analytic    = 2701/162
+  palinstrophy_analytic = 126277/1458
 
   prob = TwoDNavierStokes.Problem(dev; nx, Lx, ny, Ly, stepper="ForwardEuler")
 
@@ -278,14 +305,16 @@ function test_twodnavierstokes_energyenstrophy(dev::Device=CPU())
   TwoDNavierStokes.set_ζ!(prob, ζ₀)
   TwoDNavierStokes.updatevars!(prob)
 
-  energyζ₀ = TwoDNavierStokes.energy(prob)
-  enstrophyζ₀ = TwoDNavierStokes.enstrophy(prob)
+  energy₀       = TwoDNavierStokes.energy(prob)
+  enstrophy₀    = TwoDNavierStokes.enstrophy(prob)
+  palinstrophy₀ = TwoDNavierStokes.palinstrophy(prob)
 
   params = TwoDNavierStokes.Params(p.ν, p.nν)
 
-  (isapprox(energyζ₀, energy_calc, rtol=rtol_twodnavierstokes) &&
-   isapprox(enstrophyζ₀, enstrophy_calc, rtol=rtol_twodnavierstokes) &&
-   TwoDNavierStokes.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing && p == params)
+  return (isapprox(energy₀, energy_analytic, rtol=rtol_twodnavierstokes) &&
+          isapprox(enstrophy₀, enstrophy_analytic, rtol=rtol_twodnavierstokes) &&
+          isapprox(palinstrophy₀, palinstrophy_analytic, rtol=rtol_twodnavierstokes) &&
+          TwoDNavierStokes.addforcing!(prob.timestepper.N, sol, cl.t, cl, v, p, g)==nothing && p == params)
 end
 
 function test_twodnavierstokes_problemtype(dev, T)