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 <[email protected]>
CliMA · Mar 7, 2024 · 3d9668b · 3d9668b
1 parent ce3b755
commit 3d9668b
Show file tree

Hide file tree

Showing 4 changed files with 35 additions and 29 deletions.
diff --git 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
@@ -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
@@ -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
@@ -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