Implementing another SGS closure model #4559

anpa5279 · 2025-05-29T00:58:06Z

anpa5279
May 29, 2025

Hi developers,

I am hoping to add a new SGS model into Oceananigans developed by Peter Sullivan et al., 1994. However, from what I can understand, Oceananigans only manipulates the viscosity, not the entire term. In this paper, the SGS closure is defined as follows:

γ = isotropy factor (accounts for theoretical SGS constant due to anisotropy of the mean flow and controls the transition between SGS and averaged turbulence parameterizations)
ν_t = fluctuating eddy viscosity
ν_T = mean-field eddy viscosity
S_ij = resolved strain rate tensor
⟨S_ij⟩ = horizontal average of the resolved strain rate tensor

How should I implement this model? This can be a modified version of Smagorinsky or a one equation model like Deardorff (NCAR LES lets you chose between these two for this SGS). If we just think about it as a modified version of Smagorinsky, I think that the first term can be incorporated into the current Smagorinsky model like it is for the Smagorinsky-Lily model and then the SGS model adds the second term as a forcing function to the velocities and scalars.

On my own, I have been applying the classic Smagorinsky closure model as a forcing function to get a better understanding of how I can apply Sullivan's version in the future. It runs on GPUs, but have not confirmed if I did it correctly yet. Here is my work so far. I tried to mimic the closure model as much as possible for the forcing functions (∂ⱼ_τ₁ⱼ, ∂ⱼ_τ₂ⱼ, ∂ⱼ_τ₃ⱼ, and ∇_dot_qᶜ assuming Pr = 1):

using KernelAbstractions: @kernel, @index

using Oceananigans.Operators
using Oceananigans.TurbulenceClosures: Σ₁₁, Σ₂₂, Σ₃₃, Σ₁₂, Σ₁₃, Σ₂₃         
using Oceananigans.TurbulenceClosures: tr_Σ², Σ₁₂², Σ₁₃², Σ₂₃² 
using Oceananigans.Operators: Δy_qᶠᶜᶜ, Δx_qᶜᶠᶜ, Δx_qᶠᶜᶜ
using Oceananigans.Operators: volume
using Oceananigans.AbstractOperations: KernelFunctionOperation

# strain tensor squared and summed 
@inline function ΣᵢⱼΣᵢⱼᶜᶜᶜ(i, j, k, grid, u, v, w)
    s = tr_Σ²(i, j, k, grid, u, v, w)
    s .+ 2 * ℑxyᶜᶜᵃ(i, j, k, grid, Σ₁₂², u, v, w)
    s .+ 2 * ℑxzᶜᵃᶜ(i, j, k, grid, Σ₁₃², u, v, w)
    s .+ 2 * ℑyzᵃᶜᶜ(i, j, k, grid, Σ₂₃², u, v, w)
    return s
end

# viscosity
@kernel function smagorinsky_visc!(grid, velocities, νₑ)
    i, j, k = @index(Global, NTuple)

    u = velocities.u
    v = velocities.v
    w = velocities.w
    # Strain tensor dot product
    Σ² = ΣᵢⱼΣᵢⱼᶜᶜᶜ(i, j, k, grid, u, v, w)
    # Filter width
    Δ³ = Δxᶜᶜᶜ(i, j, k, grid) * Δyᶜᶜᶜ(i, j, k, grid) * Δzᶜᶜᶜ(i, j, k, grid)
    Δᶠ = cbrt(Δ³)
    C = 0.1
    cˢ² = C^2

    @inbounds νₑ[i, j, k] = cˢ² * Δᶠ^2 * sqrt(2Σ²)
end

# Horizontal viscous fluxes for isotropic diffusivities

@inline function viscous_flux_ux(i, j, k, grid, ν, u)
    return @inbounds -2 * ν[i, j, k] * Σ₁₁(i, j, k, grid, u)
end
@inline function viscous_flux_vx(i, j, k, grid, ν, u, v)
    return @inbounds -2 * ν[i, j, k] * Σ₁₂(i, j, k, grid, u, v)
end
@inline function viscous_flux_wx(i, j, k, grid, ν, u, w)
    return @inbounds -2 * ν[i, j, k] * Σ₁₃(i, j, k, grid, u, w)
end
@inline function viscous_flux_uy(i, j, k, grid, ν, u, v)
    return @inbounds -2 * ν[i, j, k] * Σ₁₂(i, j, k, grid, u, v)
end
@inline function viscous_flux_vy(i, j, k, grid, ν, v)
    return @inbounds -2 * ν[i, j, k] * Σ₂₂(i, j, k, grid, v)
end
@inline function viscous_flux_wy(i, j, k, grid, ν, v, w)
    return @inbounds -2 * ν[i, j, k] * Σ₂₃(i, j, k, grid, v, w)
end

# Vertical viscous fluxes for isotropic diffusivities
@inline function viscous_flux_uz(i, j, k, grid, ν, u, w)
    return @inbounds -2 * ν[i, j, k] * Σ₁₃(i, j, k, grid, u, w)
end
@inline function viscous_flux_vz(i, j, k, grid, ν, v, w)
    return @inbounds -2 * ν[i, j, k] * Σ₂₃(i, j, k, grid, v, w)
end
@inline function viscous_flux_wz(i, j, k, grid, ν, w)
    return @inbounds -2 * ν[i, j, k] * Σ₃₃(i, j, k, grid, w)
end

#diffusivity
@inline function diffusive_flux_x(i, j, k, grid, ν, c)
    return @inbounds - ν[i, j, k]  * ∂xᶠᶜᶜ(i, j, k, grid, c)
end 
@inline function diffusive_flux_y(i, j, k, grid, ν, c)
    return @inbounds - ν[i, j, k]  * ∂yᶜᶠᶜ(i, j, k, grid, c)
end
@inline function diffusive_flux_z(i, j, k, grid, ν, c)
    return @inbounds - ν[i, j, k]  * ∂zᶜᶜᶠ(i, j, k, grid, c)
end

#these are the discrete forcing functions
@inline function ∂ⱼ_τ₁ⱼ(i, j, k, grid, clock, model_fields)
    u = model_fields.u 
    v = model_fields.v
    w = model_fields.w
    ν = model_fields.νₑ
    return -1 / Vᶠᶜᶜ(i, j, k, grid) * (δxᶠᵃᵃ(i, j, k, grid, Ax_qᶜᶜᶜ, viscous_flux_ux, ν, u) +
                                      δyᵃᶜᵃ(i, j, k, grid, Ay_qᶠᶠᶜ, viscous_flux_uy, ν, u, v) +
                                      δzᵃᵃᶜ(i, j, k, grid, Az_qᶠᶜᶠ, viscous_flux_uz, ν, u, w))
end

@inline function ∂ⱼ_τ₂ⱼ(i, j, k, grid, clock, model_fields)
    u = model_fields.u 
    v = model_fields.v
    w = model_fields.w
    ν = model_fields.νₑ
    return -1 / Vᶜᶠᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶠᶜ, viscous_flux_vx, ν, u, v) +
                                      δyᵃᶠᵃ(i, j, k, grid, Ay_qᶜᶜᶜ, viscous_flux_vy, ν, v) +
                                      δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶠᶠ, viscous_flux_vz, ν, v, w))
end

@inline function ∂ⱼ_τ₃ⱼ(i, j, k, grid, clock, model_fields)
    u = model_fields.u 
    v = model_fields.v
    w = model_fields.w
    ν = model_fields.νₑ
    return -1 / Vᶜᶜᶠ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶠ, viscous_flux_wx, ν, u, w) +
                                      δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶠ, viscous_flux_wy, ν, v, w) +
                                      δzᵃᵃᶠ(i, j, k, grid, Az_qᶜᶜᶜ, viscous_flux_wz, ν, w))
end

@inline function ∇_dot_qᶜ(i, j, k, grid, clock, model_fields)
    scalar = model_fields.T
    ν = model_fields.νₑ
    return -1/Vᶜᶜᶜ(i, j, k, grid) * (δxᶜᵃᵃ(i, j, k, grid, Ax_qᶠᶜᶜ, diffusive_flux_x, ν, scalar) +
                                    δyᵃᶜᵃ(i, j, k, grid, Ay_qᶜᶠᶜ, diffusive_flux_y, ν, scalar) +
                                    δzᵃᵃᶜ(i, j, k, grid, Az_qᶜᶜᶠ, diffusive_flux_z, ν, scalar))
end

In main script:

u_SGS = Forcing(∂ⱼ_τ₁ⱼ, discrete_form=true)
v_SGS = Forcing(∂ⱼ_τ₂ⱼ, discrete_form=true)
w_SGS = Forcing(∂ⱼ_τ₃ⱼ, discrete_form=true)
T_SGS = Forcing(∇_dot_qᶜ, discrete_form=true)

νₑ = Field{Center, Center, Center}(grid)

model = NonhydrostaticModel(; grid, buoyancy, #coriolis,
                            advection = WENO(),
                            timestepper = :RungeKutta3,
                            tracers = (:T),
                            closure = nothing, #closure = Smagorinsky(coefficient=0.1)
                            stokes_drift = UniformStokesDrift(∂z_uˢ=new_dUSDdz),
                            boundary_conditions = (u=u_bcs, T=T_bcs),
                            forcing = (u=u_SGS, v = v_SGS, w = w_SGS, T = T_SGS),
                            auxiliary_fields = (νₑ = νₑ,))

function update_viscosity(model)
    u = model.velocities.u
    v = model.velocities.v
    w = model.velocities.w
    grid = model.grid
    νₑ = model.auxiliary_fields.νₑ
    launch!(arch, grid, :xyz, smagorinsky_visc!, grid, u, v, w, νₑ)
end 
simulation.callbacks[:visc_update] = Callback(update_viscosity, IterationInterval(1), callsite=UpdateStateCallsite())

I would like advice on how to implement Sullivan's SGS model. If this is not feasible to incorporate it into the current version of Oceananigans, is there a clean way to do this with discrete forcing functions? I need help on how to incorporate this for my own research because my current idea ignores the Deardorff option. Do you have any suggestions? Thanks!

glwagner · 2025-06-03T15:31:26Z

glwagner
Jun 3, 2025
Maintainer

There is nothing "special" about a closure --- anything that can be implemented as a forcing can be implemented as a closure and vice versa. The definition of a "closure" involves specifying the tendencies, eg:

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl

Line 98 in 267cb45

    
           - ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, diffusivities, clock, closure_model_fields, buoyancy)

as well as, optionally, a function called compute_diffusivities! which gives you the option to precompute certain things --- such as an eddy viscosity, or the entire subgrid stress.

This encompasses all explicitly computed closures, including the one you've written.

In addition to the generic closure interface, which consists of tendency functions + compute_diffusivities!, there is a higher-level interface via AbstractScalarDiffusivity. This more specific interface for defining closures allows implementers to define just an eddy viscosity / diffusivity. The AbstractScalarDiffusivity interface is implemented here:

https://github.com/CliMA/Oceananigans.jl/blob/main/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl

But note this is a convenience. It's not required to use AbstractScalarDiffusivity.

Turning to the question at hand, I guess you are interested in implementing a kind of closure which involves two terms: one involving a term computed from a typical LES closure, and a second term which involves horizontal reductions. This suggests that you want to use a design that's something like

struct BridgingTurbulenceClosure
    isotropic_closure
    horizontally_averaged_closure
    isotropy_factor
end

You also need to design a corresponding object to hold the diffusivity_fields.

Then you can write your tendency kernel function as something like

function ∂ⱼ_τ₁ⱼ(i, j, k, grid, bridging::BridgingTurbulenceClosure, diffusivities, args...)
    ∂ⱼ_τ₁ⱼ_isotropic = ∂ⱼ_τ₁ⱼ(i, j, k, grid, bridging.isotropic_closure, diffusivities.isotropic_closure, args...)
    ∂ⱼ_τ₁ⱼ_averaged = ∂ⱼ_τ₁ⱼ(i, j, k, grid, bridging.horizontally_averaged_closure, diffusivities.horizontally_averaged_closure, args...)
    return ∂ⱼ_τ₁ⱼ_isotropic + ∂ⱼ_τ₁ⱼ_averaged
end

You can also probably re-use the compute_diffusivities! for each component closure, provided you scale the computed eddy viscosity and diffusivities by the isotropy factor afterwards.

Here are some additional thoughts:

While there are many options for isotropic_closure currently implemented, I am not sure what you would plan to use for horizontally_averaged_closure here. You may want to implement that first, otherwise there is no way to use BridgingTurbulenceClosure.
I have anecdotally heard of issues with closures like this (that they "do more harm than good") when trying to correct LES. More broadly we find that almost all closures are overly-diffusive and occasionally produce misleading results compared to a simple WENO(order=9) advection scheme. The main motivation to use an explicit closure seems to be if one requires explicit computation of the dissipation rate. But, you may have a case for which WENO does not work and therefore have a specific need for this kind of explicit bridging closure (that would be interesting to know about!)
Implementing the Deardorff TKE-based closure is not too difficult in principle, but may be computationally expensive and overly-diffusive compared to the alternatives, such as AnisotropicMinimumDissipation, or the simple WENO(order=9) advection scheme mentioned above.

0 replies

anpa5279 · 2025-06-11T15:05:01Z

anpa5279
Jun 11, 2025
Author

Thank you! You gave me a lot to think about. I have updated my Smagorinsky forcing function and I am now matching the closure model version. I updated my forcing function to calculate at the faces of the grid cell for the respective direction like Oceananigan's Smagorinsky does in AbstractScalarDiffusivity and made sure the halos are being filled in. I will try to implement your struct suggestion soon.

Quick clarifying question, which aspect of the SGS model are you referring to does more harm than good? From my preliminary search, I have not found much in the literature as to why. Sullivan et al., 1994 does acknowledge that closures can be too dissipative (they reference Moin and Kim, 1982) and they try to correct this in their paper. However, they do not compare to a DNS or field data, but to Monin-Obukhov similarity forms.

2 replies

glwagner Jun 11, 2025
Maintainer

Can you open a PR with your contribution?

@BrodiePearson and @qingli411 can give more information about NCAR LES closures. The specific statement I heard was from Baylor FK but I don't know if he is on github?

anpa5279 Jun 27, 2025
Author

Hi! I was implementing it in hopes to have a better "apples to apples" comparison to NCAR LES, but I decided to implement Smagorinsky into NCAR LES instead of trying to add a new closure to Oceananigans. After talking to my advisor, I no longer have outstanding questions about my previous comment. Thanks for your help! Do you still want me to do a PR? I did not complete implementation because I changed NCAR LES instead. I only have my SGS forcing function testing completed.

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Uh oh!

Implementing another SGS closure model #4559

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Replies: 2 comments 2 replies

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{editor}}'s edit

{{editor}}'s edit

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Uh oh!

{{title}}

Uh oh!

Select a reply

Uh oh!

Implementing another SGS closure model #4559

Uh oh!

Uh oh!

anpa5279 May 29, 2025

Replies: 2 comments · 2 replies

Uh oh!

Uh oh!

glwagner Jun 3, 2025 Maintainer

Uh oh!

anpa5279 Jun 11, 2025 Author

Uh oh!

glwagner Jun 11, 2025 Maintainer

Uh oh!

anpa5279 Jun 27, 2025 Author

anpa5279
May 29, 2025

Replies: 2 comments 2 replies

glwagner
Jun 3, 2025
Maintainer

anpa5279
Jun 11, 2025
Author

glwagner Jun 11, 2025
Maintainer

anpa5279 Jun 27, 2025
Author