Error with default HydrostaticFreeSurfaceModel on TripolarGrid

[glwagner]

I was actually expecting a different error. But this isn't a good one either.

julia> using Oceananigans

julia> using OrthogonalSphericalShellGrids

julia> grid = TripolarGrid(size=(4, 4, 4), z=(-1, 0))
4×4×4 OrthogonalSphericalShellGrid{Float64, Periodic, RightConnected, Bounded} on CPU with 4×4×4 halo and with precomputed metrics
├── centered at (λ, φ) = (0.0, -0.3201)
├── longitude: Periodic  extent 360.112 degrees     variably spaced with min(Δλ)=15.6283, max(Δλ)=95.2245
├── latitude:  RightConnected  extent 212.5 degrees variably spaced with min(Δφ)=23.4259, max(Δφ)=56.6667
└── z:         Bounded  z ∈ [-1.0, 0.0]             regularly spaced with Δz=0.25

julia> free_surface = SplitExplicitFreeSurface(grid; cfl=0.7)
SplitExplicitFreeSurface with Barotropic time step equal to 8.112 days


julia> model = HydrostaticFreeSurfaceModel(; grid, free_surface)
ERROR: Centered{1, Float64, Nothing, Nothing, Nothing, Nothing} is not supported with Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{Float64, Periodic, RightConnected, Bounded, OffsetArrays.OffsetMatrix{Float64, Matrix{Float64}}, OffsetArrays.OffsetVector{Float64, StepRangeLen{Float64, Base.TwicePrecision{Float64}, Base.TwicePrecision{Float64}, Int64}}, Float64, OrthogonalSphericalShellGrids.Tripolar{Int64, Int64, Int64}, CPU}
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:35
 [2] validate_momentum_advection(momentum_advection::Centered{…}, grid::Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{…})
   @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/Projects/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl:227
 [3] #construct_regionally#62
   @ ~/Projects/Oceananigans.jl/src/Utils/multi_region_transformation.jl:148 [inlined]
 [4] construct_regionally
   @ ~/Projects/Oceananigans.jl/src/Utils/multi_region_transformation.jl:143 [inlined]
 [5] macro expansion
   @ ~/Projects/Oceananigans.jl/src/Utils/multi_region_transformation.jl:221 [inlined]
 [6] HydrostaticFreeSurfaceModel(; grid::Oceananigans.Grids.ZRegOrthogonalSphericalShellGrid{…}, clock::Clock{…}, momentum_advection::Centered{…}, tracer_advection::Centered{…}, buoyancy::Nothing, coriolis::Nothing, free_surface::SplitExplicitFreeSurface{…}, tracers::Nothing, forcing::@NamedTuple{}, closure::Nothing, boundary_conditions::@NamedTuple{}, particles::Nothing, biogeochemistry::Nothing, velocities::Nothing, pressure::Nothing, diffusivity_fields::Nothing, auxiliary_fields::@NamedTuple{})
   @ Oceananigans.Models.HydrostaticFreeSurfaceModels ~/Projects/Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl:129
 [7] top-level scope
   @ REPL[18]:1
Some type information was truncated. Use `show(err)` to see complete types.

[simone-silvestri]

This looks like an error when using the Centered advection scheme for momentum (the default).
I do not think the tripolar grid spits a useful error when using the tripolar grid. We can add one.

[glwagner]

Why would there be an error when using a Centered advection scheme?

[glwagner]

I'm amazed though - do we not have a test for building HydrostaticFreeSurfaceModel in this simplest case? It seems we need quite a few more tests for the TripolarGrid. We need to test various combination of model inputs and make sure that all of the ones we intent to support are working (free surfaces, advection schemes, coriolis, closures, etc).

[simone-silvestri]

That is a problem of the OrthogonalSphericalShellGrid in general.
The error in this issue is not a bug, but an argument error:

Oceananigans.jl/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl

Line 227 in 13bf409

validate_momentum_advection(momentum_advection, grid::OrthogonalSphericalShellGrid) = error("$(typeof(momentum_advection)) is not supported with $(typeof(grid))")

We can change that if you think it's not useful.

[glwagner]

It is certainly a bug, it should work and does not.

[simone-silvestri]

Are you suggesting to change the default momentum advection for the hydrostatic model?

[simone-silvestri]

Or maybe to adapt the default based on the incoming grid?

[glwagner]

I didn't make a suggestion except the extremely Holy Principle that the default absolutely must work. Otherwise the code is just broken and honestly, when I encounter this in other packages I often decide not to use them.

It wouldn't be unreasonable to use VectorInvariant as the default for the hydrostatic model honestly.

But if we want to have a flux form default then yes it has to depend on the grid to satisfy Holy Principle 1.