Are `GradientBoundaryCondition`s set by the outward normal component or the coordinate-aligned component?

[hdrake]

The documentation about how GradientBoundaryConditions are implemented suggests that the user specifies a value $\gamma$ as the boundary-normal component of the gradient vector (with an implied outward-normal orientation), $\gamma = \hat{\mathbf{n}} \cdot \nabla c = \gamma$.

By contrast, my reading of the source code is that the GradientBoundaryConditions are used to specify $\hat{\mathbf{x}}_{i} \cdot \nabla c = \gamma$.

Because the outward unit normal vector at the right boundaries point in the same directions as the coordinate unit vectors, the two implementations give the same answer anyway: $\hat{\mathbf{x}}_{i} \cdot \nabla c = \hat{\mathbf{n}} \cdot \nabla c$.

@inline function right_gradient(i, j, k, ibg, κ, Δ, bc::VBC, c, clock, fields)
    cᵇ = getbc(bc, i, j, k, ibg, clock, fields)
    cⁱʲᵏ = @inbounds c[i, j, k]
    return 2 * (cᵇ - cⁱʲᵏ) / Δ
end

@inline function _east_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δx(index_right(i, LX), j, k, ibg, LX, LY, LZ)
    κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock)
    ∇c = right_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

The problem arises for the left boundary, where the gradient is implemented with as $\hat{\mathbf{x}}_{i} \cdot \nabla c = \gamma$, which has the opposite sign of $\hat{\mathbf{n}} \cdot \nabla c$ because the outward vector at the left boundaries point in the direction of $-\hat{\mathbf{x}}_{i}$.

@inline function left_gradient(i, j, k, ibg, κ, Δ, bc::VBC, c, clock, fields)
    cᵇ = getbc(bc, i, j, k, ibg, clock, fields)
    cⁱʲᵏ = @inbounds c[i, j, k]
    return 2 * (cⁱʲᵏ - cᵇ) / Δ
end

@inline function _west_ib_flux(i, j, k, ibg, bc::VBCorGBC, (LX, LY, LZ), c, closure::ASD, K, id, clock, fields)
    Δ = Δx(index_left(i, LX), j, k, ibg, LX, LY, LZ)
    κ = h_diffusivity(i, j, k, ibg, flip(LX), LY, LZ, closure, K, id, clock)
    ∇c = left_gradient(i, j, k, ibg, κ, Δ, bc, c, clock, fields)
    return - κ * ∇c
end

A similar subtlety applies to the sign convention of user-specified fluxes, as explained in the documentation:

and commented in a relevant part of the source code:

##### Very Important Note: For FluxBoundaryCondition,
##### we assume fluxes are directed along the "inward-facing normal".
##### For example, east_immersed_flux = - user_flux.
##### With this convention, positive fluxes _increase_ boundary-adjacent
##### cell values, and negative fluxes _decrease_ them.
#####

[hdrake]

Contributing to my confusion is that I think the documentation has a typo in that there is a sign error between equations (3) and (4) of the description of the implementation of the Gradient Boundary Condition. I think $c{i,j,0}$ and $c{i,j,1}$ should maybe be swapped in eq. (3)?