From 34a87fbe46bb3ee01f0e1058d090d300317f6f20 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 10:01:31 +0100 Subject: [PATCH 01/12] bugfix catke --- .../hydrostatic_free_surface_ab2_step.jl | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 6d53944ebd..84f374c728 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_closures = hasclosure(closure, FlavorOfCATKE) + td_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_closures && tracer_name == :e @debug "Skipping AB2 step for e" - elseif closure isa FlavorOfTD && tracer_name == :ϵ + elseif td_closures && tracer_name == :ϵ @debug "Skipping AB2 step for ϵ" - elseif closure isa FlavorOfTD && tracer_name == :e + elseif td_closures && tracer_name == :e @debug "Skipping AB2 step for e" else Gⁿ = model.timestepper.Gⁿ[tracer_name] From 7deeb38e5c48949b47c4075b7811f465c60eb1a0 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Tue, 26 Nov 2024 10:03:21 +0100 Subject: [PATCH 02/12] better comment --- .../hydrostatic_free_surface_ab2_step.jl | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl index 84f374c728..ad0dbe784a 100644 --- a/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl +++ b/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_ab2_step.jl @@ -72,17 +72,17 @@ function ab2_step_tracers!(tracers, model, Δt, χ) closure = model.closure - catke_closures = hasclosure(closure, FlavorOfCATKE) - td_closures = hasclosure(closure, FlavorOfTD) + catke_in_closures = hasclosure(closure, FlavorOfCATKE) + td_in_closures = hasclosure(closure, FlavorOfTD) # Tracer update kernels for (tracer_index, tracer_name) in enumerate(propertynames(tracers)) - if catke_closures && tracer_name == :e + if catke_in_closures && tracer_name == :e @debug "Skipping AB2 step for e" - elseif td_closures && tracer_name == :ϵ + elseif td_in_closures && tracer_name == :ϵ @debug "Skipping AB2 step for ϵ" - elseif td_closures && tracer_name == :e + elseif td_in_closures && tracer_name == :e @debug "Skipping AB2 step for e" else Gⁿ = model.timestepper.Gⁿ[tracer_name] From 1ca62e8596ffebb4482937f4b1e825214692a88c Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 11:15:48 +0100 Subject: [PATCH 03/12] add catke timestepping tests --- test/test_turbulence_closures.jl | 34 ++++++++++++++++++++++++++++++++ 1 file changed, 34 insertions(+) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 28c8d226b2..0227571905 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -205,6 +205,35 @@ function time_step_with_tupled_closure(FT, arch) return true end +function run_catke_tke_substepping_tests(arch, closure) + grid = RectilinearGrid(arch, size=(2, 2, 2), extent=(1, 2, 3)) + + model = HydrostaticFreeSurfaceModel(; grid, momentum_advection = nothing, tracer = nothing, + closure, buoyancy=BuoyancyTracer(), tracers=(:b, :e)) + + # set random velocities + 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) + + eⁿ⁺¹ = compute!(Field(eⁿ + 1.5 * G⁻ - 0.5 * G⁻)) + + # Check that eⁿ⁺¹ = eⁿ + Δt * (1.5Gⁿ.e - 0.5G⁻.e) with Δt = 1 + @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() @@ -350,23 +379,28 @@ end @info " Testing time-stepping CATKE by itself..." closure = CATKEVerticalDiffusivity() run_time_step_with_catke_tests(arch, closure) + run_catke_tke_substepping_tests(arch, closure) @info " Testing time-stepping CATKE in a 2-tuple with HorizontalScalarDiffusivity..." closure = (CATKEVerticalDiffusivity(), HorizontalScalarDiffusivity()) model = run_time_step_with_catke_tests(arch, closure) @test first(model.closure) === closure[1] + 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()) model = run_time_step_with_catke_tests(arch, closure) @test first(model.closure) === closure[2] + 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()) model = run_time_step_with_catke_tests(arch, closure) @test first(model.closure) === closure[2] + run_catke_tke_substepping_tests(arch, closure) end end From 8533e0802d7259b5fcc1cebec3e95ae8f7bb77fc Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 11:19:15 +0100 Subject: [PATCH 04/12] add a test --- test/test_turbulence_closures.jl | 20 +++++++++++++------- 1 file changed, 13 insertions(+), 7 deletions(-) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 0227571905..d08a7817cb 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -206,7 +206,9 @@ function time_step_with_tupled_closure(FT, arch) end function run_catke_tke_substepping_tests(arch, closure) - grid = RectilinearGrid(arch, size=(2, 2, 2), extent=(1, 2, 3)) + # 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 = nothing, closure, buoyancy=BuoyancyTracer(), tracers=(:b, :e)) @@ -377,29 +379,33 @@ 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) - run_catke_tke_substepping_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 From 07045e689b1e67c26e309f83dfef10bc4dc4dc81 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 11:24:49 +0100 Subject: [PATCH 05/12] test fix --- test/test_turbulence_closures.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index d08a7817cb..44e2e774a8 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -227,8 +227,9 @@ function run_catke_tke_substepping_tests(arch, closure) # time step the model again time_step!(model, 1) + Gⁿ = model.timestepper.Gⁿ.e - eⁿ⁺¹ = compute!(Field(eⁿ + 1.5 * G⁻ - 0.5 * G⁻)) + eⁿ⁺¹ = compute!(Field(eⁿ + 1.5 * Gⁿ - 0.5 * G⁻)) # Check that eⁿ⁺¹ = eⁿ + Δt * (1.5Gⁿ.e - 0.5G⁻.e) with Δt = 1 @test model.tracers.e ≈ eⁿ⁺¹ From 729a113004f6ac0f7c456c19d3c6f188d7345991 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 12:17:22 +0100 Subject: [PATCH 06/12] fix test --- test/test_turbulence_closures.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 44e2e774a8..6d1a6be814 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -210,7 +210,7 @@ function run_catke_tke_substepping_tests(arch, closure) # 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 = nothing, + model = HydrostaticFreeSurfaceModel(; grid, momentum_advection = nothing, tracer_advection = nothing, closure, buoyancy=BuoyancyTracer(), tracers=(:b, :e)) # set random velocities From 530d79e99b277f6e5de22670104d253754452491 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Wed, 27 Nov 2024 18:45:34 +0100 Subject: [PATCH 07/12] do not rewrite tendencies --- src/TimeSteppers/store_tendencies.jl | 25 ++++++++++++++++--- src/TurbulenceClosures/TurbulenceClosures.jl | 7 ++++++ .../time_step_catke_equation.jl | 1 + test/test_turbulence_closures.jl | 14 ++++++----- 4 files changed, 37 insertions(+), 10 deletions(-) diff --git a/src/TimeSteppers/store_tendencies.jl b/src/TimeSteppers/store_tendencies.jl index 2eae58dc2c..7b12de464f 100644 --- a/src/TimeSteppers/store_tendencies.jl +++ b/src/TimeSteppers/store_tendencies.jl @@ -8,14 +8,31 @@ using Oceananigans.Utils: launch! @inbounds G⁻[i, j, k] = G⁰[i, j, k] end +has_catke_closure() = false +has_td_closure() = false + """ Store previous source terms before updating them. """ function store_tendencies!(model) model_fields = prognostic_fields(model) - for field_name in keys(model_fields) - launch!(model.architecture, model.grid, :xyz, store_field_tendencies!, - model.timestepper.G⁻[field_name], - model.timestepper.Gⁿ[field_name]) + closure = model.closure + + catke_in_closures = has_catke_closure(closure) + td_in_closures = has_td_closure(closure) + + # Tracer update kernels + for field_name in keys(model_fields) + if catke_in_closures && field_name == :e + @debug "Skipping AB2 step for e" + elseif td_in_closures && field_name == :ϵ + @debug "Skipping AB2 step for ϵ" + elseif td_in_closures && field_name == :e + @debug "Skipping AB2 step for e" + else + launch!(model.architecture, model.grid, :xyz, store_field_tendencies!, + model.timestepper.G⁻[field_name], + model.timestepper.Gⁿ[field_name]) + end end return nothing diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl index 2e988816a4..3ed52001a1 100644 --- a/src/TurbulenceClosures/TurbulenceClosures.jl +++ b/src/TurbulenceClosures/TurbulenceClosures.jl @@ -58,6 +58,7 @@ using Oceananigans.ImmersedBoundaries: AbstractGridFittedBottom import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z import Oceananigans.Architectures: on_architecture +import Oceananigans.TimeSteppers: has_catke_closure, has_td_closure const VerticallyBoundedGrid{FT} = AbstractGrid{FT, <:Any, <:Any, <:Bounded} @@ -196,6 +197,12 @@ using .Smagorinskys: LillyCoefficient, DynamicCoefficient, LagrangianAveraging include("diffusivity_fields.jl") include("turbulence_closure_diagnostics.jl") +has_catke_closure(closure) = closure isa CATKEVerticalDiffusivity +has_catke_closure(closure::Tuple) = any(hasclosure(c, CATKEVerticalDiffusivity) for c in closure) + +has_td_closure(closure) = closure isa TKEVerticalDiffusivity +has_td_closure(closure::Tuple) = any(hasclosure(c, TKEVerticalDiffusivity) for c in closure) + ##### ##### Some value judgements here ##### 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..6d76fcf4d1 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 @@ -46,6 +46,7 @@ function time_step_catke_equation!(model) FT = eltype(grid) for m = 1:M # substep + @show m if m == 1 && M != 1 χ = convert(FT, -0.5) # Euler step for the first substep else diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 6d1a6be814..e264028e7c 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, @@ -211,25 +212,26 @@ function run_catke_tke_substepping_tests(arch, closure) 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)) + closure = explicit_catke, 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 + @test model.tracers.e ≈ model.timestepper.G⁻.e - eⁿ = deepcopy(model.tracers.e) - G⁻ = deepcopy(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 + G⁻ = model.timestepper.G⁻.e - eⁿ⁺¹ = compute!(Field(eⁿ + 1.5 * Gⁿ - 0.5 * G⁻)) + eⁿ⁺¹ = compute!(Field(eⁿ + 1.5 * G⁻ - 0.5 * G⁻⁻)) # Check that eⁿ⁺¹ = eⁿ + Δt * (1.5Gⁿ.e - 0.5G⁻.e) with Δt = 1 @test model.tracers.e ≈ eⁿ⁺¹ From 8855d2227c8803b49273fdc4b2616585b9406332 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 09:59:26 +0100 Subject: [PATCH 08/12] tested and everything --- test/test_turbulence_closures.jl | 10 +++++++--- 1 file changed, 7 insertions(+), 3 deletions(-) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index e264028e7c..3e575b6cc5 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -212,7 +212,7 @@ function run_catke_tke_substepping_tests(arch, closure) grid = RectilinearGrid(arch, size=(2, 2, 2), extent=(100, 200, 300)) model = HydrostaticFreeSurfaceModel(; grid, momentum_advection = nothing, tracer_advection = nothing, - closure = explicit_catke, buoyancy=BuoyancyTracer(), tracers=(:b, :e)) + closure, buoyancy=BuoyancyTracer(), tracers=(:b, :e)) # set random velocities Random.seed!(1234) @@ -231,9 +231,13 @@ function run_catke_tke_substepping_tests(arch, closure) time_step!(model, 1) G⁻ = model.timestepper.G⁻.e - eⁿ⁺¹ = compute!(Field(eⁿ + 1.5 * G⁻ - 0.5 * G⁻⁻)) + C₁ = 1.5 + model.timestepper.χ + C₂ = 0.5 + model.timestepper.χ - # Check that eⁿ⁺¹ = eⁿ + Δt * (1.5Gⁿ.e - 0.5G⁻.e) with Δt = 1 + eⁿ⁺¹ = compute!(Field(eⁿ + C₁ * G⁻ - C₂ * G⁻⁻)) + + # Check that eⁿ⁺¹ = eⁿ + Δt * Gⁿ.e + # (with CATKE we always use an Euler Step at the first substep) @test model.tracers.e ≈ eⁿ⁺¹ return model From c5d5928a4b0612c21c9a4daaca65494f06049bc1 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 09:59:34 +0100 Subject: [PATCH 09/12] more testing --- ...ore_hydrostatic_free_surface_tendencies.jl | 8 ++++--- src/TimeSteppers/store_tendencies.jl | 22 +++---------------- .../time_step_catke_equation.jl | 3 +-- 3 files changed, 9 insertions(+), 24 deletions(-) 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/TimeSteppers/store_tendencies.jl b/src/TimeSteppers/store_tendencies.jl index 7b12de464f..3585d5d768 100644 --- a/src/TimeSteppers/store_tendencies.jl +++ b/src/TimeSteppers/store_tendencies.jl @@ -8,31 +8,15 @@ using Oceananigans.Utils: launch! @inbounds G⁻[i, j, k] = G⁰[i, j, k] end -has_catke_closure() = false -has_td_closure() = false - """ Store previous source terms before updating them. """ function store_tendencies!(model) model_fields = prognostic_fields(model) - closure = model.closure - - catke_in_closures = has_catke_closure(closure) - td_in_closures = has_td_closure(closure) - # Tracer update kernels for field_name in keys(model_fields) - if catke_in_closures && field_name == :e - @debug "Skipping AB2 step for e" - elseif td_in_closures && field_name == :ϵ - @debug "Skipping AB2 step for ϵ" - elseif td_in_closures && field_name == :e - @debug "Skipping AB2 step for e" - else - launch!(model.architecture, model.grid, :xyz, store_field_tendencies!, - model.timestepper.G⁻[field_name], - model.timestepper.Gⁿ[field_name]) - end + launch!(model.architecture, model.grid, :xyz, store_field_tendencies!, + model.timestepper.G⁻[field_name], + model.timestepper.Gⁿ[field_name]) end return nothing 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 6d76fcf4d1..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 @@ -46,7 +46,6 @@ function time_step_catke_equation!(model) FT = eltype(grid) for m = 1:M # substep - @show m if m == 1 && M != 1 χ = convert(FT, -0.5) # Euler step for the first substep else @@ -164,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]) From c1c9240e3ed06972475c8046b24648302826b4e4 Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 10:09:58 +0100 Subject: [PATCH 10/12] remove vestigial code --- src/TurbulenceClosures/TurbulenceClosures.jl | 7 ------- 1 file changed, 7 deletions(-) diff --git a/src/TurbulenceClosures/TurbulenceClosures.jl b/src/TurbulenceClosures/TurbulenceClosures.jl index 3ed52001a1..2e988816a4 100644 --- a/src/TurbulenceClosures/TurbulenceClosures.jl +++ b/src/TurbulenceClosures/TurbulenceClosures.jl @@ -58,7 +58,6 @@ using Oceananigans.ImmersedBoundaries: AbstractGridFittedBottom import Oceananigans.Grids: required_halo_size_x, required_halo_size_y, required_halo_size_z import Oceananigans.Architectures: on_architecture -import Oceananigans.TimeSteppers: has_catke_closure, has_td_closure const VerticallyBoundedGrid{FT} = AbstractGrid{FT, <:Any, <:Any, <:Bounded} @@ -197,12 +196,6 @@ using .Smagorinskys: LillyCoefficient, DynamicCoefficient, LagrangianAveraging include("diffusivity_fields.jl") include("turbulence_closure_diagnostics.jl") -has_catke_closure(closure) = closure isa CATKEVerticalDiffusivity -has_catke_closure(closure::Tuple) = any(hasclosure(c, CATKEVerticalDiffusivity) for c in closure) - -has_td_closure(closure) = closure isa TKEVerticalDiffusivity -has_td_closure(closure::Tuple) = any(hasclosure(c, TKEVerticalDiffusivity) for c in closure) - ##### ##### Some value judgements here ##### From a7c2ae08f45c4ab36d944de0040946984c7acaed Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 10:10:39 +0100 Subject: [PATCH 11/12] remove vestigial code 2 --- src/TimeSteppers/store_tendencies.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/TimeSteppers/store_tendencies.jl b/src/TimeSteppers/store_tendencies.jl index 3585d5d768..2eae58dc2c 100644 --- a/src/TimeSteppers/store_tendencies.jl +++ b/src/TimeSteppers/store_tendencies.jl @@ -12,8 +12,7 @@ end function store_tendencies!(model) model_fields = prognostic_fields(model) - # Tracer update kernels - for field_name in keys(model_fields) + for field_name in keys(model_fields) launch!(model.architecture, model.grid, :xyz, store_field_tendencies!, model.timestepper.G⁻[field_name], model.timestepper.Gⁿ[field_name]) From dd15de3c69e66585a7c548c71989ad0a0cd7f38b Mon Sep 17 00:00:00 2001 From: Simone Silvestri Date: Thu, 28 Nov 2024 10:11:56 +0100 Subject: [PATCH 12/12] correct comment --- test/test_turbulence_closures.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/test/test_turbulence_closures.jl b/test/test_turbulence_closures.jl index 3e575b6cc5..aacb08a7dd 100644 --- a/test/test_turbulence_closures.jl +++ b/test/test_turbulence_closures.jl @@ -236,8 +236,7 @@ function run_catke_tke_substepping_tests(arch, closure) eⁿ⁺¹ = compute!(Field(eⁿ + C₁ * G⁻ - C₂ * G⁻⁻)) - # Check that eⁿ⁺¹ = eⁿ + Δt * Gⁿ.e - # (with CATKE we always use an Euler Step at the first substep) + # Check that eⁿ⁺¹ = eⁿ + Δt * (C₁ Gⁿ.e - C₂ G⁻.e) @test model.tracers.e ≈ eⁿ⁺¹ return model