Supporting non-zero or time-dependent wall-normal velocities

[glwagner]

Wall-normal velocities can depend on model_fields:

Oceananigans.jl/src/BoundaryConditions/fill_halo_regions_normal_flow.jl

Lines 15 to 18 in 5aafe8e

	@kernel function set_top_bottom_w_velocity!(w, k_boundary, bc, grid, clock, model_fields)
	i, j = @index(Global, NTuple)
	@inbounds w[i, j, k_boundary] = getbc(bc, i, j, grid, clock, model_fields)
	end

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

Oceananigans.jl/src/TimeSteppers/pressure_correction.jl

Lines 6 to 10 in 5aafe8e

	function calculate_pressure_correction!(model, Δt)
	fill_halo_regions!(model.velocities, model.architecture, model.clock, fields(model))
	solve_for_pressure!(model.pressures.pNHS, model.pressure_solver, model.architecture, model.grid, Δt, model.velocities)

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:

Oceananigans.jl/src/TimeSteppers/runge_kutta_3.jl

Line 124 in 5aafe8e

substep_velocities_kernel! = rk3_substep_velocities!(device(model.architecture), workgroup, worksize)

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:

Oceananigans.jl/src/TimeSteppers/runge_kutta_3.jl

Lines 178 to 186 in 5aafe8e

	@kernel function rk3_substep_velocities!(U, Δt, γⁿ, ζⁿ, Gⁿ, G⁻)
	i, j, k = @index(Global, NTuple)
	@inbounds begin
	U.u[i, j, k] += Δt * (γⁿ * Gⁿ.u[i, j, k] + ζⁿ * G⁻.u[i, j, k])
	U.v[i, j, k] += Δt * (γⁿ * Gⁿ.v[i, j, k] + ζⁿ * G⁻.v[i, j, k])
	U.w[i, j, k] += Δt * (γⁿ * Gⁿ.w[i, j, k] + ζⁿ * G⁻.w[i, j, k])
	end
	end

[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)

[glwagner]

I changed the name of this issue; the real point here is that we can't support flow through a Bounded direction. If that flow is time-dependent, there are some details we need to fix regarding the pressure solve.

[glwagner]

This is a hard problem so I'm going to close this since it's not a current priority.