tomchor / Oceanostics.jl

Diagnostics for Oceananigans
https://tomchor.github.io/Oceanostics.jl/
MIT License
24 stars 8 forks source link

Some questions about calculating the fluctuation variables in TKE budget equation #103

Closed zwei961120 closed 1 year ago

zwei961120 commented 1 year ago

In the TKE budget equation ∂_t E = - ∂_z ⟨w′e′ + w′p′⟩ - ⟨w′u′⟩ ∂_z U - ⟨w′v′⟩ ∂_z V + ⟨w′b′⟩ - ϵ, for example, when calculating the shear production term in x-direction, in my opinion, one need to calculate the horizontally averaged variables (U, V and W) and the fluctuation variables (u'= u-U, v'= v-V, w'= w-W) first, then calculate the shear production term [-(u’u’∂xU + v’u’∂xV + w’u’∂xW)], but in the script TKEBudgetTerms.jl, it seems like the model.velocities (u, v and w) is directly used for calculating the shear production terms:

@inline function shear_production_x_ccc(i, j, k, grid, u, v, w, U, V, W) u_int = ℑxᶜᵃᵃ(i, j, k, grid, u) # F, C, C → C, C, C

∂xU = ∂xᶜᶜᶜ(i, j, k, grid, U) # F, C, C  → C, C, C
uu = ℑxᶜᵃᵃ(i, j, k, grid, ψ², u)
uu∂xU = uu * ∂xU

∂xV = ℑxyᶜᶜᵃ(i, j, k, grid, ∂xᶠᶠᶜ, V) # C, F, C  → F, F, C  → C, C, C
vu = ℑyᵃᶜᵃ(i, j, k, grid, v) * u_int
vu∂xV = vu * ∂xV

∂xW = ℑxzᶜᵃᶜ(i, j, k, grid, ∂xᶠᶜᶠ, W) # C, C, F  → F, C, F  → C, C, C
wu = ℑzᵃᵃᶜ(i, j, k, grid, w) * u_int
wu∂xW = wu * ∂xW

return -(uu∂xU + vu∂xV + wu∂xW)

end

Also, when calculating the pressure redistribution term, the pressure p is directly used without any operations. I don't know if my understanding is correct, and if not, are there any special considerations for doing in that way?

tomchor commented 1 year ago

That's a good point. When I've used this script in the past I've used it in the form:

U = Average(model.velocities.u, dims=(1,2))
V = Average(model.velocities.v, dims=(1,2))
ZSP = ZShearProduction(model, u-U, v-V, w, U, V, 0)

Which passes what you'd call $u'$ directly and would therefore compute the correct quantity. It calls this method:

https://github.com/tomchor/Oceanostics.jl/blob/0c23460c2d3f4b4bd4fd39af808849f34595ba5c/src/TKEBudgetTerms.jl#L271-L275

But you're right that if a user uses ths other method:

https://github.com/tomchor/Oceanostics.jl/blob/0c23460c2d3f4b4bd4fd39af808849f34595ba5c/src/TKEBudgetTerms.jl#L277-L278

then the calculation isn't correct.

I'll open a PR about this!

zwei961120 commented 1 year ago

It seems that the ZpressureRedistribution function has not been adjusted in Oceanostics.jl repo yet, right? So I tested on GPUs in the following way. The model setup is the same as the Langmuir Turbulence example setup in: https://clima.github.io/OceananigansDocumentation/stable/generated/langmuir_turbulence/.

The test codes and outputs in Jupyternotebook are attached here: https://github.com/zwei961120/TKE_test/blob/main/TKE_TEST.ipynb, and everything seems to be going well.

Besides, I am little confused, when calculating the Z-direction pressure redistribution term, why performing velocity and pressure averages in Y and Z directions like W = Field(Average(w, dims=(2, 3))) rather than performing the horizontally averages like W = Field(Average(w, dims=(1, 2)))?

Feel free to tell me your thoughts. Thanks!

tomchor commented 1 year ago

It seems that the ZpressureRedistribution function has not been adjusted in Oceanostics.jl repo yet, right? So I tested on GPUs in the following way. The model setup is the same as the Langmuir Turbulence example setup in: https://clima.github.io/OceananigansDocumentation/stable/generated/langmuir_turbulence/.

The test codes and outputs in Jupyternotebook are attached here: https://github.com/zwei961120/TKE_test/blob/main/TKE_TEST.ipynb, and everything seems to be going well.

Just FYI the urls you passed aren't properly linked. The proper way to do it is:

[some text](https://github.com/zwei961120/TKE_test/blob/main/TKE_TEST.ipynb)

hope this helps in the future :)

Besides, I am little confused, when calculating the Z-direction pressure redistribution term, why performing velocity and pressure averages in Y and Z directions like W = Field(Average(w, dims=(2, 3))) rather than performing the horizontally averages like W = Field(Average(w, dims=(1, 2)))?

Feel free to tell me your thoughts. Thanks!

Let's continue this discussion at https://github.com/tomchor/Oceanostics.jl/pull/104 since we're discussing the code there!

I'll answer you over there