From 8e36421f4b3e4380ef5d6cd3cade40d497858588 Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 30 Oct 2024 12:30:24 -0600 Subject: [PATCH 1/4] Add parameter to epsilon equation so buoyancy flux term may always be positive --- .../tke_dissipation_equations.jl | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl index d24bb8d84b..ab10f29cbe 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/tke_dissipation_equations.jl @@ -10,7 +10,8 @@ using CUDA Base.@kwdef struct TKEDissipationEquations{FT} Cᵋϵ :: FT = 1.92 Cᴾϵ :: FT = 1.44 - Cᵇϵ :: FT = -0.65 + Cᵇϵ⁺ :: FT = -0.65 + Cᵇϵ⁻ :: FT = -0.65 Cᵂu★ :: FT = 0.0 CᵂwΔ :: FT = 0.0 Cᵂα :: FT = 0.11 # Charnock parameter @@ -133,7 +134,11 @@ end # Patankar trick for ϵ-equation Cᵋϵ = closure_ij.tke_dissipation_equations.Cᵋϵ - Cᵇϵ = closure_ij.tke_dissipation_equations.Cᵇϵ + Cᵇϵ⁺ = closure_ij.tke_dissipation_equations.Cᵇϵ⁺ + Cᵇϵ⁻ = closure_ij.tke_dissipation_equations.Cᵇϵ⁻ + + N² = ℑzᵃᵃᶜ(i, j, k, grid, ∂z_b, buoyancy, tracers) + Cᵇϵ = ifelse(N² ≥ 0, Cᵇϵ⁺, Cᵇϵ⁻) Cᵇϵ_wb⁻ = min(Cᵇϵ * wb, zero(grid)) Cᵇϵ_wb⁺ = max(Cᵇϵ * wb, zero(grid)) From 4cc8d5daaa600ad6d479bc20ea4bb99454e35b9f Mon Sep 17 00:00:00 2001 From: glwagner Date: Fri, 1 Nov 2024 02:03:33 -0400 Subject: [PATCH 2/4] Fix on_architecture for ensemble column grids --- src/Grids/rectilinear_grid.jl | 30 ++++++++++++++++--- .../single_column_model_mode.jl | 14 +-------- 2 files changed, 27 insertions(+), 17 deletions(-) diff --git a/src/Grids/rectilinear_grid.jl b/src/Grids/rectilinear_grid.jl index 25a30e7cbc..ec5f7234d1 100644 --- a/src/Grids/rectilinear_grid.jl +++ b/src/Grids/rectilinear_grid.jl @@ -284,6 +284,8 @@ function RectilinearGrid(architecture::AbstractArchitecture = CPU(), Δzᵃᵃᶠ, Δzᵃᵃᶜ, zᵃᵃᶠ, zᵃᵃᶜ) end + + """ Validate user input arguments to the `RectilinearGrid` constructor. """ function validate_rectilinear_grid_args(topology, size, halo, FT, extent, x, y, z) TX, TY, TZ = topology = validate_topology(topology) @@ -343,6 +345,20 @@ function Base.show(io::IO, grid::RectilinearGrid, withsummary=true) "└── ", z_summary) end +##### +##### For "column ensemble models" +##### + +struct ColumnEnsembleSize{C<:Tuple{Int, Int}} + ensemble :: C + Nz :: Int + Hz :: Int +end + +ColumnEnsembleSize(; Nz, ensemble=(0, 0), Hz=1) = ColumnEnsembleSize(ensemble, Nz, Hz) +validate_size(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(e.ensemble[1], e.ensemble[2], e.Nz) +validate_halo(TX, TY, TZ, size, e::ColumnEnsembleSize) = tuple(0, 0, e.Hz) + ##### ##### Utilities ##### @@ -378,10 +394,16 @@ function constructor_arguments(grid::RectilinearGrid) # Kwargs topo = topology(grid) - size = (grid.Nx, grid.Ny, grid.Nz) - halo = (grid.Hx, grid.Hy, grid.Hz) - size = pop_flat_elements(size, topo) - halo = pop_flat_elements(halo, topo) + + if (topo[1] == Flat && grid.Nx > 1) || + (topo[2] == Flat && grid.Ny > 1) + size = halo = ColumnEnsembleSize(Nz=grid.Nz, Hz=grid.Nz, ensemble=(grid.Nx, grid.Ny)) + else + size = (grid.Nx, grid.Ny, grid.Nz) + halo = (grid.Hx, grid.Hy, grid.Hz) + size = pop_flat_elements(size, topo) + halo = pop_flat_elements(halo, topo) + end kwargs = Dict(:size => size, :halo => halo, diff --git a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl index d4d3dec8d8..2b02b2d73d 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/single_column_model_mode.jl @@ -2,13 +2,12 @@ using CUDA: @allowscalar using Oceananigans: UpdateStateCallsite using Oceananigans.Advection: AbstractAdvectionScheme -using Oceananigans.Grids: Flat, Bounded +using Oceananigans.Grids: Flat, Bounded, ColumnEnsembleSize using Oceananigans.Fields: ZeroField using Oceananigans.Coriolis: AbstractRotation using Oceananigans.TurbulenceClosures: AbstractTurbulenceClosure using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities: CATKEVDArray -import Oceananigans.Grids: validate_size, validate_halo import Oceananigans.Models: validate_tracer_advection import Oceananigans.BoundaryConditions: fill_halo_regions! import Oceananigans.TurbulenceClosures: time_discretization, compute_diffusivities! @@ -102,17 +101,6 @@ end return ∇_dot_qᶜ(i, j, k, grid, closure, c, tracer_index, args...) end -struct ColumnEnsembleSize{C<:Tuple{Int, Int}} - ensemble :: C - Nz :: Int - Hz :: Int -end - -ColumnEnsembleSize(; Nz, ensemble=(0, 0), Hz=1) = ColumnEnsembleSize(ensemble, Nz, Hz) - -validate_size(TX, TY, TZ, e::ColumnEnsembleSize) = tuple(e.ensemble[1], e.ensemble[2], e.Nz) -validate_halo(TX, TY, TZ, size, e::ColumnEnsembleSize) = tuple(0, 0, e.Hz) - @inline function time_discretization(closure_array::AbstractArray) first_closure = @allowscalar first(closure_array) # assumes all closures have same time-discretization return time_discretization(first_closure) From 805a7e69e99d8c5c425fe1caa8f95125eeb4b6f3 Mon Sep 17 00:00:00 2001 From: glwagner Date: Wed, 6 Nov 2024 18:25:48 -0500 Subject: [PATCH 3/4] Bugfix in ColumnEnsembleSize in on_architecture --- src/Grids/rectilinear_grid.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Grids/rectilinear_grid.jl b/src/Grids/rectilinear_grid.jl index ec5f7234d1..b86ca0ed6c 100644 --- a/src/Grids/rectilinear_grid.jl +++ b/src/Grids/rectilinear_grid.jl @@ -397,7 +397,7 @@ function constructor_arguments(grid::RectilinearGrid) if (topo[1] == Flat && grid.Nx > 1) || (topo[2] == Flat && grid.Ny > 1) - size = halo = ColumnEnsembleSize(Nz=grid.Nz, Hz=grid.Nz, ensemble=(grid.Nx, grid.Ny)) + size = halo = ColumnEnsembleSize(Nz=grid.Nz, Hz=grid.Hz, ensemble=(grid.Nx, grid.Ny)) else size = (grid.Nx, grid.Ny, grid.Nz) halo = (grid.Hx, grid.Hy, grid.Hz) From 7fdb4af560fd8d302ca46e48a86a53cf796f462e Mon Sep 17 00:00:00 2001 From: glwagner Date: Wed, 6 Nov 2024 18:26:12 -0500 Subject: [PATCH 4/4] Fix formatting --- src/OutputReaders/field_time_series.jl | 21 +++++++++++---------- 1 file changed, 11 insertions(+), 10 deletions(-) diff --git a/src/OutputReaders/field_time_series.jl b/src/OutputReaders/field_time_series.jl index a07377e22c..4f772496a9 100644 --- a/src/OutputReaders/field_time_series.jl +++ b/src/OutputReaders/field_time_series.jl @@ -214,16 +214,16 @@ Base.length(backend::PartlyInMemory) = backend.length ##### mutable struct FieldTimeSeries{LX, LY, LZ, TI, K, I, D, G, ET, B, χ, P, N, KW} <: AbstractField{LX, LY, LZ, G, ET, 4} - data :: D - grid :: G - backend :: K + data :: D + grid :: G + backend :: K boundary_conditions :: B - indices :: I - times :: χ - path :: P - name :: N - time_indexing :: TI - reader_kw :: KW + indices :: I + times :: χ + path :: P + name :: N + time_indexing :: TI + reader_kw :: KW function FieldTimeSeries{LX, LY, LZ}(data::D, grid::G, @@ -350,7 +350,8 @@ new_data(FT, grid, loc, indices, ::Nothing) = nothing # Apparently, not explicitly specifying Int64 in here makes this function # fail on x86 processors where `Int` is implied to be `Int32` -# see ClimaOcean commit 3c47d887659d81e0caed6c9df41b7438e1f1cd52 at https://github.com/CliMA/ClimaOcean.jl/actions/runs/8804916198/job/24166354095) +# see ClimaOcean commit 3c47d887659d81e0caed6c9df41b7438e1f1cd52 at +# https://github.com/CliMA/ClimaOcean.jl/actions/runs/8804916198/job/24166354095) function new_data(FT, grid, loc, indices, Nt::Union{Int, Int64}) space_size = total_size(grid, loc, indices) underlying_data = zeros(FT, architecture(grid), space_size..., Nt)