From 5661358a5aa54909e258edddd018a25d04a4ca37 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Thu, 20 Jun 2024 16:08:59 -0600 Subject: [PATCH] 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 * Update src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl Co-authored-by: Navid C. Constantinou * Update src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl Co-authored-by: Navid C. Constantinou * 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 * Fix docstring * Add reference * Better docstrings and comments --------- Co-authored-by: Navid C. Constantinou Co-authored-by: Simone Silvestri --- ...ute_hydrostatic_free_surface_tendencies.jl | 2 +- .../hydrostatic_free_surface_ab2_step.jl | 4 + .../single_column_model_mode.jl | 2 +- .../slice_ensemble_model_mode.jl | 4 +- ...ore_hydrostatic_free_surface_tendencies.jl | 4 + .../nonhydrostatic_model.jl | 2 +- src/TurbulenceClosures/TurbulenceClosures.jl | 4 +- .../TKEBasedVerticalDiffusivities.jl | 186 +++++++++ .../catke_equation.jl | 167 +++++++++ .../catke_mixing_length.jl} | 26 +- .../catke_vertical_diffusivity.jl | 337 +++++++++++++++++ .../time_step_catke_equation.jl} | 22 +- .../tke_dissipation_equations.jl | 299 +++++++++++++++ .../tke_dissipation_stability_functions.jl | 292 +++++++++++++++ .../tke_dissipation_vertical_diffusivity.jl | 352 ++++++++++++++++++ .../tke_top_boundary_condition.jl | 91 +++++ .../ri_based_vertical_diffusivity.jl | 34 +- test/test_time_stepping.jl | 2 +- .../column_windy_convection.jl | 59 +-- .../compare_catke_k_epsilon.jl | 114 ++++++ .../k_epsilon_boundary_layer.jl | 123 ++++++ .../look_at_stability_functions.jl | 67 ++++ 22 files changed, 2112 insertions(+), 81 deletions(-) create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_equation.jl rename src/TurbulenceClosures/turbulence_closure_implementations/{CATKEVerticalDiffusivities/mixing_length.jl => TKEBasedVerticalDiffusivities/catke_mixing_length.jl} (94%) create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl rename src/TurbulenceClosures/turbulence_closure_implementations/{CATKEVerticalDiffusivities/time_step_turbulent_kinetic_energy.jl => TKEBasedVerticalDiffusivities/time_step_catke_equation.jl} (90%) create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_vertical_diffusivity.jl create mode 100644 src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_top_boundary_condition.jl create mode 100644 validation/vertical_mixing_closures/compare_catke_k_epsilon.jl create mode 100644 validation/vertical_mixing_closures/k_epsilon_boundary_layer.jl create mode 100644 validation/vertical_mixing_closures/look_at_stability_functions.jl diff --git a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl index 0c9f6b17b5..e47ba53cba 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/compute_hydrostatic_free_surface_tendencies.jl @@ -8,7 +8,7 @@ using Oceananigans.Utils: work_layout, KernelParameters using Oceananigans.Grids: halo_size using Oceananigans.Fields: immersed_boundary_condition using Oceananigans.Biogeochemistry: update_tendencies! -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: FlavorOfCATKE +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: FlavorOfCATKE, FlavorOfTD using Oceananigans.ImmersedBoundaries: active_interior_map, ActiveCellsIBG, InteriorMap, active_linear_index_to_tuple diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 65ee037f2e..1f220ea216 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl @@ -78,6 +78,10 @@ function ab2_step_tracers!(tracers, model, Δt, χ) # TODO: do better than this silly criteria, also need to check closure tuples if closure isa FlavorOfCATKE && tracer_name == :e @debug "Skipping AB2 step for e" + elseif closure isa FlavorOfTD && tracer_name == :ϵ + @debug "Skipping AB2 step for ϵ" + elseif closure isa FlavorOfTD && tracer_name == :e + @debug "Skipping AB2 step for e" else Gⁿ = model.timestepper.Gⁿ[tracer_name] G⁻ = model.timestepper.G⁻[tracer_name] diff --git a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl index ea5f343720..8429887ea7 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl @@ -6,7 +6,7 @@ using Oceananigans.Grids: Flat, Bounded using Oceananigans.Fields: ZeroField using Oceananigans.Coriolis: AbstractRotation using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVDArray +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVDArray import Oceananigans.Grids: validate_size, validate_halo import Oceananigans.Models: validate_tracer_advection diff --git a/src/Models/HydrostaticFreeSurfaceModels/slice_ensemble_model_mode.jl b/src/Models/HydrostaticFreeSurfaceModels/slice_ensemble_model_mode.jl index d9fda68ef1..e2963471b3 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/slice_ensemble_model_mode.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/slice_ensemble_model_mode.jl @@ -1,11 +1,11 @@ using Oceananigans.Grids: Flat, Bounded, y_domain using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: _top_tke_flux, CATKEVDArray +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: _top_tke_flux, CATKEVDArray import Oceananigans.Grids: validate_size, validate_halo, XYRegularRG import Oceananigans.TurbulenceClosures: time_discretization, compute_diffusivities!, with_tracers import Oceananigans.TurbulenceClosures: ∂ⱼ_τ₁ⱼ, ∂ⱼ_τ₂ⱼ, ∂ⱼ_τ₃ⱼ, ∇_dot_qᶜ -import Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: top_tke_flux +import Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: top_tke_flux import Oceananigans.Coriolis: x_f_cross_U, y_f_cross_U, z_f_cross_U ##### diff --git a/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl index 85f83bf932..3b0ceadd13 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl @@ -37,6 +37,10 @@ function store_tendencies!(model::HydrostaticFreeSurfaceModel) if closure isa FlavorOfCATKE && field_name == :e @debug "Skipping store tendencies for e" + elseif closure isa FlavorOfTD && field_name == :ϵ + @debug "Skipping store tendencies for ϵ" + elseif closure isa FlavorOfTD && field_name == :e + @debug "Skipping store tendencies for e" else launch!(model.architecture, model.grid, :xyz, store_field_tendencies!, diff --git a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl index b5a22163c3..44ecdce89b 100644 --- a/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl +++ b/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl @@ -15,7 +15,7 @@ using Oceananigans.Models: AbstractModel, NaNChecker, extract_boundary_condition using Oceananigans.Solvers: FFTBasedPoissonSolver using Oceananigans.TimeSteppers: Clock, TimeStepper, update_state!, AbstractLagrangianParticles using Oceananigans.TurbulenceClosures: validate_closure, with_tracers, DiffusivityFields, time_discretization, implicit_diffusion_solver -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: FlavorOfCATKE +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: FlavorOfCATKE using Oceananigans.Utils: tupleit using Oceananigans.Grids: topology diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl index f401a8eab3..f810ce5e2f 100644 --- a/src/TurbulenceClosures/TurbulenceClosures.jl +++ b/src/TurbulenceClosures/TurbulenceClosures.jl @@ -165,7 +165,7 @@ include("turbulence_closure_implementations/scalar_biharmonic_diffusivity.jl") include("turbulence_closure_implementations/smagorinsky_lilly.jl") include("turbulence_closure_implementations/anisotropic_minimum_dissipation.jl") include("turbulence_closure_implementations/convective_adjustment_vertical_diffusivity.jl") -include("turbulence_closure_implementations/CATKEVerticalDiffusivities/CATKEVerticalDiffusivities.jl") +include("turbulence_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl") include("turbulence_closure_implementations/ri_based_vertical_diffusivity.jl") # Special non-abstracted diffusivities: @@ -173,7 +173,7 @@ include("turbulence_closure_implementations/ri_based_vertical_diffusivity.jl") include("turbulence_closure_implementations/isopycnal_skew_symmetric_diffusivity.jl") include("turbulence_closure_implementations/leith_enstrophy_diffusivity.jl") -using .CATKEVerticalDiffusivities: CATKEVerticalDiffusivity +using .TKEBasedVerticalDiffusivities: CATKEVerticalDiffusivity, TKEDissipationVerticalDiffusivity # Miscellaneous utilities include("diffusivity_fields.jl") diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl new file mode 100644 index 0000000000..77d0d465fb --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/TKEBasedVerticalDiffusivities.jl @@ -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 + diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_equation.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_equation.jl new file mode 100644 index 0000000000..80275730f3 --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_equation.jl @@ -0,0 +1,167 @@ +""" + struct CATKEEquation{FT} + +Parameters for the evolution of oceanic turbulent kinetic energy at the O(1 m) scales associated with +isotropic turbulence and diapycnal mixing. +""" +Base.@kwdef struct CATKEEquation{FT} + CʰⁱD :: FT = 0.357 # Dissipation length scale shear coefficient for high Ri + CˡᵒD :: FT = 0.926 # Dissipation length scale shear coefficient for low Ri + CᵘⁿD :: FT = 1.437 # Dissipation length scale shear coefficient for high Ri + CᶜD :: FT = 2.556 # Dissipation length scale convecting layer coefficient + CᵉD :: FT = 0.0 # Dissipation length scale penetration layer coefficient + Cᵂu★ :: FT = 0.405 # Surface shear-driven TKE flux coefficient + CᵂwΔ :: FT = 0.873 # Surface convective TKE flux coefficient + Cᵂϵ :: FT = 1.0 # Dissipative near-bottom TKE flux coefficient +end + +##### +##### Terms in the turbulent kinetic energy equation, all at cell centers +##### + +#= +@inline buoyancy_flux(i, j, k, grid, closure::FlavorOfCATKE, velocities, tracers, buoyancy, diffusivities) = + explicit_buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities) + +@inline function buoyancy_flux(i, j, k, grid, closure::FlavorOfCATKE{<:VITD}, velocities, tracers, buoyancy, diffusivities) + wb = explicit_buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities) + + # "Patankar trick" for buoyancy production (cf Patankar 1980 or Burchard et al. 2003) + # If buoyancy flux is a _sink_ of TKE, we treat it implicitly, and return zero here for + # the explicit buoyancy flux. + return max(zero(grid), wb) +end +=# + +@inline dissipation(i, j, k, grid, closure::FlavorOfCATKE{<:VITD}, args...) = zero(grid) + +@inline function dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure::FlavorOfCATKE, velocities, tracers, + buoyancy, surface_buoyancy_flux) + + # Convective dissipation length + Cᶜ = closure.turbulent_kinetic_energy_equation.CᶜD + Cᵉ = closure.turbulent_kinetic_energy_equation.CᵉD + Cˢᵖ = closure.mixing_length.Cˢᵖ + Jᵇ = surface_buoyancy_flux + ℓʰ = convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, velocities, tracers, buoyancy, Jᵇ) + + # "Stable" dissipation length + Cˡᵒ = closure.turbulent_kinetic_energy_equation.CˡᵒD + Cʰⁱ = closure.turbulent_kinetic_energy_equation.CʰⁱD + Cᵘⁿ = closure.turbulent_kinetic_energy_equation.CᵘⁿD + σᴰ = stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cᵘⁿ, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) + ℓ★ = stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy) + ℓ★ = ℓ★ / σᴰ + + # Dissipation length + ℓʰ = ifelse(isnan(ℓʰ), zero(grid), ℓʰ) + ℓ★ = ifelse(isnan(ℓ★), zero(grid), ℓ★) + ℓᴰ = max(ℓ★, ℓʰ) + + H = total_depthᶜᶜᵃ(i, j, grid) + return min(H, ℓᴰ) +end + +@inline function dissipation_rate(i, j, k, grid, closure::FlavorOfCATKE, + velocities, tracers, buoyancy, diffusivities) + + ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities.Jᵇ) + e = tracers.e + FT = eltype(grid) + eᵢ = @inbounds e[i, j, k] + + # Note: + # Because ∂t e + ⋯ = ⋯ + L e = ⋯ - ϵ, + # + # then L e = - ϵ + # = - Cᴰ e³² / ℓ + # + # and thus L = - Cᴰ √e / ℓ . + + ω_numerical = 1 / closure.negative_tke_damping_time_scale + ω_physical = sqrt(abs(eᵢ)) / ℓᴰ + + return ifelse(eᵢ < 0, ω_numerical, ω_physical) +end + +# Fallbacks for explicit time discretization +@inline function dissipation(i, j, k, grid, closure::FlavorOfCATKE, velocities, tracers, args...) + eᵢ = @inbounds tracers.e[i, j, k] + ω = dissipation_rate(i, j, k, grid, closure, velocities, tracers, args...) + return ω * eᵢ +end + +##### +##### TKE top boundary condition +##### + +@inline function top_tke_flux(i, j, grid, clock, fields, parameters, closure::FlavorOfCATKE, buoyancy) + closure = getclosure(i, j, closure) + + top_tracer_bcs = parameters.top_tracer_boundary_conditions + top_velocity_bcs = parameters.top_velocity_boundary_conditions + tke_parameters = closure.turbulent_kinetic_energy_equation + + return _top_tke_flux(i, j, grid, clock, fields, tke_parameters, closure, + buoyancy, top_tracer_bcs, top_velocity_bcs) +end + +@inline function _top_tke_flux(i, j, grid, clock, fields, + tke::CATKEEquation, closure::CATKEVD, + buoyancy, top_tracer_bcs, top_velocity_bcs) + + wΔ³ = top_convective_turbulent_velocity_cubed(i, j, grid, clock, fields, buoyancy, top_tracer_bcs) + u★ = friction_velocity(i, j, grid, clock, fields, top_velocity_bcs) + + Cᵂu★ = tke.Cᵂu★ + CᵂwΔ = tke.CᵂwΔ + + return - Cᵂu★ * u★^3 - CᵂwΔ * wΔ³ +end + +##### +##### Utilities for model constructors +##### + +""" Add TKE boundary conditions specific to `CATKEVerticalDiffusivity`. """ +function add_closure_specific_boundary_conditions(closure::FlavorOfCATKE, + user_bcs, + grid, + tracer_names, + buoyancy) + + top_tracer_bcs = top_tracer_boundary_conditions(grid, tracer_names, user_bcs) + top_velocity_bcs = top_velocity_boundary_conditions(grid, user_bcs) + parameters = TKETopBoundaryConditionParameters(top_tracer_bcs, top_velocity_bcs) + top_tke_bc = FluxBoundaryCondition(top_tke_flux, discrete_form=true, parameters=parameters) + + if :e ∈ keys(user_bcs) + e_bcs = user_bcs[:e] + + tke_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), + top = top_tke_bc, + bottom = e_bcs.bottom, + north = e_bcs.north, + south = e_bcs.south, + east = e_bcs.east, + west = e_bcs.west) + else + tke_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), top=top_tke_bc) + end + + new_boundary_conditions = merge(user_bcs, (; e = tke_bcs)) + + return new_boundary_conditions +end + +Base.summary(::CATKEEquation) = "TKEBasedVerticalDiffusivities.CATKEEquation" +Base.show(io::IO, tke::CATKEEquation) = + print(io, "TKEBasedVerticalDiffusivities.CATKEEquation parameters:", '\n', + "├── CʰⁱD: ", tke.CʰⁱD, '\n', + "├── CˡᵒD: ", tke.CˡᵒD, '\n', + "├── CᵘⁿD: ", tke.CᵘⁿD, '\n', + "├── CᶜD: ", tke.CᶜD, '\n', + "├── CᵉD: ", tke.CᵉD, '\n', + "├── Cᵂu★: ", tke.Cᵂu★, '\n', + "└── CᵂwΔ: ", tke.CᵂwΔ) + diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/mixing_length.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_mixing_length.jl similarity index 94% rename from src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/mixing_length.jl rename to src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_mixing_length.jl index 792166a6ec..856b2e681d 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/mixing_length.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_mixing_length.jl @@ -8,11 +8,11 @@ using ..TurbulenceClosures: total_depthᶜᶜᵃ """ - struct MixingLength{FT} + struct CATKEMixingLength{FT} Contains mixing length parameters for CATKE vertical diffusivity. """ -Base.@kwdef struct MixingLength{FT} +Base.@kwdef struct CATKEMixingLength{FT} Cˢ :: FT = 1.131 # Surface distance coefficient for shear length scale Cᵇ :: FT = Inf # Bottom distance coefficient for shear length scale Cˢᵖ :: FT = 0.505 # Sheared convective plume coefficient @@ -39,22 +39,6 @@ end ##### Mixing length ##### -@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 function stratification_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, e, tracers, buoyancy) FT = eltype(grid) N² = ∂z_b(i, j, k, grid, buoyancy, tracers) @@ -292,10 +276,10 @@ end return min(H, ℓe) end -Base.summary(::MixingLength) = "CATKEVerticalDiffusivities.MixingLength" +Base.summary(::CATKEMixingLength) = "TKEBasedVerticalDiffusivities.CATKEMixingLength" -Base.show(io::IO, ml::MixingLength) = - print(io, "CATKEVerticalDiffusivities.MixingLength parameters:", '\n', +Base.show(io::IO, ml::CATKEMixingLength) = + print(io, "TKEBasedVerticalDiffusivities.CATKEMixingLength parameters:", '\n', " ├── Cˢ: ", ml.Cˢ, '\n', " ├── Cᵇ: ", ml.Cᵇ, '\n', " ├── Cʰⁱu: ", ml.Cʰⁱu, '\n', diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl new file mode 100644 index 0000000000..2adf895c82 --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/catke_vertical_diffusivity.jl @@ -0,0 +1,337 @@ + +struct CATKEVerticalDiffusivity{TD, CL, FT, DT, TKE} <: AbstractScalarDiffusivity{TD, VerticalFormulation, 2} + mixing_length :: CL + turbulent_kinetic_energy_equation :: TKE + maximum_tracer_diffusivity :: FT + maximum_tke_diffusivity :: FT + maximum_viscosity :: FT + minimum_tke :: FT + minimum_convective_buoyancy_flux :: FT + negative_tke_damping_time_scale :: FT + tke_time_step :: DT +end + +function CATKEVerticalDiffusivity{TD}(mixing_length::CL, + turbulent_kinetic_energy_equation::TKE, + maximum_tracer_diffusivity::FT, + maximum_tke_diffusivity::FT, + maximum_viscosity::FT, + minimum_tke::FT, + minimum_convective_buoyancy_flux::FT, + negative_tke_damping_time_scale::FT, + tke_time_step::DT) where {TD, CL, FT, DT, TKE} + + return CATKEVerticalDiffusivity{TD, CL, FT, DT, TKE}(mixing_length, + turbulent_kinetic_energy_equation, + maximum_tracer_diffusivity, + maximum_tke_diffusivity, + maximum_viscosity, + minimum_tke, + minimum_convective_buoyancy_flux, + negative_tke_damping_time_scale, + tke_time_step) +end + +CATKEVerticalDiffusivity(FT::DataType; kw...) = + CATKEVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), FT; kw...) + +const CATKEVD{TD} = CATKEVerticalDiffusivity{TD} where TD +const CATKEVDArray{TD} = AbstractArray{<:CATKEVD{TD}} where TD +const FlavorOfCATKE{TD} = Union{CATKEVD{TD}, CATKEVDArray{TD}} where TD + +""" + CATKEVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(), + FT = Float64;] + mixing_length = CATKEMixingLength(), + turbulent_kinetic_energy_equation = CATKEEquation(), + maximum_tracer_diffusivity = Inf, + maximum_tke_diffusivity = Inf, + maximum_viscosity = Inf, + minimum_tke = 1e-9, + minimum_convective_buoyancy_flux = 1e-11, + negative_tke_damping_time_scale = 1minute, + tke_time_step = nothing) + +Return the `CATKEVerticalDiffusivity` turbulence closure for vertical mixing by +small-scale ocean turbulence based on the prognostic evolution of subgrid +Turbulent Kinetic Energy (TKE). + +!!! note "CATKE vertical diffusivity" + `CATKEVerticalDiffusivity` is new turbulence closure diffusivity. The default + values for its free parameters are obtained from calibration against large eddy + simulations. For more details please refer to [Wagner23catke](@cite). + + Use with caution and report any issues with the physics at https://github.com/CliMA/Oceananigans.jl/issues. + +Arguments +========= + +- `time_discretization`: Either `ExplicitTimeDiscretization()` or `VerticallyImplicitTimeDiscretization()`; + default `VerticallyImplicitTimeDiscretization()`. + +- `FT`: Float type; default `Float64`. + + +Keyword arguments +================= + +- `maximum_diffusivity`: Maximum value for tracer, momentum, and TKE diffusivities. + Used to clip the diffusivity when/if CATKE predicts + diffusivities that are too large. + Default: `Inf`. + +- `minimum_tke`: Minimum value for the turbulent kinetic energy. + Can be used to model the presence "background" TKE + levels due to, for example, mixing by breaking internal waves. + Default: 0. + +- `negative_tke_damping_time_scale`: Damping time-scale for spurious negative values of TKE, + typically generated by oscillatory errors associated + with TKE advection. + Default: 1 minute. + +Note that for numerical stability, it is recommended to either have a relative short +`negative_tke_damping_time_scale` or a reasonable +`minimum_tke`, or both. +""" +function CATKEVerticalDiffusivity(time_discretization::TD = VerticallyImplicitTimeDiscretization(), + FT = Float64; + mixing_length = CATKEMixingLength(), + turbulent_kinetic_energy_equation = CATKEEquation(), + maximum_tracer_diffusivity = Inf, + maximum_tke_diffusivity = Inf, + maximum_viscosity = Inf, + minimum_tke = 1e-9, + minimum_convective_buoyancy_flux = 1e-11, + negative_tke_damping_time_scale = 1minute, + tke_time_step = nothing) where TD + + mixing_length = convert_eltype(FT, mixing_length) + turbulent_kinetic_energy_equation = convert_eltype(FT, turbulent_kinetic_energy_equation) + + return CATKEVerticalDiffusivity{TD}(mixing_length, + turbulent_kinetic_energy_equation, + convert(FT, maximum_tracer_diffusivity), + convert(FT, maximum_tke_diffusivity), + convert(FT, maximum_viscosity), + convert(FT, minimum_tke), + convert(FT, minimum_convective_buoyancy_flux), + convert(FT, negative_tke_damping_time_scale), + tke_time_step) +end + +function with_tracers(tracer_names, closure::FlavorOfCATKE) + :e ∈ tracer_names || + throw(ArgumentError("Tracers must contain :e to represent turbulent kinetic energy " * + "for `CATKEVerticalDiffusivity`.")) + + return closure +end + +# For tuples of closures, we need to know _which_ closure is CATKE. +# Here we take a "simple" approach that sorts the tuple so CATKE is first. +# This is not sustainable though if multiple closures require this. +# The two other possibilities are: +# 1. Recursion to find which closure is CATKE in a compiler-inferrable way +# 2. Store the "CATKE index" inside CATKE via validate_closure. +validate_closure(closure_tuple::Tuple) = Tuple(sort(collect(closure_tuple), lt=catke_first)) + +catke_first(closure1, catke::FlavorOfCATKE) = false +catke_first(catke::FlavorOfCATKE, closure2) = true +catke_first(closure1, closure2) = false +catke_first(catke1::FlavorOfCATKE, catke2::FlavorOfCATKE) = error("Can't have two CATKEs in one closure tuple.") + +##### +##### Diffusivities and diffusivity fields utilities +##### + +function DiffusivityFields(grid, tracer_names, bcs, closure::FlavorOfCATKE) + + default_diffusivity_bcs = (κu = FieldBoundaryConditions(grid, (Center, Center, Face)), + κc = FieldBoundaryConditions(grid, (Center, Center, Face)), + κe = FieldBoundaryConditions(grid, (Center, Center, Face))) + + bcs = merge(default_diffusivity_bcs, bcs) + + κu = ZFaceField(grid, boundary_conditions=bcs.κu) + κc = ZFaceField(grid, boundary_conditions=bcs.κc) + κe = ZFaceField(grid, boundary_conditions=bcs.κe) + Le = CenterField(grid) + Jᵇ = Field{Center, Center, Nothing}(grid) + previous_compute_time = Ref(zero(grid)) + + # Note: we may be able to avoid using the "previous velocities" in favor of a "fully implicit" + # discretization of shear production + u⁻ = XFaceField(grid) + v⁻ = YFaceField(grid) + previous_velocities = (; u=u⁻, v=v⁻) + + # Secret tuple for getting tracer diffusivities with tuple[tracer_index] + _tupled_tracer_diffusivities = NamedTuple(name => name === :e ? κe : κc for name in tracer_names) + _tupled_implicit_linear_coefficients = NamedTuple(name => name === :e ? Le : ZeroField() for name in tracer_names) + + return (; κu, κc, κe, Le, Jᵇ, + previous_compute_time, previous_velocities, + _tupled_tracer_diffusivities, _tupled_implicit_linear_coefficients) +end + +@inline viscosity_location(::FlavorOfCATKE) = (c, c, f) +@inline diffusivity_location(::FlavorOfCATKE) = (c, c, f) + +function compute_diffusivities!(diffusivities, closure::FlavorOfCATKE, model; parameters = :xyz) + + arch = model.architecture + grid = model.grid + velocities = model.velocities + tracers = model.tracers + buoyancy = model.buoyancy + clock = model.clock + top_tracer_bcs = NamedTuple(c => tracers[c].boundary_conditions.top for c in propertynames(tracers)) + Δt = model.clock.time - diffusivities.previous_compute_time[] + diffusivities.previous_compute_time[] = model.clock.time + + if isfinite(model.clock.last_Δt) # Check that we have taken a valid time-step first. + # Compute e at the current time: + # * update tendency Gⁿ using current and previous velocity field + # * use tridiagonal solve to take an implicit step + time_step_catke_equation!(model) + end + + # Update "previous velocities" + u, v, w = model.velocities + u⁻, v⁻ = diffusivities.previous_velocities + parent(u⁻) .= parent(u) + parent(v⁻) .= parent(v) + + launch!(arch, grid, :xy, + compute_average_surface_buoyancy_flux!, + diffusivities.Jᵇ, grid, closure, velocities, tracers, buoyancy, top_tracer_bcs, clock, Δt) + + launch!(arch, grid, parameters, + compute_CATKE_diffusivities!, + diffusivities, grid, closure, velocities, tracers, buoyancy) + + return nothing +end + +@kernel function compute_average_surface_buoyancy_flux!(Jᵇ, grid, closure, velocities, tracers, + buoyancy, top_tracer_bcs, clock, Δt) + i, j = @index(Global, NTuple) + k = grid.Nz + + closure = getclosure(i, j, closure) + + model_fields = merge(velocities, tracers) + Jᵇ★ = top_buoyancy_flux(i, j, grid, buoyancy, top_tracer_bcs, clock, model_fields) + ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, Jᵇ) + + Jᵇᵋ = closure.minimum_convective_buoyancy_flux + Jᵇᵢⱼ = @inbounds Jᵇ[i, j, 1] + Jᵇ⁺ = max(Jᵇᵋ, Jᵇᵢⱼ, Jᵇ★) # selects fastest (dominant) time-scale + t★ = (ℓᴰ^2 / Jᵇ⁺)^(1/3) + ϵ = Δt / t★ + + @inbounds Jᵇ[i, j, 1] = (Jᵇᵢⱼ + ϵ * Jᵇ★) / (1 + ϵ) +end + +@kernel function compute_CATKE_diffusivities!(diffusivities, grid, closure::FlavorOfCATKE, velocities, tracers, buoyancy) + i, j, k = @index(Global, NTuple) + + # Ensure this works with "ensembles" of closures, in addition to ordinary single closures + closure_ij = getclosure(i, j, closure) + Jᵇ = diffusivities.Jᵇ + + # Note: we also compute the TKE diffusivity here for diagnostic purposes, even though it + # is recomputed in time_step_turbulent_kinetic_energy. + κu★ = κuᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, Jᵇ) + κc★ = κcᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, Jᵇ) + κe★ = κeᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, Jᵇ) + + κu★ = mask_diffusivity(i, j, k, grid, κu★) + κc★ = mask_diffusivity(i, j, k, grid, κc★) + κe★ = mask_diffusivity(i, j, k, grid, κe★) + + @inbounds begin + diffusivities.κu[i, j, k] = κu★ + diffusivities.κc[i, j, k] = κc★ + diffusivities.κe[i, j, k] = κe★ + end +end + +@inline function κuᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) + w★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocityᶜᶜᶜ, closure, tracers.e) + ℓu = momentum_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) + κu = ℓu * w★ + κu_max = closure.maximum_viscosity + return min(κu, κu_max) +end + +@inline function κcᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) + w★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocityᶜᶜᶜ, closure, tracers.e) + ℓc = tracer_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) + κc = ℓc * w★ + κc_max = closure.maximum_tracer_diffusivity + return min(κc, κc_max) +end + +@inline function κeᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) + w★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocityᶜᶜᶜ, closure, tracers.e) + ℓe = TKE_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) + κe = ℓe * w★ + κe_max = closure.maximum_tke_diffusivity + return min(κe, κe_max) +end + +@inline viscosity(::FlavorOfCATKE, diffusivities) = diffusivities.κu +@inline diffusivity(::FlavorOfCATKE, diffusivities, ::Val{id}) where id = diffusivities._tupled_tracer_diffusivities[id] + +##### +##### Show +##### + +function Base.summary(closure::CATKEVD) + TD = nameof(typeof(time_discretization(closure))) + return string("CATKEVerticalDiffusivity{$TD}") +end + +function Base.show(io::IO, clo::CATKEVD) + print(io, summary(clo)) + print(io, '\n') + print(io, "├── maximum_tracer_diffusivity: ", prettysummary(clo.maximum_tracer_diffusivity), '\n', + "├── maximum_tke_diffusivity: ", prettysummary(clo.maximum_tke_diffusivity), '\n', + "├── maximum_viscosity: ", prettysummary(clo.maximum_viscosity), '\n', + "├── minimum_tke: ", prettysummary(clo.minimum_tke), '\n', + "├── negative_tke_time_scale: ", prettysummary(clo.negative_tke_damping_time_scale), '\n', + "├── minimum_convective_buoyancy_flux: ", prettysummary(clo.minimum_convective_buoyancy_flux), '\n', + "├── tke_time_step: ", prettysummary(clo.tke_time_step), '\n', + "├── mixing_length: ", prettysummary(clo.mixing_length), '\n', + "│ ├── Cˢ: ", prettysummary(clo.mixing_length.Cˢ), '\n', + "│ ├── Cᵇ: ", prettysummary(clo.mixing_length.Cᵇ), '\n', + "│ ├── Cʰⁱu: ", prettysummary(clo.mixing_length.Cʰⁱu), '\n', + "│ ├── Cʰⁱc: ", prettysummary(clo.mixing_length.Cʰⁱc), '\n', + "│ ├── Cʰⁱe: ", prettysummary(clo.mixing_length.Cʰⁱe), '\n', + "│ ├── Cˡᵒu: ", prettysummary(clo.mixing_length.Cˡᵒu), '\n', + "│ ├── Cˡᵒc: ", prettysummary(clo.mixing_length.Cˡᵒc), '\n', + "│ ├── Cˡᵒe: ", prettysummary(clo.mixing_length.Cˡᵒe), '\n', + "│ ├── Cᵘⁿu: ", prettysummary(clo.mixing_length.Cᵘⁿu), '\n', + "│ ├── Cᵘⁿc: ", prettysummary(clo.mixing_length.Cᵘⁿc), '\n', + "│ ├── Cᵘⁿe: ", prettysummary(clo.mixing_length.Cᵘⁿe), '\n', + "│ ├── Cᶜu: ", prettysummary(clo.mixing_length.Cᶜu), '\n', + "│ ├── Cᶜc: ", prettysummary(clo.mixing_length.Cᶜc), '\n', + "│ ├── Cᶜe: ", prettysummary(clo.mixing_length.Cᶜe), '\n', + "│ ├── Cᵉc: ", prettysummary(clo.mixing_length.Cᵉc), '\n', + "│ ├── Cᵉe: ", prettysummary(clo.mixing_length.Cᵉe), '\n', + "│ ├── Cˢᵖ: ", prettysummary(clo.mixing_length.Cˢᵖ), '\n', + "│ ├── CRiᵟ: ", prettysummary(clo.mixing_length.CRiᵟ), '\n', + "│ └── CRi⁰: ", prettysummary(clo.mixing_length.CRi⁰), '\n', + "└── turbulent_kinetic_energy_equation: ", prettysummary(clo.turbulent_kinetic_energy_equation), '\n', + " ├── CʰⁱD: ", prettysummary(clo.turbulent_kinetic_energy_equation.CʰⁱD), '\n', + " ├── CˡᵒD: ", prettysummary(clo.turbulent_kinetic_energy_equation.CˡᵒD), '\n', + " ├── CᵘⁿD: ", prettysummary(clo.turbulent_kinetic_energy_equation.CᵘⁿD), '\n', + " ├── CᶜD: ", prettysummary(clo.turbulent_kinetic_energy_equation.CᶜD), '\n', + " ├── CᵉD: ", prettysummary(clo.turbulent_kinetic_energy_equation.CᵉD), '\n', + " ├── Cᵂu★: ", prettysummary(clo.turbulent_kinetic_energy_equation.Cᵂu★), '\n', + " ├── CᵂwΔ: ", prettysummary(clo.turbulent_kinetic_energy_equation.CᵂwΔ), '\n', + " └── Cᵂϵ: ", prettysummary(clo.turbulent_kinetic_energy_equation.Cᵂϵ)) +end + diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/time_step_turbulent_kinetic_energy.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl similarity index 90% rename from src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/time_step_turbulent_kinetic_energy.jl rename to src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl index a52168e4b4..3c5130b320 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/time_step_turbulent_kinetic_energy.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl @@ -7,15 +7,9 @@ using Oceananigans.TimeSteppers: store_field_tendencies!, ab2_step_field!, impli using Oceananigans.TurbulenceClosures: ∇_dot_qᶜ, immersed_∇_dot_qᶜ, hydrostatic_turbulent_kinetic_energy_tendency using CUDA -tke_time_step(closure::CATKEVerticalDiffusivity) = closure.turbulent_kinetic_energy_time_step +get_time_step(closure::CATKEVerticalDiffusivity) = closure.tke_time_step -function tke_time_step(closure_array::AbstractArray) - # assume they are all the same - closure = CUDA.@allowscalar closure_array[1, 1] - return tke_time_step(closure) -end - -function time_step_turbulent_kinetic_energy!(model) +function time_step_catke_equation!(model) # TODO: properly handle closure tuples closure = model.closure @@ -33,13 +27,13 @@ function time_step_turbulent_kinetic_energy!(model) tracer_index = findfirst(k -> k == :e, keys(model.tracers)) implicit_solver = model.timestepper.implicit_solver - Δt = model.clock.last_Δt # simulation time-step (for velocities, tracers, etc) - Δτ = tke_time_step(closure) # special time-step for turbulent kinetic energy + Δt = model.clock.last_Δt + Δτ = get_time_step(closure) if isnothing(Δτ) Δτ = Δt M = 1 - else # limit TKE time-step Δτ by Δt. + else M = ceil(Int, Δt / Δτ) # number of substeps Δτ = Δt / M end @@ -98,7 +92,7 @@ end wb⁺ = max(zero(grid), wb) eⁱʲᵏ = @inbounds e[i, j, k] - eᵐⁱⁿ = closure_ij.minimum_turbulent_kinetic_energy + eᵐⁱⁿ = closure_ij.minimum_tke wb⁻_e = wb⁻ / eⁱʲᵏ * (eⁱʲᵏ > eᵐⁱⁿ) # Treat the divergence of TKE flux at solid bottoms implicitly. @@ -122,7 +116,7 @@ end on_bottom = !inactive_cell(i, j, k, grid) & inactive_cell(i, j, k-1, grid) Δz = Δzᶜᶜᶜ(i, j, k, grid) Cᵂϵ = closure_ij.turbulent_kinetic_energy_equation.Cᵂϵ - e⁺ = clip(eⁱʲᵏ) # ensure that eⁱʲᵏ > 0 + e⁺ = clip(eⁱʲᵏ) w★ = sqrt(e⁺) div_Jᵉ_e = - on_bottom * Cᵂϵ * w★ / Δz @@ -178,7 +172,7 @@ end end #= -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: FlavorOfCATKE +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: FlavorOfCATKE @inline tracer_tendency_kernel_function(model::HFSM, name, c, K) = compute_hydrostatic_free_surface_Gc!, c, K @inline tracer_tendency_kernel_function(model::HFSM, ::Val{:e}, c::FlavorOfCATKE, K) = compute_hydrostatic_free_surface_Ge!, c, K diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl new file mode 100644 index 0000000000..f142db2133 --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl @@ -0,0 +1,299 @@ +using Oceananigans: fields +using Oceananigans.Advection: div_Uc, U_dot_∇u, U_dot_∇v +using Oceananigans.Fields: immersed_boundary_condition +using Oceananigans.Grids: active_interior_map +using Oceananigans.BoundaryConditions: apply_x_bcs!, apply_y_bcs!, apply_z_bcs! +using Oceananigans.TimeSteppers: store_field_tendencies!, ab2_step_field!, implicit_step! +using Oceananigans.TurbulenceClosures: ∇_dot_qᶜ, immersed_∇_dot_qᶜ, hydrostatic_turbulent_kinetic_energy_tendency +using CUDA + +Base.@kwdef struct TKEDissipationEquations{FT} + Cᵋϵ :: FT = 1.92 + Cᴾϵ :: FT = 1.44 + Cᵇϵ :: FT = -0.65 + Cᵂu★ :: FT = 0.0 + CᵂwΔ :: FT = 0.0 + Cᵂα :: FT = 0.11 # Charnock parameter + gravitational_acceleration :: FT = 9.8065 + minimum_roughness_length :: FT = 1e-4 +end + +get_time_step(closure::TKEDissipationVerticalDiffusivity) = closure.tke_dissipation_time_step + +function time_step_tke_dissipation_equations!(model) + + # TODO: properly handle closure tuples + closure = model.closure + + e = model.tracers.e + ϵ = model.tracers.ϵ + arch = model.architecture + grid = model.grid + Gⁿe = model.timestepper.Gⁿ.e + G⁻e = model.timestepper.G⁻.e + Gⁿϵ = model.timestepper.Gⁿ.ϵ + G⁻ϵ = model.timestepper.G⁻.ϵ + + diffusivity_fields = model.diffusivity_fields + κe = diffusivity_fields.κe + κϵ = diffusivity_fields.κϵ + Le = diffusivity_fields.Le + Lϵ = diffusivity_fields.Lϵ + previous_velocities = diffusivity_fields.previous_velocities + e_index = findfirst(k -> k == :e, keys(model.tracers)) + ϵ_index = findfirst(k -> k == :ϵ, keys(model.tracers)) + implicit_solver = model.timestepper.implicit_solver + + Δt = model.clock.last_Δt + Δτ = get_time_step(closure) + + if isnothing(Δτ) + Δτ = Δt + M = 1 + else + M = ceil(Int, Δt / Δτ) # number of substeps + Δτ = Δt / M + end + + FT = eltype(grid) + + for m = 1:M # substep + if m == 1 && M != 1 + χ = convert(FT, -0.5) # Euler step for the first substep + else + χ = model.timestepper.χ + end + + # Compute the linear implicit component of the RHS (diffusivities, L) + # and step forward + launch!(arch, grid, :xyz, + substep_tke_dissipation!, + κe, κϵ, Le, Lϵ, + grid, closure, + model.velocities, previous_velocities, # try this soon: model.velocities, model.velocities, + model.tracers, model.buoyancy, diffusivity_fields, + Δτ, χ, Gⁿe, G⁻e, Gⁿϵ, G⁻ϵ) + + implicit_step!(e, implicit_solver, closure, + model.diffusivity_fields, Val(e_index), + model.clock, Δτ) + + implicit_step!(ϵ, implicit_solver, closure, + model.diffusivity_fields, Val(ϵ_index), + model.clock, Δτ) + end + + return nothing +end + +@kernel function substep_tke_dissipation!(κe, κϵ, Le, Lϵ, + grid, closure, + next_velocities, previous_velocities, + tracers, buoyancy, diffusivities, + Δτ, χ, slow_Gⁿe, G⁻e, slow_Gⁿϵ, G⁻ϵ) + + i, j, k = @index(Global, NTuple) + + e = tracers.e + ϵ = tracers.ϵ + + closure_ij = getclosure(i, j, closure) + + # Compute TKE and dissipation diffusivities + κe★ = κeᶜᶜᶠ(i, j, k, grid, closure_ij, next_velocities, tracers, buoyancy) + κϵ★ = κϵᶜᶜᶠ(i, j, k, grid, closure_ij, next_velocities, tracers, buoyancy) + + κe★ = mask_diffusivity(i, j, k, grid, κe★) + κϵ★ = mask_diffusivity(i, j, k, grid, κϵ★) + + @inbounds κe[i, j, k] = κe★ + @inbounds κϵ[i, j, k] = κϵ★ + + # Compute TKE and dissipation tendencies + ϵ★ = dissipationᶜᶜᶜ(i, j, k, grid, closure_ij, tracers, buoyancy) + e★ = turbulent_kinetic_energyᶜᶜᶜ(i, j, k, grid, closure_ij, tracers) + eⁱʲᵏ = @inbounds e[i, j, k] + ϵⁱʲᵏ = @inbounds ϵ[i, j, k] + + # Different destruction time-scales for TKE vs dissipation for numerical reasons + ω★ = ϵ★ / e★ # target / physical dissipation time scale + ωe⁻ = closure_ij.negative_tke_damping_time_scale + ωe = ifelse(eⁱʲᵏ < 0, ωe⁻, ω★) + ωϵ = ϵⁱʲᵏ / e★ + + # Compute additional diagonal component of the linear TKE operator + wb = explicit_buoyancy_flux(i, j, k, grid, closure_ij, next_velocities, tracers, buoyancy, diffusivities) + + # Patankar trick for TKE equation + wb⁻ = min(wb, zero(grid)) + wb⁺ = max(wb, zero(grid)) + + eᵐⁱⁿ = closure_ij.minimum_tke + wb⁻_e = wb⁻ / eⁱʲᵏ * (eⁱʲᵏ > eᵐⁱⁿ) + + # Patankar trick for ϵ-equation + Cᵋϵ = closure_ij.tke_dissipation_equations.Cᵋϵ + Cᵇϵ = closure_ij.tke_dissipation_equations.Cᵇϵ + + Cᵇϵ_wb⁻ = min(Cᵇϵ * wb, zero(grid)) + Cᵇϵ_wb⁺ = max(Cᵇϵ * wb, zero(grid)) + + # ∂t e = Lⁱ e + ⋯, + @inbounds Le[i, j, k] = wb⁻_e - ωe + @inbounds Lϵ[i, j, k] = Cᵇϵ_wb⁻ / e★ - Cᵋϵ * ωϵ + + # Compute fast TKE and dissipation RHSs + u⁺ = next_velocities.u + v⁺ = next_velocities.v + uⁿ = previous_velocities.u + vⁿ = previous_velocities.v + κu = diffusivities.κu + Cᴾϵ = closure_ij.tke_dissipation_equations.Cᴾϵ + + # TODO: correctly handle closure / diffusivity tuples + # TODO: the shear_production is actually a slow term so we _could_ precompute. + P = shear_production(i, j, k, grid, κu, uⁿ, u⁺, vⁿ, v⁺) + + @inbounds begin + fast_Gⁿe = P + wb⁺ # - ϵ (no implicit time stepping for now) + fast_Gⁿϵ = ωϵ * (Cᴾϵ * P + Cᵇϵ_wb⁺) + end + + # Advance TKE and store tendency + FT = eltype(χ) + Δτ = convert(FT, Δτ) + + # See below. + α = convert(FT, 1.5) + χ + β = convert(FT, 0.5) + χ + + @inbounds begin + total_Gⁿe = slow_Gⁿe[i, j, k] + fast_Gⁿe + total_Gⁿϵ = slow_Gⁿϵ[i, j, k] + fast_Gⁿϵ + + e[i, j, k] += Δτ * (α * total_Gⁿe - β * G⁻e[i, j, k]) + ϵ[i, j, k] += Δτ * (α * total_Gⁿϵ - β * G⁻ϵ[i, j, k]) + + G⁻e[i, j, k] = total_Gⁿe + G⁻ϵ[i, j, k] = total_Gⁿϵ + end +end + +@inline function implicit_linear_coefficient(i, j, k, grid, closure::FlavorOfTD{<:VITD}, K, ::Val{id}, args...) where id + L = K._tupled_implicit_linear_coefficients[id] + return @inbounds L[i, j, k] +end + +##### +##### TKE top boundary condition +##### + +@inline function top_tke_flux(i, j, grid, clock, fields, parameters, closure::FlavorOfTD, buoyancy) + closure = getclosure(i, j, closure) + + top_tracer_bcs = parameters.top_tracer_boundary_conditions + top_velocity_bcs = parameters.top_velocity_boundary_conditions + tke_dissipation_parameters = closure.tke_dissipation_equations + + return _top_tke_flux(i, j, grid, clock, fields, tke_dissipation_parameters, closure, + buoyancy, top_tracer_bcs, top_velocity_bcs) +end + +@inline function _top_tke_flux(i, j, grid, clock, fields, + parameters::TKEDissipationEquations, closure::TDVD, + buoyancy, top_tracer_bcs, top_velocity_bcs) + + wΔ³ = top_convective_turbulent_velocity_cubed(i, j, grid, clock, fields, buoyancy, top_tracer_bcs) + u★ = friction_velocity(i, j, grid, clock, fields, top_velocity_bcs) + + Cᵂu★ = parameters.Cᵂu★ + CᵂwΔ = parameters.CᵂwΔ + + return - Cᵂu★ * u★^3 #- CᵂwΔ * wΔ³ +end + +@inline function top_dissipation_flux(i, j, grid, clock, fields, parameters, closure::FlavorOfTD, buoyancy) + closure = getclosure(i, j, closure) + + top_tracer_bcs = parameters.top_tracer_boundary_conditions + top_velocity_bcs = parameters.top_velocity_boundary_conditions + tke_dissipation_parameters = closure.tke_dissipation_equations + + return _top_dissipation_flux(i, j, grid, clock, fields, tke_dissipation_parameters, closure, + buoyancy, top_tracer_bcs, top_velocity_bcs) +end + +@inline function _top_dissipation_flux(i, j, grid, clock, fields, parameters::TKEDissipationEquations, + closure::TDVD, buoyancy, top_tracer_bcs, top_velocity_bcs) + + 𝕊u₀ = closure.stability_functions.𝕊u₀ + σϵ = closure.stability_functions.Cσϵ + + u★ = friction_velocity(i, j, grid, clock, fields, top_velocity_bcs) + α = parameters.Cᵂα + g = parameters.gravitational_acceleration + ℓ_charnock = α * u★^2 / g + + ℓmin = parameters.minimum_roughness_length + ℓᵣ = max(ℓmin, ℓ_charnock) + + k = grid.Nz + e★ = turbulent_kinetic_energyᶜᶜᶜ(i, j, k, grid, closure, fields) + z = znode(i, j, k, grid, c, c, c) + d = - z + + return - 𝕊u₀^4 / σϵ * e★^2 / (d + ℓᵣ) +end + +##### +##### Utilities for model constructors +##### + +""" Add TKE boundary conditions specific to `TKEDissipationVerticalDiffusivity`. """ +function add_closure_specific_boundary_conditions(closure::FlavorOfTD, + user_bcs, + grid, + tracer_names, + buoyancy) + + top_tracer_bcs = top_tracer_boundary_conditions(grid, tracer_names, user_bcs) + top_velocity_bcs = top_velocity_boundary_conditions(grid, user_bcs) + parameters = TKETopBoundaryConditionParameters(top_tracer_bcs, top_velocity_bcs) + top_tke_bc = FluxBoundaryCondition(top_tke_flux, discrete_form=true, parameters=parameters) + + top_dissipation_bc = FluxBoundaryCondition(top_dissipation_flux, discrete_form=true, parameters=parameters) + + + if :e ∈ keys(user_bcs) + e_bcs = user_bcs[:e] + + tke_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), + top = top_tke_bc, + bottom = e_bcs.bottom, + north = e_bcs.north, + south = e_bcs.south, + east = e_bcs.east, + west = e_bcs.west) + else + tke_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), top=top_tke_bc) + end + + if :ϵ ∈ keys(user_bcs) + ϵ_bcs = user_bcs[:ϵ] + + dissipation_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), + top = top_dissipation_bc, + bottom = e_bcs.bottom, + north = e_bcs.north, + south = e_bcs.south, + east = e_bcs.east, + west = e_bcs.west) + else + dissipation_bcs = FieldBoundaryConditions(grid, (Center, Center, Center), top=top_dissipation_bc) + end + + new_boundary_conditions = merge(user_bcs, (e=tke_bcs, ϵ=dissipation_bcs)) + + return new_boundary_conditions +end + diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl new file mode 100644 index 0000000000..ad130a2c42 --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_stability_functions.jl @@ -0,0 +1,292 @@ +abstract type AbstractConstantSchmidtStabilityFunctions end + +const ConstantSchmidtStabilityTDVD = TKEDissipationVerticalDiffusivity{<:Any, <:Any, <:AbstractConstantSchmidtStabilityFunctions} + +@inline function tke_stability_functionᶜᶜᶠ(i, j, k, grid, closure::ConstantSchmidtStabilityTDVD, args...) + Cσe = closure.stability_functions.Cσe + 𝕊u = momentum_stability_functionᶜᶜᶠ(i, j, k, grid, closure, args...) + return 𝕊u / Cσe +end + +@inline function dissipation_stability_functionᶜᶜᶠ(i, j, k, grid, closure::ConstantSchmidtStabilityTDVD, args...) + Cσϵ = closure.stability_functions.Cσϵ + 𝕊u = momentum_stability_functionᶜᶜᶠ(i, j, k, grid, closure, args...) + return 𝕊u / Cσϵ +end + +Base.@kwdef struct ConstantStabilityFunctions{FT} <: AbstractConstantSchmidtStabilityFunctions + Cσe :: FT = 1.0 + Cσϵ :: FT = 1.2 + Cu₀ :: FT = 0.53 # √3 + Cc₀ :: FT = 0.53 # √3 + 𝕊u₀ :: FT = 0.53 # √3 +end + +Base.summary(s::ConstantStabilityFunctions{FT}) where FT = "ConstantStabilityFunctions{$FT}" + +summarize_stability_functions(s::ConstantStabilityFunctions{FT}, prefix="", sep="│ ") where FT = + string(prefix, "ConstantStabilityFunctions{$FT}:", '\n', + " ├── 𝕊u₀: ", prettysummary(s.𝕊u₀), '\n', + " ├── Cσe: ", prettysummary(s.Cσe), '\n', + " ├── Cσϵ: ", prettysummary(s.Cσϵ), '\n', + " ├── Cu₀: ", prettysummary(s.Cu₀), '\n', + " └── Cc₀: ", prettysummary(s.Cc₀)) + +const ConstantStabilityTDVD = TKEDissipationVerticalDiffusivity{<:Any, <:Any, <:ConstantStabilityFunctions} + +@inline momentum_stability_functionᶜᶜᶠ(i, j, k, grid, c::ConstantStabilityTDVD, args...) = c.stability_functions.Cu₀ +@inline tracer_stability_functionᶜᶜᶠ(i, j, k, grid, c::ConstantStabilityTDVD, args...) = c.stability_functions.Cc₀ + +struct VariableStabilityFunctions{FT} <: AbstractConstantSchmidtStabilityFunctions + Cσe :: FT + Cσϵ :: FT + Cu₀ :: FT + Cu₁ :: FT + Cu₂ :: FT + Cc₀ :: FT + Cc₁ :: FT + Cc₂ :: FT + Cd₀ :: FT + Cd₁ :: FT + Cd₂ :: FT + Cd₃ :: FT + Cd₄ :: FT + Cd₅ :: FT + 𝕊u₀ :: FT +end + +function VariableStabilityFunctions(FT=Float64; + Cσe = 1.0, + Cσϵ = 1.2, + Cu₀ = 0.1067, + Cu₁ = 0.0173, + Cu₂ = -0.0001205, + Cc₀ = 0.1120, + Cc₁ = 0.003766, + Cc₂ = 0.0008871, + Cd₀ = 1.0, + Cd₁ = 0.2398, + Cd₂ = 0.02872, + Cd₃ = 0.005154, + Cd₄ = 0.006930, + Cd₅ = -0.0003372, + 𝕊u₀ = nothing) + + if isnothing(𝕊u₀) + # Compute 𝕊u₀ for the logarithmic boundary layer where production + # balances dissipation. For more information see the discussion + # surrounding equation (13) in Umlauf and Burchard (2003). + a = Cd₅ - Cu₂ + b = Cd₂ - Cu₀ + c = Cd₀ + 𝕊u₀ = (2a / (-b - sqrt(b^2 - 4a * c)))^(1/4) + end + + return VariableStabilityFunctions(convert(FT, Cσe), + convert(FT, Cσϵ), + convert(FT, Cu₀), + convert(FT, Cu₁), + convert(FT, Cu₂), + convert(FT, Cc₀), + convert(FT, Cc₁), + convert(FT, Cc₂), + convert(FT, Cd₀), + convert(FT, Cd₁), + convert(FT, Cd₂), + convert(FT, Cd₃), + convert(FT, Cd₄), + convert(FT, Cd₅), + convert(FT, 𝕊u₀)) +end + +Base.summary(s::VariableStabilityFunctions{FT}) where FT = "VariableStabilityFunctions{$FT}" + +summarize_stability_functions(s::VariableStabilityFunctions{FT}, prefix="", sep="") where FT = + string("VariableStabilityFunctions{$FT}:", '\n', + " ├── Cσe: ", prettysummary(s.Cσe), '\n', + " ├── Cσϵ: ", prettysummary(s.Cσϵ), '\n', + " ├── Cu₀: ", prettysummary(s.Cu₀), '\n', + " ├── Cu₁: ", prettysummary(s.Cu₁), '\n', + " ├── Cu₂: ", prettysummary(s.Cu₂), '\n', + " ├── Cc₀: ", prettysummary(s.Cc₀), '\n', + " ├── Cc₁: ", prettysummary(s.Cc₁), '\n', + " ├── Cc₂: ", prettysummary(s.Cc₂), '\n', + " ├── Cd₀: ", prettysummary(s.Cd₀), '\n', + " ├── Cd₁: ", prettysummary(s.Cd₁), '\n', + " ├── Cd₂: ", prettysummary(s.Cd₂), '\n', + " ├── Cd₃: ", prettysummary(s.Cd₃), '\n', + " ├── Cd₄: ", prettysummary(s.Cd₄), '\n', + " └── Cd₅: ", prettysummary(s.Cd₅)) + +@inline function square_time_scaleᶜᶜᶜ(i, j, k, grid, closure, tracers, buoyancy) + e★ = turbulent_kinetic_energyᶜᶜᶜ(i, j, k, grid, closure, tracers) + ϵ★ = dissipationᶜᶜᶜ(i, j, k, grid, closure, tracers, buoyancy) + return e★^2 / ϵ★^2 +end + +@inline function shear_numberᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy) + τ² = ℑzᵃᵃᶠ(i, j, k, grid, square_time_scaleᶜᶜᶜ, closure, tracers, buoyancy) + S² = shearᶜᶜᶠ(i, j, k, grid, velocities.u, velocities.v) + return τ² * S² +end + +@inline function stratification_numberᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) + τ² = ℑzᵃᵃᶠ(i, j, k, grid, square_time_scaleᶜᶜᶜ, closure, tracers, buoyancy) + N² = ∂z_b(i, j, k, grid, buoyancy, tracers) + return τ² * N² +end + +@inline maximum_stratification_number(closure) = 1e10 # ? + +""" +Based on an argument for free convection, assuming a balance between +buoyancy production and dissipation. + +See Umlauf and Burchard (2005) equation A.22. + +Note that _another_ condition could arise depending on the time discretization, +as discussed in the text surrounding equation 45-46 in Umlauf and Buchard (2005). +""" +@inline function minimum_stratification_number(closure) + m₀ = closure.stability_functions.Cc₀ + m₁ = closure.stability_functions.Cc₁ + m₂ = closure.stability_functions.Cc₂ + + d₀ = closure.stability_functions.Cd₀ + d₁ = closure.stability_functions.Cd₁ + d₂ = closure.stability_functions.Cd₂ + d₃ = closure.stability_functions.Cd₃ + d₄ = closure.stability_functions.Cd₄ + d₅ = closure.stability_functions.Cd₅ + + a = d₄ + m₁ + b = d₁ + m₀ + c = d₀ + + αᴺmin = (- b + sqrt(b^2 - 4a*c)) / 2a + + # Reduce by the "safety factor" + ϵ = closure.minimum_stratification_number_safety_factor + αᴺmin *= ϵ + + return αᴺmin +end + +@inline minimum_shear_number(closure::FlavorOfTD) = zero(eltype(closure)) + +""" +Based on the condition that shear aniostropy must increase. + +See Umlauf and Burchard (2005) equation 44. +""" +@inline function maximum_shear_number(closure, αᴺ) + n₀ = closure.stability_functions.Cu₀ + n₁ = closure.stability_functions.Cu₁ + n₂ = closure.stability_functions.Cu₂ + + d₀ = closure.stability_functions.Cd₀ + d₁ = closure.stability_functions.Cd₁ + d₂ = closure.stability_functions.Cd₂ + d₃ = closure.stability_functions.Cd₃ + d₄ = closure.stability_functions.Cd₄ + d₅ = closure.stability_functions.Cd₅ + + ϵ₀ = d₀ * n₀ + ϵ₁ = d₀ * n₁ + d₁ * n₀ + ϵ₂ = d₁ * n₁ + d₄ * n₀ + ϵ₃ = d₄ * n₁ + ϵ₄ = d₂ * n₀ + ϵ₅ = d₂ * n₁ + d₃ * n₀ + ϵ₆ = d₃ * n₁ + + num = ϵ₀ + ϵ₁ * αᴺ + ϵ₂ * αᴺ^2 + ϵ₃ * αᴺ^3 + den = ϵ₄ + ϵ₅ * αᴺ + ϵ₆ * αᴺ^2 + + return num / den +end + +const VariableStabilityTDVD = TKEDissipationVerticalDiffusivity{<:Any, <:Any, <:VariableStabilityFunctions} + +@inline function momentum_stability_functionᶜᶜᶠ(i, j, k, grid, closure::VariableStabilityTDVD, velocities, tracers, buoyancy) + αᴺ = stratification_numberᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) + αᴹ = shear_numberᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy) + + αᴺmin = minimum_stratification_number(closure) + αᴺmax = maximum_stratification_number(closure) + αᴺ = clamp(αᴺ, αᴺmin, αᴺmax) + + αᴹmin = minimum_shear_number(closure) + αᴹmax = maximum_shear_number(closure, αᴺ) + αᴹ = clamp(αᴹ, αᴹmin, αᴹmax) + + 𝕊u = momentum_stability_function(closure, αᴺ, αᴹ) + return 𝕊u +end + +@inline function momentum_stability_function(closure::VariableStabilityTDVD, αᴺ::Number, αᴹ::Number) + Cu₀ = closure.stability_functions.Cu₀ + Cu₁ = closure.stability_functions.Cu₁ + Cu₂ = closure.stability_functions.Cu₂ + + Cd₀ = closure.stability_functions.Cd₀ + Cd₁ = closure.stability_functions.Cd₁ + Cd₂ = closure.stability_functions.Cd₂ + Cd₃ = closure.stability_functions.Cd₃ + Cd₄ = closure.stability_functions.Cd₄ + Cd₅ = closure.stability_functions.Cd₅ + + num = Cu₀ + + Cu₁ * αᴺ + + Cu₂ * αᴹ + + den = Cd₀ + Cd₁ * αᴺ + + Cd₂ * αᴹ + + Cd₃ * αᴺ * αᴹ + + Cd₄ * αᴺ^2 + + Cd₅ * αᴹ^2 + + return num / den +end + +@inline function tracer_stability_functionᶜᶜᶠ(i, j, k, grid, closure::VariableStabilityTDVD, velocities, tracers, buoyancy) + αᴺ = stratification_numberᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) + αᴹ = shear_numberᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy) + + αᴺmin = minimum_stratification_number(closure) + αᴺmax = maximum_stratification_number(closure) + αᴺ = clamp(αᴺ, αᴺmin, αᴺmax) + + αᴹmin = minimum_shear_number(closure) + αᴹmax = maximum_shear_number(closure, αᴺ) + αᴹ = clamp(αᴹ, αᴹmin, αᴹmax) + + 𝕊c = tracer_stability_function(closure, αᴺ, αᴹ) + return 𝕊c +end + +@inline function tracer_stability_function(closure::VariableStabilityTDVD, αᴺ::Number, αᴹ::Number) + Cc₀ = closure.stability_functions.Cc₀ + Cc₁ = closure.stability_functions.Cc₁ + Cc₂ = closure.stability_functions.Cc₂ + + Cd₀ = closure.stability_functions.Cd₀ + Cd₁ = closure.stability_functions.Cd₁ + Cd₂ = closure.stability_functions.Cd₂ + Cd₃ = closure.stability_functions.Cd₃ + Cd₄ = closure.stability_functions.Cd₄ + Cd₅ = closure.stability_functions.Cd₅ + + num = Cc₀ + + Cc₁ * αᴺ + + Cc₂ * αᴹ + + den = Cd₀ + + Cd₁ * αᴺ + + Cd₂ * αᴹ + + Cd₃ * αᴺ * αᴹ + + Cd₄ * αᴺ^2 + + Cd₅ * αᴹ^2 + + return num / den +end + diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_vertical_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_vertical_diffusivity.jl new file mode 100644 index 0000000000..d94035e507 --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_vertical_diffusivity.jl @@ -0,0 +1,352 @@ +struct TKEDissipationVerticalDiffusivity{TD, KE, ST, LMIN, FT, DT} <: AbstractScalarDiffusivity{TD, VerticalFormulation, 2} + tke_dissipation_equations :: KE + stability_functions :: ST + minimum_length_scale :: LMIN + maximum_tracer_diffusivity :: FT + maximum_tke_diffusivity :: FT + maximum_dissipation_diffusivity :: FT + maximum_viscosity :: FT + minimum_tke :: FT + minimum_stratification_number_safety_factor :: FT + negative_tke_damping_time_scale :: FT + tke_dissipation_time_step :: DT +end + +function TKEDissipationVerticalDiffusivity{TD}(tke_dissipation_equations::KE, + stability_functions::ST, + minimum_length_scale :: LMIN, + maximum_tracer_diffusivity::FT, + maximum_tke_diffusivity::FT, + maximum_dissipation_diffusivity::FT, + maximum_viscosity::FT, + minimum_tke::FT, + minimum_stratification_number_safety_factor::FT, + negative_tke_damping_time_scale::FT, + tke_dissipation_time_step::DT) where {TD, KE, ST, LMIN, FT, DT} + + return TKEDissipationVerticalDiffusivity{TD, KE, ST, LMIN, FT, DT}(tke_dissipation_equations, + stability_functions, + minimum_length_scale, + maximum_tracer_diffusivity, + maximum_tke_diffusivity, + maximum_dissipation_diffusivity, + maximum_viscosity, + minimum_tke, + minimum_stratification_number_safety_factor, + negative_tke_damping_time_scale, + tke_dissipation_time_step) +end + +TKEDissipationVerticalDiffusivity(FT::DataType; kw...) = + TKEDissipationVerticalDiffusivity(VerticallyImplicitTimeDiscretization(), FT; kw...) + +const TDVD{TD} = TKEDissipationVerticalDiffusivity{TD} where TD +const TDVDArray{TD} = AbstractArray{<:TDVD{TD}} where TD +const FlavorOfTD{TD} = Union{TDVD{TD}, TDVDArray{TD}} where TD + +@inline Base.eltype(::TKEDissipationVerticalDiffusivity{<:Any, <:Any, <:Any, <:Any, FT}) where FT = FT + +""" + TKEDissipationVerticalDiffusivity([time_discretization = VerticallyImplicitTimeDiscretization(), + FT = Float64;] + tke_dissipation_equations = TKEDissipationEquations(), + stability_functions = VariableStabilityFunctions(), + minimum_length_scale = StratifiedDisplacementScale(), + maximum_tracer_diffusivity = Inf, + maximum_tke_diffusivity = Inf, + maximum_dissipation_diffusivity = Inf, + maximum_viscosity = Inf, + minimum_tke = 1e-6, + minimum_stratification_number_safety_factor = 0.73, + negative_tke_damping_time_scale = 1minute, + tke_dissipation_time_step = nothing) + +Return the `TKEDissipationVerticalDiffusivity` turbulence closure for vertical mixing by +microscale ocean turbulence based on the prognostic evolution of two variables: the +turbulent kinetic energy (TKE), and the turbulent kinetic energy dissipation. +Elsewhere this is referred to as "k-ϵ". For more information about k-ϵ, see +Burchard and Bolding (2001), Umlauf and Buchard (2003), and Umlauf and Burchard (2005). + +Arguments +========= + +- `time_discretization`: Either `ExplicitTimeDiscretization()` or `VerticallyImplicitTimeDiscretization()`; + default `VerticallyImplicitTimeDiscretization()`. + +- `FT`: Float type; default `Float64`. + + +Keyword arguments +================= + +- `maximum_diffusivity`: Maximum value for tracer, momentum, and TKE diffusivities. + Used to clip the diffusivity when/if + TKEDissipationVerticalDiffusivity predicts diffusivities + that are too large. + Default: `Inf`. + +- `minimum_tke`: Minimum value for the turbulent kinetic energy. + Can be used to model the presence "background" TKE + levels due to, for example, mixing by breaking internal waves. + Default: 1e-9. + +- `negative_tke_damping_time_scale`: Damping time-scale for spurious negative values of TKE, + typically generated by oscillatory errors associated + with TKE advection. + Default: 1 minute. + +Note that for numerical stability, it is recommended to either have a relative short +`negative_turbulent_kinetic_energy_damping_time_scale` or a reasonable +`minimum_turbulent_kinetic_energy`, or both. +""" +function TKEDissipationVerticalDiffusivity(time_discretization::TD = VerticallyImplicitTimeDiscretization(), + FT = Float64; + tke_dissipation_equations = TKEDissipationEquations(), + stability_functions = VariableStabilityFunctions(), + minimum_length_scale = StratifiedDisplacementScale(), + maximum_tracer_diffusivity = Inf, + maximum_tke_diffusivity = Inf, + maximum_dissipation_diffusivity = Inf, + maximum_viscosity = Inf, + minimum_tke = 1e-6, + minimum_stratification_number_safety_factor = 0.73, + negative_tke_damping_time_scale = 1minute, + tke_dissipation_time_step = nothing) where TD + + stability_functions = convert_eltype(FT, stability_functions) + + return TKEDissipationVerticalDiffusivity{TD}(tke_dissipation_equations, + stability_functions, + minimum_length_scale, + convert(FT, maximum_tracer_diffusivity), + convert(FT, maximum_tke_diffusivity), + convert(FT, maximum_dissipation_diffusivity), + convert(FT, maximum_viscosity), + convert(FT, minimum_tke), + convert(FT, minimum_stratification_number_safety_factor), + convert(FT, negative_tke_damping_time_scale), + tke_dissipation_time_step) +end + +function with_tracers(tracer_names, closure::FlavorOfTD) + :e ∈ tracer_names && :ϵ ∈ tracer_names || + throw(ArgumentError("Tracers must contain :e and :ϵ to represent turbulent kinetic energy " * + "for `TKEDissipationVerticalDiffusivity`.")) + + return closure +end + +##### +##### Stratified displacement length scale limiter +##### + +Base.@kwdef struct StratifiedDisplacementScale{FT} + Cᴺ :: FT = 0.75 + minimum_buoyancy_frequency :: FT = 1e-14 +end + +##### +##### Diffusivities and diffusivity fields utilities +##### + +function DiffusivityFields(grid, tracer_names, bcs, closure::FlavorOfTD) + + default_diffusivity_bcs = (κu = FieldBoundaryConditions(grid, (Center, Center, Face)), + κc = FieldBoundaryConditions(grid, (Center, Center, Face)), + κe = FieldBoundaryConditions(grid, (Center, Center, Face)), + κϵ = FieldBoundaryConditions(grid, (Center, Center, Face))) + + bcs = merge(default_diffusivity_bcs, bcs) + + κu = ZFaceField(grid, boundary_conditions=bcs.κu) + κc = ZFaceField(grid, boundary_conditions=bcs.κc) + κe = ZFaceField(grid, boundary_conditions=bcs.κe) + κϵ = ZFaceField(grid, boundary_conditions=bcs.κϵ) + Le = CenterField(grid) + Lϵ = CenterField(grid) + + # Note: we may be able to avoid using the "previous velocities" in favor of a "fully implicit" + # discretization of shear production + u⁻ = XFaceField(grid) + v⁻ = YFaceField(grid) + previous_velocities = (; u=u⁻, v=v⁻) + + # Secret tuple for getting tracer diffusivities with tuple[tracer_index] + _tupled_tracer_diffusivities = Dict{Symbol, Any}(name => κc for name in tracer_names) + _tupled_tracer_diffusivities[:e] = κe + _tupled_tracer_diffusivities[:ϵ] = κϵ + _tupled_tracer_diffusivities = NamedTuple(name => _tupled_tracer_diffusivities[name] + for name in tracer_names) + + _tupled_implicit_linear_coefficients = Dict{Symbol, Any}(name => ZeroField() for name in tracer_names) + _tupled_implicit_linear_coefficients[:e] = Le + _tupled_implicit_linear_coefficients[:ϵ] = Lϵ + _tupled_implicit_linear_coefficients = NamedTuple(name => _tupled_implicit_linear_coefficients[name] + for name in tracer_names) + + return (; κu, κc, κe, κϵ, Le, Lϵ, previous_velocities, + _tupled_tracer_diffusivities, _tupled_implicit_linear_coefficients) +end + +@inline viscosity_location(::FlavorOfTD) = (c, c, f) +@inline diffusivity_location(::FlavorOfTD) = (c, c, f) + +function compute_diffusivities!(diffusivities, closure::FlavorOfTD, model; parameters = :xyz) + + arch = model.architecture + grid = model.grid + velocities = model.velocities + tracers = model.tracers + buoyancy = model.buoyancy + clock = model.clock + top_tracer_bcs = NamedTuple(c => tracers[c].boundary_conditions.top for c in propertynames(tracers)) + + if isfinite(model.clock.last_Δt) # Check that we have taken a valid time-step first. + # Compute e at the current time: + # * update tendency Gⁿ using current and previous velocity field + # * use tridiagonal solve to take an implicit step + time_step_tke_dissipation_equations!(model) + end + + # Update "previous velocities" + u, v, w = model.velocities + u⁻, v⁻ = diffusivities.previous_velocities + parent(u⁻) .= parent(u) + parent(v⁻) .= parent(v) + + launch!(arch, grid, parameters, + compute_TKEDissipation_diffusivities!, + diffusivities, grid, closure, velocities, tracers, buoyancy) + + return nothing +end + +@kernel function compute_TKEDissipation_diffusivities!(diffusivities, grid, closure::FlavorOfTD, + velocities, tracers, buoyancy) + i, j, k = @index(Global, NTuple) + + # Ensure this works with "ensembles" of closures, in addition to ordinary single closures + closure_ij = getclosure(i, j, closure) + + # Note: we also compute the TKE diffusivity here for diagnostic purposes, even though it + # is recomputed in time_step_turbulent_kinetic_energy. + κu★ = κuᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy) + κc★ = κcᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy) + κe★ = κeᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy) + κϵ★ = κϵᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy) + + κu★ = mask_diffusivity(i, j, k, grid, κu★) + κc★ = mask_diffusivity(i, j, k, grid, κc★) + κe★ = mask_diffusivity(i, j, k, grid, κe★) + κϵ★ = mask_diffusivity(i, j, k, grid, κϵ★) + + @inbounds begin + diffusivities.κu[i, j, k] = κu★ + diffusivities.κc[i, j, k] = κc★ + diffusivities.κe[i, j, k] = κe★ + diffusivities.κϵ[i, j, k] = κϵ★ + end +end + +@inline function turbulent_kinetic_energyᶜᶜᶜ(i, j, k, grid, closure, tracers) + eᵐⁱⁿ = closure.minimum_tke + eⁱʲᵏ = @inbounds tracers.e[i, j, k] + return max(eᵐⁱⁿ, eⁱʲᵏ) +end + +@inline max_a_b(i, j, k, grid, a::Number, b, args...) = max(a, b(i, j, k, grid, args...)) + +@inline maximum_dissipation(i, j, k, grid, closure, tracers, buoyancy) = convert(eltype(grid), Inf) + +@inline function minimum_dissipation(i, j, k, grid, closure, tracers, buoyancy) + FT = eltype(grid) + + N²min = closure.minimum_length_scale.minimum_buoyancy_frequency + N²⁺ = ℑbzᵃᵃᶜ(i, j, k, grid, max_a_b, N²min, ∂z_b, buoyancy, tracers) + + Cᴺ = closure.minimum_length_scale.Cᴺ + e★ = turbulent_kinetic_energyᶜᶜᶜ(i, j, k, grid, closure, tracers) + ℓst = Cᴺ * sqrt(e★ / N²⁺) + + 𝕊u₀ = closure.stability_functions.𝕊u₀ + ℓmin = min(grid.Lz, ℓst) + ϵmin = 𝕊u₀^3 * sqrt(e★)^3 / ℓmin + + another_ϵmin = convert(FT, 1e-12) + return max(another_ϵmin, ϵmin) +end + +@inline function dissipationᶜᶜᶜ(i, j, k, grid, closure, tracers, buoyancy) + ϵᵐⁱⁿ = minimum_dissipation(i, j, k, grid, closure, tracers, buoyancy) + ϵᵐᵃˣ = maximum_dissipation(i, j, k, grid, closure, tracers, buoyancy) + ϵⁱʲᵏ = @inbounds tracers.ϵ[i, j, k] + return clamp(ϵⁱʲᵏ, ϵᵐⁱⁿ, ϵᵐᵃˣ) +end + +@inline function κuᶜᶜᶠ(i, j, k, grid, closure::TDVD, velocities, tracers, buoyancy) + e² = ℑzᵃᵃᶠ(i, j, k, grid, ϕ², turbulent_kinetic_energyᶜᶜᶜ, closure, tracers) + ϵ = ℑzᵃᵃᶠ(i, j, k, grid, dissipationᶜᶜᶜ, closure, tracers, buoyancy) + 𝕊u = momentum_stability_functionᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy) + κu = 𝕊u * e² / ϵ + κu_max = closure.maximum_viscosity + return min(κu, κu_max) +end + +@inline function κcᶜᶜᶠ(i, j, k, grid, closure::TDVD, velocities, tracers, buoyancy) + e² = ℑzᵃᵃᶠ(i, j, k, grid, ϕ², turbulent_kinetic_energyᶜᶜᶜ, closure, tracers) + ϵ = ℑzᵃᵃᶠ(i, j, k, grid, dissipationᶜᶜᶜ, closure, tracers, buoyancy) + 𝕊c = tracer_stability_functionᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy) + κc = 𝕊c * e² / ϵ + κc_max = closure.maximum_tracer_diffusivity + return min(κc, κc_max) +end + +@inline function κeᶜᶜᶠ(i, j, k, grid, closure::TDVD, velocities, tracers, buoyancy) + e² = ℑzᵃᵃᶠ(i, j, k, grid, ϕ², turbulent_kinetic_energyᶜᶜᶜ, closure, tracers) + ϵ = ℑzᵃᵃᶠ(i, j, k, grid, dissipationᶜᶜᶜ, closure, tracers, buoyancy) + 𝕊e = tke_stability_functionᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy) + κe = 𝕊e * e² / ϵ + κe_max = closure.maximum_tke_diffusivity + return min(κe, κe_max) +end + +@inline function κϵᶜᶜᶠ(i, j, k, grid, closure::TDVD, velocities, tracers, buoyancy) + e² = ℑzᵃᵃᶠ(i, j, k, grid, ϕ², turbulent_kinetic_energyᶜᶜᶜ, closure, tracers) + ϵ = ℑzᵃᵃᶠ(i, j, k, grid, dissipationᶜᶜᶜ, closure, tracers, buoyancy) + 𝕊ϵ = dissipation_stability_functionᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy) + κϵ = 𝕊ϵ * e² / ϵ + κϵ_max = closure.maximum_dissipation_diffusivity + return min(κϵ, κϵ_max) +end + +@inline viscosity(::FlavorOfTD, diffusivities) = diffusivities.κu +@inline diffusivity(::FlavorOfTD, diffusivities, ::Val{id}) where id = diffusivities._tupled_tracer_diffusivities[id] + +##### +##### Show +##### + +function Base.summary(closure::TDVD) + TD = nameof(typeof(time_discretization(closure))) + return string("TKEDissipationVerticalDiffusivity{$TD}") +end + +function Base.show(io::IO, clo::TDVD) + print(io, summary(clo)) + print(io, '\n') + print(io, "├── maximum_tracer_diffusivity: ", prettysummary(clo.maximum_tracer_diffusivity), '\n', + "├── maximum_tke_diffusivity: ", prettysummary(clo.maximum_tke_diffusivity), '\n', + "├── maximum_dissipation_diffusivity: ", prettysummary(clo.maximum_dissipation_diffusivity), '\n', + "├── maximum_viscosity: ", prettysummary(clo.maximum_viscosity), '\n', + "├── minimum_tke: ", prettysummary(clo.minimum_tke), '\n', + "├── negative_tke_damping_time_scale: ", prettysummary(clo.negative_tke_damping_time_scale), '\n', + "├── tke_dissipation_time_step: ", prettysummary(clo.tke_dissipation_time_step), '\n', + "├── tke_dissipation_equations: ", prettysummary(clo.tke_dissipation_equations), '\n', + "│ ├── Cᵋϵ: ", prettysummary(clo.tke_dissipation_equations.Cᵋϵ), '\n', + "│ ├── Cᴾϵ: ", prettysummary(clo.tke_dissipation_equations.Cᴾϵ), '\n', + "│ ├── Cᵇϵ: ", prettysummary(clo.tke_dissipation_equations.Cᵇϵ), '\n', + "│ ├── Cᵂu★: ", prettysummary(clo.tke_dissipation_equations.Cᵂu★), '\n', + "│ └── CᵂwΔ: ", prettysummary(clo.tke_dissipation_equations.CᵂwΔ), '\n') + print(io, "└── ", summarize_stability_functions(clo.stability_functions), "", " ") +end + diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_top_boundary_condition.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_top_boundary_condition.jl new file mode 100644 index 0000000000..f67db36c3e --- /dev/null +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_top_boundary_condition.jl @@ -0,0 +1,91 @@ +using Oceananigans.BoundaryConditions: DiscreteBoundaryFunction, BoundaryCondition, Flux + +struct TKETopBoundaryConditionParameters{C, U} + top_tracer_boundary_conditions :: C + top_velocity_boundary_conditions :: U +end + +const TKEBoundaryFunction = DiscreteBoundaryFunction{<:TKETopBoundaryConditionParameters} +const TKEBoundaryCondition = BoundaryCondition{<:Flux, <:TKEBoundaryFunction} + +@inline Adapt.adapt_structure(to, p::TKETopBoundaryConditionParameters) = + TKETopBoundaryConditionParameters(adapt(to, p.top_tracer_boundary_conditions), + adapt(to, p.top_velocity_boundary_conditions)) + +@inline on_architecture(to, p::TKETopBoundaryConditionParameters) = + TKETopBoundaryConditionParameters(on_architecture(to, p.top_tracer_boundary_conditions), + on_architecture(to, p.top_velocity_boundary_conditions)) + +@inline getbc(bc::TKEBoundaryCondition, i::Integer, j::Integer, grid::AbstractGrid, clock, fields, clo, buoyancy) = + bc.condition.func(i, j, grid, clock, fields, bc.condition.parameters, clo, buoyancy) + +@inline getbc(bc::TKEBoundaryCondition, i::Integer, j::Integer, k::Integer, grid::AbstractGrid, clock, fields, clo, buoyancy) = + bc.condition.func(i, j, k, grid, clock, fields, bc.condition.parameters, clo, buoyancy) + +""" + top_tke_flux(i, j, grid, clock, fields, parameters, closure, buoyancy) + +Compute the flux of TKE through the surface / top boundary. +Designed to be used with TKETopBoundaryConditionParameters in a FluxBoundaryCondition, eg: + +``` +top_tracer_bcs = top_tracer_boundary_conditions(grid, tracer_names, user_bcs) +top_velocity_bcs = top_velocity_boundary_conditions(grid, user_bcs) +parameters = TKETopBoundaryConditionParameters(top_tracer_bcs, top_velocity_bcs) +top_tke_bc = FluxBoundaryCondition(top_tke_flux, discrete_form=true, parameters=parameters) +``` + +See the implementation in catke_equation.jl. +""" +@inline top_tke_flux(i, j, grid, clock, fields, parameters, closure, buoyancy) = zero(grid) + +##### +##### For model constructors +##### + +""" Infer tracer boundary conditions from user_bcs and tracer_names. """ +function top_tracer_boundary_conditions(grid, tracer_names, user_bcs) + default_tracer_bcs = NamedTuple(c => FieldBoundaryConditions(grid, (Center, Center, Center)) for c in tracer_names) + bcs = merge(default_tracer_bcs, user_bcs) + return NamedTuple(c => bcs[c].top for c in tracer_names) +end + +""" Infer velocity boundary conditions from `user_bcs` and `tracer_names`. """ +function top_velocity_boundary_conditions(grid, user_bcs) + default_top_bc = default_prognostic_bc(topology(grid, 3)(), Center(), DefaultBoundaryCondition()) + + user_bc_names = keys(user_bcs) + u_top_bc = :u ∈ user_bc_names ? user_bcs.u.top : default_top_bc + v_top_bc = :v ∈ user_bc_names ? user_bcs.v.top : default_top_bc + + return (u=u_top_bc, v=v_top_bc) +end + +""" Computes the friction velocity u★ based on fluxes of u and v. """ +@inline function friction_velocity(i, j, grid, clock, fields, velocity_bcs) + FT = eltype(grid) + τx = getbc(velocity_bcs.u, i, j, grid, clock, fields) + τy = getbc(velocity_bcs.v, i, j, grid, clock, fields) + return sqrt(sqrt(τx^2 + τy^2)) +end + +""" Computes the convective velocity w★. """ +@inline function top_convective_turbulent_velocity_cubed(i, j, grid, clock, fields, buoyancy, tracer_bcs) + Jᵇ = top_buoyancy_flux(i, j, grid, buoyancy, tracer_bcs, clock, fields) + Δz = Δzᶜᶜᶜ(i, j, grid.Nz, grid) + return clip(Jᵇ) * Δz +end + +@inline top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple::Tuple{<:Any}, buoyancy) = + top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple[1], buoyancy) + +@inline top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple::Tuple{<:Any, <:Any}, buoyancy) = + top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple[1], buoyancy) + + top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple[2], buoyancy) + +@inline top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple::Tuple{<:Any, <:Any, <:Any}, buoyancy) = + top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple[1], buoyancy) + + top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple[2], buoyancy) + + top_tke_flux(i, j, grid, clock, fields, parameters, closure_tuple[3], buoyancy) + + diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/ri_based_vertical_diffusivity.jl b/src/TurbulenceClosures/turbulence_closure_implementations/ri_based_vertical_diffusivity.jl index d06d7df1b1..2d4b7bb2eb 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/ri_based_vertical_diffusivity.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/ri_based_vertical_diffusivity.jl @@ -177,17 +177,17 @@ const f = Face() @inline viscosity_location(::FlavorOfRBVD) = (c, c, f) @inline diffusivity_location(::FlavorOfRBVD) = (c, c, f) -@inline viscosity(::FlavorOfRBVD, diffusivities) = diffusivities.κᵘ -@inline diffusivity(::FlavorOfRBVD, diffusivities, id) = diffusivities.κᶜ +@inline viscosity(::FlavorOfRBVD, diffusivities) = diffusivities.κu +@inline diffusivity(::FlavorOfRBVD, diffusivities, id) = diffusivities.κc with_tracers(tracers, closure::FlavorOfRBVD) = closure # Note: computing diffusivities at cell centers for now. function DiffusivityFields(grid, tracer_names, bcs, closure::FlavorOfRBVD) - κᶜ = Field((Center, Center, Face), grid) - κᵘ = Field((Center, Center, Face), grid) + κc = Field((Center, Center, Face), grid) + κu = Field((Center, Center, Face), grid) Ri = Field((Center, Center, Face), grid) - return (; κᶜ, κᵘ, Ri) + return (; κc, κu, Ri) end function compute_diffusivities!(diffusivities, closure::FlavorOfRBVD, model; parameters = :xyz) @@ -312,30 +312,30 @@ end # Shear mixing diffusivity and viscosity τ = taper(tapering, Ri, Ri₀, Riᵟ) - κᶜ★ = κ₀ * τ - κᵘ★ = ν₀ * τ + κc★ = κ₀ * τ + κu★ = ν₀ * τ # Previous diffusivities - κᶜ = diffusivities.κᶜ - κᵘ = diffusivities.κᵘ + κc = diffusivities.κc + κu = diffusivities.κu # New diffusivities - κᶜ⁺ = κᶜᵃ + κᵉⁿ + κᶜ★ - κᵘ⁺ = κᵘ★ + κc⁺ = κᶜᵃ + κᵉⁿ + κc★ + κu⁺ = κu★ # Limit by specified maximum - κᶜ⁺ = min(κᶜ⁺, closure_ij.maximum_diffusivity) - κᵘ⁺ = min(κᵘ⁺, closure_ij.maximum_viscosity) + κc⁺ = min(κc⁺, closure_ij.maximum_diffusivity) + κu⁺ = min(κu⁺, closure_ij.maximum_viscosity) # Set to zero on periphery and NaN within inactive region on_periphery = peripheral_node(i, j, k, grid, c, c, f) within_inactive = inactive_node(i, j, k, grid, c, c, f) - κᶜ⁺ = ifelse(on_periphery, zero(grid), ifelse(within_inactive, NaN, κᶜ⁺)) - κᵘ⁺ = ifelse(on_periphery, zero(grid), ifelse(within_inactive, NaN, κᵘ⁺)) + κc⁺ = ifelse(on_periphery, zero(grid), ifelse(within_inactive, NaN, κc⁺)) + κu⁺ = ifelse(on_periphery, zero(grid), ifelse(within_inactive, NaN, κu⁺)) # Update by averaging in time - @inbounds κᶜ[i, j, k] = (Cᵃᵛ * κᶜ[i, j, k] + κᶜ⁺) / (1 + Cᵃᵛ) - @inbounds κᵘ[i, j, k] = (Cᵃᵛ * κᵘ[i, j, k] + κᵘ⁺) / (1 + Cᵃᵛ) + @inbounds κc[i, j, k] = (Cᵃᵛ * κc[i, j, k] + κc⁺) / (1 + Cᵃᵛ) + @inbounds κu[i, j, k] = (Cᵃᵛ * κu[i, j, k] + κu⁺) / (1 + Cᵃᵛ) return nothing end diff --git a/test/test_time_stepping.jl b/test/test_time_stepping.jl index 9b07247248..5e6c51b993 100644 --- a/test/test_time_stepping.jl +++ b/test/test_time_stepping.jl @@ -4,7 +4,7 @@ using TimesDates: TimeDate using Oceananigans.Grids: topological_tuple_length, total_size using Oceananigans.Fields: BackgroundField using Oceananigans.TimeSteppers: Clock -using Oceananigans.TurbulenceClosures.CATKEVerticalDiffusivities: CATKEVerticalDiffusivity +using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity function time_stepping_works_with_flat_dimensions(arch, topology) size = Tuple(1 for i = 1:topological_tuple_length(topology...)) diff --git a/validation/vertical_mixing_closures/column_windy_convection.jl b/validation/vertical_mixing_closures/column_windy_convection.jl index 16a5b60446..2dc96f17ba 100644 --- a/validation/vertical_mixing_closures/column_windy_convection.jl +++ b/validation/vertical_mixing_closures/column_windy_convection.jl @@ -5,12 +5,11 @@ using Printf using Oceananigans.TurbulenceClosures: RiBasedVerticalDiffusivity, - CATKEVerticalDiffusivity, - ConvectiveAdjustmentVerticalDiffusivity, - ExplicitTimeDiscretization + TKEDissipationVerticalDiffusivity, + CATKEVerticalDiffusivity # Parameters -Δz = 1 # Vertical resolution +Δz = 2 # Vertical resolution Lz = 256 # Extent of vertical domain Nz = Int(Lz/Δz) # Vertical resolution f₀ = 1e-4 # Coriolis parameter (s⁻¹) @@ -19,7 +18,7 @@ Jᵇ = +1e-7 # Surface buoyancy flux (m² s⁻³) τˣ = -2e-3 # Surface kinematic momentum flux (m s⁻¹) stop_time = 4days -convective_adjustment = ConvectiveAdjustmentVerticalDiffusivity(convective_κz=0.1, convective_νz=0.01) +tke_dissipation = TKEDissipationVerticalDiffusivity() catke = CATKEVerticalDiffusivity() ri_based = RiBasedVerticalDiffusivity() @@ -29,19 +28,36 @@ grid = RectilinearGrid(size=Nz, z=(-Lz, 0), topology=(Flat, Flat, Bounded)) coriolis = FPlane(f=f₀) b_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(Jᵇ)) u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(τˣ)) -closures_to_run = [catke] #, ri_based, convective_adjustment] +closures_to_run = [tke_dissipation, catke] +closure = tke_dissipation + +function progress(sim) + model = sim.model + + msg = @sprintf("Iter: % 5d, time: % 16s, max(b): %6.2e", + iteration(sim), prettytime(sim), maximum(model.tracers.b)) + + msg *= @sprintf(", max(κ): %6.2e, max(e): %6.2e, max(ϵ): %6.2e", + maximum(model.diffusivity_fields.κc), + maximum(model.tracers.e), + maximum(model.tracers.ϵ)) + + @info msg + + return nothing +end for closure in closures_to_run model = HydrostaticFreeSurfaceModel(; grid, closure, coriolis, - tracers = (:b, :e), + tracers = (:b, :e, :ϵ), buoyancy = BuoyancyTracer(), boundary_conditions = (; b=b_bcs, u=u_bcs)) bᵢ(z) = N² * z set!(model, b=bᵢ, e=1e-6) - simulation = Simulation(model; Δt=20minutes, stop_time) + simulation = Simulation(model; Δt=10.0, stop_time) closurename = string(nameof(typeof(closure))) @@ -50,16 +66,18 @@ for closure in closures_to_run outputs = merge(model.velocities, model.tracers, diffusivities) - simulation.output_writers[:fields] = JLD2OutputWriter(model, outputs, - schedule = TimeInterval(20minutes), - filename = "windy_convection_" * closurename, - overwrite_existing = true) + output_writer = JLD2OutputWriter(model, outputs, + schedule = TimeInterval(20minutes), + filename = "windy_convection_" * closurename, + overwrite_existing = true) - progress(sim) = @info string("Iter: ", iteration(sim), " t: ", prettytime(sim), - ", max(b): ", maximum(model.tracers.b)) - simulation.callbacks[:progress] = Callback(progress, IterationInterval(100)) + simulation.output_writers[:fields] = output_writer + + + add_callback!(simulation, progress, IterationInterval(10)) - @info "Running a simulation of $model..." + @info "Running a simulation of " + @info "$model" run!(simulation) end @@ -92,13 +110,12 @@ end b1 = first(b_ts) e1 = first(e_ts) κ1 = first(κc_ts) -@show maximum(e1) zc = znodes(b1) zf = znodes(κ1) Nt = length(b1.times) -fig = Figure(size=(1800, 600)) +fig = Figure(size=(1400, 600)) slider = Slider(fig[2, 1:4], range=1:Nt, startvalue=1) n = slider.value @@ -114,9 +131,9 @@ axe = Axis(fig[1, 3], xlabel=TKE_label, ylabel="z (m)") axκ = Axis(fig[1, 4], xlabel=diffusivities_label, ylabel="z (m)") xlims!(axb, -grid.Lz * N², 0) -xlims!(axu, -0.1, 0.1) -xlims!(axe, -1e-4, 2e-4) -xlims!(axκ, -1e-1, 5e-1) +xlims!(axu, -0.2, 0.2) +xlims!(axe, -1e-4, 1e-2) +xlims!(axκ, -1e-1, 1e1) colors = [:black, :blue, :red, :orange] diff --git a/validation/vertical_mixing_closures/compare_catke_k_epsilon.jl b/validation/vertical_mixing_closures/compare_catke_k_epsilon.jl new file mode 100644 index 0000000000..d54033ea8f --- /dev/null +++ b/validation/vertical_mixing_closures/compare_catke_k_epsilon.jl @@ -0,0 +1,114 @@ +using Oceananigans +using Oceananigans.Units + +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: + TKEDissipationVerticalDiffusivity, + CATKEVerticalDiffusivity, + VariableStabilityFunctions, + minimum_stratification_number, + maximum_shear_number, + ConstantStabilityFunctions, + stratification_numberᶜᶜᶠ, + shear_numberᶜᶜᶠ + +using GLMakie +using Printf + +grid = RectilinearGrid(size=128, z=(-128, 0), topology=(Flat, Flat, Bounded)) + +f = 0 +N² = 1e-5 +τˣ = -1e-4 +Jᵇ = 1e-7 + +u_top_bc = FluxBoundaryCondition(τˣ) +u_bcs = FieldBoundaryConditions(top=u_top_bc) + +b_top_bc = FluxBoundaryCondition(Jᵇ) +b_bcs = FieldBoundaryConditions(top=b_top_bc) + +coriolis = FPlane(; f) + +const_stability_functions = ConstantStabilityFunctions() +k_epsilon = TKEDissipationVerticalDiffusivity() +k_epsilon_const_stability = TKEDissipationVerticalDiffusivity(stability_functions=const_stability_functions) +catke = CATKEVerticalDiffusivity() + +bᵢ(z) = N² * z + +bn = [] +un = [] +vn = [] +en = [] +ϵn = [] + +#stratification_numberᶜᶜᶠ(i, j, k, grid, closure, tracers, buoyancy) + +for closure in (k_epsilon, catke) #, k_epsilon_const_stability) + global model + + model = HydrostaticFreeSurfaceModel(; grid, closure, coriolis, + tracers = (:b, :e, :ϵ), + buoyancy = BuoyancyTracer(), + boundary_conditions=(u=u_bcs, b=b_bcs)) + + set!(model, b=bᵢ) + + simulation = Simulation(model, Δt=60.0, stop_time=1day) + + local b, u, v, w, e, ϵ + + u, v, w = model.velocities + e = model.tracers.e + ϵ = model.tracers.ϵ + b = model.tracers.b + + progress(sim) = @info @sprintf("Iter: % 4d, time: % 24s, max(e): %6.2e, max(ϵ): %6.2e", + iteration(sim), prettytime(sim), maximum(e), maximum(ϵ)) + + add_callback!(simulation, progress, IterationInterval(100)) + + run!(simulation) + + push!(bn, deepcopy(interior(b, 1, 1, :))) + push!(un, deepcopy(interior(u, 1, 1, :))) + push!(vn, deepcopy(interior(v, 1, 1, :))) + push!(en, deepcopy(interior(e, 1, 1, :))) + push!(ϵn, deepcopy(interior(ϵ, 1, 1, :))) +end + +z = znodes(model.tracers.e) + +colors = [:black, :blue, :tomato] + +fig = Figure(size=(800, 400)) + +axb = Axis(fig[1, 1], title="Velocities") +axu = Axis(fig[1, 2], title="Velocities") +axe = Axis(fig[1, 3], title="TKE") +axϵ = Axis(fig[1, 4], title="Epsilon") + +lines!(axb, bn[1], z, color=colors[1]) +lines!(axu, un[1], z, color=colors[1], label="u") +lines!(axu, vn[1], z, color=colors[1], linestyle=:dash, label="v") +lines!(axe, en[1], z, color=colors[1], label="k-ϵ") +lines!(axϵ, ϵn[1], z, color=colors[1]) + +lines!(axb, bn[2], z, color=colors[2]) +lines!(axu, un[2], z, color=colors[2], label="u") +lines!(axu, vn[2], z, color=colors[2], linestyle=:dash, label="v") +lines!(axe, en[2], z, color=colors[2], label="CATKE") +lines!(axϵ, ϵn[2], z, color=colors[2]) + +#= +lines!(axb, bn[3], z, color=colors[3]) +lines!(axu, un[3], z, color=colors[3], label="u") +lines!(axu, vn[3], z, color=colors[3], linestyle=:dash, label="v") +lines!(axe, en[3], z, color=colors[3], label="k-ϵ (constant stability functions)") +lines!(axϵ, ϵn[3], z, color=colors[3]) +=# + +Legend(fig[0, 1:4], axe, nbanks=3, framevisible=false, tellheight=true) + +fig + diff --git a/validation/vertical_mixing_closures/k_epsilon_boundary_layer.jl b/validation/vertical_mixing_closures/k_epsilon_boundary_layer.jl new file mode 100644 index 0000000000..cb3506a0c2 --- /dev/null +++ b/validation/vertical_mixing_closures/k_epsilon_boundary_layer.jl @@ -0,0 +1,123 @@ +using Oceananigans +using Oceananigans.Units + +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: + TKEDissipationVerticalDiffusivity, + CATKEVerticalDiffusivity, + VariableStabilityFunctions, + minimum_stratification_number, + maximum_shear_number, + ConstantStabilityFunctions, + stratification_numberᶜᶜᶠ, + shear_numberᶜᶜᶠ, + momentum_stability_functionᶜᶜᶠ, + tracer_stability_functionᶜᶜᶠ + +using GLMakie +using Printf + +grid = RectilinearGrid(size=256, z=(-128, 0), topology=(Flat, Flat, Bounded)) + +f = 1e-4 +N² = 1e-5 +τˣ = -1e-4 +Jᵇ = 5e-8 + +u_top_bc = FluxBoundaryCondition(τˣ) +u_bcs = FieldBoundaryConditions(top=u_top_bc) + +b_top_bc = FluxBoundaryCondition(Jᵇ) +b_bcs = FieldBoundaryConditions(top=b_top_bc) + +coriolis = FPlane(; f) + +closure = TKEDissipationVerticalDiffusivity() +#closure = CATKEVerticalDiffusivity() + +model = HydrostaticFreeSurfaceModel(; grid, closure, coriolis, + tracers = (:b, :e, :ϵ), + buoyancy = BuoyancyTracer(), + boundary_conditions=(u=u_bcs, b=b_bcs)) + +bᵢ(z) = N² * z +set!(model, b=bᵢ) + +simulation = Simulation(model, Δt=60.0, stop_time=1day) + +u, v, w = model.velocities +e = model.tracers.e +ϵ = model.tracers.ϵ +b = model.tracers.b + +tracers = model.tracers +buoyancy = model.buoyancy +velocities = model.velocities + +κc = model.diffusivity_fields.κc +κe = model.diffusivity_fields.κe +κϵ = model.diffusivity_fields.κϵ + +progress(sim) = @info @sprintf("Iter: % 4d, time: % 24s, max(e): %6.2e, extrema(ϵ): (%6.2e, %6.2e)", + iteration(sim), prettytime(sim), maximum(e), minimum(ϵ), maximum(ϵ)) + +add_callback!(simulation, progress, IterationInterval(100)) + +run!(simulation) + +κcn = interior(κc, 1, 1, :) +bn = interior(b, 1, 1, :) +un = interior(u, 1, 1, :) +vn = interior(v, 1, 1, :) +en = interior(e, 1, 1, :) +ϵn = interior(ϵ, 1, 1, :) + +zc = znodes(model.tracers.e) +zf = znodes(κc) + +fig = Figure(size=(1600, 400)) + +axb = Axis(fig[1, 1], title="Buoyancy") +axu = Axis(fig[1, 2], title="Velocities") +axe = Axis(fig[1, 3], title="TKE") +axϵ = Axis(fig[1, 4], title="Epsilon") +axκ = Axis(fig[1, 5], title="Diffusivity") +axα = Axis(fig[1, 6], title="αᴺ, αᴹ") +axs = Axis(fig[1, 7], title="Stability functions") + +scatterlines!(axb, bn, zc) +lines!(axu, un, zc, label="u") +lines!(axu, vn, zc, label="v") +axislegend(axu, position=:rb) + +lines!(axe, en, zc, label="k-ϵ") +lines!(axϵ, ϵn, zc) +lines!(axκ, κcn, zf) + +𝕊ᵘ_op = KernelFunctionOperation{Center, Center, Face}(momentum_stability_functionᶜᶜᶠ, grid, closure, velocities, tracers, buoyancy) +𝕊ᶜ_op = KernelFunctionOperation{Center, Center, Face}(tracer_stability_functionᶜᶜᶠ, grid, closure, velocities, tracers, buoyancy) +𝕊ᵘ = Field(𝕊ᵘ_op) +𝕊ᶜ = Field(𝕊ᶜ_op) +compute!(𝕊ᵘ) +compute!(𝕊ᶜ) +𝕊ᵘn = interior(𝕊ᵘ, 1, 1, :) +𝕊ᶜn = interior(𝕊ᶜ, 1, 1, :) + +αᴺ_op = KernelFunctionOperation{Center, Center, Face}(stratification_numberᶜᶜᶠ, grid, closure, tracers, buoyancy) +αᴹ_op = KernelFunctionOperation{Center, Center, Face}(shear_numberᶜᶜᶠ, grid, closure, velocities, tracers, buoyancy) +αᴺ = Field(αᴺ_op) +αᴹ = Field(αᴹ_op) +compute!(αᴺ) +compute!(αᴹ) +αᴹn = interior(αᴹ, 1, 1, :) +αᴺn = interior(αᴺ, 1, 1, :) +scatterlines!(axα, αᴺn, zf, label="αᴺ") +scatterlines!(axα, αᴹn, zf, label="αᴹ") +axislegend(axα) +xlims!(axα, -30, 50) + +lines!(axs, 𝕊ᵘn, zf, label="𝕊ᵘ") +lines!(axs, 𝕊ᶜn, zf, label="𝕊ᶜ") +axislegend(axs, position=:rb) + +fig + diff --git a/validation/vertical_mixing_closures/look_at_stability_functions.jl b/validation/vertical_mixing_closures/look_at_stability_functions.jl new file mode 100644 index 0000000000..88f8afe025 --- /dev/null +++ b/validation/vertical_mixing_closures/look_at_stability_functions.jl @@ -0,0 +1,67 @@ +using Oceananigans + +using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: + TKEDissipationVerticalDiffusivity, + momentum_stability_function, + maximum_stratification_number, + minimum_stratification_number, + maximum_shear_number, + minimum_shear_number, + tracer_stability_function + +using GLMakie +using Printf + +closure = TKEDissipationVerticalDiffusivity() + +αᴺmin = minimum_stratification_number(closure) +αᴺmax = 30.0 +αᴺ = αᴺmin:0.01:αᴺmax + +αᴹmin = minimum_shear_number(closure) +αᴹmax = maximum_shear_number.(Ref(closure), αᴺ) +αᴹ = αᴹmin:0.01:maximum(αᴹmax) + +fig = Figure(size=(1600, 600)) +ax1 = Axis(fig[1, 1], title="Stability functions", xlabel="αᴺ", ylabel="𝕊c") +ax2 = Axis(fig[2, 1], title="Prandtl number", xlabel="αᴺ", ylabel="𝕊c") +ax3 = Axis(fig[1:2, 2], title="Tracer stability functions", xlabel="αᴺ", ylabel="αᴹ") +ax4 = Axis(fig[1:2, 3], title="Momentum stability functions", xlabel="αᴺ", ylabel="αᴹ") +ax5 = Axis(fig[1:2, 4], title="Prandtl number", xlabel="αᴺ", ylabel="αᴹ") + +𝕊c_max_αᴹ = tracer_stability_function.(Ref(closure), αᴺ, αᴹmax) +𝕊u_max_αᴹ = momentum_stability_function.(Ref(closure), αᴺ, αᴹmax) + +𝕊c_min_αᴹ = tracer_stability_function.(Ref(closure), αᴺ, αᴹmin) +𝕊u_min_αᴹ = momentum_stability_function.(Ref(closure), αᴺ, αᴹmin) + +NN = length(αᴺ) +NM = length(αᴹ) +𝕊c = tracer_stability_function.(Ref(closure), reshape(αᴺ, NN, 1), reshape(αᴹ, 1, NM)) +𝕊u = momentum_stability_function.(Ref(closure), reshape(αᴺ, NN, 1), reshape(αᴹ, 1, NM)) +Pr = 𝕊u ./ 𝕊c + +Pr_max = maximum(Pr, dims=2)[:] +Pr_min = minimum(Pr, dims=2)[:] + +lines!(ax1, αᴺ, 𝕊c_max_αᴹ, label="max(αᴹ)", color=:blue) +lines!(ax1, αᴺ, 𝕊c_min_αᴹ, label="min(αᴹ)", linestyle=:dash, color=:blue) + +lines!(ax1, αᴺ, 𝕊u_max_αᴹ, color=:red) +lines!(ax1, αᴺ, 𝕊u_min_αᴹ, linestyle=:dash, color=:red) + +band!(ax1, αᴺ, 𝕊c_min_αᴹ, 𝕊c_max_αᴹ, label="𝕊c", color=(:blue, 0.5)) +band!(ax1, αᴺ, 𝕊u_min_αᴹ, 𝕊u_max_αᴹ, label="𝕊u", color=(:red, 0.5)) +axislegend(ax1) + +band!(ax2, αᴺ, Pr_min, Pr_max, label="Pr", color=(:blue, 0.5)) + +cf = contourf!(ax3, αᴺ, αᴹ, 𝕊c, levels=0.01:0.03:0.2, colorrrange=(0.01, 0.2)) +cf = contourf!(ax4, αᴺ, αᴹ, 𝕊u, levels=0.01:0.03:0.2, colorrrange=(0.01, 0.2)) +Colorbar(fig[3, 2:3], cf, vertical=false, tellwidth=false, label="Stability functions", flipaxis=false) + +cf = contourf!(ax6, αᴺ, αᴹ, Pr, levels=0.3:0.1:3.0, colorrrange=(0.35, 2.8), colormap=:solar) +Colorbar(fig[3, 4], cf, vertical=false, tellwidth=false, label="Prandtl number", flipaxis=false) + +display(fig) +