HydrostaticFreeSurface model on a grid with Flat z

[navidcy]

Does this make sense (physically)?

               _
   _       _ _(_)_     |  Documentation: https://docs.julialang.org
  (_)     | (_) (_)    |
   _ _   _| |_  __ _   |  Type "?" for help, "]?" for Pkg help.
  | | | | | | |/ _` |  |
  | | |_| | | | (_| |  |  Version 1.6.5 (2021-12-19)
 _/ |\__'_|_|_|\__'_|  |
|__/                   |

julia> using Oceananigans
[ Info: Oceananigans will use 12 threads

julia> grid = RectilinearGrid(topology=(Periodic, Bounded, Flat), size=(1, 1), extent=(1, 2))
1×1×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 1×1×0 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=1.0
├── Bounded  y ∈ [0.0, 2.0] regularly spaced with Δy=2.0
└── Flat z

julia> model = HydrostaticFreeSurfaceModel(grid=grid)
HydrostaticFreeSurfaceModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 1×1×0 halo
├── tracers: (:T, :S)
├── closure: Nothing
├── buoyancy: Buoyancy{SeawaterBuoyancy{Float64, LinearEquationOfState{Float64}, Nothing, Nothing}, Oceananigans.Grids.ZDirection}
└── coriolis: Nothing

julia> keys(fields(model))
(:u, :v, :η, :T, :S)

I guess not. Perhaps the HydrostaticFreeSurfaceModel should print a warning to users that try to do that?

[francispoulin]

Good idea. Maybe even suggest the user try ShallowWaterModel, as that is a better way to try this configuration? ;)

[glwagner]

Hmm... we don't support Flat z right now with ExplicitFreeSurface because we require two grid points in the vertical velocity (one grid point at the "bottom", one at the "top" coincident with the location of the free surface).

I'm not sure about ImplicitFreeSurface. Can you check?

We might be able to refactor the algorithm to generally support hydrostatic free surface simulations with Flat z. Maybe we just need to generalize the free surface tendency for ExplicitFreeSurface?

On a regular RectilinearGrid, NonhydrostaticModel with Flat z is hydrostatic and has equivalent physics (because the free surface physics are linearized in HydrostaticFreeSurfaceModel).

Before generalizing anything we could add an error message to the HydrostaticFreeSurfaceModel constructor to point this out (and maybe link to this issue).

ShallowWaterModel has a nonlinear free surface and thus finite "Mach number" (ratio between gravity wave speed and flow speed). So I agree with @francispoulin
that ShallowWaterModel is more general and might be a nice choice for some applications. I think it still lacks a diffusive flux implementation though?

[francispoulin]

Yes, I agree that ShallowWaterModel needs diffusive flux. I think it can borrow a lot from hydrostatic model, but there's still an issue as to how to write them properly. I'll think about it and try and start something in the next week. When we decide exactly what to do, it should not be difficult to do. Thanks for the remindier @glwagner !

[navidcy]

Hmm... we don't support Flat z right now with ExplicitFreeSurface because we require two grid points in the vertical velocity (one grid point at the "bottom", one at the "top" coincident with the location of the free surface).

well, using a314f1f I get

julia> using Oceananigans

julia> grid
1×1×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 1×1×0 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=1.0
├── Bounded  y ∈ [0.0, 2.0] regularly spaced with Δy=2.0
└── Flat z

julia> model = HydrostaticFreeSurfaceModel(grid=grid)
HydrostaticFreeSurfaceModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 1×1×0 halo
├── tracers: (:T, :S)
├── closure: Nothing
├── buoyancy: Buoyancy{SeawaterBuoyancy{Float64, LinearEquationOfState{Float64}, Nothing, Nothing}, Oceananigans.Grids.ZDirection}
└── coriolis: Nothing

julia> simulation = Simulation(model, Δt=1.0, stop_iteration=2)
Simulation{typename(HydrostaticFreeSurfaceModel){typename(CPU), Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step: 1 second
├── Elapsed wall time: 0 seconds
├── Stop time: Inf years
├── Stop iteration : 2.0
├── Wall time limit: Inf
├── Callbacks: typename(OrderedCollections.OrderedDict) with 4 entries:
│   ├── stop_time_exceeded => typename(Callback)
│   ├── stop_iteration_exceeded => typename(Callback)
│   ├── wall_time_limit_exceeded => typename(Callback)
│   └── nan_checker => typename(Callback)
├── Output writers: typename(OrderedCollections.OrderedDict) with no entries
└── Diagnostics: typename(OrderedCollections.OrderedDict) with no entries

julia> run!(simulation)
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (169.437 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (23.728 seconds).
[ Info: Simulation is stopping. Model iteration 2 has hit or exceeded simulation stop iteration 2.

so something's running and not complaining!

I'm not sure about ImplicitFreeSurface. Can you check?

Same! Runs again no complaining :)

julia> using Oceananigans

julia> grid
1×1×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 1×1×0 halo
├── Periodic x ∈ [0.0, 1.0) regularly spaced with Δx=1.0
├── Bounded  y ∈ [0.0, 2.0] regularly spaced with Δy=2.0
└── Flat z

julia> model = HydrostaticFreeSurfaceModel(grid=grid, free_surface = ImplicitFreeSurface())
HydrostaticFreeSurfaceModel{CPU, Float64}(time = 0 seconds, iteration = 0)
├── grid: 1×1×1 RectilinearGrid{Float64, Periodic, Bounded, Flat} on CPU with 1×1×0 halo
├── tracers: (:T, :S)
├── closure: Nothing
├── buoyancy: Buoyancy{SeawaterBuoyancy{Float64, LinearEquationOfState{Float64}, Nothing, Nothing}, Oceananigans.Grids.ZDirection}
└── coriolis: Nothing

julia> simulation = Simulation(model, Δt=1.0, stop_iteration=2)
Simulation{typename(HydrostaticFreeSurfaceModel){typename(CPU), Float64}}
├── Model clock: time = 0 seconds, iteration = 0
├── Next time step: 1 second
├── Elapsed wall time: 0 seconds
├── Stop time: Inf years
├── Stop iteration : 2.0
├── Wall time limit: Inf
├── Callbacks: typename(OrderedCollections.OrderedDict) with 4 entries:
│   ├── stop_time_exceeded => typename(Callback)
│   ├── stop_iteration_exceeded => typename(Callback)
│   ├── wall_time_limit_exceeded => typename(Callback)
│   └── nan_checker => typename(Callback)
├── Output writers: typename(OrderedCollections.OrderedDict) with no entries
└── Diagnostics: typename(OrderedCollections.OrderedDict) with no entries

julia> run!(simulation)
[ Info: Initializing simulation...
[ Info:     ... simulation initialization complete (198.852 ms)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.013 minutes).
[ Info: Simulation is stopping. Model iteration 2 has hit or exceeded simulation stop iteration 2.

Before generalizing anything we could add an error message to the HydrostaticFreeSurfaceModel constructor to point this out (and maybe link to this issue).

yeah!

[glwagner]

Is this closed?

[glwagner]

Wait nevermind, I know that its closed.