Nonhydrostatic pressure correction broken?

[jagoosw]

Hi all,

I've just come across a new problem where the velocity fields in non-hydrostatic models is abruptly becoming the average of the whole field (I think, all positions become one value anyway). This appears to be occurring in the pressure correction step.

For example:

julia> using Oceananigans
julia> grid = RectilinearGrid(size = (128, 128, 128), extent = (1000, 1000, 1000));
julia> model = NonhydrostaticModel(; grid);
julia> uᵢ(x, y, z) = (1-tanh((x-500)/10));
julia> set!(model, u = uᵢ)

julia> model.velocities.u
128×128×128 Field{Face, Center, Center} on RectilinearGrid on CPU
├── grid: 128×128×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 134×134×134 OffsetArray(::Array{Float64, 3}, -2:131, -2:131, -2:131) with eltype Float64 with indices -2:131×-2:131×-2:131
    └── max=1.00781, min=1.00781, mean=1.00781

julia> set!(model, u = uᵢ, enforce_incompressibility=false)

julia> model.velocities.u
128×128×128 Field{Face, Center, Center} on RectilinearGrid on CPU
├── grid: 128×128×128 RectilinearGrid{Float64, Periodic, Periodic, Bounded} on CPU with 3×3×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Periodic, north: Periodic, bottom: ZeroFlux, top: ZeroFlux, immersed: ZeroFlux
└── data: 134×134×134 OffsetArray(::Array{Float64, 3}, -2:131, -2:131, -2:131) with eltype Float64 with indices -2:131×-2:131×-2:131
    └── max=2.0, min=0.0, mean=1.00781

I initially spotted this issue after time-stepping a model and the same thing occurs, presumably after the pressure correction, so this isn't isolated to set!.

This is occurring in all of my julia projects where Oceananigans is in various versions and I know Nonhydrostatic models have worked correctly in before, downgrading is also not working so perhaps this is similar to #3381 with a broken dependency.

I will try updating to 1.10 to see if that fixes the issue.

[jagoosw]

This is not fixed by julia-1.10 in either 0.90.1 or 0.85

[jagoosw]

Also just to confirm this does not happen with the hydrostatic free surface model

[glwagner]

This velocity field:

uᵢ(x, y, z) = (1-tanh((x-500)/10))

is compressible, because uᵢ depends on x. When we project out the compressible component (the part that depends on x), we are left with just u=1.

What do we expect to happen in this case?

[glwagner]

Also just to confirm this does not happen with the hydrostatic free surface model

For the hydrostatic free surface model, w should be calculated using continuity; in other words you should find that dw/dz = - du/dx after using set! (or just more generally, there will be a non-zero vertical velocity).

[johnryantaylor]

@glwagner is right about this, for me it works as expected if you use an incompressible initial condition, e.g.
uᵢ(x, y, z) = (1-tanh((y-500)/10))
@jagoosw, I would be interested to know more about the context where you saw this during timestepping.

[jagoosw]

Oh sorry I thought I must have missed something! Thank you