Topology-aware operators

[simone-silvestri]

these operators encode impenetrability boundary conditions for face fields and no-flux boundary conditions for tracers on Bounded grids and periodic boundary conditions on Periodic grids, without requiring a fill_halo procedure.

Previously only used for the split explicit free surface, we move them to the Operators module to allow more widespread use.

[glwagner]

How about

δxUᶜᵃᵃ

instead of

δxᶜᵃᵃ_U

so that the new operators keep the compact look of the current operators?

I would also support a new notation that involves just one letter for all of them, ie something like

δxTᶜᵃᵃ

where T is for "topology".

[glwagner]

@simone-silvestri what do you think about my notation suggestion? Personally I like keeping letters like "U" out of the operators so they are more physics agnostic

[simone-silvestri]

I was thinking about it. The problem is that they do in fact encode some physical statements, like, for example, impenetrability for U (XFaceFields) and V (YFaceFields) and no flux for center-fields (C).
Periodic topologies act according to periodic BC.
We probably need a way to make this clear. I am afraid a compact notation would not convey this entirely, maybe we can be a bit more verbose like

δyᵃᶜᵃ_nopenetration
δxᶠᵃᵃ_noflux

?

[glwagner]

I was thinking about it. The problem is that they do in fact encode some physical statements, like, for example, impenetrability for U (XFaceFields) and V (YFaceFields) and no flux for center-fields (C). Periodic topologies act according to periodic BC. We probably need a way to make this clear. I am afraid a compact notation would not convey this entirely, maybe we can be a bit more verbose like

δyᵃᶜᵃ_nopenetration
δxᶠᵃᵃ_noflux

?

Two thoughts:

• "penetration" and "flux" are physical concepts. More precisely the operator returns zero when differenced across a Bounded direction.
• The operator names are "mathematical" rather than English, ie we use symbols like δx, rather than x_difference. I think we should try to strive to keep names in one category or another rather than blending them

[glwagner]

For the compact notation, the symbol T merely denotes a difference with operators that do not use T. The meaning has to be looked up

[glwagner]

How is no penetration enforced? Do we actually need that one? Presumably, for faces we need interpolation operators (but not differences)

[simone-silvestri]

No penetration is enforced in this way:

Oceananigans.jl/src/Operators/topology_aware_operators.jl

Lines 47 to 48 in 5c7fafc

	@inline δxUᶜᵃᵃ(i, j, k, grid::AGXB, U★::Function, args...) = ifelse(i == grid.Nx, - U★(i, j, k, grid, args...),
	ifelse(i == 1, U★(2, j, k, grid, args...), δxᶜᵃᵃ(i, j, k, grid, U★, args...)))

Probably these operators should not typically live in the operator module as they encore "more" than the actual topology, they assume a boundary conditions. It would probably better to call them boundary_aware_operators.jl and include them in the Field module?

[glwagner]

is there another topology-aware operator of the form

δxᶜᵃᵃ

?

[simone-silvestri]

Thinking of revamping this, I think T is a good notation; we do not need to specify U and C since probably (in the future) we will need to refactor these a little to pass through boundary conditions.

[glwagner]

Thinking of revamping this, I think T is a good notation; we do not need to specify U and C since probably (in the future) we will need to refactor these a little to pass through boundary conditions.

Hopefully we don't have to pass boundary conditions 🥺

Not all complexity is justified by the performance gains...

[ali-ramadhan]

Is this PR looking for an adopter? Happy to try to complete it since it seems like the first step towards open BCs for hydrostatic models.

Thinking of revamping this, I think T is a good notation; we do not need to specify U and C since probably (in the future) we will need to refactor these a little to pass through boundary conditions.

Hopefully we don't have to pass boundary conditions 🥺

Not all complexity is justified by the performance gains...

Ah I'm actually not sure where the performance gains would come from, but I thought that the operators need to be aware of the boundary condition in order to correctly implement open BCs (at least based on @simone-silvestri's comments in #3628 (comment)).

[simone-silvestri]

If you want to adopt this PR, that would be great. I would leave the boundary conditions outside the operator for the moment, then we can look at the influence of those in performance in another PR. To test the performance we should push PR #3596 that implements split explicit with fill halos at every substep. Using that implementation for split-explicit will make boundary conditions for barotropic variables quite trivial to implement, but probably extremely slow. However, we can use it as a benchmark

[glwagner]

I don't support putting boundary conditions inside operators.

[simone-silvestri]

I am on the fence. In the cubed sphere aquaplanet (not even MPI but just on one GPU) the gain of performance is a factor 5 by using the efficient split explicit rather than filling the halos at each substep. We do not have to tackle this problem here or now, but we have keep in mind that fill halo is very inefficient and probably not the way to go. Maybe @ali-ramadhan finds a better way to have a gpu-compatible code that does not require adding boundary conditions in the operators.

[simone-silvestri]

I am on the fence. In the cubed sphere aquaplanet (not even MPI but just on one GPU) the gain of performance is a factor 5 by using the efficient split explicit rather than filling the halos at each substep. We do not have to tackle this problem here or now, but we have keep in mind that fill halo is very inefficient and probably not the way to go. Maybe @ali-ramadhan finds a better way to have a gpu-compatible code that does not require adding boundary conditions in the operators.

At least that is what @siddharthabishnu told me, I am not sure about the timings, can you confirm @siddharthabishnu?

[simone-silvestri]

I am on the fence. In the cubed sphere aquaplanet (not even MPI but just on one GPU) the gain of performance is a factor 5 by using the efficient split explicit rather than filling the halos at each substep. We do not have to tackle this problem here or now, but we have keep in mind that fill halo is very inefficient and probably not the way to go. Maybe @ali-ramadhan finds a better way to have a gpu-compatible code that does not require adding boundary conditions in the operators.

For example, if the boundary is fixed and not dependent on the interior variables, it might be possible just to fix the BC once before the sub stepping and just iterate without needing to use any update

[glwagner]

the gain of performance is a factor 5 by using the efficient split explicit rather than filling the halos at each substep

What does this have to do with putting boundary conditions in operators?

[ali-ramadhan]

To get this PR merged, it sounds like we just want to agree on a naming convention and move the topologically-aware operators to the Operators module?

I like the T for topology suggestion, e.g. δxTᶜᵃᵃ. These operators are not exported so they don't need to have user-friendly names, only developer-friendly names. I'm happy to move forward with this suggestion and try to get this PR review-ready.

I also don't think operators should depend on boundary conditions. But I can write down some thoughts in a more appropriate issue.

[ali-ramadhan]

I notice that here we're dispatching on ::Function whereas elsewhere for operators we've dispatched on ::F where F<:Function. Is there a difference?

@inline δxᶜᵃᵃ(i, j, k, grid, f::F, args...) where F<:Function = f(i+1, j, k, grid, args...) - f(i,   j, k, grid, args...)

vs.

@inline δxCᶠᵃᵃ(i, j, k, grid, c★::Function, args...) = δxᶠᵃᵃ(i, j, k, grid, c★, args...)

[glwagner]

Long ago it was argued that F<:Function encourages the compiler to specialize these operators for the specific type of F which is needed for GPU compilation

[ali-ramadhan]

Did some refactoring and CPU and GPU hydrostatic regression tests pass locally so I think this PR is ready for review!

@simone-silvestri Requesting a review from you on your own PR again lol but I won't approve unless you do.

[ali-ramadhan]

Just bumping my review request 🥺

[simone-silvestri]

I can't approve it, being the creator of the PR, but I approve here via text. I guess this PR does not introduce all the possible operators (for example, there are no δzTᵃᵃᶠ) but only the ones we need. I think it is okay at the moment. Then, we can introduce the other operators when/if we need them

[ali-ramadhan]

Thanks @simone-silvestri, approved! And yeah hopefully it'll be easy/straightforward to add any other operators we'll need in the future.