From ab1f6cf2cb6a05409fcf5159628be94f1a6718d2 Mon Sep 17 00:00:00 2001 From: "Gregory L. Wagner" Date: Wed, 6 Mar 2024 09:20:13 -0700 Subject: [PATCH 1/7] Change default scheme for `HydrostaticSphericalCoriolis` The default should be `ActiveCellEnstrophyConserving` according to @simone-silvestri --- src/Coriolis/hydrostatic_spherical_coriolis.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Coriolis/hydrostatic_spherical_coriolis.jl b/src/Coriolis/hydrostatic_spherical_coriolis.jl index 6dc98c6496..117d7b9fdd 100644 --- a/src/Coriolis/hydrostatic_spherical_coriolis.jl +++ b/src/Coriolis/hydrostatic_spherical_coriolis.jl @@ -26,7 +26,7 @@ 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. From 622bf6dd0f85dcb929c339b74079a95ee3150b7e Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Wed, 6 Mar 2024 18:42:30 +0200 Subject: [PATCH 2/7] change default EnergyConserving -> ActiveCellEnstrophyConserving --- src/Coriolis/hydrostatic_spherical_coriolis.jl | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/src/Coriolis/hydrostatic_spherical_coriolis.jl b/src/Coriolis/hydrostatic_spherical_coriolis.jl index 117d7b9fdd..ed96a01760 100644 --- a/src/Coriolis/hydrostatic_spherical_coriolis.jl +++ b/src/Coriolis/hydrostatic_spherical_coriolis.jl @@ -29,16 +29,15 @@ end 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()`. +- `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 From fd990ba435d82b4833b1f44662ec59e35fc5537f Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 6 Mar 2024 22:02:39 -0700 Subject: [PATCH 3/7] Fix unit test --- test/test_coriolis.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/test_coriolis.jl b/test/test_coriolis.jl index 143c18beb1..113bbd08b6 100644 --- a/test/test_coriolis.jl +++ b/test/test_coriolis.jl @@ -1,7 +1,7 @@ include("dependencies_for_runtests.jl") using Oceananigans.Coriolis: Ω_Earth -using Oceananigans.Advection: EnergyConserving, EnstrophyConserving +using Oceananigans.Advection: EnergyConserving, EnstrophyConserving, ActiveCellEnstrophyConserving function instantiate_fplane_1(FT) coriolis = FPlane(FT, f=π) @@ -68,7 +68,7 @@ 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 From 35dc07cde1e4776f15927eef042c5f5f4174d40d Mon Sep 17 00:00:00 2001 From: Gregory Wagner Date: Wed, 6 Mar 2024 22:05:23 -0700 Subject: [PATCH 4/7] Restore coriolis formulation for regression tests --- .../hydrostatic_free_turbulence_regression_test.jl | 13 ++++++++----- 1 file changed, 8 insertions(+), 5 deletions(-) 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: From e4e6fed2cd75de5ad78ac69911258cb2c49a6e5c Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 08:32:11 +0200 Subject: [PATCH 5/7] fix import of ActiveCellEnstrophyConserving --- test/test_coriolis.jl | 22 +++++++++++----------- 1 file changed, 11 insertions(+), 11 deletions(-) diff --git a/test/test_coriolis.jl b/test/test_coriolis.jl index 113bbd08b6..dcd1718dab 100644 --- a/test/test_coriolis.jl +++ b/test/test_coriolis.jl @@ -1,7 +1,7 @@ include("dependencies_for_runtests.jl") -using Oceananigans.Coriolis: Ω_Earth -using Oceananigans.Advection: EnergyConserving, EnstrophyConserving, ActiveCellEnstrophyConserving +using Oceananigans.Coriolis: Ω_Earth, ActiveCellEnstrophyConserving +using Oceananigans.Advection: EnergyConserving, EnstrophyConserving function instantiate_fplane_1(FT) coriolis = FPlane(FT, f=π) @@ -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,19 +49,19 @@ 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 From 0bda9411f5b7207d920a9d1b3f580125dbc8d346 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 08:52:49 +0200 Subject: [PATCH 6/7] add rotation_rate kwarg in doc --- src/Coriolis/Coriolis.jl | 4 +++- src/Coriolis/hydrostatic_spherical_coriolis.jl | 3 ++- 2 files changed, 5 insertions(+), 2 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 ed96a01760..803184e4c5 100644 --- a/src/Coriolis/hydrostatic_spherical_coriolis.jl +++ b/src/Coriolis/hydrostatic_spherical_coriolis.jl @@ -33,7 +33,8 @@ Return a parameter object for Coriolis forces on a sphere rotating at `rotation_ Keyword arguments ================= -- `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `ActiveCellEnstrophyConserving()` (default). +- `rotation_rate`: Sphere's rotation rate; default: [`Ω_Earth`](@ref). +- `scheme`: Either `EnergyConserving()`, `EnstrophyConserving()`, or `ActiveCellEnstrophyConserving()` (default). """ function HydrostaticSphericalCoriolis(FT::DataType=Float64; rotation_rate = Ω_Earth, From 8e54ba98db53419d7945dcc92109ef7de7e7f752 Mon Sep 17 00:00:00 2001 From: "Navid C. Constantinou" Date: Thu, 7 Mar 2024 08:58:24 +0200 Subject: [PATCH 7/7] remove some unnecessary inbounds --- src/Coriolis/hydrostatic_spherical_coriolis.jl | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/src/Coriolis/hydrostatic_spherical_coriolis.jl b/src/Coriolis/hydrostatic_spherical_coriolis.jl index 803184e4c5..294e05487c 100644 --- a/src/Coriolis/hydrostatic_spherical_coriolis.jl +++ b/src/Coriolis/hydrostatic_spherical_coriolis.jl @@ -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 -