Use stress divergence and diffusive flux divergence from Turbulence Closures #245
Conversation
…ion, cleans up code and adds comments
…dary transport coefficients with tensor diffusivity abstraction
Merge branch 'master' into integrate-turbulence-closures
| @@ -152,29 +154,29 @@ function calculate_interior_source_terms!(grid::Grid, constants, eos, cfg, u, v, | |||
| @inbounds Gu[i, j, k] = (-u∇u(grid, u, v, w, i, j, k) | |||
| + Gu_cori(grid, v, fCor, i, j, k) | |||
| - δx_c2f(grid, pHY′, i, j, k) / (Δx * ρ₀) | |||
| + 𝜈∇²u(grid, u, 𝜈h, 𝜈v, i, j, k) | |||
| + ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, closure, eos, grav, u, v, w, T, S) | |||
There was a problem hiding this comment.
Dumb question but is Σ the strain-rate tensor?
There was a problem hiding this comment.
quite! I will add some comments to turbulence_closures.jl
|
Cool! Looks pretty neat considering how much tedious math is buried underneath. Thanks for cleaning up the
+1. Do you want to close #120 once this is merged? |
| κ_bottom = κ(i, j, grid.Nz, grid, closure, eos, g, u, v, w, T, S) | ||
|
|
||
| apply_z_top_bc!(top_bc, i, j, grid, ϕ, Gϕ, κ_top, t, iteration, u, v, w, T, S) | ||
| apply_z_bottom_bc!(bottom_bc, i, j, grid, ϕ, Gϕ, κ_bottom, t, iteration, u, v, w, T, S) |
There was a problem hiding this comment.
It looks like if you want to impose a z boundary condition that does not depend on κ, you still have to calculate κ using the full closure which can be expensive if using an LES closure. Not sure how to get around this, probably some clever multiple dispatch?
This is probably fine for now as constant Smagorinsky isn't integrated yet, and the performance hit probably isn't big enough to worry about right now.
There was a problem hiding this comment.
Hmm yes that is possible. It doesn't have to clever :-P --- we will just pass the function to calculate the transport coefficient into apply_z_top_bc!, which dispatches on the type of its first argument (the boundary condition now). That's possibly a better design anyways.
But actually, it might get compiled away if not used anyways --- the type of both boundary conditions is known to the function apply_z_bcs!.
I didn't put a significant amount of effort into this problem because I'm expecting a lot to simplify once we have halo regions.
|
Shouldn't the default be anisotropic diffusion with the anisotrop scaled by
delx, delz. This will be by far the most common model configuration.
…On Mon, May 27, 2019, 4:04 PM Ali Ramadhan ***@***.***> wrote:
***@***.**** commented on this pull request.
------------------------------
In src/time_steppers.jl
<#245 (comment)>
:
> @loop for j in (1:grid.Ny; (blockIdx().y - 1) * blockDim().y + threadIdx().y)
@loop for i in (1:grid.Nx; (blockIdx().x - 1) * blockDim().x + threadIdx().x)
- apply_z_top_bc!(top_bc, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S)
- apply_z_bottom_bc!(bottom_bc, i, j, grid, ϕ, Gϕ, κ, t, iteration, u, v, w, T, S)
+
+ κ_top = κ(i, j, 1, grid, closure, eos, g, u, v, w, T, S)
+ κ_bottom = κ(i, j, grid.Nz, grid, closure, eos, g, u, v, w, T, S)
+
+ apply_z_top_bc!(top_bc, i, j, grid, ϕ, Gϕ, κ_top, t, iteration, u, v, w, T, S)
+ apply_z_bottom_bc!(bottom_bc, i, j, grid, ϕ, Gϕ, κ_bottom, t, iteration, u, v, w, T, S)
It looks like if you want to impose a z boundary condition that does not
depend on κ, you still have to calculate κ using the full closure which
can be expensive if using an LES closure. Not sure how to get around this,
probably some clever multiple dispatch?
This is probably fine for now as constant Smagorinsky isn't integrated
yet, and the performance hit probably isn't big enough to worry about right
now.
—
You are receiving this because you are subscribed to this thread.
Reply to this email directly, view it on GitHub
<#245>,
or mute the thread
<https://github.com/notifications/unsubscribe-auth/AKXUEQXS3APE4PYPCYWURM3PXQ5FPANCNFSM4HP6AWMA>
.
|
|
@johncmarshall54 we can make that the default. I wouldn't expect any constant diffusivity closure to be common for LES application considering the ready availability of better subgrid-scale models. The purpose of making ConstantIsotropicDiffusivity the default has a purpose: I want users to choose certain closures (and realize that any specification which is not molecular diffusivity is, in fact, a sub grid model for turbulence), which will motivate them to think about the consequences of their choice. |
Let's wait --- I'm working on test cases for |
Codecov Report
@@ Coverage Diff @@
## master #245 +/- ##
==========================================
- Coverage 75.37% 70.59% -4.78%
==========================================
Files 23 23
Lines 877 857 -20
==========================================
- Hits 661 605 -56
- Misses 216 252 +36
Continue to review full report at Codecov.
|
Codecov Report
@@ Coverage Diff @@
## master #245 +/- ##
==========================================
- Coverage 74.71% 70.59% -4.12%
==========================================
Files 23 23
Lines 866 857 -9
==========================================
- Hits 647 605 -42
- Misses 219 252 +33
Continue to review full report at Codecov.
|
|
any suggestions for tests I might add to get the coverage up? |
|
Another q: should I shower the code in a hailstorm of |
|
@ali-ramadhan do you mind running some benchmarks to test for performance regression under this PR? |
|
Lastly, I am thinking that all the doc strings in |
Hmmm, I think for now it's sufficient that the regression tests pass as this PR should preserve existing functionality. If you're going to implement more rigorous/high-level LES tests in the future then the coverage will go up. And it'll probably become clearer which unit tests are needed.
I kind of agree, but with the docstrings we can integrate them into the documentation, and if the docstrings have LaTeX then we can view the operators alongside the math in the docs. I guess it's readable documentation vs. more readable code? Good practice says we should probably keep them, but maybe we can separate them somehow? I guess right now we only read the code but maybe in the future we'll mainly be reading the docs and not the code. |
Looks really bad, still waiting on GPU results but CPU is ~100x slower! We can't merge this as is. Something must be wrong somewhere. Memory allocations went from I benchmarked master yesterday so I know it's still performing well (maybe 2-3% slower because of the extra stuff and fixes we recently introduced). |
|
I'm getting an error when trying to compile constant Smagorinsky: ERROR: LoadError: InvalidIRError: compiling #12(RegularCartesianGrid{Float64,StepRangeLen{Float64,Base.TwicePrecision{Flo
at64},Base.TwicePrecision{Float64}}}, PlanetaryConstants{Float64}, LinearEquationOfState{Float64}, ConstantSmagorinsky{Fl
oat64}, CUDAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global
}, CUDAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CU
DAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnat
ive.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnative.C
uDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnative.CuDeviceArray{Float64,3,CUDAnative.AS.Global}, CUDAnative.CuDevi
ceArray{Float64,3,CUDAnative.AS.Global}, Forcing{typeof(Oceananigans.zero_func),typeof(Oceananigans.zero_func),typeof(Oce
ananigans.zero_func),typeof(Oceananigans.zero_func),typeof(Oceananigans.zero_func)}) resulted in invalid LLVM IR
Reason: unsupported call to the Julia runtime (call to jl_f__apply)
Stacktrace:
[1] overdub at /data5/glwagner/.julia/packages/Cassette/xggAf/src/context.jl:260
[2] ν_ccc at /data5/glwagner/Projects/Oceananigans.jl/src/closures/constant_smagorinsky.jl:109
[3] ν_Σᵢⱼ at /data5/glwagner/Projects/Oceananigans.jl/src/closures/closure_operators.jl:405
[4] ∂x_faa at /data5/glwagner/Projects/Oceananigans.jl/src/closures/closure_operators.jl:64
[5] ∂x_2ν_Σ₁₁ at /data5/glwagner/Projects/Oceananigans.jl/src/closures/closure_operators.jl:409
[6] ∂ⱼ_2ν_Σ₁ⱼ at /data5/glwagner/Projects/Oceananigans.jl/src/closures/closure_operators.jl:432
[7] calculate_interior_source_terms! at /data5/glwagner/Projects/Oceananigans.jl/src/time_steppers.jl:152
[8] #12 at /data5/glwagner/.julia/packages/GPUifyLoops/hBRid/src/context.jl:136
Reason: unsupported dynamic function invocation (call to Cassette.overdub)I think this is specific to the package upgrades. I also got this error when running the |
|
@ali-ramadhan are those regressions on the CPU or GPU? On the CPU, we might need Edit: just realize you said its CPU. Ok, well we might need to shower the code in And of course allocation could also be a separate issue; we will see. |
|
Which version of GPUifyLoops are you using? Maybe post |
Just CPU right now. They're taking too long so I'm just gonna kill it and run the GPU tests by themselves. I'd be pretty surprised if |
|
@glwagner You might be right about the Before: Now: |
|
which closure is this using? |
The model is constructed using model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), arch=arch, float_type=float_type)so looks like it should be using the constructor defaults ν = 1.05e-6, νh=ν, νv=ν,
κ = 1.43e-7, κh=κ, κv=κ,
closure = ConstantAnisotropicDiffusivity(νh=νh, νv=νv, κh=κh, κv=κv) |
|
Hmm... with a vanilla closure, the change is completely encapsulated in the addition of two layers of abstraction (we are just calling a simple diffusion operator). So, let's figure out how to make the abstractions fast. I think the slow down for vanilla closures should be The 'abstraction slowdown' causes much larger problems with the complicated closures, so we need to solve that problem anyways. Edit: I see your post, so what I said above holds. |
|
Hmmm, since this is a CPU issue might be a good time to just profile the code. |
…nstants and eos to type of model
|
@ali-ramadhan do you mind running the benchmarks again? There are still some performance issues, especially with |
|
Some comments:
|
|
Huh --- on my GPU I was getting 2x speed up for The slow down has to do with the abstractions I have introduced. 30% is a huge slow down for one function, indicative of a major problem --- probably a type inference issue? I think that once this problem is solved the code may become faster because of the disambiguation this PR lends to the innermost kernels. This problem becomes catastrophic for the closures, which make heavy use of the abstraction. So solving this problem is imperative. We can restore the performance of the default closure by simply pasting the old operators into Unfortunately, I'm in I'm wondering whether these problems will vanish once we eliminate branches from the inmost functions... |
|
There are still some type problems --- the Adams-Bashforth parameter is explicitly |
Use stress divergence and diffusive flux divergence from Turbulence Closures Former-commit-id: 0b36067
This PR integrates the
TurbulenceClosuresmodule into time stepping and boundary conditions.The need to abstractly deal with anisotropic transport coefficients for arbitrary boundaries introduces considerable complexity. This problem is solved by exporting
NamedTuplesthat collect functions to calculate the diagonal components of viscosity and diffusivity at the necessary locations. The consequence of this implementation is viewed in theapply_bcs!function.Only the
ConstantAnisotropicDiffusivityclosure (corresponding to the former default) is currently tested.In the future, we should probably make
ConstantIsotropicDiffusivitythe default, and remove the option to set the horizontal and vertical diffusion tensor components in theModelconstructor.