-
Notifications
You must be signed in to change notification settings - Fork 191
Commit
This commit does not belong to any branch on this repository, and may belong to a fork outside of the repository.
TKEDissipationVerticalDiffusivity (aka k-epsilon) turbulence closure …
…(take 2) (#3605) * Generalize definition of Ri * Non-negative Ri * Use quasi-conservative reconstruction of shear production and buoyancy flux * Add convective mixing length for momentum * Update minimum TKE * Major reorganizations to get CATKE time-stepping to work * Working on comparison * Some updates * Heroics to get things to wokrk * Update tke time stepping * Bit of cleanup * Figure out the right temporal stencil * Updated default convective buoyancy flux and delete ccc methods * Add back convective mixing length for ccc * Grr also add back stability function for ccc * Better error message when setting fails for field time series * Better notation for surface buoyancy flux; * Implement split-explicit time-stepping for CATKE * Optimize CATKE substepping a bit more * Performance enhancements * Add third method to CATKE comparison * Minor improvements * Clean up * Update * Move around some inbounds * Delete unused code * Clean up time stepping * Little things * Improve show and consistently use surface-ignoring reconstruction * Update definition of Ri and try a new stability funciton * Delete show * Queue up change to sheared convection * Update scale * Another tweak * Add four more parameters for shear turbulence in unstable stratification (#3600) * Add special negative Ri branch to the stability function * Update parameters * Improve show for CATKE * Skip test for tuples-with-CATKE for now * Try to fix checkpointing / regression test issues * Its a special test so we have to comment it out * Make room for TKEDissipationVerticalDiffusivity * Things are more or less wokring * Trying to get k-epsilon to work * Success? * Wow its working * Updates to permit ensembles of closures * Change TimeInterval to avoid roundoff error issue plus some cleanup * fix typo * Import fix * Potentially address overflow problem * Minor improvements * Bug fix for safety factor * Add a script to inspect stability functions * Make schedule work even with fake models * Fix tests for TimeInterval * Fix JLD2 output writer tests * Fix NetCDF writer test? * Rm reference to CATKEVerticalDiffusivities * fix typo Co-authored-by: Simone Silvestri <[email protected]> * Update src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl Co-authored-by: Navid C. Constantinou <[email protected]> * Update src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl Co-authored-by: Navid C. Constantinou <[email protected]> * Restore convective adjustment notation * Address a few review comments * Get rid of background diffusivities * Update src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_vertical_diffusivity.jl Co-authored-by: Navid C. Constantinou <[email protected]> * Fix docstring * Add reference * Better docstrings and comments --------- Co-authored-by: Navid C. Constantinou <[email protected]> Co-authored-by: Simone Silvestri <[email protected]>
- Loading branch information
3 people
committed
Jun 20, 2024
1 parent
56c6dd9
commit 5661358
Showing
22 changed files
with
2,112 additions
and
81 deletions.
There are no files selected for viewing
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
4 changes: 2 additions & 2 deletions
4
src/Models/HydrostaticFreeSurfaceModels/slice_ensemble_model_mode.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
186 changes: 186 additions & 0 deletions
186
...ce_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| Original file line number | Diff line number | Diff line change |
|---|---|---|
| @@ -0,0 +1,186 @@ | ||
| module TKEBasedVerticalDiffusivities | ||
|
|
||
| using Adapt | ||
| using CUDA | ||
| using KernelAbstractions: @kernel, @index | ||
|
|
||
| using Oceananigans.Architectures | ||
| using Oceananigans.Grids | ||
| using Oceananigans.Utils | ||
| using Oceananigans.Units | ||
| using Oceananigans.Fields | ||
| using Oceananigans.Operators | ||
|
|
||
| using Oceananigans.Utils: prettysummary | ||
| using Oceananigans.Grids: peripheral_node, inactive_node, inactive_cell | ||
| using Oceananigans.Fields: ZeroField | ||
| using Oceananigans.BoundaryConditions: default_prognostic_bc, DefaultBoundaryCondition | ||
| using Oceananigans.BoundaryConditions: BoundaryCondition, FieldBoundaryConditions | ||
| using Oceananigans.BoundaryConditions: DiscreteBoundaryFunction, FluxBoundaryCondition | ||
| using Oceananigans.BuoyancyModels: ∂z_b, top_buoyancy_flux | ||
| using Oceananigans.Grids: inactive_cell | ||
|
|
||
| using Oceananigans.TurbulenceClosures: | ||
| getclosure, | ||
| time_discretization, | ||
| AbstractScalarDiffusivity, | ||
| VerticallyImplicitTimeDiscretization, | ||
| VerticalFormulation | ||
|
|
||
| import Oceananigans.BoundaryConditions: getbc | ||
| import Oceananigans.Utils: with_tracers | ||
| import Oceananigans.TurbulenceClosures: | ||
| validate_closure, | ||
| shear_production, | ||
| buoyancy_flux, | ||
| dissipation, | ||
| add_closure_specific_boundary_conditions, | ||
| compute_diffusivities!, | ||
| DiffusivityFields, | ||
| implicit_linear_coefficient, | ||
| viscosity, | ||
| diffusivity, | ||
| viscosity_location, | ||
| diffusivity_location, | ||
| diffusive_flux_x, | ||
| diffusive_flux_y, | ||
| diffusive_flux_z | ||
|
|
||
| const c = Center() | ||
| const f = Face() | ||
| const VITD = VerticallyImplicitTimeDiscretization | ||
|
|
||
| @inline ϕ²(i, j, k, grid, ϕ, args...) = ϕ(i, j, k, grid, args...)^2 | ||
|
|
||
| @inline function shearᶜᶜᶠ(i, j, k, grid, u, v) | ||
| ∂z_u² = ℑxᶜᵃᵃ(i, j, k, grid, ϕ², ∂zᶠᶜᶠ, u) | ||
| ∂z_v² = ℑyᵃᶜᵃ(i, j, k, grid, ϕ², ∂zᶜᶠᶠ, v) | ||
| S² = ∂z_u² + ∂z_v² | ||
| return S² | ||
| end | ||
|
|
||
| @inline function shearᶜᶜᶜ(i, j, k, grid, u, v) | ||
| ∂z_u² = ℑxᶜᵃᵃ(i, j, k, grid, ℑbzᵃᵃᶜ, ϕ², ∂zᶠᶜᶠ, u) | ||
| ∂z_v² = ℑyᵃᶜᵃ(i, j, k, grid, ℑbzᵃᵃᶜ, ϕ², ∂zᶜᶠᶠ, v) | ||
| S² = ∂z_u² + ∂z_v² | ||
| return S² | ||
| end | ||
|
|
||
| @inline Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy) = | ||
| ℑbzᵃᵃᶜ(i, j, k, grid, Riᶜᶜᶠ, velocities, tracers, buoyancy) | ||
|
|
||
| @inline function Riᶜᶜᶠ(i, j, k, grid, velocities, tracers, buoyancy) | ||
| u = velocities.u | ||
| v = velocities.v | ||
| S² = shearᶜᶜᶠ(i, j, k, grid, u, v) | ||
| N² = ∂z_b(i, j, k, grid, buoyancy, tracers) | ||
| Ri = N² / S² | ||
| return ifelse(N² == 0, zero(grid), Ri) | ||
| end | ||
|
|
||
| # @inline ℑbzᵃᵃᶜ(i, j, k, grid, fᵃᵃᶠ, args...) = ℑzᵃᵃᶜ(i, j, k, grid, fᵃᵃᶠ, args...) | ||
|
|
||
| # A particular kind of reconstruction that ignores peripheral nodes | ||
| @inline function ℑbzᵃᵃᶜ(i, j, k, grid, fᵃᵃᶠ, args...) | ||
| k⁺ = k + 1 | ||
| k⁻ = k | ||
|
|
||
| f⁺ = fᵃᵃᶠ(i, j, k⁺, grid, args...) | ||
| f⁻ = fᵃᵃᶠ(i, j, k⁻, grid, args...) | ||
|
|
||
| p⁺ = peripheral_node(i, j, k⁺, grid, c, c, f) | ||
| p⁻ = peripheral_node(i, j, k⁻, grid, c, c, f) | ||
|
|
||
| f⁺ = ifelse(p⁺, f⁻, f⁺) | ||
| f⁻ = ifelse(p⁻, f⁺, f⁻) | ||
|
|
||
| return (f⁺ + f⁻) / 2 | ||
| end | ||
|
|
||
| # To reconstruct buoyancy flux "conservatively" (ie approximately corresponding to production/destruction | ||
| # of mean potential energy): | ||
| @inline function buoyancy_fluxᶜᶜᶠ(i, j, k, grid, tracers, buoyancy, diffusivities) | ||
| κc = @inbounds diffusivities.κc[i, j, k] | ||
| N² = ∂z_b(i, j, k, grid, buoyancy, tracers) | ||
| return - κc * N² | ||
| end | ||
|
|
||
| @inline explicit_buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities) = | ||
| ℑbzᵃᵃᶜ(i, j, k, grid, buoyancy_fluxᶜᶜᶠ, tracers, buoyancy, diffusivities) | ||
|
|
||
| # Note special attention paid to averaging the vertical grid spacing correctly | ||
| @inline Δz_νₑ_az_bzᶠᶜᶠ(i, j, k, grid, νₑ, a, b) = ℑxᶠᵃᵃ(i, j, k, grid, νₑ) * ∂zᶠᶜᶠ(i, j, k, grid, a) * | ||
| Δzᶠᶜᶠ(i, j, k, grid) * ∂zᶠᶜᶠ(i, j, k, grid, b) | ||
|
|
||
| @inline Δz_νₑ_az_bzᶜᶠᶠ(i, j, k, grid, νₑ, a, b) = ℑyᵃᶠᵃ(i, j, k, grid, νₑ) * ∂zᶜᶠᶠ(i, j, k, grid, a) * | ||
| Δzᶜᶠᶠ(i, j, k, grid) * ∂zᶜᶠᶠ(i, j, k, grid, b) | ||
|
|
||
| @inline function shear_production_xᶠᶜᶜ(i, j, k, grid, νₑ, uⁿ, u⁺) | ||
| Δz_Pxⁿ = ℑbzᵃᵃᶜ(i, j, k, grid, Δz_νₑ_az_bzᶠᶜᶠ, νₑ, uⁿ, u⁺) | ||
| Δz_Px⁺ = ℑbzᵃᵃᶜ(i, j, k, grid, Δz_νₑ_az_bzᶠᶜᶠ, νₑ, u⁺, u⁺) | ||
| return (Δz_Pxⁿ + Δz_Px⁺) / (2 * Δzᶠᶜᶜ(i, j, k, grid)) | ||
| end | ||
|
|
||
| @inline function shear_production_yᶜᶠᶜ(i, j, k, grid, νₑ, vⁿ, v⁺) | ||
| Δz_Pyⁿ = ℑbzᵃᵃᶜ(i, j, k, grid, Δz_νₑ_az_bzᶜᶠᶠ, νₑ, vⁿ, v⁺) | ||
| Δz_Py⁺ = ℑbzᵃᵃᶜ(i, j, k, grid, Δz_νₑ_az_bzᶜᶠᶠ, νₑ, v⁺, v⁺) | ||
| return (Δz_Pyⁿ + Δz_Py⁺) / (2 * Δzᶜᶠᶜ(i, j, k, grid)) | ||
| end | ||
|
|
||
| @inline function shear_production(i, j, k, grid, νₑ, uⁿ, u⁺, vⁿ, v⁺) | ||
| # Reconstruct the shear production term in an "approximately conservative" manner | ||
| # (ie respecting the spatial discretization and using a stencil commensurate with the | ||
| # loss of mean kinetic energy due to shear production --- but _not_ respecting the | ||
| # the temporal discretization. Note that also respecting the temporal discretization, would | ||
| # require storing the velocity field at n and n+1): | ||
|
|
||
| return ℑxᶜᵃᵃ(i, j, k, grid, shear_production_xᶠᶜᶜ, νₑ, uⁿ, u⁺) + | ||
| ℑyᵃᶜᵃ(i, j, k, grid, shear_production_yᶜᶠᶜ, νₑ, vⁿ, v⁺) | ||
| end | ||
|
|
||
| @inline function turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, e) | ||
| eᵢ = @inbounds e[i, j, k] | ||
| eᵐⁱⁿ = closure.minimum_tke | ||
| return sqrt(max(eᵐⁱⁿ, eᵢ)) | ||
| end | ||
|
|
||
| @inline function mask_diffusivity(i, j, k, grid, κ★) | ||
| on_periphery = peripheral_node(i, j, k, grid, c, c, f) | ||
| within_inactive = inactive_node(i, j, k, grid, c, c, f) | ||
| nan = convert(eltype(grid), NaN) | ||
| return ifelse(on_periphery, zero(grid), ifelse(within_inactive, nan, κ★)) | ||
| end | ||
|
|
||
| @inline clip(x) = max(zero(x), x) | ||
|
|
||
| function get_time_step(closure_array::AbstractArray) | ||
| # assume they are all the same | ||
| closure = CUDA.@allowscalar closure_array[1, 1] | ||
| return get_time_step(closure) | ||
| end | ||
|
|
||
| include("tke_top_boundary_condition.jl") | ||
|
|
||
| include("catke_vertical_diffusivity.jl") | ||
| include("catke_mixing_length.jl") | ||
| include("catke_equation.jl") | ||
| include("time_step_catke_equation.jl") | ||
|
|
||
| include("tke_dissipation_vertical_diffusivity.jl") | ||
| include("tke_dissipation_stability_functions.jl") | ||
| include("tke_dissipation_equations.jl") | ||
|
|
||
| for S in (:CATKEMixingLength, | ||
| :CATKEEquation, | ||
| :StratifiedDisplacementScale, | ||
| :ConstantStabilityFunctions, | ||
| :VariableStabilityFunctions) | ||
|
|
||
| @eval @inline convert_eltype(::Type{FT}, s::$S) where FT = | ||
| $S{FT}(; Dict(p => getproperty(s, p) for p in propertynames(s))...) | ||
|
|
||
| @eval @inline convert_eltype(::Type{FT}, s::$S{FT}) where FT = s | ||
| end | ||
|
|
||
| end # module | ||
|
|
Oops, something went wrong.