diff --git a/Project.toml b/Project.toml index cec2f85661..8451e115c0 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "Oceananigans" uuid = "9e8cae18-63c1-5223-a75c-80ca9d6e9a09" authors = ["Climate Modeling Alliance and contributors"] -version = "0.94.3" +version = "0.94.4" [deps] Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e" diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 6d53944ebd..ad0dbe784a 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl @@ -63,21 +63,26 @@ end const EmptyNamedTuple = NamedTuple{(),Tuple{}} +hasclosure(closure, ClosureType) = closure isa ClosureType +hasclosure(closure_tuple::Tuple, ClosureType) = any(hasclosure(c, ClosureType) for c in closure_tuple) + ab2_step_tracers!(::EmptyNamedTuple, model, Δt, χ) = nothing function ab2_step_tracers!(tracers, model, Δt, χ) closure = model.closure + catke_in_closures = hasclosure(closure, FlavorOfCATKE) + td_in_closures = hasclosure(closure, FlavorOfTD) + # Tracer update kernels for (tracer_index, tracer_name) in enumerate(propertynames(tracers)) - # TODO: do better than this silly criteria, also need to check closure tuples - if closure isa FlavorOfCATKE && tracer_name == :e + if catke_in_closures && tracer_name == :e @debug "Skipping AB2 step for e" - elseif closure isa FlavorOfTD && tracer_name == :ϵ + elseif td_in_closures && tracer_name == :ϵ @debug "Skipping AB2 step for ϵ" - elseif closure isa FlavorOfTD && tracer_name == :e + elseif td_in_closures && tracer_name == :e @debug "Skipping AB2 step for e" else Gⁿ = model.timestepper.Gⁿ[tracer_name] diff --git a/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl b/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl index a84d03e13f..6c6fe5395c 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/store_hydrostatic_free_surface_tendencies.jl @@ -32,14 +32,16 @@ function store_tendencies!(model::HydrostaticFreeSurfaceModel) three_dimensional_prognostic_field_names = filter(name -> name != :η, prognostic_field_names) closure = model.closure + catke_in_closures = hasclosure(closure, FlavorOfCATKE) + td_in_closures = hasclosure(closure, FlavorOfTD) for field_name in three_dimensional_prognostic_field_names - if closure isa FlavorOfCATKE && field_name == :e + if catke_in_closures && field_name == :e @debug "Skipping store tendencies for e" - elseif closure isa FlavorOfTD && field_name == :ϵ + elseif td_in_closures && field_name == :ϵ @debug "Skipping store tendencies for ϵ" - elseif closure isa FlavorOfTD && field_name == :e + elseif td_in_closures && field_name == :e @debug "Skipping store tendencies for e" else launch!(model.architecture, model.grid, :xyz, diff --git a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl index 892a8c041c..9cf74fcc5b 100644 --- a/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl +++ b/src/TurbulenceClosures/turbulence_closure_implementations/TKEBasedVerticalDiffusivities/time_step_catke_equation.jl @@ -163,7 +163,7 @@ end # See below. α = convert(FT, 1.5) + χ β = convert(FT, 0.5) + χ - + @inbounds begin total_Gⁿe = slow_Gⁿe[i, j, k] + fast_Gⁿe e[i, j, k] += Δτ * (α * total_Gⁿe - β * G⁻e[i, j, k]) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 28c8d226b2..98af6d56b0 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -1,5 +1,6 @@ include("dependencies_for_runtests.jl") +using Random using Oceananigans.TurbulenceClosures: CATKEVerticalDiffusivity, RiBasedVerticalDiffusivity, DiscreteDiffusionFunction using Oceananigans.TurbulenceClosures: viscosity_location, diffusivity_location, @@ -205,6 +206,42 @@ function time_step_with_tupled_closure(FT, arch) return true end +function run_catke_tke_substepping_tests(arch, closure) + # A large domain to make sure we do not have viscous CFL problems + # with the explicit CATKE time-stepping necessary for this test + grid = RectilinearGrid(arch, size=(2, 2, 2), extent=(100, 200, 300)) + + model = HydrostaticFreeSurfaceModel(; grid, momentum_advection = nothing, tracer_advection = nothing, + closure, buoyancy=BuoyancyTracer(), tracers=(:b, :e)) + + # set random velocities + Random.seed!(1234) + set!(model, u = (x, y, z) -> rand(), v = (x, y, z) -> rand()) + + # time step the model + time_step!(model, 1) + + # Check that eⁿ⁺¹ == Δt * Gⁿ.e with Δt = 1 (euler step) + @test model.tracers.e ≈ model.timestepper.G⁻.e + + eⁿ = deepcopy(model.tracers.e) + G⁻⁻ = deepcopy(model.timestepper.G⁻.e) + + # time step the model again + time_step!(model, 1) + G⁻ = model.timestepper.G⁻.e + + C₁ = 1.5 + model.timestepper.χ + C₂ = 0.5 + model.timestepper.χ + + eⁿ⁺¹ = compute!(Field(eⁿ + C₁ * G⁻ - C₂ * G⁻⁻)) + + # Check that eⁿ⁺¹ == eⁿ + Δt * (C₁ Gⁿ.e - C₂ G⁻.e) + @test model.tracers.e ≈ eⁿ⁺¹ + + return model +end + function run_time_step_with_catke_tests(arch, closure) grid = RectilinearGrid(arch, size=(2, 2, 2), extent=(1, 2, 3)) buoyancy = BuoyancyTracer() @@ -348,25 +385,34 @@ end @info " Testing time-stepping with CATKE closure and closure tuples with CATKE..." for arch in archs @info " Testing time-stepping CATKE by itself..." - closure = CATKEVerticalDiffusivity() - run_time_step_with_catke_tests(arch, closure) + catke = CATKEVerticalDiffusivity() + explicit_catke = CATKEVerticalDiffusivity(ExplicitTimeDiscretization()) + run_time_step_with_catke_tests(arch, catke) + run_catke_tke_substepping_tests(arch, explicit_catke) @info " Testing time-stepping CATKE in a 2-tuple with HorizontalScalarDiffusivity..." - closure = (CATKEVerticalDiffusivity(), HorizontalScalarDiffusivity()) + closure = (catke, HorizontalScalarDiffusivity()) model = run_time_step_with_catke_tests(arch, closure) @test first(model.closure) === closure[1] + closure = (explicit_catke, HorizontalScalarDiffusivity()) + run_catke_tke_substepping_tests(arch, closure) + # Test that closure tuples with CATKE are correctly reordered @info " Testing time-stepping CATKE in a 2-tuple with HorizontalScalarDiffusivity..." - closure = (HorizontalScalarDiffusivity(), CATKEVerticalDiffusivity()) + closure = (HorizontalScalarDiffusivity(), catke) model = run_time_step_with_catke_tests(arch, closure) @test first(model.closure) === closure[2] + closure = (HorizontalScalarDiffusivity(), explicit_catke) + run_catke_tke_substepping_tests(arch, closure) # These are slow to compile... @info " Testing time-stepping CATKE in a 3-tuple..." - closure = (HorizontalScalarDiffusivity(), CATKEVerticalDiffusivity(), VerticalScalarDiffusivity()) + closure = (HorizontalScalarDiffusivity(), catke, VerticalScalarDiffusivity()) model = run_time_step_with_catke_tests(arch, closure) @test first(model.closure) === closure[2] + closure = (HorizontalScalarDiffusivity(), explicit_catke, VerticalScalarDiffusivity()) + run_catke_tke_substepping_tests(arch, closure) end end