Allow users to mask output

[navidcy]

Resolves #3061

Questions:

Does it make sense to mask whenwith_halos = true? With with_halos = true implies that the users can use the output to compute derivative fields etc; will masking obscure this?

[navidcy]

mask_immersed_field! was renamed to mask_immersed!; this is a breaking change, right?

[navidcy]

I was trying to run validation/immersed_boundaries/flow_over_hills.jl and all simulations except the ref one are NaNing pretty quickly. On main it works!

[glwagner]

Resolves #3061

Questions:

Does it make sense to mask whenwith_halos = true? With with_halos = true implies that the users can use the output to compute derivative fields etc; will masking obscure this?

No, because derivatives are "boundary aware" and thus don't touch the regions that are masked.

[navidcy]

I'm trying to finish up this.

@glwagner, is there a reason you here put a function to return the kernel? Usually the kernel is launched and then the function returns nothing.

https://github.com/CliMA/Oceananigans.jl/blame/d7d146bc59f966d65eb66a875c16679bef365130/src/ImmersedBoundaries/mask_immersed_field.jl#L44

(cc @simone-silvestri)

[navidcy]

@glwagner I think this won't work because what happens is that, e.g., if u is the output to be saved then model.velocities.u is masked and then it includes NaNs at the next step so the whole state NaN...

At least that's what I found happening and I boil it down to this line here masking and then returning the parent of the field. Am I right?

[glwagner]

Ok that makes sense. But the real problem is how we define masking. Consider: the mask we apply to a field for output should not affect time-stepping. If it does, there's a problem. So let's fix that before hacking around it.

The issue is velocity components I think. The "mask" is defined to cover the points that are on the boundary, not just inside it. Thus masking the velocity to 0 is how we ensure impenetrability.

The confusion is therefore: are velocity nodes on the boundary immersed or not?

I think we should probably define this points as not immersed. Therefore to fix the problem we can:

1. Implement a new function that satisfies impenetrability. We can start by doing this just for immersed grids.
2. Define "mask_immersed" as something that only affects cells that are fully immersed. This has no effect on fully centered fields, but for staggered fields this implies a change: cells are only masked if they are not adjacent to active cells (I think the distinction between inactive_node and peripheral_node can be used for this).

[navidcy]

Yes, I was thinking of point 2 you make and modified masked_immersed! to use inactive_node (just to test the hypothesis) but that didn't work either.

[simone-silvestri]

There might also be a problem with the implicit solver to fix. Our strategy right now is to zero out elements of the matrix that correspond to immersed cells. Unfortunately 0 * NaN = NaN so any NaN in the immersed boundary will be propagated immediately in the whole column

[navidcy]

Good point.

Can't we just apply the mask to the saved output?

[glwagner]

Oh that "implicit solver".

I guess that is the code here:

Oceananigans.jl/src/ImmersedBoundaries/grid_fitted_immersed_boundaries.jl

Lines 126 to 148 in a1f4f4e

	@inline immersed_ivd_peripheral_node(i, j, k, ibg, LX, LY, ::Center) = immersed_peripheral_node(i, j, k+1, ibg, LX, LY, Face())
	@inline immersed_ivd_peripheral_node(i, j, k, ibg, LX, LY, ::Face) = immersed_peripheral_node(i, j, k, ibg, LX, LY, Center())
	# Extend the upper and lower diagonal functions of the batched tridiagonal solver
	for location in (:upper_, :lower_)
	immersed_func = Symbol(:immersed_ivd_, location, :diagonal)
	ordinary_func = Symbol(:ivd_ , location, :diagonal)
	@eval begin
	# Disambiguation
	@inline $ordinary_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz::Face, clock, Δt, κz) =
	$immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)
	@inline $ordinary_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz::Center, clock, Δt, κz) =
	$immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)
	@inline function $immersed_func(i, j, k, ibg::GFIBG, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz)
	return ifelse(immersed_ivd_peripheral_node(i, j, k, ibg, ℓx, ℓy, ℓz),
	zero(eltype(ibg.underlying_grid)),
	$ordinary_func(i, j, k, ibg.underlying_grid, closure, K, id, ℓx, ℓy, ℓz, clock, Δt, κz))
	end
	end
	end

Seems like we would need

Oceananigans.jl/src/Solvers/batched_tridiagonal_solver.jl

Line 179 in a1f4f4e

f₁ = get_coefficient(f, i, j, 1, grid, p, tridiagonal_direction, args...)

to return 0 for immersed_cell.

To accomplish this without explicit masking we might have to wrap field in a ConditionalOperation:

Oceananigans.jl/src/TurbulenceClosures/vertically_implicit_diffusion_solver.jl

Lines 197 to 199 in a1f4f4e

	return solve!(field, implicit_solver, field,
	# ivd_*_diagonal gets called with these args after (i, j, k, grid):
	vi_closure, vi_diffusivity_fields, tracer_index, map(ℓ -> ℓ(), loc)..., clock, Δt, κz)

hmm...

Another possibility is to move mask_immersed! from update_state!; for example if we masked fields before ab2_step! then I think we'd be ok:

Oceananigans.jl/src/TimeSteppers/quasi_adams_bashforth_2.jl

Line 93 in a1f4f4e

ab2_step!(model, Δt, χ) # full step for tracers, fractional step for velocities.

Not sure if that would affect non-masked output. Hopefully not.

[navidcy]

Can't we just apply the mask to the saved output?

What do you mean?

I mean only mask the data written in the jld2/nc file, not the actual data in the u, v, w, fields, etc.

[glwagner]

Data has to be written from an existing array. So, masking data being saved to disk without touching the "actual" fields requires allocating a temporary array. Note, in the GPU case we already do this on the CPU; however we would then have to mask on the CPU which could slow things down for large fields.

[navidcy]

OK, I see!

[glwagner]

But maybe what you're suggesting could be achieved by digging into JLD2's array-saving code. I'm not sure about NetCDF. I took a quick look and couldn't figure it out but such masking during i/o might be possible...