Advection schemes should depend on the grid for spacings and float type inference

[ali-ramadhan]

I'm trying to understand how advection schemes are constructed, and I see that you can pass a float type FT and a grid to the constructors.

Makes perfect sense that a scheme could depend on FT to store any coefficients using the correct data type. It also makes perfect sense for a scheme to depend on grid if the coefficients vary spatially.

But some constructors take both. For example, Centered takes both but doesn't really use the grid. It only uses the grid , if passed, to overwrite FT. Well, looks like compute_reconstruction_coefficients is commented out.

Oceananigans.jl/src/Advection/centered_reconstruction.jl

Lines 39 to 59 in 82099ac

	function Centered(FT::DataType = Float64; grid = nothing, order = 2)
	if !(grid isa Nothing)
	FT = eltype(grid)
	end
	mod(order, 2) != 0 && throw(ArgumentError("Centered reconstruction scheme is defined only for even orders"))
	N = Int(order ÷ 2)
	if N > 1
	coefficients = Tuple(nothing for i in 1:6)
	# Stretched coefficient seem to be more unstable that constant spacing ones for centered reconstruction.
	# Some tests are needed to verify why this is the case (and if it is expected). We keep constant coefficients for the moment
	# coefficients = compute_reconstruction_coefficients(grid, FT, :Centered; order)
	buffer_scheme = Centered(FT; grid, order = order - 2)
	else
	coefficients = Tuple(nothing for i in 1:6)
	buffer_scheme = nothing
	end
	return Centered{N, FT}(coefficients..., buffer_scheme)
	end

Whereas WENO uses the grid to compute reconstruction coefficients:

Oceananigans.jl/src/Advection/weno_reconstruction.jl

Lines 96 to 119 in 82099ac

	function WENO(FT::DataType=Float64;
	order = 5,
	grid = nothing,
	bounds = nothing)
	if !(grid isa Nothing)
	FT = eltype(grid)
	end
	mod(order, 2) == 0 && throw(ArgumentError("WENO reconstruction scheme is defined only for odd orders"))
	if order < 3
	# WENO(order = 1) is equivalent to UpwindBiased(order = 1)
	return UpwindBiased(FT; order = 1)
	else
	N = Int((order + 1) ÷ 2)
	weno_coefficients = compute_reconstruction_coefficients(grid, FT, :WENO; order = N)
	advecting_velocity_scheme = Centered(FT; grid, order = order - 1)
	buffer_scheme = WENO(FT; grid, order = order - 2, bounds)
	end
	return WENO{N, FT}(weno_coefficients..., bounds, buffer_scheme, advecting_velocity_scheme)
	end

---

In both cases grid is an optional argument. But won't the advection scheme be wrong if grid is not passed then? If so, then shouldn't advection schemes just depend on the grid (which can be used to infer FT)?

I know some interfaces take FT while others take grid (X-Ref: #3800) so I'm not sure of the right approach here, but it seems like some advection schemes actually depend on the grid spacings.

X-Ref: #3703

[simone-silvestri]

I've been wanting to clean up this interface a bit but haven't gotten to do it.
The grid argument is used to compute coefficients for stretched advection; if it is not passed, we use uniform spacing coefficients.
Some time ago, we disabled using stretched coefficients for Centered and Upwind because it seemed they would make the simulation less stable, but there was no rigorous study to check this.

The grid is a kwarg to allow doing stuff like

WENO()

I think the idea is to move to an interface where every element requires the grid as a positional argument so if this is the idea, we could introduce a keyword argument like stretched_coefficients which would make the interface clearer

[glwagner]

grid should be a positional argument (note that this is no barrier to the syntax WENO(), we can have optional positional arguments, so I'm not sure I understand @simone-silvestri's point).

I'd favor devoting this issue to solving this problem.

PS for pure questions, I think it's best to use discussions, and reserve "issues" for things that require a PR to fix.

[ali-ramadhan]

Thanks for the explanation @simone-silvestri!

@glwagner Yeah I think this issue started off as a question (thus the title) then turned into more of a UI suggestion. I've updated the title to reflect this.

Totally support making grid a positional argument. I wonder if making it a required argument would that help with #3800 by making sure users won't forget to pass the grid to set the float type.

[simone-silvestri]

The grid is also a positional argument for WENO at the moment

Oceananigans.jl/src/Advection/weno_reconstruction.jl

Line 121 in 03b8acf

WENO(grid, FT::DataType=Float64; kwargs...) = WENO(FT; grid, kwargs...)

The question is whether we want a better keyword argument to specify if we use stretched or uniform coefficients. For that, the grid would probably be a required keyword argument and not an optional one.

[glwagner]

What's the point of having grid as both a positional and keyword argument?

[simone-silvestri]

I think there is not much thought behind it. I support removing the grid from the keyword arguments.
If we decide to make the grid a required type, we can also add a stretched_coefficients = true/false keyword argument.

[glwagner]

I thought we were going to discontinue support for stretched coefficients?