Smagorinsky closure with a damping function

[j-ostiguy]

Hello,

I would like to use the Smagorinsky closure, but where the length scale includes a van Driest damping function (i.e. l = CΔ[1-exp(-y+/A+)] ). Is there a way to do this?

Thank you!

[glwagner]

The length scale (which we call the "filter width") is implemented here:

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

Lines 114 to 116 in f3996c2

	# Filter width
	Δ³ = Δxᶜᶜᶜ(i, j, k, grid) * Δyᶜᶜᶜ(i, j, k, grid) * Δzᶜᶜᶜ(i, j, k, grid)
	Δᶠ = cbrt(Δ³)

To generalize this requires a bit of source code development (not much...)

Another route is to generalize the "model constant", which multiplies the length-scale, eg

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

Lines 117 to 119 in f3996c2

	C = closure.C # free parameter
	@inbounds νₑ[i, j, k] = ς * (C * Δᶠ)^2 * sqrt(2Σ²)

We've provided such a feature for AMD:

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl

Lines 140 to 144 in f3996c2

	# Dispatch on the type of the user-provided AMD model constant.
	# Only numbers, arrays, and functions supported now.
	@inline Cᴾᵒⁱⁿ(i, j, k, grid, C::Number) = C
	@inline Cᴾᵒⁱⁿ(i, j, k, grid, C::AbstractArray) = @inbounds C[i, j, k]
	@inline Cᴾᵒⁱⁿ(i, j, k, grid, C::Function) = C(xnode(i, grid, Center()), ynode(j, grid, Center()), znode(k, grid, Center()))

What are Δ, C, y+, and A+ in your expression? Personally, I think it makes sense to support this kind of thing by generalizing the model constant calculation as we've done for AMD.

[glwagner]

I'll do something simple to start to illustrate how these kinds of things are done. But people who use this feature should feel free to implement a more general and sophisticated interface if they want that.

[glwagner]

Let's move the discussion to #3370

[j-ostiguy]

What is δ_ν? Viscosity has units of length and time, so we'd need a time unit to convert viscosity into a length scale.

δ_ν = ν*sqrt(ρ/ τ_0) = ν/u_τ is the viscous length scale (where τ_0 is the stress at the wall and u_τ is the friction velocity). I'm not sure I understand why viscosity has to be converted to a length scale.

Distance-to-wall has to be defined by the user I think. It's a little annoying for irregular boundaries. We can get to that.

y+, the distance to the wall, can be written as y+ = u_τ*y/ν. So, it's a function of the shear stress at the boundary - can you get that from the model? Or I could calculate u_τ from the momentum flux at the boundary.

[glwagner]

δ_ν = ν*sqrt(ρ/ τ_0) = ν/u_τ is the viscous length scale (where τ_0 is the stress at the wall and u_τ is the friction velocity). I'm not sure I understand why viscosity has to be converted to a length scale.

Maybe I miscommunicated, but u_τ is precisely the scale (here a velocity scale) that mediates between viscosity and the length scale δ_ν. So you provided what I asked.

If u_τ is not prescribed by forcing, then we would have to estimate the friction velocity dynamically based on the near-wall velocity gradient, right? If we need to do that, things are a little tricker (we probably want to compute the friction velocity within a callback). Let me know if you'd like me to show you how to do that. I provided an example in the PR that assumes the near-wall length scale is known a priori.

[glwagner]

So, it's a function of the shear stress at the boundary - can you get that from the model? Or I could calculate u_τ from the momentum flux at the boundary.

Indeed, if the shear stress at the wall is associated with a no-slip boundary condition, then it is calculated on-the-fly within the boundary tendency calculation, and is not stored anywhere. In order to use it for a subsequent computation, we'd have to re-compute it within a callback.

It also seems that one might want to, for example, average the momentum flux in horizontal / homogeneous directions before using it for this model constant calculation. But I'm not sure.

Either way, such a model is probably best implemented by pre-allocating an updatable object for the friction velocity (eg using Ref) and, using that updatable object within the FunctionField that computes the model constant, and then updating the object each time-step using a Callback within simulation.