Implementation of the scale-dependent Lagrangian-averaged dynamic Smagorinsky #3638

xkykai · 2024-06-24T19:08:00Z

This is an implementation of the scale-dependent Lagrangian-averaged dynamic Smagorinsky model as outlined in section C III. in Bou-Zeid et al. 2005.

@simone-silvestri @tomchor

xkykai · 2024-06-24T19:10:55Z

src/TurbulenceClosures/turbulence_closure_implementations/lagrangian_dynamic_smagorinsky.jl

@@ -225,6 +225,8 @@ end
    u̅w̅ = ℑxyzᶜᶜᶜ(i, j, k, grid, ℑxyzᶠᶠᶠ, ℑxyzᶜᶜᶜ, ℑxyzᶠᶠᶠ, uᵢuⱼ, ℑxᶜᵃᵃ, ℑzᵃᵃᶜ, u, w)
    v̅w̅ = ℑxyzᶜᶜᶜ(i, j, k, grid, ℑxyzᶠᶠᶠ, ℑxyzᶜᶜᶜ, ℑxyzᶠᶠᶠ, uᵢuⱼ, ℑyᵃᶜᵃ, ℑzᵃᵃᶜ, v, w)

+    S̅ₘ = ℑxyzᶜᶜᶜ(i, j, k, grid, ℑxyzᶠᶠᶠ, ℑxyzᶜᶜᶜ, ℑxyzᶠᶠᶠ, ϕ⁰⁵, ΣᵢⱼΣᵢⱼᶜᶜᶜ, u, v, w)


Currently the code goes into StackOverflowError, likely due to the way the filter averaging is done as we might have function compositions which are too deep.

Does it work if we use a simpler filter?

we can maybe use a box filter here instead of all this combination of interpolations.
We can add a filter function to the struct

Hmm. But this should not be an overflow error, right? Are we hitting some heuristic?

I would try combining steps to see what happens. For example you could define

ϕ⁰⁵_ΣᵢⱼΣᵢⱼᶜᶜᶜ(i, j, k, grid, u, v, w) = ϕ⁰⁵(i, j, k, grid, ΣᵢⱼΣᵢⱼᶜᶜᶜ, u, v, w)

You could also simply write the whole thing manually. With composition this deep, it might also make this code easier to understand.

Does it work if we use a simpler filter?

What is complicated about this?

tomchor · 2024-06-24T19:16:10Z

We should "Close #3637" with this one I think

tomchor · 2024-06-24T19:18:03Z

src/TurbulenceClosures/turbulence_closure_implementations/lagrangian_dynamic_smagorinsky.jl

+    M₁₁ = 2Δᶠ * (S̅S̅₁₁ - α^2 * β^2 * S̅ₘ * S̅₁₁)
+    M₂₂ = 2Δᶠ * (S̅S̅₂₂ - α^2 * β^2 * S̅ₘ * S̅₂₂)
+    M₃₃ = 2Δᶠ * (S̅S̅₃₃ - α^2 * β^2 * S̅ₘ * S̅₃₃)
+    M₁₂ = 2Δᶠ * (S̅S̅₁₂ - α^2 * β^2 * S̅ₘ * S̅₁₂)
+    M₁₃ = 2Δᶠ * (S̅S̅₁₃ - α^2 * β^2 * S̅ₘ * S̅₁₃)
+    M₂₃ = 2Δᶠ * (S̅S̅₂₃ - α^2 * β^2 * S̅ₘ * S̅₂₃)


Should this be

Suggested change

M₁₁ = 2Δᶠ * (S̅S̅₁₁ - α^2 * β^2 * S̅ₘ * S̅₁₁)

M₂₂ = 2Δᶠ * (S̅S̅₂₂ - α^2 * β^2 * S̅ₘ * S̅₂₂)

M₃₃ = 2Δᶠ * (S̅S̅₃₃ - α^2 * β^2 * S̅ₘ * S̅₃₃)

M₁₂ = 2Δᶠ * (S̅S̅₁₂ - α^2 * β^2 * S̅ₘ * S̅₁₂)

M₁₃ = 2Δᶠ * (S̅S̅₁₃ - α^2 * β^2 * S̅ₘ * S̅₁₃)

M₂₃ = 2Δᶠ * (S̅S̅₂₃ - α^2 * β^2 * S̅ₘ * S̅₂₃)

M₁₁ = 2Δᶠ * (S̅S̅₁₁ - α^2 * β * S̅ₘ * S̅₁₁)

M₂₂ = 2Δᶠ * (S̅S̅₂₂ - α^2 * β * S̅ₘ * S̅₂₂)

M₃₃ = 2Δᶠ * (S̅S̅₃₃ - α^2 * β * S̅ₘ * S̅₃₃)

M₁₂ = 2Δᶠ * (S̅S̅₁₂ - α^2 * β * S̅ₘ * S̅₁₂)

M₁₃ = 2Δᶠ * (S̅S̅₁₃ - α^2 * β * S̅ₘ * S̅₁₃)

M₂₃ = 2Δᶠ * (S̅S̅₂₃ - α^2 * β * S̅ₘ * S̅₂₃)

?

Can you paste in the reference that shows this?

Taking a second look at things I now think there may be more typos than this one. So I suggest @xkykai and @simone-silvestri double-check this, but if I understand correctly, what's called 𝒥ᴿᴺ in this PR is called $\mathcal J^{QN}$in Bou-Zeid et al. (2005), and here's its advection equation:

so using the language of this PR the forcing here should be

1 / T4Δ * (QᵢⱼNᵢⱼ - 𝒥ᴿᴺ)

instead of what's written below: 1 / TΔ * (LᵢⱼMᵢⱼ - 𝒥ᴿᴺ).

And the functions $Q_{ij}$ and $N_{ij}$ are:

Similar comments are valid for the function below (𝒥ᴺᴺ_forcing_function()).

For the record, my original comment had to do with $M_{ij}$ being implemented as M₁₁ = 2Δᶠ * (S̅S̅₁₁ - α^2 * β^2 * S̅ₘ * S̅₁₁) when it should be

but now (per my previous comment) I realize that $M_{ij}$ isn't supposed to be in this forcing function at all.

The difference between Q and L (and M and N) is the filter size. The stresses for the 𝒥ᴿᴺ forcings are averaged twice as much as the 𝒥ᴸᴹ forcings (if you look at how the averages are calculated). However, I would say that to push this PR forward we need to:

define an averaging operators (to be held in the closure struct) for the M - L fluxes and for the R - N fluxes

remove all this combination of interpolations in favor of those averaging operators

The difference between Q and L (and M and N) is the filter size.

Agreed. The $N$ and $M$ nomenclature difference makes it clearer though. That said, the combination $\alpha^2 \beta^2$ shouldn't appear anywhere, right? It's either $\alpha^4 \beta^2$ or $\alpha^2 \beta$ (unless you hard-coded some sutff I guess)

What do the tilde, hat, and bar mean?

Tilde: LES-filtered velocity. In our case it's simply the velocities that come out of model.velocities

Bar: First test filter

Hat: Second test filter (at a larger spatial scale than the first)

tomchor · 2024-06-24T19:19:06Z

src/TurbulenceClosures/turbulence_closure_implementations/lagrangian_dynamic_smagorinsky.jl

+    M₁₁ = 2Δᶠ * (S̅S̅₁₁ - α^2 * β^2 * S̅ₘ * S̅₁₁)
+    M₂₂ = 2Δᶠ * (S̅S̅₂₂ - α^2 * β^2 * S̅ₘ * S̅₂₂)
+    M₃₃ = 2Δᶠ * (S̅S̅₃₃ - α^2 * β^2 * S̅ₘ * S̅₃₃)
+    M₁₂ = 2Δᶠ * (S̅S̅₁₂ - α^2 * β^2 * S̅ₘ * S̅₁₂)
+    M₁₃ = 2Δᶠ * (S̅S̅₁₃ - α^2 * β^2 * S̅ₘ * S̅₁₃)
+    M₂₃ = 2Δᶠ * (S̅S̅₂₃ - α^2 * β^2 * S̅ₘ * S̅₂₃)


Same comment here

Suggested change

M₁₁ = 2Δᶠ * (S̅S̅₁₁ - α^2 * β^2 * S̅ₘ * S̅₁₁)

M₂₂ = 2Δᶠ * (S̅S̅₂₂ - α^2 * β^2 * S̅ₘ * S̅₂₂)

M₃₃ = 2Δᶠ * (S̅S̅₃₃ - α^2 * β^2 * S̅ₘ * S̅₃₃)

M₁₂ = 2Δᶠ * (S̅S̅₁₂ - α^2 * β^2 * S̅ₘ * S̅₁₂)

M₁₃ = 2Δᶠ * (S̅S̅₁₃ - α^2 * β^2 * S̅ₘ * S̅₁₃)

M₂₃ = 2Δᶠ * (S̅S̅₂₃ - α^2 * β^2 * S̅ₘ * S̅₂₃)

M₁₁ = 2Δᶠ * (S̅S̅₁₁ - α^2 * β * S̅ₘ * S̅₁₁)

M₂₂ = 2Δᶠ * (S̅S̅₂₂ - α^2 * β * S̅ₘ * S̅₂₂)

M₃₃ = 2Δᶠ * (S̅S̅₃₃ - α^2 * β * S̅ₘ * S̅₃₃)

M₁₂ = 2Δᶠ * (S̅S̅₁₂ - α^2 * β * S̅ₘ * S̅₁₂)

M₁₃ = 2Δᶠ * (S̅S̅₁₃ - α^2 * β * S̅ₘ * S̅₁₃)

M₂₃ = 2Δᶠ * (S̅S̅₂₃ - α^2 * β * S̅ₘ * S̅₂₃)

tomchor · 2024-06-24T19:22:30Z

I'm curious as to why you picked α=4. It makes for a pretty big stencil and α=2 works pretty well already. Does it have to do with the filtering being two consecutive interpolations?

simone-silvestri · 2024-06-24T19:24:11Z

This is also not super clean because the forcing functions for the additional tracers are added manually to the model.
It would be better to have a way to include them directly when prescribing this closure.

xkykai · 2024-06-24T19:28:39Z

I believe what we are doing here is using two filter scales α=2 and α=4 to compute β, which is ratio of coefficients at different filter scales Δ, αΔ, and α²Δ, assuming that β is scale-invariant.

tomchor · 2024-06-24T19:34:08Z

I believe what we are doing here is using two filter scales α=2 and α=2 to compute β, which is ratio of coefficients at different filter scales Δ, αΔ, and α²Δ, assuming that β is scale-invariant.

Hmm, then I'm confused. When you assume β to be scale-invariant you only need two filters: the grid (Δ) and a test filter (αΔ). afaik the second test filter is only needed when assuming that β can vary between scales.

xkykai · 2024-06-24T19:45:05Z

afaik the second test filter is only needed when assuming that β can vary between scales.

I think you're correct here! Sorry for the confusion here, but I suppose by scale-invariance I mean β has a power law dependence on the scale size. And since we were trying to implement the scale-dependent version, β would vary across different scales? Hence we need the α=4 filter as well. Am I understanding Equation (26) in Bou-Zeid et al. 2005 correctly here?

tomchor · 2024-06-24T19:52:37Z

afaik the second test filter is only needed when assuming that β can vary between scales.

I think you're correct here! Sorry for the confusion here, but I suppose by scale-invariance I mean β has a power law dependence on the scale size. And since we were trying to implement the scale-dependent version, β would vary across different scales? Hence we need the α=4 filter as well. Am I understanding Equation (26) in Bou-Zeid et al. 2005 correctly here?

Ah, I see, I was under the impression you were trying to implement the scale-invariant version, but it does make sense to if you're implementing the scale-dependent version.

glwagner · 2024-06-25T17:59:32Z

validation/stratified_couette_flow/stratified_couette_flow.jl

+Pr=0.7
+Ri = 0.5
+Ni=10
+end_time=1000


What is going on here? Isn't there already a script for invoking the utilities in this script?

glwagner · 2024-06-25T18:01:27Z

validation/stratified_couette_flow/stratified_couette_flow.jl

@@ -113,7 +125,7 @@ function simulate_stratified_couette_flow(; Nxy, Nz, arch=GPU(), h=1, U_wall=1,
    ##### Impose boundary conditions
    #####

-    grid = RectilinearGrid(arch, size = (Nxy, Nxy, Nz), extent = (4π*h, 2π*h, 2h))
+    grid = RectilinearGrid(arch, size = (Nxy, Nxy, Nz), extent = (4π*h, 2π*h, 2h), halo = (6, 6, 6))


Suggested change

grid = RectilinearGrid(arch, size = (Nxy, Nxy, Nz), extent = (4π*h, 2π*h, 2h), halo = (6, 6, 6))

grid = RectilinearGrid(arch; size = (Nxy, Nxy, Nz), extent = (4π*h, 2π*h, 2h), halo)

make this cleaner by adding halo as an argument to simulate_stratified_couette_flow.

Also, it doesn't look like the closure is an argument to simulate_stratified_couette_flow. You'll need that too. (You can also infer the halo size.)

Just as a general tip, I think you will save time if you generalize these functions first. Now, you are in the position of having to 1) walk back the changes you made and then 2) implement an interface that allows different closures. This is actually more work then if you had simply done (2) to begin with.

glwagner · 2024-07-07T23:19:22Z

This is also not super clean because the forcing functions for the additional tracers are added manually to the model. It would be better to have a way to include them directly when prescribing this closure.

I suggest implementing the tracer-specific forcing in a similar manner that we do for CATKE and k-epsilon

simone-silvestri added 11 commits March 1, 2024 15:49

bugfix

6bac83d

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

a077af4

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

d0215e9

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

b9fbdbc

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

fea4bed

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

c4bfcd8

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

02bfe37

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

e1be27f

Merge branch 'main' of github.com:CliMA/Oceananigans.jl

bd4fca1

first try

921b333

first try goes in stack overflow

0b4825e

xkykai commented Jun 24, 2024

View reviewed changes

xkykai added numerics 🧮 So things don't blow up and boil the lobsters alive turbulence closures 🎐 labels Jun 24, 2024

xkykai mentioned this pull request Jun 24, 2024

More / more flexible Smagorinsky models #3637

Open

tomchor reviewed Jun 24, 2024

View reviewed changes

glwagner reviewed Jun 25, 2024

View reviewed changes

tomchor mentioned this pull request Jun 28, 2024

Implement directionally-averaged, scale-invariant dynamic Smagorinsky closure #3642

Draft

5 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Implementation of the scale-dependent Lagrangian-averaged dynamic Smagorinsky #3638

Implementation of the scale-dependent Lagrangian-averaged dynamic Smagorinsky #3638

xkykai commented Jun 24, 2024 •

edited

Loading

xkykai Jun 24, 2024

tomchor Jun 24, 2024

simone-silvestri Jun 24, 2024

glwagner Jun 25, 2024

glwagner Jun 25, 2024

tomchor commented Jun 24, 2024

tomchor Jun 24, 2024

glwagner Jun 27, 2024

tomchor Jul 5, 2024

tomchor Jul 5, 2024

simone-silvestri Jul 5, 2024

tomchor Jul 6, 2024

glwagner Jul 7, 2024

tomchor Jul 8, 2024

tomchor Jun 24, 2024 •

edited

Loading

tomchor commented Jun 24, 2024

simone-silvestri commented Jun 24, 2024

xkykai commented Jun 24, 2024 •

edited

Loading

tomchor commented Jun 24, 2024

xkykai commented Jun 24, 2024 •

edited

Loading

tomchor commented Jun 24, 2024

glwagner Jun 25, 2024

glwagner Jun 25, 2024 •

edited

Loading

glwagner commented Jul 7, 2024

	grid = RectilinearGrid(arch, size = (Nxy, Nxy, Nz), extent = (4πh, 2πh, 2h), halo = (6, 6, 6))
	grid = RectilinearGrid(arch; size = (Nxy, Nxy, Nz), extent = (4πh, 2πh, 2h), halo)

Implementation of the scale-dependent Lagrangian-averaged dynamic Smagorinsky #3638

Are you sure you want to change the base?

Implementation of the scale-dependent Lagrangian-averaged dynamic Smagorinsky #3638

Conversation

xkykai commented Jun 24, 2024 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tomchor commented Jun 24, 2024

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

tomchor Jun 24, 2024 • edited Loading

Choose a reason for hiding this comment

tomchor commented Jun 24, 2024

simone-silvestri commented Jun 24, 2024

xkykai commented Jun 24, 2024 • edited Loading

tomchor commented Jun 24, 2024

xkykai commented Jun 24, 2024 • edited Loading

tomchor commented Jun 24, 2024

Choose a reason for hiding this comment

glwagner Jun 25, 2024 • edited Loading

Choose a reason for hiding this comment

glwagner commented Jul 7, 2024

xkykai commented Jun 24, 2024 •

edited

Loading

tomchor Jun 24, 2024 •

edited

Loading

xkykai commented Jun 24, 2024 •

edited

Loading

xkykai commented Jun 24, 2024 •

edited

Loading

glwagner Jun 25, 2024 •

edited

Loading