Problems at the periodic boundary

[andrewooish]

I am trying to set up Oceananigans to run buoyant flow in a vertical channel. The left wall (y=-0.5) is hot and the right wall (y=+0.5) is cold. I am implementing the no-slip boundary condition at y=-05 and y=0.5. I have specified buoyancy to be constant value os positive +0.5 at the south wall and -0.5 at the north wall. Due to buoyancy, the flow should go up one wall and come down the other wall. The flow in the streamwise (z) direction (the same direction as gravity g) is periodic, so I have specified periodic boundary condition in the z direction. At this stage, I am only trying to calculate 2d flow (y and z coordinates, x is flat).

The problem I am seeing is that the flow seem to be messed up at the periodic boundaries (z=0 and z=8). I have run for many (100s) of time units. I have chacked the boundary conditions and they all seem to be correct. The flow away from the periodic boundaries (say at z=4) seem to be OK, but close to the periodic boundaries(z=0 and z=8)., they just do not look right. I am attaching code and results (contours of w) with this post.

Any help appreciated. Thank you.

---

grid = RectilinearGrid(CPU(),size = (Ny, Nz),
y = y_faces,
z = (0,Lz),
topology = (Flat,Bounded,Periodic))

Pr = 1 # Prandtl number
Ra = 1e5 # Rayleigh number
b★ = 1
ν = sqrt(Pr * b★ * H^3 / Ra) # Laplacian viscosity
κ = ν * Pr # Laplacian diffusivity

g̃ = (0.0,0.0,-1.0)

no_slip_field_bcs = FieldBoundaryConditions(north=ValueBoundaryCondition(0.0),south=ValueBoundaryCondition(0.0));
buoyancy_bc=FieldBoundaryConditions(north=ValueBoundaryCondition(-0.5),south=ValueBoundaryCondition(0.5))

model = NonhydrostaticModel(; grid,
boundary_conditions=(w=no_slip_field_bcs,b=buoyancy_bc),
advection=Centered(order=2),
timestepper = :RungeKutta3,
tracers = :b,
buoyancy = BuoyancyForce(BuoyancyTracer(), gravity_unit_vector=g̃),
closure = ScalarDiffusivity(; ν, κ))

initial_buoyancy(y,z)=0.5-(y+0.5)
width=0.1
initial_w(y,z)=(1-tanh((y-0)/width))/2.0-0.5
initial_v(y,z)=0.0

u, v, w = model.velocities
set!(model, b=initial_buoyancy, w=initial_w, v=initial_v)

Checking boundary conditions

model.velocities.w.boundary_conditions
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: Nothing
├── east: Nothing
├── south: ValueBoundaryCondition: 0.0
├── north: ValueBoundaryCondition: 0.0
├── bottom: PeriodicBoundaryCondition
├── top: PeriodicBoundaryCondition
└── immersed: FluxBoundaryCondition: Nothing

model.velocities.v.boundary_conditions
Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: Nothing
├── east: Nothing
├── south: OpenBoundaryCondition{Nothing}: Nothing
├── north: OpenBoundaryCondition{Nothing}: Nothing
├── bottom: PeriodicBoundaryCondition
├── top: PeriodicBoundaryCondition
└── immersed: FluxBoundaryCondition: Nothing

Oceananigans.FieldBoundaryConditions, with boundary conditions
├── west: Nothing
├── east: Nothing
├── south: ValueBoundaryCondition: 0.5
├── north: ValueBoundaryCondition: -0.5
├── bottom: PeriodicBoundaryCondition
├── top: PeriodicBoundaryCondition
└── immersed: FluxBoundaryCondition: Nothing

---

[glwagner]

I suspect the vertical pressure integral does not correctly implement the Periodic z-direction, so basically, there is a bug. Here's two quick fixes:

1. rotate the problem so that x or y is Periodic by changing gravity_unit_vector (which you have already defined)
2. eliminate the vertical pressure integral by passing hydrostatic_pressure_anomaly=nothing to NonhydrostaticModel

I'm not 100% option 2 will work so it will be good to here.

Something actionable is to throw an error or warning for this configuration.

[andrewooish]

Thank you for your respnose. Both items 1 and 2 seem to have worked for version 0.95.5

I had another suggestion by a colleague to try version 0.91.12. In this version, periodic in the z direction (without the hydrostatic_pressure_anomaly=nothing) seem to work. So it looks like periodic in the z direction is a bug in version 0.95.5

[glwagner]

Isn't that merely because of a change to the default setting for hydrostatic_pressure_anomaly? Or is there more to it?