Change default advection scheme and halo size for grids, and add utilities for inferring needed halo sizes?

[glwagner]

There are a few applications that require building fields before IncompressibleModel, in order to form AbstractOperations and ComputedFields that need be computed during a model time step (eg PR #1091). I think this use case will only become more and more important in the future.

Currently this functionality is possible but plagued by a huge useability issue: the default grid has a halo size of 1, while most applications benefit from higher-order advection schemes. To hide the need to specify halo region sizes from users, we currently "inflate" halos inside the constructor for IncompressibleModel:

Oceananigans.jl/src/Models/IncompressibleModels/incompressible_model.jl

Lines 111 to 112 in 5ddace1

	Hx, Hy, Hz = inflate_halo_size(grid.Hx, grid.Hy, grid.Hz, advection, closure)
	grid = with_halo((Hx, Hy, Hz), grid)

This means that users who want to build fields before IncompressibleModel do, in fact, have to know the halo size they need to specify for their chosen advection scheme. This isn't well-documented right now...

Somehow, we have to figure out how to smooth this whole process out. One huge help will simply be to choose a default advection scheme that is useful for science: either UpwindBiasedFifthOrder or WENO5, and to set the default halo size for the grid to 5. Having these default will mitigate the problem greatly I think.

But we also probably need utilities (or documentation at the very least) that explains this issue and how to choose the halo size if one needed to build the grid outside the model.

Or, perhaps there are even better solutions to this issue. Basically the point is that the grid depends "circularly" on aspects of IncompressibleModel, which becomes an issue when things like VelocityFields (which also depend on the grid) need to be constructed prior to IncompressibleModel. We don't want users to have to replicate the IncompressibleModel constructor in their scripts...

[glwagner]

After discussing with @ali-ramadhan, it seems more clear that a simple solution is just to use a default halo size of 3. For most models that use the highest order advection scheme we offer this has no effect. For models that use a lower-order advection scheme but don't change the halo size, the memory foot print of the model is ever-so-slightly larger than it needs to be. But this slightly-larger footprint probably isn't noticeable for most problems.

So in summary, minimal halo sizes are a minor optimization that has little effect on most problems. Auto-optimizing the halo size has major downsides for usability, so I think the trade-off leans towards big default halos.

[navidcy]

Sure. But let's document this somewhere in the docs or in the appropriate docstrings so that users who want the extra halo size (which I guess is required for achieving the scheme's convergence accuracy) know how to do it themselves.

[glwagner]

@navidcy if halos are too small then the code often NaNs or seg faults. The issue here is that we auto adjust halos in the model constructor:

Oceananigans.jl/src/Models/IncompressibleModels/incompressible_model.jl

Lines 111 to 112 in ebd7858

	Hx, Hy, Hz = inflate_halo_size(grid.Hx, grid.Hy, grid.Hz, advection, closure)
	grid = with_halo((Hx, Hy, Hz), grid)

this means that with a higher-order advection scheme, model.grid ends up being different from what the user passes into the constructor.

This is a usability issue, because it means that VelocityFields(arch, grid) built before the model is wrong for the default grid in common scenarios with high-order advection. Thus if users want to build VelocityFields(arch, grid) before building a model, they need to know about halos. We think this is undesirable.

Neatly summed the issue is mainly that we probably can't rely exclusively on the model constructor to infer halo sizes, because the grid is a crucial object that often needs to have correct halo sizes before the model is constructed. Thus we need another solution to this usability issue. I'm proposing that we make the default halos large enough to accommodate almost all use cases to solve this problem here.

[glwagner]

We've done a lot of this.