Sharing turbulence closures between compressible and incompressible models

[ali-ramadhan]

@thabbott and I have been thinking about potential ways we can reuse LES/turbulence closures between incompressible and compressible models.

One approach we came up is for all turbulence closures to define functions for the SGS stress tensor like ∂ⱼτᵢⱼ where τᵢⱼ = νˢᵍˢSᵢⱼ and where νˢᵍˢ will dispatch on the type of closure while the strain-rate tensor Sᵢⱼ will dispatch on whether the model is compressible or incompressible. For compressible models, Sᵢⱼ will instead compute and return the traceless strain-rate tensor Sᵢⱼ - ⅓Sₖₖδᵢⱼ.

A turbulence closure can define ∂ⱼτᵢⱼ directly to elide or short circuit computations or in case the SGS stress tensor is not of the form ∂ⱼτᵢⱼ = νˢᵍˢSᵢⱼ.

Still learning more about this topic but it seems that as a start, just using the traceless strain-rate tensor for compressible models allows us to reuse a closure from an incompressible model.

@thabbott suggested that νˢᵍˢ might have to dispatch on the thermodynamic variable as well. Is this true for the viscosity or only the diffusivities κˢᵍˢ? I.e. do we have to treat SGS fluxes of moisture and energy in a special way?

Couple of useful references for compressible LES (mostly with dynamic Smagorinsky):

1. Moin et al. (1991), A dynamic subgrid‐scale model for compressible turbulence and scalar transport: https://doi.org/10.1063/1.858164
2. Chai & Mahesh (2012), Dynamic k-equation model for large-eddy simulation of compressible flows: https://doi.org/10.1017/jfm.2012.115

---

Some notational and misc. changes I'd like to make:

1. Rename ∂ⱼ_2ν_Σ₁ⱼ to ∂ⱼτ₁ⱼ as τᵢⱼ = 2νSᵢⱼ is would no longer be true for CompressibleModel, it's τᵢⱼ = 2νSᵢⱼ + λSₘₘδᵢⱼ where λ is the second viscosity, or τᵢⱼ = 2μ(Sᵢⱼ - ⅓Sₖₖδᵢⱼ) + μᵥSₖₖδᵢⱼ where μᵥ is the bulk or volume viscosity.

2. For constant isotropic viscosity: rename ∂ⱼνᵢⱼ∂ᵢu to ν∇²u. For constant isotropic diffusivity: rename ∂ⱼκᵢⱼ∂ᵢc to κ∇²c.

3. For constant anistropic viscosity: rename ∂ⱼνᵢⱼ∂ᵢu to νⱼ∂ⱼ²u. For constant anistropic diffusivity: rename ∂ⱼκᵢⱼ∂ᵢc to κⱼ∂ⱼ²c. Potentially a misuse of summation notation but should be a valid and understandable use of it.

4. For ConstantAnisotropicDiffusivity and AnisotropicBiharmonicDiffusivity should we generalize it to be truly anisotropic with three diagonal components?

5. Use Sᵢⱼ instead of Σᵢⱼ to denote the strain-rate tensor, following a common convention in most books, papers, and Wikipedia. I believe Σ was used initially to not conflict with salinity S but I think the use of salinity has is now mainly restricted to Oceananigans.Buoyancy so there is a very small chance for confusion in Oceananigans.TurbulenceClosures.

[glwagner]

My decision to use Sigma rather than S was not to avoid confusion in the code, but to avoid confusion in documentation and papers, where it will not be possible to use S for both salinity and strain rate.

[glwagner]

I support changes 1-4. 1 may also be important for future turbulence closures in IncompressibleModel.

[glwagner]

One approach we came up is for all turbulence closures to define functions for the SGS stress tensor like ∂ⱼτᵢⱼ

With you here.

where τᵢⱼ = νˢᵍˢSᵢⱼ and where νˢᵍˢ will dispatch on the type of closure while the strain-rate tensor Sᵢⱼ will dispatch on whether the model is compressible or incompressible.

I don't think this needs to be pre-ordained somehow. Each turbulence closure can act as needed. What's being pointed out is that it's easy to use type information with the turbulence closure so that it can be used as a valid closure for either incompressible or compressible systems. For this I think we should add type information to the turbulence closure itself, eg:

struct VerstappenAnisotropicMinimumDissipation{FT, PK, PN, K} <: AbstractAnisotropicMinimumDissipation{FT, S}
    Cν :: PN
    Cκ :: PK
    Cb :: FT
     ν :: FT
     κ :: K
strain_rate :: S
    function VerstappenAnisotropicMinimumDissipation{FT}(Cν, Cκ, Cb, ν, κ, strain_rate) where FT
        return new{FT, typeof(Cκ), typeof(Cν), typeof(κ)}(Cν, Cκ, Cb, ν, convert_diffusivity(FT, κ), strain_rate)
    end
end

where strain_rate = Traceless() would inform the turbulence closure to use the traceless rate of strain.

Stress divergences are computed here:

https://github.com/climate-machine/Oceananigans.jl/blob/042f8cf12928b53b36cfa2a61a6b6a8ec7dc71bc/src/TurbulenceClosures/closure_operators.jl#L38

For example, one component of the stress divergence is:

@inline ∂x_2ν_Σ₁₁(i, j, k, grid, closure, U, diffusivities) =
    2 * ∂xᶠᵃᵃ(i, j, k, grid, ν_σᶜᶜᶜ, diffusivities.νₑ, Σ₁₁, U.u, U.v, U.w)

We can thus extend functionality via dispatch:

@inline ∂x_2ν_Σ₁₁(i, j, k, grid, ::TracelessStrainRateClosure, U, diffusivities) =
    2 * ∂xᶠᵃᵃ(i, j, k, grid, ν_σᶜᶜᶜ, diffusivities.νₑ, traceless_Σ₁₁, U.u, U.v, U.w)

where TracelessStrainRateClosure is an alias for a closure that requires the traceless strain rate.

A turbulence closure can define ∂ⱼτᵢⱼ directly to elide or short circuit computations or in case the SGS stress tensor is not of the form ∂ⱼτᵢⱼ = νˢᵍˢSᵢⱼ.

This is how the system already works for closures that don't have eddy diffusivities. Is something different from the current system being proposed?

@thabbott suggested that νˢᵍˢ might have to dispatch on the thermodynamic variable as well. Is this true for the viscosity or only the diffusivities κˢᵍˢ? I.e. do we have to treat SGS fluxes of moisture and energy in a special way?

Currently, the calculate_diffusivities! functions receives buoyancy, velocities, and tracers, as arguments:

https://github.com/climate-machine/Oceananigans.jl/blob/042f8cf12928b53b36cfa2a61a6b6a8ec7dc71bc/src/TurbulenceClosures/turbulence_closure_implementations/verstappen_anisotropic_minimum_dissipation.jl#L176

If a different abstraction is needed we can change the calculate_diffusivities function appropriately.

As for ensuring that compressible models always use the correct strain rate, this can be handled within the CompressibleModel constructor.

[glwagner]

We don't have a compressible model so I'm closing this.