Oceananigans.jl: Refilling halo regions prior to pressure solve can produce incorrect wall-normal velocities

[glwagner]

Wall-normal velocities can depend on model_fields:

https://github.com/CliMA/Oceananigans.jl/blob/5aafe8ee1d3e49a53906e7225a01fc18f1a5f165/src/BoundaryConditions/fill_halo_regions_normal_flow.jl#L15-L18

and wall-normal velocities are updated _after_ an RK3 substep, but _before_ the pressure solve:

https://github.com/CliMA/Oceananigans.jl/blob/5aafe8ee1d3e49a53906e7225a01fc18f1a5f165/src/TimeSteppers/pressure_correction.jl#L6-L10

Thus for some problems the wall-normal velocity fields are updated based on the predictor model fields (both the predictor velocity and the updated tracer fields) that result from an RK3 substep.

This devious bug can be avoided simply by _not updating wall-normal velocity components on the boundary_ in the RK3 substep by changing the indexing in the rk3 substep as well as the worksize here:

https://github.com/CliMA/Oceananigans.jl/blob/5aafe8ee1d3e49a53906e7225a01fc18f1a5f165/src/TimeSteppers/runge_kutta_3.jl#L124

Then we don't have to fill halo regions before performing the pressure correction. The resulting algorithm is both more correct and computationally less expensive.

Note that doing this could require a bit of gymnastics to get the indexing right in the rk3 substep kernel:

https://github.com/CliMA/Oceananigans.jl/blob/5aafe8ee1d3e49a53906e7225a01fc18f1a5f165/src/TimeSteppers/runge_kutta_3.jl#L178-L186

[glwagner]

Some context: @kburns and I encountered this issue when trying to use an Ekman pumping boundary condition, which prescribes

w = k * om

at the bottom boundary, where om is the vertical vorticity and k is a constant:

""" Returns the vertical component of vorticity. """
@inline function ζᶠᶠᶜ(i, j, k, grid, u, v)
    ∂x_v = δxᶠᵃᵃ(i, j, k, grid, v) / grid.Δx
    ∂y_u = δyᵃᶠᵃ(i, j, k, grid, u) / grid.Δy
    return ∂x_v - ∂y_u
end

@inline pumping_velocity(ζ, κᵖ) = κᵖ * ζ 

@inline function pumping_velocity(i, j, grid, clock, model_fields, κᵖ) 
    ζʷ = ℑxyᶜᶜᵃ(i, j, 1, grid, ζᶠᶠᶜ, model_fields.u, model_fields.v)
    return pumping_velocity(ζʷ, κᵖ) 
end

pumping_bc = BoundaryCondition(NormalFlow, pumping_velocity, discrete_form=true, parameters=kp)