Halo inflation + restriction on halo size makes it impossible to run some problems

[glwagner]

For example,

using Oceananigans

grid = RectilinearGrid(size = (16, 1, 16),
                       x = (0, 1),
                       y = (0, 1),
                       z = (0, 1),
                       topology = (Periodic, Periodic, Bounded))

model = NonhydrostaticModel(; grid, advection = WENO(order=5))

errors with

julia> include("problem.jl")
┌ Warning: Inflating model grid halo size to (3, 3, 3) and recreating grid. Note that an ImmersedBoundaryGrid requires an extra halo point in all non-flat directions compared to a non-immersed boundary grid.
└ @ Oceananigans.Models.NonhydrostaticModels ~/.julia/packages/Oceananigans/OHYQj/src/Models/NonhydrostaticModels/nonhydrostatic_model.jl:223
ERROR: LoadError: ArgumentError: halo must be ≤ size for coordinate y

Sad, because I didn't actually ask for a halo (3, 3, 3)! Manually setting things doesn't help:

using Oceananigans

grid = RectilinearGrid(size = (16, 1, 16),
                       halo = (3, 1, 3),
                       x = (0, 1),
                       y = (0, 1),
                       z = (0, 1),
                       topology = (Periodic, Periodic, Bounded))

model = NonhydrostaticModel(; grid, advection = WENO(order=5))

I guess the restriction on the halo size is supposed to be an optimization, but its not working in this case because I can't even set up my model. @navidcy

[navidcy]

okie -- there is definitely something wrong there

[glwagner]

bump ran into this again today

[simone-silvestri]

It is not an optimization; rather, to fill in periodic halos, a field must have a grid size larger than the halos (for other boundary conditions, only one is required).
@navidcy, you could probably look at how the periodic halos are filled and make sure that they work also if the dimension is smaller than the number of halos.
It might entail writing some ifelse in the fill periodic halo kernels.

For example:

@kernel function fill_periodic_west_and_east_halo!(c, ::Val{H}, N) where H
    j, k = @index(Global, NTuple)
    @unroll for i = 1:H
        @inbounds begin
            # Put some ifelse condition here to make sure that halos pick from the interior
            # and not other halos
            c[i, j, k]     = c[N+i, j, k] # west
            c[N+H+i, j, k] = c[H+i, j, k] # east
        end
    end
end

[glwagner]

@simone-silvestri we just need to be able to write

using Oceananigans

grid = RectilinearGrid(size = (16, 1, 16),
                       halo = (3, 1, 3),
                       x = (0, 1),
                       y = (0, 1),
                       z = (0, 1),
                       topology = (Periodic, Periodic, Bounded))

[simone-silvestri]

One option is the one above, other options that could achieve that results are:
2) Limit the advection scheme order based on the grid size
3) we write something similar to TracerAdvection for momentum to allow different advection schemes in different directions and use nothing in the periodic direction with one grid point

[simone-silvestri]

I like option 3. I could open a PR to tackle this issue if there is consensus in the solution

[glwagner]

But then what would you do for 2 grid points?

[glwagner]

Why not implement both 3 and 2?

[simone-silvestri]

Thinking about it again, I would prefer to stay away from 2 because it results in very implicit behavior. It is nice to be explicit about the scheme used.

[glwagner]

We already do 2 near boundaries

[glwagner]

3 alone doesn't solve the problem, because we still don't know what to do with 2 grid points.

[glwagner]

The question is just what you want to do with things like

using Oceananigans

grid = RectilinearGrid(size = (16, 1, 16),
                       halo = (3, 1, 3),
                       x = (0, 1),
                       y = (0, 1),
                       z = (0, 1),
                       topology = (Periodic, Periodic, Bounded))

model = NonhydrostaticModel(; grid, advection = WENO(order=5))

and

using Oceananigans

grid = RectilinearGrid(size = (16, 2, 16),
                       halo = (3, 2, 3),
                       x = (0, 1),
                       y = (0, 1),
                       z = (0, 1),
                       topology = (Periodic, Periodic, Bounded))

model = NonhydrostaticModel(; grid, advection = WENO(order=5))

etc

[glwagner]

I think we want to handle this automatically because for example WENO() has a default order. I feel this falls within our philosophy of "friendly" (not sure if @simone-silvestri agrees with this philosophy though 😆 )

[glwagner]

The risk of doing something automatic is that in theory someone could think "oh, I only have 2 grid poits in the y direction, but I still expect to observe WENO advection in that direction". Will someone think that? It doesn't really make sense, with 2 grid points in a direction there won't be advection at all anyways. I don't think this is a big risk. Big risks are like, simulations are dead wrong when they present an illusion that makes them seem ok.