More / more flexible Smagorinsky models

[tomchor]

I have been thinking that we should expand a bit our Smagorinsky models. Currently the only Smagorinsky-type model we have is the Smagorinsky-Lilly, which works well in some cases, notably the canonical horizontally-periodic turbulence, but has some downsides for more complex flows. For example when turbulence is not heterogeneous, or when there are persistent mean gradients in the flow. AMD is a good option in some cases, but often there are numerical artifacts in stratified, laminar regions in the flow (I don't think it's clear exactly what causes this, but it happens in many of my simulations, and I've herd similar things from others). So I think having at least one option to run a dynamic Smag would be good.

One possibility is that, as a start, we allow for user-defined Smagorinsky coefficients that could depend on the velocities. Something like

cₛ(i, j, k, grid, C::Function) = C(i, j, k, grid, u, v, w)

where C(i, j, k, grid, u, v, w) is provided by the user. That would allow the user a bit more control, and they could implement local dynamic Smag models, such as the dynamic gradient Smagorinsky. Something similar is currently done for AMD, although I think it's a bit more restrictive than what I'm proposing.

Soon, if there's interest, I'm planning on implementing a scale-independent Smagorinsky that's averaged in the symmetric flow directions (generally this is referred to as Planar-Averaged Scale-Independent Smag, but we wouldn't be limited to a planar average in our case). From then on, the extension to a Langrangian-Averaged Scale-Independent Smag would be straightforward.

Eventually it would be nice to extend these to a Scale-dependent model, but I think it's wise to start with the simpler cases.

How do people feel about this?

[glwagner]

You can start by allowing the coefficient to be an array as in AMD:

https://github.com/CliMA/Oceananigans.jl/blob/7f2d512ceaa2ce3fec3132cf0b08ab39344270d9/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl#L143

Then with a callback you can implement any model, including dynamic Smagorinsky. So this provides unlimited flexibility in developing a new closure.

I think its possible you will find that its more efficient to precompute (eg in a callback) the coefficient rather than compute it on the fly with a function, especially because the computation of the different filters is a bit non-local, eg you have to average over a few grid points, which could involve relatively expensive memory fetches.

A function-based interface could be nice too, but note that they are relatively complicated. Do we expect a lot of users to experiment with custom coefficient functions? @tomchor if this feature is just for you to experiment with dynamic Smagorinsky implementations then I think experimenting in the source code is an even quicker and easier option than designing a function-based interface...

@simone-silvestri and @xkykai have talked about dynamic Smagorinsky too but I think their main conclusion is that it is not a gimme.

[xkykai]

Here's a first attempt at Lagrangian-averaged dynamic Smagorinsky: https://github.com/CliMA/Oceananigans.jl/tree/ss-xk/dynamic-smagorinsky It doesn't run due to StackOverflowError, likely due to the way we are taking the averages. I think if we write custom functions to take each of the averages it would run, just haven't got to it yet.

[tomchor]

You can start by allowing the coefficient to be an array as in AMD:

https://github.com/CliMA/Oceananigans.jl/blob/7f2d512ceaa2ce3fec3132cf0b08ab39344270d9/src/TurbulenceClosures/turbulence_closure_implementations/anisotropic_minimum_dissipation.jl#L143

Then with a callback you can implement any model, including dynamic Smagorinsky. So this provides unlimited flexibility in developing a new closure.

I think its possible you will find that its more efficient to precompute (eg in a callback) the coefficient rather than compute it on the fly with a function, especially because the computation of the different filters is a bit non-local, eg you have to average over a few grid points, which could involve relatively expensive memory fetches.

I see what you mean. That's a good point. I think I'll start with that PR.

A function-based interface could be nice too, but note that they are relatively complicated. Do we expect a lot of users to experiment with custom coefficient functions? @tomchor if this feature is just for you to experiment with dynamic Smagorinsky implementations then I think experimenting in the source code is an even quicker and easier option than designing a function-based interface...

That's not my intention at all. I have plans to try to implement a dynamic option with or without the above flexibility I proposed. I do think there are legitimate reasons to allow for a user-defined coefficient. Like manually decreasing the coefficient in regions where it's overly diffusive, or when running wall-resolved LES. Basically, since the necessary change is a couple lines of code, I tend to think that flexibility is worth it. That said, the Array-based method you proposed can also cover these cases.

@simone-silvestri and @xkykai have talked about dynamic Smagorinsky too but I think their main conclusion is that it is not a gimme.

[glwagner]

That's not my intention at all. I have plans to try to implement a dynamic option with or without the above flexibility I proposed. I do think there are legitimate reasons to allow for a user-defined coefficient. Like manually decreasing the coefficient in regions where it's overly diffusive, or when running wall-resolved LES. Basically, since the necessary change is a couple lines of code, I tend to think that flexibility is worth it. That said, the Array-based method you proposed can also cover these cases.

Ok, but you proposed an interface that would include the velocity field. What are the use cases for that?

[tomchor]

Here's a first attempt at Lagrangian-averaged dynamic Smagorinsky: https://github.com/CliMA/Oceananigans.jl/tree/ss-xk/dynamic-smagorinsky It doesn't run due to StackOverflorError, likely due to the way we are taking the averages. I think if we write custom functions to take each of the averages it would run, just haven't got to it yet.

Nice, I wasn't aware of that effort. I see you're starting with the scale-independent version, correct?

I think starting with the lagrangian-averaged version is bold though. From what I've seen the lagrangian vs directional averaged versions of the model are very close in results for most LES we run for the ocean (i.e. doubly periodic), but the lagrangian average can be significantly more expensive. So I think for most cases we'd actually prefer the directionally-averaged one in order to save on computational costs. The use cases where lagrangian average really shines are the heterogeneous flows (like with immersed boundaries) and laminar-turbulent transitions.

Also, if you open a PR I can help you guys with that. I think I caught a couple of typos in the 𝒥ᴿᴺ and 𝒥ᴺᴺ forcing functions for example.

[xkykai]

Also, if you open a PR I can help you guys with that. I think I caught a couple of typos in the 𝒥ᴿᴺ and 𝒥ᴺᴺ forcing functions for example.

Done! https://github.com/CliMA/Oceananigans.jl/pull/3638

[glwagner]

Here's a first attempt at Lagrangian-averaged dynamic Smagorinsky: https://github.com/CliMA/Oceananigans.jl/tree/ss-xk/dynamic-smagorinsky It doesn't run due to StackOverflorError, likely due to the way we are taking the averages. I think if we write custom functions to take each of the averages it would run, just haven't got to it yet.

Nice, I wasn't aware of that effort. I see you're starting with the scale-independent version, correct?

I think starting with the lagrangian-averaged version is bold though. From what I've seen the lagrangian vs directional averaged versions of the model are very close in results for most LES we run for the ocean (i.e. doubly periodic), but the lagrangian average can be significantly more expensive. So I think for most cases we'd actually prefer the directionally-averaged one in order to save on computational costs. The use cases where lagrangian average really shines are the heterogeneous flows (like with immersed boundaries) and laminar-turbulent transitions.

Also, if you open a PR I can help you guys with that. I think I caught a couple of typos in the 𝒥ᴿᴺ and 𝒥ᴺᴺ forcing functions for example.

Why don't you implement DynamicSmagorinsky and abstract the notion of the Averaging to LagrangianAveraging and DirectionalAveraging.

Note, you should expect to spend your time validating the closures by running simulations --- not software design and abstraction. If you're having difficulties with the software design I can help.

[tomchor]

Why don't you implement DynamicSmagorinsky and abstract the notion of the Averaging to LagrangianAveraging and DirectionalAveraging.

Note, you should expect to spend your time validating the closures by running simulations --- not software design and abstraction. If you're having difficulties with the software design I can help.

I like this idea!