Choice of constant for the Smagorinsky model

[tomchor]

I have a feeling that I've asked this before but I looked and couldn't find anything. So forgive me if this is double-posting.

Usually the Smagorisnky eddy viscosity is defined as:

and with C_s=0.165.

Looking at the code I can see that the formulation is a bit different:

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

Line 94 in 32c5c5a

@inline νₑ_deardorff(ς, C, Δᶠ, Σ²) = ς * (C*Δᶠ)^2 * sqrt(2Σ²)

Unless I'm missing something, this means that the coefficient C of the closure is different from the coefficient C_s as it is usually defined (and also as it is defined in the docs). I'm curious as to why this choice was made.

Also, assuming ς=1, the default value of C=0.23 doesn't seem to match the default value in the literature for C_s=0.165.

Am I missing something here? Could someone please clarify?

[glwagner]

I'm not sure what part of the formula you're pointing out is specifically different from the image you've pasted, but it does look like something is fishy. Here's what I found:

From Pressel et al. 2015 --- which we should cite in our code since I believe it was used as a reference ---

The formula in Oceananigans is

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

Line 94 in 32c5c5a

@inline νₑ_deardorff(ς, C, Δᶠ, Σ²) = ς * (C*Δᶠ)^2 * sqrt(2Σ²)

which is called here:

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

Line 102 in 32c5c5a

return νₑ_deardorff(ς, clo.C, Δᶠ, Σ²) + clo.ν

It looks like Δᶠ is the same that's used in Pressel et al. 2015:

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_utils.jl

Lines 23 to 24 in 9248437

	@inline geo_mean_Δᶠ(i, j, k, grid::AbstractGrid) =
	cbrt(Δxᶜᶜᵃ(i, j, k, grid) * Δyᶜᶜᵃ(i, j, k, grid) * Δzᵃᵃᶜ(i, j, k, grid))

and, apparently, Σ² = S_ij S_ij (if we believe we understand the notation here):

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

Line 97 in 32c5c5a

Σ² = ΣᵢⱼΣᵢⱼᶜᶜᶜ(i, j, k, grid, U.u, U.v, U.w)

This means that our formula is identical to Pressel et al (2015), but the constant is wrong: Pressel proposes Cs=0.17, and we have a default of Cs=0.23:

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

Lines 63 to 64 in 32c5c5a

	SmagorinskyLilly(FT=Float64; C=0.23, Cb=1.0, Pr=1.0, ν=0, κ=0,
	time_discretization::TD=ExplicitTimeDiscretization()) where TD =

Amusingly,

julia> 0.165 * sqrt(2)
0.23334523779156072

which is probably not a coincidence. However, I'm not sure how an extra factor of sqrt(2) snuck into our constant.

As for clarification, I don't have much to offer. Perhaps the constant was taken from some reference that used a different formulation than either us or Pressel et al. 2015. Nobody has submitted a validation test for this closure so I don't think we know how it performs.

As a historical note, the paper cited by both Pressel et al (2015) and us is Lilly (1962), which does indeed use the same formulation:

where

@tomchor I can't tell if the formula you've pasted is actually different from ours (or what the definition of |S| is). Where does it come from?

It'd be fine to change the constant because there's no validation test. So the best we can do is theorize, and theorization on this issue suggests changing the default to C=0.17. If one wanted to set up a validation test, it could be nice to reproduce Compte-Bellot and Corrsin (1964) (this wasn't possible when the closures were written because we didn't support triply periodic domains, but is possible now). This could also be used to validate Anisotropic Minimum Dissipation. Here's a figure with such a comparison from Rozema et al 2015:

I believe for shear flows it has also been found that the constant needs to be as small as C=0.1 or smaller? It'd be nice to mention and cite these in the docstring as well, if we can find those references.

Hope that helps!

[tomchor]

@tomchor I can't tell if the formula you've pasted is actually different from ours (or what the definition of |S| is). Where does it come from?

That was my bad, forgot to cite the source. It comes from Chamecki et al. 2019. But you're right that the definition is identical. I misinterpreted the code.

It'd be fine to change the constant because there's no validation test

I indeed think we should change constant default to match the most common value used in the literature. I'd also like to change the notation used in the code from ς to \Upsilon, so as to match the docs (just to make the code easier for contributors).

I'll open a PR.

[glwagner]

Sounds good to update the code! Perhaps we should also add a warning or note of some kind about the lack of validation (we've taken to doing this recently and I think it's a good thing).

@navidcy would be happy not to use ς!

\Upsilon is disturbingly similar to "Y":

julia> Υ = 1
1

maybe there's a better choice? We could also change the docs.

[tomchor]

Sounds good to update the code! Perhaps we should also add a warning or note of some kind about the lack of validation (we've taken to doing this recently and I think it's a good thing).

@navidcy would be happy not to use ς!

\Upsilon is disturbingly similar to "Y":

julia> Υ = 1
1

maybe there's a better choice? We could also change the docs.

Yeah, and I ended up realizing the same thing. I changed the docs instead.