Inconsistency between mask_immersed_field and condition_operand

[glwagner]

Here's an MWE:

using Oceananigans
using Oceananigans.ImmersedBoundaries: mask_immersed_field!

grid = RectilinearGrid(size=(2, 2), x = (0, 4), z=(0, 4), topology = (Periodic, Flat, Bounded))
mask = CenterField(grid)
mask[1, 1, 1] = 1
grid = ImmersedBoundaryGrid(grid, GridFittedBoundary(mask))

c = CenterField(grid)
w = ZFaceField(grid)
set!(c, 1)
set!(w, 1)

@show maximum(c) maximum(w)

as might be expected this returns

maximum(c) = 1.0
maximum(w) = 1.0

but then writing

mask_immersed_field!(c, NaN)
mask_immersed_field!(w, NaN)

@show maximum(c) maximum(w)

leads to

maximum(c) = 1.0
maximum(w) = NaN

The reason is because mask_immersed_field will touch nodes that maximum does not ignore. But we want consistency between the two for many reasons, including analysis.

I believe @simone-silvestri's suggestion is to fix condition_operand to ignore all peripheral_nodes --- not just inactive and immersed_peripheral_node. The difference between the two are points that lie on non-immersed boundaries, and only affect fields with at least one Face location (like w above).