From 3d9668b7957d0770bfb7f375bd9c848efdaa9bd8 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Thu, 7 Mar 2024 07:11:46 -0700 Subject: [PATCH] Change default scheme for `HydrostaticSphericalCoriolis` (#3500) * Change default scheme for `HydrostaticSphericalCoriolis` The default should be `ActiveCellEnstrophyConserving` according to @simone-silvestri * change default EnergyConserving -> ActiveCellEnstrophyConserving * Fix unit test * Restore coriolis formulation for regression tests * fix import of ActiveCellEnstrophyConserving * add rotation_rate kwarg in doc * remove some unnecessary inbounds --------- Co-authored-by: Navid C. Constantinou --- src/Coriolis/Coriolis.jl | 4 ++- .../hydrostatic_spherical_coriolis.jl | 25 ++++++++++--------- ...ostatic_free_turbulence_regression_test.jl | 13 ++++++---- test/test_coriolis.jl | 22 ++++++++-------- 4 files changed, 35 insertions(+), 29 deletions(-) diff --git a/src/Coriolis/Coriolis.jl b/src/Coriolis/Coriolis.jl index af1687704a..0df9aed11a 100644 --- a/src/Coriolis/Coriolis.jl +++ b/src/Coriolis/Coriolis.jl @@ -12,7 +12,9 @@ using Oceananigans.Operators # Physical constants for constructors. using Oceananigans.Grids: R_Earth -const Ω_Earth = 7.292115e-5 # [s⁻¹] https://en.wikipedia.org/wiki/Earth%27s_rotation#Angular_speed + +"Earth's rotation rate [s⁻¹]; see https://en.wikipedia.org/wiki/Earth%27s_rotation#Angular_speed" +const Ω_Earth = 7.292115e-5 """ AbstractRotation diff --git a/src/Coriolis/hydrostatic_spherical_coriolis.jl b/src/Coriolis/hydrostatic_spherical_coriolis.jl index 6dc98c6496..294e05487c 100644 --- a/src/Coriolis/hydrostatic_spherical_coriolis.jl +++ b/src/Coriolis/hydrostatic_spherical_coriolis.jl @@ -26,19 +26,19 @@ end """ HydrostaticSphericalCoriolis([FT=Float64;] rotation_rate = Ω_Earth, - scheme = EnergyConserving()) + scheme = ActiveCellEnstrophyConserving()) Return a parameter object for Coriolis forces on a sphere rotating at `rotation_rate`. -By default, `rotation_rate` is assumed to be Earth's. Keyword arguments ================= -- `scheme`: Either `EnergyConserving()` (default), `EnstrophyConserving()`, or `ActiveCellEnstrophyConserving()`. +- `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref). +- `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `ActiveCellEnstrophyConserving()` (default). """ function HydrostaticSphericalCoriolis(FT::DataType=Float64; rotation_rate = Ω_Earth, - scheme :: S = EnergyConserving(FT)) where S + scheme :: S = ActiveCellEnstrophyConserving()) where S return HydrostaticSphericalCoriolis{S, FT}(rotation_rate, scheme) end @@ -68,24 +68,24 @@ const CoriolisActiveCellEnstrophyConserving = HydrostaticSphericalCoriolis{<:Act @inline not_peripheral_node(args...) = !peripheral_node(args...) @inline function mask_inactive_points_ℑxyᶠᶜᵃ(i, j, k, grid, f::Function, args...) - neighboring_active_nodes = @inbounds ℑxyᶠᶜᵃ(i, j, k, grid, not_peripheral_node, Center(), Face(), Center()) + neighboring_active_nodes = ℑxyᶠᶜᵃ(i, j, k, grid, not_peripheral_node, Center(), Face(), Center()) return ifelse(neighboring_active_nodes == 0, zero(grid), - @inbounds ℑxyᶠᶜᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes) + ℑxyᶠᶜᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes) end @inline function mask_inactive_points_ℑxyᶜᶠᵃ(i, j, k, grid, f::Function, args...) neighboring_active_nodes = @inbounds ℑxyᶜᶠᵃ(i, j, k, grid, not_peripheral_node, Face(), Center(), Center()) return ifelse(neighboring_active_nodes == 0, zero(grid), - @inbounds ℑxyᶜᶠᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes) + ℑxyᶜᶠᵃ(i, j, k, grid, f, args...) / neighboring_active_nodes) end @inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisActiveCellEnstrophyConserving, U) = @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * - mask_inactive_points_ℑxyᶠᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) / Δxᶠᶜᶜ(i, j, k, grid) + mask_inactive_points_ℑxyᶠᶜᵃ(i, j, k, grid, Δx_qᶜᶠᶜ, U[2]) / Δxᶠᶜᶜ(i, j, k, grid) @inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisActiveCellEnstrophyConserving, U) = @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * - mask_inactive_points_ℑxyᶜᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) / Δyᶜᶠᶜ(i, j, k, grid) + mask_inactive_points_ℑxyᶜᶠᵃ(i, j, k, grid, Δy_qᶠᶜᶜ, U[1]) / Δyᶜᶠᶜ(i, j, k, grid) ##### ##### Enstrophy-conserving scheme @@ -94,10 +94,12 @@ end const CoriolisEnstrophyConserving = HydrostaticSphericalCoriolis{<:EnstrophyConserving} @inline x_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) = - @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δx_qᶜᶠᶜ, U[2]) / Δxᶠᶜᶜ(i, j, k, grid) + @inbounds - ℑyᵃᶜᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * + ℑxᶠᵃᵃ(i, j, k, grid, ℑyᵃᶜᵃ, Δx_qᶜᶠᶜ, U[2]) / Δxᶠᶜᶜ(i, j, k, grid) @inline y_f_cross_U(i, j, k, grid, coriolis::CoriolisEnstrophyConserving, U) = - @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δy_qᶠᶜᶜ, U[1]) / Δyᶜᶠᶜ(i, j, k, grid) + @inbounds + ℑxᶜᵃᵃ(i, j, k, grid, fᶠᶠᵃ, coriolis) * + ℑyᵃᶠᵃ(i, j, k, grid, ℑxᶜᵃᵃ, Δy_qᶠᶜᶜ, U[1]) / Δyᶜᶠᶜ(i, j, k, grid) ##### ##### Energy-conserving scheme @@ -128,4 +130,3 @@ function Base.show(io::IO, hydrostatic_spherical_coriolis::HydrostaticSphericalC "├─ rotation rate: ", rotation_rate_str, '\n', "└─ scheme: ", summary(coriolis_scheme)) end - diff --git a/test/regression_tests/hydrostatic_free_turbulence_regression_test.jl b/test/regression_tests/hydrostatic_free_turbulence_regression_test.jl index add78c4d11..7b52aaadfc 100644 --- a/test/regression_tests/hydrostatic_free_turbulence_regression_test.jl +++ b/test/regression_tests/hydrostatic_free_turbulence_regression_test.jl @@ -1,6 +1,7 @@ using Oceananigans.Fields: FunctionField using Oceananigans.Grids: architecture using Oceananigans.BoundaryConditions: fill_halo_regions! +using Oceananigans.Coriolis: EnergyConserving using Oceananigans.Models.HydrostaticFreeSurfaceModels: HydrostaticFreeSurfaceModel, VectorInvariant using Oceananigans.TurbulenceClosures: HorizontalScalarDiffusivity @@ -27,11 +28,13 @@ function run_hydrostatic_free_turbulence_regression_test(grid, free_surface; reg ##### Constructing Grid and model ##### - model = HydrostaticFreeSurfaceModel(grid = grid, - momentum_advection = VectorInvariant(), - free_surface = free_surface, - coriolis = HydrostaticSphericalCoriolis(), - closure = HorizontalScalarDiffusivity(ν=1e+5, κ=1e+4)) + # This coriolis scheme was used to generated the regression test data + coriolis = HydrostaticSphericalCoriolis(scheme=EnergyConserving()) + + model = HydrostaticFreeSurfaceModel(; grid, coriolis, + momentum_advection = VectorInvariant(), + free_surface = free_surface, + closure = HorizontalScalarDiffusivity(ν=1e+5, κ=1e+4)) ##### ##### Imposing initial conditions: diff --git a/test/test_coriolis.jl b/test/test_coriolis.jl index 143c18beb1..dcd1718dab 100644 --- a/test/test_coriolis.jl +++ b/test/test_coriolis.jl @@ -1,6 +1,6 @@ include("dependencies_for_runtests.jl") -using Oceananigans.Coriolis: Ω_Earth +using Oceananigans.Coriolis: Ω_Earth, ActiveCellEnstrophyConserving using Oceananigans.Advection: EnergyConserving, EnstrophyConserving function instantiate_fplane_1(FT) @@ -21,9 +21,9 @@ end function instantiate_constant_coriolis_2(FT) coriolis = ConstantCartesianCoriolis(FT, f=10, rotation_axis=[√(1/3),√(1/3),√(1/3)]) - @test coriolis.fx ≈ FT(10*√(1/3)) - @test coriolis.fy ≈ FT(10*√(1/3)) - @test coriolis.fz ≈ FT(10*√(1/3)) + @test coriolis.fx ≈ FT(10 * √(1/3)) + @test coriolis.fy ≈ FT(10 * √(1/3)) + @test coriolis.fz ≈ FT(10 * √(1/3)) end function instantiate_betaplane_1(FT) @@ -49,26 +49,26 @@ end function instantiate_ntbetaplane_2(FT) Ω, φ, R = π, 17, ℯ coriolis = NonTraditionalBetaPlane(FT, rotation_rate=Ω, latitude=φ, radius=R) - @test coriolis.fz ≈ FT(+2Ω*sind(φ)) - @test coriolis.fy ≈ FT(+2Ω*cosd(φ)) - @test coriolis.β ≈ FT(+2Ω*cosd(φ)/R) - @test coriolis.γ ≈ FT(-4Ω*sind(φ)/R) + @test coriolis.fz ≈ FT(+ 2Ω * sind(φ)) + @test coriolis.fy ≈ FT(+ 2Ω * cosd(φ)) + @test coriolis.β ≈ FT(+ 2Ω * cosd(φ) / R) + @test coriolis.γ ≈ FT(- 4Ω * sind(φ) / R) end function instantiate_hydrostatic_spherical_coriolis1(FT) coriolis = HydrostaticSphericalCoriolis(FT, scheme=EnergyConserving()) - @test coriolis.rotation_rate == FT(Ω_Earth) + @test coriolis.rotation_rate == FT(Ω_Earth) # default @test coriolis.scheme isa EnergyConserving coriolis = HydrostaticSphericalCoriolis(FT, scheme=EnstrophyConserving()) - @test coriolis.rotation_rate == FT(Ω_Earth) + @test coriolis.rotation_rate == FT(Ω_Earth) # default @test coriolis.scheme isa EnstrophyConserving end function instantiate_hydrostatic_spherical_coriolis2(FT) coriolis = HydrostaticSphericalCoriolis(FT, rotation_rate=π) @test coriolis.rotation_rate == FT(π) - @test coriolis.scheme isa EnergyConserving # default + @test coriolis.scheme isa ActiveCellEnstrophyConserving # default end @testset "Coriolis" begin