One degree simulation in the examples

Project.toml

@@ -42,9 +42,9 @@ JLD2 = "0.4, 0.5"
 KernelAbstractions = "0.9"
 MPI = "0.20"
 NCDatasets = "0.12, 0.13, 0.14"
-Oceananigans = "0.94.3 - 0.99" # This needs to be 0.94.4 when PolarBoundaryCondition is implemented
+Oceananigans = "0.95 - 0.99"
 OffsetArrays = "1.14"
-OrthogonalSphericalShellGrids = "0.1.9"
+OrthogonalSphericalShellGrids = "0.2.0"
 Scratch = "1"
 SeawaterPolynomials = "0.3.4"
 StaticArrays = "1"

experiments/one_degree_simulation/one_degree_simulation.jl

@@ -46,7 +46,7 @@ grid = ImmersedBoundaryGrid(underlying_grid, GridFittedBottom(tampered_bottom_he
 
 gm = Oceananigans.TurbulenceClosures.IsopycnalSkewSymmetricDiffusivity(κ_skew=1000, κ_symmetric=1000)
 catke = ClimaOcean.OceanSimulations.default_ocean_closure()
-viscous_closure = Oceananigans.TurbulenceClosures.HorizontalScalarDiffusivity(ν=2000)
+viscous_closure = HorizontalScalarBiharmonicDiffusivity(ν = 1e11)
 
 closure = (gm, catke, viscous_closure)
 
@@ -95,7 +95,7 @@ atmosphere = JRA55_prescribed_atmosphere(arch; backend=JRA55NetCDFBackend(20))
 #####
 
 coupled_model = OceanSeaIceModel(ocean; atmosphere, radiation) 
-simulation = Simulation(coupled_model; Δt=15minutes, stop_time=2*365days)
+simulation = Simulation(coupled_model; Δt=2minutes, stop_time=30days)
 
 #####
 ##### Run it!
@@ -127,3 +127,7 @@ end
 add_callback!(simulation, progress, IterationInterval(10))
 
 run!(simulation)
+
+simulation.Δt = 25minutes
+
+run!(simulation)

src/OceanSimulations/OceanSimulations.jl

@@ -6,7 +6,7 @@ using Oceananigans
 using Oceananigans.Units
 using Oceananigans.Utils: with_tracers
 using Oceananigans.Advection: FluxFormAdvection
-using Oceananigans.Coriolis: ActiveCellEnstrophyConserving
+using Oceananigans.DistributedComputations: DistributedGrid, all_reduce
 using Oceananigans.ImmersedBoundaries: immersed_peripheral_node, inactive_node
 using OrthogonalSphericalShellGrids
 
@@ -16,6 +16,7 @@ using Oceananigans.TurbulenceClosures.TKEBasedVerticalDiffusivities:
     CATKEEquation
 
 using SeawaterPolynomials.TEOS10: TEOS10EquationOfState
+using Statistics: mean
 
 using Oceananigans.BuoyancyModels: g_Earth
 using Oceananigans.Coriolis: Ω_Earth
@@ -42,10 +43,43 @@ default_or_override(override, alternative_default=nothing) = override
 # Some defaults
 default_free_surface(grid) = SplitExplicitFreeSurface(grid; cfl=0.7)
 
-# 70 substeps is a safe rule of thumb for an ocean at 1/4 - 1/10th of a degree
-# TODO: pass the cfl and a given Δt to calculate the number of substeps?
+function compute_maximum_Δt(grid)
+    Δx = mean(xspacings(grid))
+    Δy = mean(yspacings(grid))
+    Δθ = rad2deg(mean([Δx, Δy])) / grid.radius
+
+    # The maximum Δt is roughly 40minutes / Δθ, giving:
+    # - 40 minutes for a 1 degree ocean
+    # - 20 minutes for a 1/4 degree ocean
+    # - 10 minutes for a 1/8 degree ocean
+    # - 5 minutes for a 1/16 degree ocean
+    # - 2.5 minutes for a 1/32 degree ocean
+
+    return 40minutes / Δθ
+end
+
 const TripolarOfSomeKind = Union{TripolarGrid, ImmersedBoundaryGrid{<:Any, <:Any, <:Any, <:Any, <:TripolarGrid}}
-default_free_surface(grid::TripolarOfSomeKind) = SplitExplicitFreeSurface(grid; substeps=70)
+
+function default_free_surface(grid::TripolarOfSomeKind; 
+                              fixed_Δt = compute_maximum_Δt(grid),
+                              cfl = 0.7) 
+    free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt)
+    @info "Using a $(free_surface)"
+    return free_surface
+end
+
+function default_free_surface(grid::DistributedGrid; 
+                              fixed_Δt = compute_maximum_Δt(grid),
+                              cfl = 0.7) 
+
+    free_surface = SplitExplicitFreeSurface(grid; cfl, fixed_Δt)
+    substeps = length(free_surface.substepping.averaging_weights)
+    substeps = all_reduce(max, substeps, architecture(grid))
+    free_surface = SplitExplicitFreeSurface(grid; substeps)
+    @info "Using a $(free_surface)"
+    return free_surface
+end
+
 
 function default_ocean_closure()
     mixing_length = CATKEMixingLength(Cᵇ=0.01)
@@ -75,7 +109,7 @@ default_tracer_advection() = FluxFormAdvection(WENO(order=7),
 # TODO: Specify the grid to a grid on the sphere; otherwise we can provide a different
 # function that requires latitude and longitude etc for computing coriolis=FPlane...
 function ocean_simulation(grid;
-                          Δt = 5minutes,
+                          Δt = compute_maximum_Δt(grid),
                           closure = default_ocean_closure(),
                           tracers = (:T, :S),
                           free_surface = default_free_surface(grid),