Having a domain that's periodic in the "gravity" direction doesn't work for `gravity_unit_vector = NegativeZDirection()`

[tomchor]

@zhazorken and I were trying to set up a toy model where the direction where gravity is acting is periodic and there's some sort of buoyancy source in the domain, leading eventually to a plume that's infinite in the direction of gravity.

The code below presents such a model. Here I'm leaving the gravity as is default (i.e. pointing to gravity_unit_vector = NegativeZDirection()), I'm setting up a small buoyancy source at the west wall, and I'm including a bit of noise in the same west wall to kick off a plume.

using Oceananigans
using Oceananigans.Units

grid = RectilinearGrid(topology = (Bounded, Flat, Periodic), size = (16, 16),
                       x = (0, 1), z = (0, 1))

buoyancy = BuoyancyTracer()
model = NonhydrostaticModel(; grid, tracers = (:b),
                            buoyancy = buoyancy,
                            boundary_conditions = (; b = FieldBoundaryConditions(west=FluxBoundaryCondition(5e-9))),);

noise(x, y, z) = 1e-3 * randn() * exp(-(10x)^2/grid.Lx^2)
set!(model, u=noise, w=noise)

simulation = Simulation(model, Δt=1, stop_time=10minutes);

wizard = TimeStepWizard(cfl=0.8, min_Δt=0.001seconds)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(2))

outputs = (; model.tracers.b, model.velocities.u, model.velocities.w,
           pHY = model.pressures.pHY′, pNH = model.pressures.pNHS,
           p = sum(model.pressures),)
output_filename = "2d-zper"
simulation.output_writers[:fields] = NetCDFOutputWriter(model, outputs, 
                                                        schedule = TimeInterval(10seconds),
                                                        filename = output_filename,
                                                        overwrite_existing = true)
run!(simulation)

@info "Begin plotting"
using Rasters
filename = simulation.output_writers[:fields].filepath
ds = RasterStack(filename)

using GLMakie

set_theme!(Theme(fontsize = 20))
fig = Figure()

kwargs = (xlabel="x", ylabel="z", height=150, width=250)
ax1 = Axis(fig[2, 1]; title = "w", kwargs...);
ax2 = Axis(fig[2, 2]; title = "b", kwargs...);
ax3 = Axis(fig[2, 3]; title = "pHY", kwargs...);
ax4 = Axis(fig[2, 4]; title = "pNH", kwargs...);
ax5 = Axis(fig[2, 5]; title = "total pressure", kwargs...);

# Next we use `Observable`s to lift the values and plot heatmaps and their colorbars

n = Observable(1)

speed_magnitude = 2e-2
using Statistics
temp_magnitude = max(std(ds.b), 1e-12)

wₙ = @lift ds.w[Ti=$n, yC=Near(0)]
hm1 = heatmap!(ax1, wₙ; colormap = :balance, colorrange=(-speed_magnitude, speed_magnitude))
Colorbar(fig[3, 1], hm1, vertical=false, height=8, ticklabelsize=14)

bₙ = @lift ds.b[Ti=$n, yC=Near(0)]
hm2 = heatmap!(ax2, bₙ; colormap = :inferno)
Colorbar(fig[3, 2], hm2, vertical=false, height=8, ticklabelsize=14);

pₙ = @lift ds.pHY[Ti=$n, yC=Near(0)]
hm3 = heatmap!(ax3, pₙ; colormap = :viridis)
Colorbar(fig[3, 3], hm3, vertical=false, height=8, ticklabelsize=14);

pₙ = @lift ds.pNH[Ti=$n, yC=Near(0)]
hm4 = heatmap!(ax4, pₙ; colormap = :viridis)
Colorbar(fig[3, 4], hm4, vertical=false, height=8, ticklabelsize=14);

pₙ = @lift ds.p[Ti=$n, yC=Near(0)]
hm5 = heatmap!(ax5, pₙ; colormap = :viridis)
Colorbar(fig[3, 5], hm5, vertical=false, height=8, ticklabelsize=14);

times = dims(ds, :Ti)
title = @lift "Time = " * string(prettytime(times[$n]))
fig[1, 1:5] = Label(fig, title, fontsize=24, tellwidth=false);

resize_to_layout!(fig)
@info "Animating..."
record(fig, filename * ".mp4", 1:length(times), framerate=10) do i
       n[] = i
end

This simulation runs for 10 minutes and then plots this (on a 32x32 grid):

https://github.com/CliMA/Oceananigans.jl/assets/13205162/e4fc9a13-f14f-4c6d-bda6-0d4aae79ac4c

So basically there's a discontinuity in the vertical direction's pressure, which causes an artificial vertical pressure gradient, causing a spurious flow at the top and bottom, even those aren't supposed to be boundaries since the grid is periodic in z. It seems like this comes from the fact that the hydrostatic pressure isn't periodic (since it always comes from a vertical b integral I think), but the nonhydrostatic pressure is periodic. This in turn leads the total pressure to have a discontinuity. In fact if we run the exact same configuration but make both the gravitational direction and the periodic direction the x direction (code here) we have a plume looks and behaves as expected and has no discontinuities:

https://github.com/CliMA/Oceananigans.jl/assets/13205162/97ab7df9-4b64-4f65-9433-4eb9a5e7fe22

I'm not sure if this is something that needs fixing or if it's just a consequence of the hydrostatic decomposition, which always assumes gravity is in the NegativeZDirection() and that z is Bounded. Hence I'm not sure this is better left as an issue or a discussion, so please lmk if I should move this to a discussion.

Also, if I run the example above (where gravity is in the z direction) using https://github.com/CliMA/Oceananigans.jl/pull/3080, which gets rid of the hydrostatic pressure separation, then I get what I believe to be the correct behavior (all pressures in the animation below are actually total pressure):

https://github.com/CliMA/Oceananigans.jl/assets/13205162/60ef470d-22d7-429a-9029-adb7db5687ba

Which also points to the hydrostatic decomposition being the issue here.

[glwagner]

The movies don't play for me. But it seems like a vertically-integrated component (a constant) is missing from the hydrostatic pressure anomaly in the case of a vertically-periodic domain?

https://github.com/CliMA/Oceananigans.jl/blob/e16cdc6cfb67703df9e29368017868331f68b1c0/src/Models/NonhydrostaticModels/update_hydrostatic_pressure.jl#L12-L20

Actually, maybe it's the other way around -- the vertically-integrated component needs to be subtracted when the domain is vertically periodic? This is effectively what occurs here for example:

https://github.com/CliMA/Oceananigans.jl/blob/e16cdc6cfb67703df9e29368017868331f68b1c0/src/Solvers/fft_based_poisson_solver.jl#L111

However this is not enforced for the hydrostatic pressure.

The tridiagonal solvers are also incorrect for vertically-periodic domains, I think (well, now that we have x- and y- tridiagonal solvers I believe they are also incorrect for x- and y- periodic if using an x-tridiagonal or y-tridiagonal solver).

Either way it does seem like the simplest solution is to eliminate the pressure decomposition.

Interested what @simone-silvestri and @xkykai think.

[tomchor]

The movies don't play for me.

Weird, they're playing for me on two different browsers... not sure what to do about that.

The tridiagonal solvers are also incorrect for vertically-periodic domains, I think (well, now that we have x- and y- tridiagonal solvers I believe they are also incorrect for x- and y- periodic if using an x-tridiagonal or y-tridiagonal solver).

True, but just to be clear, these simulations don't use any stretched grid direction, so I believe they don't use the tridiagonal solver, correct?

[glwagner]

The movies don't play for me.

Weird, they're playing for me on two different browsers... not sure what to do about that.

There's nothing for you to do, I was on plane internet! I can see it now. I was just letting you know...

The tridiagonal solvers are also incorrect for vertically-periodic domains, I think (well, now that we have x- and y- tridiagonal solvers I believe they are also incorrect for x- and y- periodic if using an x-tridiagonal or y-tridiagonal solver).

True, but just to be clear, these simulations don't use any stretched grid direction, so I believe they don't use the tridiagonal solver, correct?

Correct --- I just wanted to issue that warning in case there was more interest in vertically-periodic simulations (we basically don't test that situation, but it wouldn't be unreasonable to work on that).