Is HydrostaticFreeSurfaceModel supposed to support a rigid lid?

[ali-ramadhan]

I noticed a file called rigid_lid.jl exists (https://github.com/CliMA/Oceananigans.jl/blob/b3be950ef6a1f2139c2f65e083f20f1f7958ea55/src/Models/HydrostaticFreeSurfaceModels/rigid_lid.jl) so I was curious whether HydrostaticFreeSurfaceModel can run with a free surface, presumably by passing free_surface = nothing.

It does error (see MWE below) but with a couple of extra function definitions (see far below) it seems to time step. Are these fixes enough for a working rigid lid? If not, does it make sense to remove the rigid_lid.jl file?

Note: The documentation also mentions a rigid lid here: https://clima.github.io/OceananigansDocumentation/stable/numerical_implementation/elliptic_solvers/#Rigid-lid-pressure-operator

---

Minimal working example to reproduce the error:

using Oceananigans

grid = LatitudeLongitudeGrid(size=(10, 10, 10), longitude=(0, 1), latitude=(0, 1), z=(-1, 0))
model = HydrostaticFreeSurfaceModel(; grid, free_surface=nothing)
time_step!(model, 1)

Error:

ERROR: MethodError: no method matching materialize_free_surface(::Nothing, ::@NamedTuple{…}, ::LatitudeLongitudeGrid{…})

Closest candidates are:
  materialize_free_surface(::ExplicitFreeSurface{Nothing}, ::Any, ::Any)
   @ Oceananigans ~/.julia/packages/Oceananigans/Hkk5J/src/Models/HydrostaticFreeSurfaceModels/explicit_free_surface.jl:32
  materialize_free_surface(::SplitExplicitFreeSurface, ::Any, ::Any)
   @ Oceananigans ~/.julia/packages/Oceananigans/Hkk5J/src/Models/HydrostaticFreeSurfaceModels/split_explicit_free_surface.jl:114
  materialize_free_surface(::ImplicitFreeSurface{Nothing}, ::Any, ::Any)
   @ Oceananigans ~/.julia/packages/Oceananigans/Hkk5J/src/Models/HydrostaticFreeSurfaceModels/implicit_free_surface.jl:93
  ...

Stacktrace:
 [1] HydrostaticFreeSurfaceModel(; grid::LatitudeLongitudeGrid{…}, clock::Clock{…}, momentum_advection::Centered{…}, tracer_advection::Centered{…}, buoyancy::Nothing, coriolis::Nothing, free_surface::Nothing, 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 ~/.julia/packages/Oceananigans/Hkk5J/src/Models/HydrostaticFreeSurfaceModels/hydrostatic_free_surface_model.jl:177
 [2] top-level scope
   @ REPL[3]:1
Some type information was truncated. Use `show(err)` to see complete types.

---

"Fixes":

using Oceananigans

import Oceananigans.Models.HydrostaticFreeSurfaceModels: materialize_free_surface, compute_free_surface_tendency!

materialize_free_surface(free_surface::Nothing, velocities, grid) = nothing

compute_free_surface_tendency!(grid, model::HydrostaticFreeSurfaceModel{<:Any, <:Any, <:Any, Nothing}, kernel_parameters) = nothing

grid = LatitudeLongitudeGrid(size=(10, 10, 10), longitude=(0, 1), latitude=(0, 1), z=(-1, 0))
model = HydrostaticFreeSurfaceModel(; grid, free_surface=nothing)
time_step!(model, 1)

---

Environment: Oceananigans.jl v0.91.11 and

julia> versioninfo()
Julia Version 1.10.4
Commit 48d4fd48430 (2024-06-04 10:41 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 48 × AMD Ryzen Threadripper 7960X 24-Cores
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-15.0.7 (ORCJIT, znver3)
Threads: 1 default, 0 interactive, 1 GC (on 48 virtual cores)

[simone-silvestri]

Nice. I think, these fixes are enough for a rigid lid to work because pressure_correct_velocities! does nothing by default.
Maybe we can also add a

pressure_correct_velocities!(model::RigidLidHFSM, Δt) = nothing

in the barotropic_pressure_correction.jl file, to be explicit about what is happening in case of a rigid lid, although it is not needed

[ali-ramadhan]

Awesome! I can start a PR to support rigid lids then.

[glwagner]

A rigid lid requires a solver for barotropic pressure and a pressure correction step.

[ali-ramadhan]

Ah is it worth doing something with #3740 then so the model doesn't error with free_surface = nothing or would that just be misleading?

I can also update the PR to get rid of rigid_lid.jl to avoid confusing users into thinking a rigid lid mode exists.

[simone-silvestri]

A rigid lid requires a solver for barotropic pressure and a pressure correction step.

Right, I blundered a bit here. Should we even allow free_surface = nothing then with non prescribed velocities?

[glwagner]

Hmmmmmmmm

Right, if we really want to do stuff with a rigid lid then we need to add another field to HydrostaticFreeSurfaceModel to hold the barotropic pressure (a 2D field in x, y) and get the algorithm working for that... it's definitely possible...

As for free_surface=nothing, this should be ok in fact if we are using PrescribedVelocityFields, but otherwise I think it should error (until we support rigid lids?)