Optimize masking xy planes

[simone-silvestri]

There is a problem at the moment with the mask_immersed_field_xy! function. The implementation leverages a scalar_mask which receives a mask optional argument that seems to create problems with inference on the CPU.
this PR removes the scalar mask function and restructures the masking function such that we remove all the extra allocation (that causes GC time).
For example within the following MWE

Random.seed!(1234)
grid = RectilinearGrid(size=(50, 50, 50), extent=(1, 1, 1))
grid = ImmersedBoundaryGrid(grid, GridFittedBottom((x, y)->-rand()))
c = CenterField(grid)
cx = sum(c, dims=1)
cy = sum(c, dims=2)
cz = sum(c, dims=3)

on main we get

julia> @time mask_immersed_field!(c)
  0.000286 seconds (20 allocations: 10.234 KiB)

julia> @time mask_immersed_field!(cx)
  0.000266 seconds (19 allocations: 10.141 KiB)

julia> @time mask_immersed_field!(cy)
  0.000257 seconds (19 allocations: 10.141 KiB)

julia> @time mask_immersed_field!(cz)
  0.001513 seconds (12.52 k allocations: 6.190 MiB)

In this PR we get

julia> @time mask_immersed_field!(c)
  0.000299 seconds (20 allocations: 10.234 KiB)

julia> @time mask_immersed_field!(cx)
  0.000249 seconds (19 allocations: 10.141 KiB)

julia> @time mask_immersed_field!(cy)
  0.000269 seconds (19 allocations: 10.141 KiB)

julia> @time mask_immersed_field!(cz)
  0.000063 seconds (20 allocations: 10.234 KiB)

This PR also unifies the interface for mask_immersed_field_xy! with the one for mask_immersed_field! by providing a mask keyword argument also for the mask_immersed_field! function

[glwagner]

I would favor not having an interface where one of the arguments is a function. It doesn't really help that much compared to the pattern we usually implement based on types and dispatch. Also I don't see the point of being able to change mask.

[simone-silvestri]

It is the previous pattern we had, this PR does not change it. We can think about a new interface. However, being able to change the mask is quite important to define the behavior of the maks_immersed_field! function.
Think inactive_node vs peripheral_node vs immersed_inactive_node etc...

[glwagner]

This PR adds mask as an argument to mask_immersed_field!. But if we are working to make the interfaces more uniform, I would favor removing mask from mask_immersed_field_xy!.

[glwagner]

Where do we use different "masking functions"?

[glwagner]

I'm really just worried that this is an example of poor design. Obviously, if we never need to develop or change this code, it doesn't matter. But it could lead to problems down the line because it is not well thought out.

The first issue is that this pattern of passing a function to a function, while valid, is atypical in our source code. It's a nice informal "trick" for getting what we want with less code. The other annoying thing to me is that this pattern effectively defines what a "masking function" is. This code will fail unless a mask can be called with i, j, k, grid. Can all "masks" be called that way? This seems to have arisen organically, but we do not formally define masks in this way. The convention is hidden in the ImmersedBoundaries module, quite remote from where the first mask inactive_node is defined. So we are hoping that all future developers will be able to take in the entirety of the code base and intuitively recognize this pattern for masking. I think it's a little sloppy as a practice. Either we should formalize the definition of a mask, or not, but we can't have it both ways just to save ourselves a tiny, negligible amount of effort.

It always seemed like a temporary hack for mask_immersed_field_xy! that hopefully we would fix in the future, but this PR goes the other direction and makes it permanent. It doesn't seem like we actually need a mask argument for mask_immersed_field!, the only point of adding it is for the sake of "uniformity"?

[glwagner]

Why don't we just start by identifying where we have to call mask_immersed_field! and mask_immersed_field_xy! with a different masking function. That basic motivation is missing for me. If we can't write this code any other way easily then we can merge and push coming up with a better design to the future.

[simone-silvestri]

The problem stems from the masking of the free surface.
The free surface is a field that lives at size(grid, 3) + 1 which, in theory, is a peripheral_node (since the z-direction is bounded). So the free surface needs to be masked with an inactive_node function instead of a peripheral_node. We could also just mask immersed_peripheral_nodes that would allow the distinction between actual boundaries and immersed boundaries and allow masking everything in the same way without touching domain boundaries

[simone-silvestri]

So let's maybe do this, we remove the masking function and use only the immersed_peripheral_node option for the moment, if we find cases where we need to provide another masking function we can think at another interface

[glwagner]

The problem stems from the masking of the free surface. The free surface is a field that lives at size(grid, 3) + 1 which, in theory, is a peripheral_node (since the z-direction is bounded). So the free surface needs to be masked with an inactive_node function instead of a peripheral_node. We could also just mask immersed_peripheral_nodes that would allow the distinction between actual boundaries and immersed boundaries and allow masking everything in the same way without touching domain boundaries

It's wrong that Nz+1 is a peripheral node when we have a free surface, right? The problem is not actually mask_immersed_field but the definition of peripheral_node. A peripheral_node is supposed to be a point that has a boundary condition for Face fields; for a model with a free surface this is not true at k = Nz + 1.

[simone-silvestri]

Right, it might be a bit difficult to enforce because the function needs to be model-dependent

[glwagner]

Right, it might be a bit difficult to enforce because the function needs to be model-dependent

I disagree. It's not that it's difficult to enforce. It's that the definition of "peripheral" is associated with the topology, and we may not have a topology to properly represent this concept.

It also seems wrong that "mask immersed field" masks all inactive nodes. It should be rather mask_inactive_nodes!. We can then have mask_peripheral_nodes! and mask_immersed_peripheral_nodes!. The code might be easier to understand and read then.

The confusion in this conversation is, I think, a symptom of the fact that the abstraction isn't quite right here. It's not holding up to logical analysis.

[glwagner]

A few more points:

• masking should only be needed for output?
• masking and boundary conditions are redundant for Face fields, I think. In other words we use "masking" to enforce boundary conditions.

[simone-silvestri]

exactly, we mask to have zero velocities on the boundary, otherwise velocities will have a value probably. So masking is required for velocities (unless we really make sure that tendencies are zero on immersed peripheral nodes but that is quite difficult to ensure for every new parameterization and sounds like it would be a bit brittle) but not for tracers. For tracers it's only required for outputs

[simone-silvestri]

It also seems wrong that "mask immersed field" masks all inactive nodes. It should be rather mask_inactive_nodes!. We can then have mask_peripheral_nodes! and mask_immersed_peripheral_nodes!. The code might be easier to understand and read then.

We could do that, it would be a little bit more code but I would be ok with it.

[glwagner]

exactly, we mask to have zero velocities on the boundary, otherwise velocities will have a value probably. So masking is required for velocities (unless we really make sure that tendencies are zero on immersed peripheral nodes but that is quite difficult to ensure for every new parameterization and sounds like it would be a bit brittle) but not for tracers. For tracers it's only required for outputs

Only the tendencies / vertical tridiagonal solve update the velocities I guess. But why do we only mask on immersed peripheral nodes? Is this again assuming that we redundantly satisfy boundary conditions? I think that algorithm is sort of overcomplicated. I can say that because I implemented as a simple thing to do on first try so long ago and we are still using it!

[simone-silvestri]

Ok, I think we can merge this to fix the inference issue and think about a better interface for masking immersed fields.