diff --git a/src/DistributedComputations/halo_communication.jl b/src/DistributedComputations/halo_communication.jl index 06264dd39c..2ddcaae60f 100644 --- a/src/DistributedComputations/halo_communication.jl +++ b/src/DistributedComputations/halo_communication.jl @@ -242,6 +242,8 @@ for (side, opposite_side) in zip([:west, :south], [:east, :north]) function $fill_both_halo!(c, bc_side::DCBCT, bc_opposite_side::DCBCT, size, offset, loc, arch::Distributed, grid::DistributedGrid, buffers, args...; only_local_halos = false, kwargs...) + sync_device!(arch) + only_local_halos && return nothing sync_device!(arch) @@ -273,6 +275,8 @@ for side in [:west, :east, :south, :north] function $fill_side_halo!(c, bc_side::DCBCT, size, offset, loc, arch::Distributed, grid::DistributedGrid, buffers, args...; only_local_halos = false, kwargs...) + sync_device!(arch) + only_local_halos && return nothing sync_device!(arch) diff --git a/src/ImmersedBoundaries/grid_fitted_bottom.jl b/src/ImmersedBoundaries/grid_fitted_bottom.jl index e693e20e2d..f1ca0f1e4b 100644 --- a/src/ImmersedBoundaries/grid_fitted_bottom.jl +++ b/src/ImmersedBoundaries/grid_fitted_bottom.jl @@ -72,12 +72,24 @@ Computes ib.bottom_height and wraps in an array. function ImmersedBoundaryGrid(grid, ib::GridFittedBottom) bottom_field = Field{Center, Center, Nothing}(grid) set!(bottom_field, ib.bottom_height) + + # Make sure that `abs(bottom_height) <= grid.Lz` to constrain the bottom + @apply_regionally launch!(architecture(grid), grid, :xy, _limit_bottom_height!, bottom_field, grid.Lz) + fill_halo_regions!(bottom_field) new_ib = GridFittedBottom(bottom_field, ib.immersed_condition) TX, TY, TZ = topology(grid) return ImmersedBoundaryGrid{TX, TY, TZ}(grid, new_ib) end +# Make sure that `abs(bottom_height) <= grid.Lz` to constrain the bottom +@kernel function _limit_bottom_height!(bottom_field, Lz) + i, j = @index(Global, NTuple) + if abs(bottom_field[i, j, 1]) > Lz + bottom_field[i, j, 1] = sign(bottom_field[i, j, 1]) * Lz + end +end + @inline function _immersed_cell(i, j, k, underlying_grid, ib::GridFittedBottom{<:Any, <:InterfaceImmersedCondition}) z = znode(i, j, k+1, underlying_grid, c, c, f) h = @inbounds ib.bottom_height[i, j, 1] diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 717bb1fd2d..10363d7ac7 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl @@ -43,7 +43,7 @@ function ab2_step_velocities!(velocities, model, Δt, χ) velocity_field = model.velocities[name] launch!(model.architecture, model.grid, :xyz, - ab2_step_field!, velocity_field, Δt, χ, Gⁿ, G⁻) + ab2_step_field!, velocity_field, Δt, χ, Gⁿ, G⁻, model.grid) # TODO: let next implicit solve depend on previous solve + explicit velocity step # Need to distinguish between solver events and tendency calculation events. @@ -78,7 +78,7 @@ function ab2_step_tracers!(tracers, model, Δt, χ) closure = model.closure launch!(model.architecture, model.grid, :xyz, - ab2_step_field!, tracer_field, Δt, χ, Gⁿ, G⁻) + ab2_step_field!, tracer_field, Δt, χ, Gⁿ, G⁻, model.grid) implicit_step!(tracer_field, model.timestepper.implicit_solver, diff --git a/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl b/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl index d32b5786d8..525eb45995 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl @@ -213,7 +213,6 @@ function SplitExplicitAuxiliaryFields(grid::AbstractGrid) fill_halo_regions!((Hᶠᶜ, Hᶜᶠ)) kernel_parameters = :xy - return SplitExplicitAuxiliaryFields(Gᵁ, Gⱽ, Hᶠᶜ, Hᶜᶠ, kernel_parameters) end diff --git a/src/TimeSteppers/quasi_adams_bashforth_2.jl b/src/TimeSteppers/quasi_adams_bashforth_2.jl index 780df779b6..5e26927234 100644 --- a/src/TimeSteppers/quasi_adams_bashforth_2.jl +++ b/src/TimeSteppers/quasi_adams_bashforth_2.jl @@ -125,7 +125,8 @@ function ab2_step!(model, Δt, χ) step_field_kernel!(field, Δt, χ, model.timestepper.Gⁿ[i], - model.timestepper.G⁻[i]) + model.timestepper.G⁻[i], + model.grid) # TODO: function tracer_index(model, field_index) = field_index - 3, etc... tracer_index = Val(i - 3) # assumption @@ -148,7 +149,7 @@ Time step velocity fields via the 2nd-order quasi Adams-Bashforth method `U^{n+1} = U^n + Δt ((3/2 + χ) * G^{n} - (1/2 + χ) G^{n-1})` """ -@kernel function ab2_step_field!(u, Δt, χ, Gⁿ, G⁻) +@kernel function ab2_step_field!(u, Δt, χ, Gⁿ, G⁻, grid) i, j, k = @index(Global, NTuple) FT = eltype(χ) @@ -158,4 +159,4 @@ Time step velocity fields via the 2nd-order quasi Adams-Bashforth method @inbounds u[i, j, k] += convert(FT, Δt) * ((one_point_five + χ) * Gⁿ[i, j, k] - (oh_point_five + χ) * G⁻[i, j, k]) end -@kernel ab2_step_field!(::FunctionField, Δt, χ, Gⁿ, G⁻) = nothing +@kernel ab2_step_field!(::FunctionField, Δt, χ, Gⁿ, G⁻, grid) = nothing diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/CATKEVerticalDiffusivities.jl b/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/CATKEVerticalDiffusivities.jl index 38aaa5bd12..7d0455557c 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/CATKEVerticalDiffusivities.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/CATKEVerticalDiffusivities.jl @@ -93,6 +93,7 @@ include("turbulent_kinetic_energy_equation.jl") maximum_viscosity = Inf, minimum_turbulent_kinetic_energy = 1e-6, minimum_convective_buoyancy_flux = 1e-8, + minimum_dissipation_length_scale = 1e-2, negative_turbulent_kinetic_energy_damping_time_scale = 1minute) Return the `CATKEVerticalDiffusivity` turbulence closure for vertical mixing by @@ -226,7 +227,11 @@ function DiffusivityFields(grid, tracer_names, bcs, closure::FlavorOfCATKE) _tupled_tracer_diffusivities = NamedTuple(name => name === :e ? κᵉ : κᶜ for name in tracer_names) _tupled_implicit_linear_coefficients = NamedTuple(name => name === :e ? Lᵉ : ZeroField() for name in tracer_names) - return (; κᵘ, κᶜ, κᵉ, Lᵉ, Qᵇ, previous_compute_time, _tupled_tracer_diffusivities, _tupled_implicit_linear_coefficients) + S² = ZFaceField(grid) + N² = ZFaceField(grid) + w★ = CenterField(grid) + + return (; κᵘ, κᶜ, κᵉ, Lᵉ, Qᵇ, S², N², w★, previous_compute_time, _tupled_tracer_diffusivities, _tupled_implicit_linear_coefficients) end const c = Center() @@ -249,47 +254,76 @@ function compute_diffusivities!(diffusivities, closure::FlavorOfCATKE, model; pa Δt = model.clock.time - diffusivities.previous_compute_time[] diffusivities.previous_compute_time[] = model.clock.time + launch!(arch, grid, parameters, + _compute_CATKE_auxiliaries!, + diffusivities, grid, closure, velocities, tracers, buoyancy) + launch!(arch, grid, :xy, - compute_average_surface_buoyancy_flux!, - diffusivities.Qᵇ, grid, closure, velocities, tracers, buoyancy, top_tracer_bcs, clock, Δt) + _compute_average_surface_buoyancy_flux!, + diffusivities.Qᵇ, grid, closure, diffusivities, velocities, tracers, buoyancy, top_tracer_bcs, clock, Δt) launch!(arch, grid, parameters, - compute_CATKE_diffusivities!, + _compute_CATKE_diffusivities!, diffusivities, grid, closure, velocities, tracers, buoyancy) return nothing end -@kernel function compute_average_surface_buoyancy_flux!(Qᵇ, grid, closure, velocities, tracers, buoyancy, top_tracer_bcs, clock, Δt) +@kernel function _compute_CATKE_auxiliaries!(diffusivities, grid, closure, velocities, tracers, buoyancy) + i, j, k = @index(Global, NTuple) + + e = tracers.e + S² = diffusivities.S² + N² = diffusivities.N² + w★ = diffusivities.w★ + u, v, _ = velocities + + @inbounds begin + S²[i, j, k] = shearᶜᶜᶠ(i, j, k, grid, u, v) + N²[i, j, k] = ∂z_b(i, j, k, grid, buoyancy, tracers) + w★[i, j, k] = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, tracers.e) + end +end + +@kernel function _compute_average_surface_buoyancy_flux!(Qᵇ, grid, closure, diffusivities, velocities, tracers, buoyancy, top_tracer_bcs, clock, Δt) i, j = @index(Global, NTuple) + S² = diffusivities.S² + N² = diffusivities.N² + w★ = diffusivities.w★ + closure = getclosure(i, j, closure) Qᵇ★ = top_buoyancy_flux(i, j, grid, buoyancy, top_tracer_bcs, clock, merge(velocities, tracers)) k = grid.Nz - ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, Qᵇ) + ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Qᵇ, S², N², w★) Qᵇᵋ = closure.minimum_convective_buoyancy_flux Qᵇᵢⱼ = @inbounds Qᵇ[i, j, 1] Qᵇ⁺ = max(Qᵇᵋ, Qᵇᵢⱼ, Qᵇ★) # selects fastest (dominant) time-scale - t★ = (ℓᴰ^2 / Qᵇ⁺)^(1/3) + t★ = (ℓᴰ^2 / Qᵇ⁺)^(1/3) ϵ = Δt / t★ - @inbounds Qᵇ[i, j, 1] = (Qᵇᵢⱼ + ϵ * Qᵇ★) / (1 + ϵ) + Qᵇ⁺ᵢⱼ = (Qᵇᵢⱼ + ϵ * Qᵇ★) / (1 + ϵ) + + @inbounds Qᵇ[i, j, 1] = ifelse(ϵ == Inf, Qᵇ★, Qᵇ⁺ᵢⱼ) # avoid problems when `ϵ == Inf` that leads to `NaN` end -@kernel function compute_CATKE_diffusivities!(diffusivities, grid, closure::FlavorOfCATKE, velocities, tracers, buoyancy) +@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) Qᵇ = diffusivities.Qᵇ + S² = diffusivities.S² + N² = diffusivities.N² + w★ = diffusivities.w★ @inbounds begin - κᵘ★ = κuᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, Qᵇ) - κᶜ★ = κcᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, Qᵇ) - κᵉ★ = κeᶜᶜᶠ(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, Qᵇ) + κᵘ★ = κuᶜᶜᶠ(i, j, k, grid, closure_ij, Qᵇ, S², N², w★) + κᶜ★ = κcᶜᶜᶠ(i, j, k, grid, closure_ij, Qᵇ, S², N², w★) + κᵉ★ = κeᶜᶜᶠ(i, j, k, grid, closure_ij, Qᵇ, S², N², w★) on_periphery = peripheral_node(i, j, k, grid, c, c, f) within_inactive = inactive_node(i, j, k, grid, c, c, f) @@ -304,7 +338,7 @@ end # "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. - wb = explicit_buoyancy_flux(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, diffusivities) + wb = explicit_buoyancy_flux(i, j, k, grid, closure, diffusivities) eⁱʲᵏ = @inbounds tracers.e[i, j, k] eᵐⁱⁿ = closure_ij.minimum_turbulent_kinetic_energy @@ -333,12 +367,10 @@ 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ⁱʲᵏ) - w★ = sqrt(e⁺) - div_Jᵉ_e = - on_bottom * Cᵂϵ * w★ / Δz + div_Jᵉ_e = - on_bottom * Cᵂϵ * w[i, j, k] / Δz # Implicit TKE dissipation - ω = dissipation_rate(i, j, k, grid, closure_ij, velocities, tracers, buoyancy, diffusivities) + ω_e = dissipation_rate(i, j, k, grid, closure_ij, tracers, diffusivities) # The interior contributions to the linear implicit term `L` are defined via # @@ -353,7 +385,7 @@ end # # where ω = ϵ / e ∼ √e / ℓ. - diffusivities.Lᵉ[i, j, k] = wb_e - ω + div_Jᵉ_e + diffusivities.Lᵉ[i, j, k] = wb_e - ω_e + div_Jᵉ_e end end @@ -368,42 +400,40 @@ end return sqrt(max(eᵐⁱⁿ, eᵢ)) 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) - ℓᵘ = momentum_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) - κᵘ = ℓᵘ * w★ +@inline function κuᶜᶜᶠ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + w★ᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, w★) + ℓᵘ = momentum_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, S², N², w★) + κᵘ = ℓᵘ * w★ᶜᶜᶠ κᵘ_max = closure.maximum_viscosity return min(κᵘ, κᵘ_max) end -@inline function κuᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) - w★ = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, tracers.e) - ℓᵘ = momentum_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) - κᵘ = ℓᵘ * w★ +@inline function κuᶜᶜᶜ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + ℓᵘ = momentum_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, S², N², w★) + κᵘ = @inbounds ℓᵘ * w★[i, j, k] κᵘ_max = closure.maximum_viscosity return min(κᵘ, κᵘ_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) - ℓᶜ = tracer_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) - κᶜ = ℓᶜ * w★ +@inline function κcᶜᶜᶠ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + w★ᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, w★) + ℓᶜ = tracer_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + κᶜ = ℓᶜ * w★ᶜᶜᶠ κᶜ_max = closure.maximum_tracer_diffusivity return min(κᶜ, κᶜ_max) end -@inline function κcᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) - w★ = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, tracers.e) - ℓᶜ = tracer_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) - κᶜ = ℓᶜ * w★ +@inline function κcᶜᶜᶜ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + ℓᶜ = tracer_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + κᶜ = @inbounds ℓᶜ * w★[i, j, k] κᶜ_max = closure.maximum_tracer_diffusivity return min(κᶜ, κᶜ_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) - ℓᵉ = TKE_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) - κᵉ = ℓᵉ * w★ +@inline function κeᶜᶜᶠ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + w★ᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, w★) + ℓᵉ = TKE_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) + κᵉ = ℓᵉ * w★ᶜᶜᶠ κᵉ_max = closure.maximum_tke_diffusivity return min(κᵉ, κᵉ_max) end diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/mixing_length.jl b/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/mixing_length.jl index 1267a025a7..1a92119672 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/mixing_length.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/mixing_length.jl @@ -50,23 +50,22 @@ end return S² end -@inline function stratification_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, e, tracers, buoyancy) +@inline function stratification_mixing_lengthᶜᶜᶠ(i, j, k, grid, N², w★) FT = eltype(grid) - N² = ∂z_b(i, j, k, grid, buoyancy, tracers) - N²⁺ = clip(N²) - w★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocityᶜᶜᶜ, closure, e) - return ifelse(N²⁺ == 0, FT(Inf), w★ / sqrt(N²⁺)) + @inbounds N²⁺ = clip(N²[i, j, k]) + w★ᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, w★) + return ifelse(N²⁺ == 0, FT(Inf), w★ᶜᶜᶠ / sqrt(N²⁺)) end -@inline function stratification_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy) +@inline function stratification_mixing_lengthᶜᶜᶜ(i, j, k, grid, N², w★) FT = eltype(grid) - N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers) - N²⁺ = clip(N²) - w★ = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, e) - return ifelse(N²⁺ == 0, FT(Inf), w★ / sqrt(N²⁺)) + N²ᶜᶜᶜ = ℑzᵃᵃᶜ(i, j, k, grid, N²) + N²⁺ = clip(N²ᶜᶜᶜ) + w★ᶜᶜᶜ = @inbounds w★[i, j, k] + return ifelse(N²⁺ == 0, FT(Inf), w★ᶜᶜᶜ / sqrt(N²⁺)) end -@inline function stable_length_scaleᶜᶜᶠ(i, j, k, grid, closure, e, velocities, tracers, buoyancy) +@inline function stable_length_scaleᶜᶜᶠ(i, j, k, grid, closure, N², w★) Cˢ = closure.mixing_length.Cˢ Cᵇ = closure.mixing_length.Cᵇ @@ -74,7 +73,7 @@ end d_down = Cᵇ * height_above_bottomᶜᶜᶠ(i, j, k, grid) d = min(d_up, d_down) - ℓᴺ = stratification_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, e, tracers, buoyancy) + ℓᴺ = stratification_mixing_lengthᶜᶜᶠ(i, j, k, grid, N², w★) ℓ = min(d, ℓᴺ) ℓ = ifelse(isnan(ℓ), d, ℓ) @@ -82,14 +81,14 @@ end return ℓ end -@inline function stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, e, velocities, tracers, buoyancy) +@inline function stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, N², w★) Cˢ = closure.mixing_length.Cˢ Cᵇ = closure.mixing_length.Cᵇ d_up = Cˢ * depthᶜᶜᶜ(i, j, k, grid) d_down = Cᵇ * height_above_bottomᶜᶜᶜ(i, j, k, grid) d = min(d_up, d_down) - ℓᴺ = stratification_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, e, tracers, buoyancy) + ℓᴺ = stratification_mixing_lengthᶜᶜᶜ(i, j, k, grid, N², w★) ℓ = min(d, ℓᴺ) ℓ = ifelse(isnan(ℓ), d, ℓ) @@ -97,34 +96,31 @@ end return ℓ end -@inline three_halves_tkeᶜᶜᶜ(i, j, k, grid, closure, e) = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, e)^3 -@inline squared_tkeᶜᶜᶜ(i, j, k, grid, closure, e) = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, e)^2 +@inline three_halves_tkeᶜᶜᶜ(i, j, k, grid, w★) = @inbounds w★[i, j, k]^3 +@inline squared_tkeᶜᶜᶜ(i, j, k, grid, w★) = @inbounds w★[i, j, k]^2 @inline function convective_length_scaleᶜᶜᶠ(i, j, k, grid, closure, Cᶜ::Number, Cᵉ::Number, Cˢᵖ::Number, - velocities, tracers, buoyancy, surface_buoyancy_flux) - - u = velocities.u - v = velocities.v + surface_buoyancy_flux, S², N², w★) Qᵇᵋ = closure.minimum_convective_buoyancy_flux Qᵇ = @inbounds surface_buoyancy_flux[i, j, 1] - w★ = ℑzᵃᵃᶠ(i, j, k, grid, turbulent_velocityᶜᶜᶜ, closure, tracers.e) - w★² = ℑzᵃᵃᶠ(i, j, k, grid, squared_tkeᶜᶜᶜ, closure, tracers.e) - w★³ = ℑzᵃᵃᶠ(i, j, k, grid, three_halves_tkeᶜᶜᶜ, closure, tracers.e) - S² = shearᶜᶜᶠ(i, j, k, grid, u, v) - N² = ∂z_b(i, j, k, grid, buoyancy, tracers) - N²_above = ∂z_b(i, j, k+1, grid, buoyancy, tracers) + w★¹ᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, w★) + w★²ᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, squared_tkeᶜᶜᶜ, w★) + w★³ᶜᶜᶠ = ℑzᵃᵃᶠ(i, j, k, grid, three_halves_tkeᶜᶜᶜ, w★) + S²ᶜᶜᶠ = @inbounds S²[i, j, k] + N²_local = @inbounds N²[i, j, k] + N²_above = @inbounds N²[i, j, k+1] # "Convective length" # ℓᶜ ∼ boundary layer depth according to Deardorff scaling - ℓᶜ = Cᶜ * w★³ / (Qᵇ + Qᵇᵋ) + ℓᶜ = Cᶜ * w★³ᶜᶜᶠ / (Qᵇ + Qᵇᵋ) ℓᶜ = ifelse(isnan(ℓᶜ), zero(grid), ℓᶜ) # Figure out which mixing length applies - convecting = (Qᵇ > Qᵇᵋ) & (N² < 0) + convecting = (Qᵇ > Qᵇᵋ) & (N²_local < 0) # Model for shear-convection interaction - Sp = sqrt(S²) * w★² / (Qᵇ + Qᵇᵋ) # Sp = "Sheared convection number" + Sp = sqrt(S²ᶜᶜᶠ) * w★²ᶜᶜᶠ / (Qᵇ + Qᵇᵋ) # Sp = "Sheared convection number" ϵˢᵖ = 1 - Cˢᵖ * Sp # ϵ = Sheared convection factor # Reduce convective and entraining mixing lengths by sheared convection factor @@ -133,10 +129,10 @@ end # "Entrainment length" # Ensures that w′b′ ~ Qᵇ at entrainment depth - ℓᵉ = Cᵉ * Qᵇ / (w★ * N² + Qᵇᵋ) + ℓᵉ = Cᵉ * Qᵇ / (w★¹ᶜᶜᶠ * N²_local + Qᵇᵋ) ℓᵉ = clip(ϵˢᵖ * ℓᵉ) - entraining = (Qᵇ > Qᵇᵋ) & (N² > 0) & (N²_above < 0) + entraining = (Qᵇ > Qᵇᵋ) & (N²_local > 0) & (N²_above < 0) ℓ = ifelse(convecting, ℓᶜ, ifelse(entraining, ℓᵉ, zero(grid))) @@ -145,19 +141,16 @@ end end @inline function convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ::Number, Cᵉ::Number, Cˢᵖ::Number, - velocities, tracers, buoyancy, surface_buoyancy_flux) - - u = velocities.u - v = velocities.v + surface_buoyancy_flux, S², N², w★) Qᵇᵋ = closure.minimum_convective_buoyancy_flux Qᵇ = @inbounds surface_buoyancy_flux[i, j, 1] - w★ = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, tracers.e) - w★² = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, tracers.e)^2 - w★³ = turbulent_velocityᶜᶜᶜ(i, j, k, grid, closure, tracers.e)^3 - S² = shearᶜᶜᶜ(i, j, k, grid, u, v) - N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers) - N²_above = ℑzᵃᵃᶜ(i, j, k+1, grid, ∂z_b, buoyancy, tracers) + w★¹ = @inbounds w★[i, j, k] + w★² = @inbounds w★[i, j, k]^2 + w★³ = @inbounds w★[i, j, k]^3 + S²ᶜᶜᶜ = ℑzᵃᵃᶜ(i, j, k, grid, S²) + N²_local = ℑzᵃᵃᶜ(i, j, k, grid, N²) + N²_above = ℑzᵃᵃᶜ(i, j, k+1, grid, N²) # "Convective length" # ℓᶜ ∼ boundary layer depth according to Deardorff scaling @@ -165,10 +158,10 @@ end ℓᶜ = ifelse(isnan(ℓᶜ), zero(grid), ℓᶜ) # Figure out which mixing length applies - convecting = (Qᵇ > Qᵇᵋ) & (N² < 0) + convecting = (Qᵇ > Qᵇᵋ) & (N²_local < 0) # Model for shear-convection interaction - Sp = sqrt(S²) * w★² / (Qᵇ + Qᵇᵋ) # Sp = "Sheared convection number" + Sp = sqrt(S²ᶜᶜᶜ) * w★² / (Qᵇ + Qᵇᵋ) # Sp = "Sheared convection number" ϵˢᵖ = 1 - Cˢᵖ * Sp # ϵ = Sheared convection factor # Reduce convective and entraining mixing lengths by sheared convection factor @@ -177,10 +170,10 @@ end # "Entrainment length" # Ensures that w′b′ ~ Qᵇ at entrainment depth - ℓᵉ = Cᵉ * Qᵇ / (w★ * N² + Qᵇᵋ) + ℓᵉ = Cᵉ * Qᵇ / (w★¹ * N²_local + Qᵇᵋ) ℓᵉ = clip(ϵˢᵖ * ℓᵉ) - entraining = (Qᵇ > Qᵇᵋ) & (N² > 0) & (N²_above < 0) + entraining = (Qᵇ > Qᵇᵋ) & (N²_local > 0) & (N²_above < 0) ℓ = ifelse(convecting, ℓᶜ, ifelse(entraining, ℓᵉ, zero(grid))) @@ -192,26 +185,30 @@ end @inline step(x, c, w) = max(zero(x), min(one(x), (x - c) / w)) @inline scale(Ri, σ⁻, σ⁺ , c, w) = σ⁻ + (σ⁺ - σ⁻) * step(Ri, c, w) -@inline function stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) - Ri = Riᶜᶜᶠ(i, j, k, grid, velocities, tracers, buoyancy) +@inline function stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) + N²ᶜᶜᶠ = @inbounds N²[i, j, k] + S²ᶜᶜᶠ = @inbounds S²[i, j, k] + Ri = ifelse(N²ᶜᶜᶠ ≤ 0, zero(grid), N²ᶜᶜᶠ / S²ᶜᶜᶠ) CRi⁰ = closure.mixing_length.CRi⁰ CRiᵟ = closure.mixing_length.CRiᵟ return scale(Ri, Cˡᵒ, Cʰⁱ, CRi⁰, CRiᵟ) end -@inline function stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) - Ri = Riᶜᶜᶜ(i, j, k, grid, velocities, tracers, buoyancy) +@inline function stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) + N²ᶜᶜᶜ = ℑzᵃᵃᶠ(i, j, k, grid, N²) + S²ᶜᶜᶜ = ℑzᵃᵃᶠ(i, j, k, grid, S²) + Ri = ifelse(N²ᶜᶜᶜ ≤ 0, zero(grid), N²ᶜᶜᶜ / S²ᶜᶜᶜ) CRi⁰ = closure.mixing_length.CRi⁰ CRiᵟ = closure.mixing_length.CRiᵟ return scale(Ri, Cˡᵒ, Cʰⁱ, CRi⁰, CRiᵟ) end -@inline function momentum_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) +@inline function momentum_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, S², N², w★) Cˡᵒ = closure.mixing_length.Cˡᵒu Cʰⁱ = closure.mixing_length.Cʰⁱu - σ = stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) + σ = stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) - ℓ★ = σ * stable_length_scaleᶜᶜᶠ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy) + ℓ★ = σ * stable_length_scaleᶜᶜᶠ(i, j, k, grid, closure, N², w★) ℓ★ = ifelse(isnan(ℓ★), zero(grid), ℓ★) H = total_depthᶜᶜᵃ(i, j, grid) @@ -219,12 +216,12 @@ end return min(H, ℓ★) end -@inline function momentum_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) +@inline function momentum_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, S², N², w★) Cˡᵒ = closure.mixing_length.Cˡᵒu Cʰⁱ = closure.mixing_length.Cʰⁱu - σ = stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) + σ = stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) - ℓ★ = σ * stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy) + ℓ★ = σ * stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, N², w★) ℓ★ = ifelse(isnan(ℓ★), zero(grid), ℓ★) H = total_depthᶜᶜᵃ(i, j, grid) @@ -232,16 +229,16 @@ end return min(H, ℓ★) end -@inline function tracer_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) +@inline function tracer_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) Cᶜ = closure.mixing_length.Cᶜc Cᵉ = closure.mixing_length.Cᵉc Cˢᵖ = closure.mixing_length.Cˢᵖ - ℓʰ = convective_length_scaleᶜᶜᶠ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, velocities, tracers, buoyancy, surface_buoyancy_flux) + ℓʰ = convective_length_scaleᶜᶜᶠ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, surface_buoyancy_flux, S², N², w★) Cˡᵒ = closure.mixing_length.Cˡᵒc Cʰⁱ = closure.mixing_length.Cʰⁱc - σ = stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) - ℓ★ = σ * stable_length_scaleᶜᶜᶠ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy) + σ = stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) + ℓ★ = σ * stable_length_scaleᶜᶜᶠ(i, j, k, grid, closure, N², w★) ℓʰ = ifelse(isnan(ℓʰ), zero(grid), ℓʰ) ℓ★ = ifelse(isnan(ℓ★), zero(grid), ℓ★) @@ -251,16 +248,16 @@ end return min(H, ℓᶜ) end -@inline function tracer_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) +@inline function tracer_mixing_lengthᶜᶜᶜ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) Cᶜ = closure.mixing_length.Cᶜc Cᵉ = closure.mixing_length.Cᵉc Cˢᵖ = closure.mixing_length.Cˢᵖ - ℓʰ = convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, velocities, tracers, buoyancy, surface_buoyancy_flux) + ℓʰ = convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, surface_buoyancy_flux, S², N², w★) Cˡᵒ = closure.mixing_length.Cˡᵒc Cʰⁱ = closure.mixing_length.Cʰⁱc - σ = stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) - ℓ★ = σ * stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy) + σ = stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) + ℓ★ = σ * stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, N², w★) ℓʰ = ifelse(isnan(ℓʰ), zero(grid), ℓʰ) ℓ★ = ifelse(isnan(ℓ★), zero(grid), ℓ★) @@ -271,16 +268,16 @@ end return min(H, ℓᶜ) end -@inline function TKE_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, velocities, tracers, buoyancy, surface_buoyancy_flux) +@inline function TKE_mixing_lengthᶜᶜᶠ(i, j, k, grid, closure, surface_buoyancy_flux, S², N², w★) Cᶜ = closure.mixing_length.Cᶜe Cᵉ = closure.mixing_length.Cᵉe Cˢᵖ = closure.mixing_length.Cˢᵖ - ℓʰ = convective_length_scaleᶜᶜᶠ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, velocities, tracers, buoyancy, surface_buoyancy_flux) + ℓʰ = convective_length_scaleᶜᶜᶠ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, surface_buoyancy_flux, S², N², w★) Cˡᵒ = closure.mixing_length.Cˡᵒe Cʰⁱ = closure.mixing_length.Cʰⁱe - σ = stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, velocities, tracers, buoyancy) - ℓ★ = σ * stable_length_scaleᶜᶜᶠ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy) + σ = stability_functionᶜᶜᶠ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) + ℓ★ = σ * stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, N², w★) ℓʰ = ifelse(isnan(ℓʰ), zero(grid), ℓʰ) ℓ★ = ifelse(isnan(ℓ★), zero(grid), ℓ★) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/turbulent_kinetic_energy_equation.jl b/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/turbulent_kinetic_energy_equation.jl index 58d82af1af..f7adf81b33 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/turbulent_kinetic_energy_equation.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/CATKEVerticalDiffusivities/turbulent_kinetic_energy_equation.jl @@ -33,8 +33,8 @@ end # Non-conservative reconstruction of shear production: closure = getclosure(i, j, closure) - κᵘ = κuᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities.Qᵇ) - S² = shearᶜᶜᶜ(i, j, k, grid, u, v) + κᵘ = @inbounds diffusivities.κᵘ[i, j, k] + S² = ℑzᵃᵃᶜ(i, j, k, grid, diffusivities.S²) return κᵘ * S² end @@ -51,20 +51,20 @@ end # ℑzᵃᵃᶜ(i, j, k, grid, buoyancy_fluxᶜᶜᶠ, tracers, buoyancy, diffusivities) # Non-conservative reconstruction of buoyancy flux: -@inline function explicit_buoyancy_flux(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities) +@inline function explicit_buoyancy_flux(i, j, k, grid, closure, diffusivities) closure = getclosure(i, j, closure) - κᶜ = κcᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities.Qᵇ) - N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers) + κᶜ = ℑzᵃᵃᶜ(i, j, k, grid, diffusivities.κᶜ) + N² = ℑzᵃᵃᶜ(i, j, k, grid, diffusivities.N²) return - κᶜ * N² end @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) + explicit_buoyancy_flux(i, j, k, grid, closure, diffusivities) const VITD = VerticallyImplicitTimeDiscretization @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) + wb = explicit_buoyancy_flux(i, j, k, grid, closure, diffusivities) eⁱʲᵏ = @inbounds tracers.e[i, j, k] dissipative_buoyancy_flux = sign(wb) * sign(eⁱʲᵏ) < 0 @@ -75,23 +75,20 @@ const VITD = VerticallyImplicitTimeDiscretization return ifelse(dissipative_buoyancy_flux, 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) +@inline function dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure::FlavorOfCATKE, surface_buoyancy_flux, S², N², w★) # Convective dissipation length Cᶜ = closure.turbulent_kinetic_energy_equation.CᶜD Cᵉ = closure.turbulent_kinetic_energy_equation.CᵉD Cˢᵖ = closure.mixing_length.Cˢᵖ Qᵇ = surface_buoyancy_flux - ℓʰ = convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, velocities, tracers, buoyancy, Qᵇ) + ℓʰ = convective_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Cᶜ, Cᵉ, Cˢᵖ, Qᵇ, S², N², w★) # "Stable" dissipation length 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ʰⁱ, velocities, tracers, buoyancy) - ℓ★ = stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, tracers.e, velocities, tracers, buoyancy) + σᴰ = stability_functionᶜᶜᶜ(i, j, k, grid, closure, Cˡᵒ, Cʰⁱ, S², N²) + ℓ★ = stable_length_scaleᶜᶜᶜ(i, j, k, grid, closure, N², w★) ℓ★ = ℓ★ / σᴰ # Dissipation length @@ -103,12 +100,15 @@ end return min(H, ℓᴰ) end -@inline function dissipation_rate(i, j, k, grid, closure::FlavorOfCATKE, - velocities, tracers, buoyancy, diffusivities) +@inline function dissipation_rate(i, j, k, grid, closure::FlavorOfCATKE, tracers, diffusivities) + + Qᵇ = diffusivities.Qᵇ + S² = diffusivities.S² + N² = diffusivities.N² + w★ = diffusivities.w★ - ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, velocities, tracers, buoyancy, diffusivities.Qᵇ) + ℓᴰ = dissipation_length_scaleᶜᶜᶜ(i, j, k, grid, closure, Qᵇ, S², N², w★) e = tracers.e - FT = eltype(grid) eᵢ = @inbounds e[i, j, k] # Note: @@ -132,6 +132,8 @@ end return ω * eᵢ end +dissipation(i, j, k, grid, closure::FlavorOfCATKE{<:VITD}, velocities, tracers, args...) = zero(grid) + ##### ##### For closure tuples... ##### diff --git a/src/Utils/kernel_launching.jl b/src/Utils/kernel_launching.jl index b0d3338669..dd0348abe5 100644 --- a/src/Utils/kernel_launching.jl +++ b/src/Utils/kernel_launching.jl @@ -167,9 +167,9 @@ function configured_kernel(arch, grid, workspec, kernel!; kwargs...) workgroup, worksize = work_layout(grid, workspec; - include_right_boundaries, - reduced_dimensions, - location) + include_right_boundaries, + reduced_dimensions, + location) offset = offsets(workspec)