Changes default AMD poincare coefficient (#4494)

* Changes default AMD poincare coefficient

* Update test_nonhydrostatic_regression.jl

* Update anisotropic_minimum_dissipation.jl

* update amd dockets

* updated docstrings + added validation case

* Update src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl

Co-authored-by: Gregory L. Wagner <[email protected]>

* Update anisotropic_minimum_dissipation.jl

---------

Co-authored-by: Gregory L. Wagner <[email protected]>
Co-authored-by: Navid C. Constantinou <[email protected]>

Author: 3 people

docs/oceananigans.bib

@@ -990,4 +990,13 @@ @article{Lan2022
 	pages = {111050},
 	year = {2022},
 	doi = {https://doi.org/10.1016/j.jcp.2022.111050},
+}
+
+@article{Verstappen14,
+    author = {Verstappen, Roel and Rozema, Wybe and Bae, H. Jane},
+	year = {2014},
+	month = {12},
+	pages = {417-426},
+	title = {Numerical scale separation in large-eddy simulation},
+	url = {https://www.researchgate.net/publication/311885631}
 }

docs/src/model_setup/turbulent_diffusivity_closures_and_les_models.md

@@ -82,8 +82,8 @@ julia> using Oceananigans.TurbulenceClosures
 
 julia> closure = AnisotropicMinimumDissipation()
 AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
-           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
-    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.08333333333333333
+           Poincaré constant for momentum eddy viscosity Cν: 0.3333333333333333
+    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.3333333333333333
                         Buoyancy modification multiplier Cb: nothing
 ```
 

src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl

@@ -24,13 +24,13 @@ const AMD = AnisotropicMinimumDissipation
 
 Base.show(io::IO, closure::AMD{TD}) where TD =
     print(io, "AnisotropicMinimumDissipation{$TD} turbulence closure with:\n",
-              "           Poincaré constant for momentum eddy viscosity Cν: ", closure.Cν, "\n",
-              "    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: ", closure.Cκ, "\n",
+              "           Poincaré constant for momentum eddy viscosity Cν: ", closure.Cν, "\n",
+              "    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: ", closure.Cκ, "\n",
               "                        Buoyancy modification multiplier Cb: ", closure.Cb)
 
 """
     AnisotropicMinimumDissipation([time_discretization = ExplicitTimeDiscretization, FT = Float64;]
-                                  C = 1/12, Cν = nothing, Cκ = nothing, Cb = nothing)
+                                  C = 1/3, Cν = nothing, Cκ = nothing, Cb = nothing)
 
 
 Return parameters of type `FT` for the `AnisotropicMinimumDissipation`
@@ -49,24 +49,29 @@ Arguments
 
 Keyword arguments
 =================
-* `C`: Poincaré constant for both eddy viscosity and eddy diffusivities. `C` is overridden
-       for eddy viscosity or eddy diffusivity if `Cν` or `Cκ` are set, respecitvely.
+* `C`: Poincaré constant for both eddy viscosity and eddy diffusivities. `C` is overridden
+       for eddy viscosity or eddy diffusivity if `Cν` or `Cκ` are set, respectively.
 
-* `Cν`: Poincaré constant for momentum eddy viscosity.
+* `Cν`: Poincaré constant for momentum eddy viscosity.
 
-* `Cκ`: Poincaré constant for tracer eddy diffusivities. If one number or function, the same
+* `Cκ`: Poincaré constant for tracer eddy diffusivities. If one number or function, the same
         number or function is applied to all tracers. If a `NamedTuple`, it must possess
-        a field specifying the Poncaré constant for every tracer.
+        a field specifying the Poncaré constant for every tracer.
 
 * `Cb`: Buoyancy modification multiplier (`Cb = nothing` turns it off, `Cb = 1` was used by [Abkar16](@citet)).
         *Note*: that we _do not_ subtract the horizontally-average component before computing this
         buoyancy modification term. This implementation differs from [Abkar16](@citet)'s proposal
         and the impact of this approximation has not been tested or validated.
 
-By default: `C = Cν = Cκ = 1/12`, which is appropriate for a finite-volume method employing a
-second-order advection scheme, and `Cb = nothing`, which turns off the buoyancy modification term.
+By default: `C = Cν = Cκ = 1/3`, and `Cb = nothing`, which turns off the buoyancy modification term.
+The default Poincaré constant is found by discretizing subgrid scale energy production, assuming a 
+second-order advection scheme. [Verstappen14](@citeo) shows that the Poincaré constant 
+should be 4 times larger than for straightforward (spectral) discretisation, resulting in `C = 1/3` 
+in our formulation. They also empirically demonstrated that this coefficient produces the correct
+discrete production-dissipation balance. We further demonstrated this in 
+https://github.com/CliMA/Oceananigans.jl/issues/4367.
 
-`Cν` or `Cκ` may be numbers, or functions of `x, y, z`.
+`C`, `Cν` and `Cκ` may be numbers, or functions of `x, y, z`.
 
 Examples
 ========
@@ -76,8 +81,8 @@ julia> using Oceananigans
 
 julia> pretty_diffusive_closure = AnisotropicMinimumDissipation(C=1/2)
 AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
-           Poincaré constant for momentum eddy viscosity Cν: 0.5
-    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.5
+           Poincaré constant for momentum eddy viscosity Cν: 0.5
+    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: 0.5
                         Buoyancy modification multiplier Cb: nothing
 ```
 
@@ -90,8 +95,8 @@ julia> surface_enhanced_tracer_C(x, y, z) = 1/12 * (1 + exp((z + Δz/2) / 8Δz))
 
 julia> fancy_closure = AnisotropicMinimumDissipation(Cκ=surface_enhanced_tracer_C)
 AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
-           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
-    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: surface_enhanced_tracer_C
+           Poincaré constant for momentum eddy viscosity Cν: 0.3333333333333333
+    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: surface_enhanced_tracer_C
                         Buoyancy modification multiplier Cb: nothing
 ```
 
@@ -100,23 +105,27 @@ julia> using Oceananigans
 
 julia> tracer_specific_closure = AnisotropicMinimumDissipation(Cκ=(c₁=1/12, c₂=1/6))
 AnisotropicMinimumDissipation{ExplicitTimeDiscretization} turbulence closure with:
-           Poincaré constant for momentum eddy viscosity Cν: 0.08333333333333333
-    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: (c₁ = 0.08333333333333333, c₂ = 0.16666666666666666)
+           Poincaré constant for momentum eddy viscosity Cν: 0.3333333333333333
+    Poincaré constant for tracer(s) eddy diffusivit(ies) Cκ: (c₁ = 0.08333333333333333, c₂ = 0.16666666666666666)
                         Buoyancy modification multiplier Cb: nothing
 ```
 
 References
 ==========
 
+Verstappen, R. & Rozema, W. & Bae, J.H. (2014), "Numerical scale separation in large-eddy 
+    simulation", Center for Turbulence ResearchProceedings of the Summer Program 2014.
+
 Vreugdenhil C., and Taylor J. (2018), "Large-eddy simulations of stratified plane Couette
     flow using the anisotropic minimum-dissipation model", Physics of Fluids 30, 085104.
 
 Verstappen, R. (2018), "How much eddy dissipation is needed to counterbalance the nonlinear
     production of small, unresolved scales in a large-eddy simulation of turbulence?",
     Computers & Fluids 176, pp. 276-284.
+}
 """
 function AnisotropicMinimumDissipation(time_disc::TD = ExplicitTimeDiscretization(), FT = Oceananigans.defaults.FloatType;
-                                       C = FT(1/12), Cν = nothing, Cκ = nothing, Cb = nothing) where TD
+                                       C = FT(1/3), Cν = nothing, Cκ = nothing, Cb = nothing) where TD
 
     Cν = Cν === nothing ? C : Cν
     Cκ = Cκ === nothing ? C : Cκ

test/test_nonhydrostatic_regression.jl

@@ -64,7 +64,7 @@ include("regression_tests/ocean_large_eddy_simulation_regression_test.jl")
                     run_thermal_bubble_regression_test(arch, grid_type)
                 end
 
-                amd_closure = (AnisotropicMinimumDissipation(), ScalarDiffusivity(ν=1.05e-6, κ=1.46e-7))
+                amd_closure = (AnisotropicMinimumDissipation(C=1/12), ScalarDiffusivity(ν=1.05e-6, κ=1.46e-7))
                 smag_closure = (SmagorinskyLilly(C=0.23, Cb=1, Pr=1), ScalarDiffusivity(ν=1.05e-6, κ=1.46e-7))
 
                 for closure in (amd_closure, smag_closure)
@@ -77,4 +77,4 @@ include("regression_tests/ocean_large_eddy_simulation_regression_test.jl")
             end
         end
     end
-end
+end

validation/les/decaying_homogeneous_isotropic_turbulence.jl

@@ -0,0 +1,85 @@
+using Oceananigans, FFTW, StatsBase, CairoMakie
+
+using  Oceananigans.Models.HydrostaticFreeSurfaceModels: compute_w_from_continuity!
+
+arch = CPU()
+
+grid = RectilinearGrid(arch, size = (32, 32, 32), extent = (1, 1, 1))
+
+closure = AnisotropicMinimumDissipation(VerticallyImplicitTimeDiscretization(), C=1/3)
+
+model = NonhydrostaticModel(; grid, closure)
+
+u, v, w = model.velocities
+
+set!(u, (args...)->randn())
+set!(v, (args...)->randn())
+
+compute_w_from_continuity!((; u, v, w), arch, grid)
+
+function extract_energy(u, v, w, grid)
+    û = fftshift(fft(@at (Center, Center, Center) u))
+    v̂ = fftshift(fft(@at (Center, Center, Center) v))
+    ŵ = fftshift(fft(@at (Center, Center, Center) w))
+
+    Nx, Ny, Nz = u.grid.Nx, u.grid.Ny, u.grid.Nz
+
+    kx = fftshift(fftfreq(Nx, 1/u.grid.Δxᶜᵃᵃ))
+    ky = fftshift(fftfreq(Ny, 1/u.grid.Δyᵃᶜᵃ))
+    kz = fftshift(fftfreq(Nz, 1/u.grid.z.Δᵃᵃᶜ))
+
+    Nx, Ny, Nz = length(kx), length(ky), length(kz)
+
+    KX = reshape(kx, Nx, 1, 1)
+    KY = reshape(ky, 1, Ny, 1)
+    KZ = reshape(kz, 1, 1, Nz)
+
+    K = sqrt.(KX.^2 .+ KY.^2 .+ KZ.^2)
+
+    E = abs.((û.^2 .+ v̂.^2 .+ ŵ.^2)./2)
+
+    return E, K
+end
+
+fig = Figure()
+
+ax = Axis(fig[1, 1], xscale=log, yscale=log, xlabel = "Wavenumber (1/m)", ylabel = "Energy density (m³/s²)")
+
+k_bins = exp.(0:0.25:4)
+
+xlims!(ax, minimum(k_bins), maximum(k_bins))
+ylims!(ax, exp(-1), exp(25))
+
+E, K = extract_energy(u, v, w, grid)
+
+E_binned = fit(Histogram, [K...], weights([E...]), k_bins).weights
+
+lines!(ax, (k_bins[1:end-1] .+ k_bins[2:end])./2, E_binned, color = 0, colorrange = (0, 10), colormap = :oslo)
+
+simulation = Simulation(model, Δt = 0.5 * minimum_xspacing(grid) / 5, stop_time = 10)
+
+add_callback!(simulation, (sim)->(@info "$(prettytime(time(sim))) in $(prettytime(sim.run_wall_time))"), IterationInterval(100))
+
+function add_line!(sim)
+    E, K = extract_energy(u, v, w, grid)
+    E_binned = fit(Histogram, [K...], weights([E...]), k_bins).weights
+    lines!(ax, (k_bins[1:end-1] .+ k_bins[2:end])./2, E_binned, color = time(sim), colorrange = (0, 10), colormap = :oslo)
+end
+
+add_callback!(simulation, add_line!, TimeInterval(0.1))
+
+run!(simulation)
+
+lines!(ax, [exp(1), exp(3)], x->(x/exp(1))^(-5/3)*exp(15), color = :red, linestyle = :dash)
+
+text!(ax, exp(0.8), exp(13); text = L"$E(k)\sim k^{-5/3}$", color = :white) 
+
+Colorbar(fig[1, 2], colormap = :oslo, colorrange = (0, 10), label = "Time (s)")
+
+k_filt = 1/sqrt(closure.Cν * 3 / (1/(2*minimum_xspacing(grid))^2 + 1/(2*minimum_yspacing(grid))^2 + 1/(2*minimum_zspacing(grid))^2))
+
+lines!(ax, ones(2) .* k_filt, [exp(0), exp(10)], color = :red, linestyle = :dash)
+
+text!(ax, k_filt, exp(10); text = "1/δ")
+
+save("energy.png", fig)