Implement an implicit free surface solver in the NonhydrostaticModel

[glwagner]

Closes #3946

[glwagner]

So here's how the algorithm changes with an implicit free surface (which is all I'd like to attempt for the time being):

1. Update velocities to obtain first predictor velocities
2. Step forward the free surface with an implicit solve
3. Update the velocities to obtain the "second" predictor velocities
4. Compute the pressure correction

I think the simplest way to implement this is therefore to add the implicit free surface solve as a precursor to the pressure correction.

[glwagner]

Okay, this script:

using Oceananigans
using Oceananigans.Models.HydrostaticFreeSurfaceModels: ImplicitFreeSurface
using GLMakie

grid = RectilinearGrid(size=(128, 32), halo=(4, 4), x=(-5, 5), z=(0, 1), topology=(Bounded, Flat, Bounded))

mountain(x) = (x - 3) / 2
grid = ImmersedBoundaryGrid(grid, GridFittedBottom(mountain))

Fu(x, z, t) = sin(t)
free_surface = ImplicitFreeSurface(gravitational_acceleration=10)
model = NonhydrostaticModel(; grid, free_surface, advection=WENO(order=5), forcing=(; u=Fu))

simulation = Simulation(model, Δt=0.1, stop_time=20*2π)
conjure_time_step_wizard!(simulation, cfl=0.7)

progress(sim) = @info string(iteration(sim), ": ", time(sim))
add_callback!(simulation, progress, IterationInterval(100))

ow = JLD2OutputWriter(model, merge(model.velocities, (; η=model.free_surface.η)),
                      filename = "nonhydrostatic_internal_tide.jld2",
                      schedule = TimeInterval(0.1),
                      overwrite_existing = true)

simulation.output_writers[:jld2] = ow

run!(simulation)

fig = Figure()

axη = Axis(fig[1, 1], xlabel="x", ylabel="Free surface \n displacement")
axw = Axis(fig[2, 1], xlabel="x", ylabel="Surface vertical velocity")
axu = Axis(fig[3, 1], xlabel="x", ylabel="z")

ut = FieldTimeSeries("nonhydrostatic_internal_tide.jld2", "u")
wt = FieldTimeSeries("nonhydrostatic_internal_tide.jld2", "w")
ηt = FieldTimeSeries("nonhydrostatic_internal_tide.jld2", "η")
Nt = length(wt)

slider = Slider(fig[4, 1], range=1:Nt, startvalue=1)
n = slider.value
Nz = size(ut.grid, 3)

u = @lift ut[$n]
η = @lift interior(ηt[$n], :, 1, 1)
w = @lift interior(wt[$n], :, 1, Nz+1)
x = xnodes(wt)

ulim = maximum(abs, ut) * 3/4

lines!(axη, x, η)
lines!(axw, x, w)
heatmap!(axu, u)

ylims!(axη, -0.1, 0.1)
ylims!(axw, -0.01, 0.01)

record(fig, "nonhydrostatic_internal_tide.mp4", 1:Nt) do nn
    @info "Drawing frame $nn of $Nt..."
    n[] = nn
end

produces this movie

nonhydrostatic_internal_tide.mp4

some weird grid artifacts in there but maybe higher resolution will help with that.

[glwagner]

Here's the script:

https://github.com/CliMA/Oceananigans.jl/blob/glw/nh-free-surface/validation/implicit_free_surface/nonhydrostatic_internal_tide.jl

[glwagner]

@shriyafruitwala let me know if this code works for the problem you are interested in.

The implementation is fairly clean, but there are a few things we could consider before merging, like tests, some validation in the constructor.

@simone-silvestri I feel this code exposes some messiness with the peformance optimization stuff regarding kernel parameters. Please check it over to make sure what I did will work and add any tests that may be missing...

[simone-silvestri]

I am unsure of the algorithm, but wouldn't this replace the w-velocity that should be prognostic?

[glwagner]

For sure. That's what is written in MITgcm docs. @jm-c can you confirm?

[simone-silvestri]

@simone-silvestri I feel this code exposes some messiness with the peformance optimization stuff regarding kernel parameters. Please check it over to make sure what I did will work and add any tests that may be missing...

it seems that you are missing required_halo_size_x and required_halo_size_y in the new file src/Models/buffer_tendency_kernel_parameters.jl. I can add it