Proposal for simplifying the Field Interpolator API

[erikvansebille]

(This Issue continues on #2059, but is a bit more specific)

The current v4-dev signature for Interpolators is

def Interpolator(
    field: Field,
    ti: int,
    position: dict[_XGRID_AXES, tuple[int, float | np.ndarray]],
    tau: np.float32 | np.float64,
    t: np.float32 | np.float64,
    z: np.float32 | np.float64,
    y: np.float32 | np.float64,
    x: np.float32 | np.float64,
) -> np.float32 | np.float64:

Now, I would propose that we simplify this to:

def Interpolator(
    positions: dict[_XGRID_AXES, tuple[int, float | np.ndarray]],
    field: Field,
) -> np.float32 | np.float64:

where the (ti, tau) tuple also becomes part of positions plural, as it is vectorized (thinking spacetime-dimensions). We also don't need the physical coordinates (t, z, y, x) as input because we can always recompute these from the computational coordinates if needed.

I also swapped the order of the two remaining arguments to mirror the Kernel signature (which is Kernel(particles, fieldset)).

Any objections to this proposal?

[fluidnumerics-joe]

I like the idea of encapsulating time indices and linear interpolation weights into the positions object. It changes things a little bit on our (the developer's) side in that the xgrid.search / uxgrid.search methods fill in the spatial component and time index searching is currently done withing field.eval; I suppose we could stitch the ti,tau information into the positions object within field.eval prior to calling interpolate. Not a big change for tidying up the presentation of the interpolation API to users.

We also don't need the physical coordinates (t, z, y, x) as input because we can always recompute these from the computational coordinates if needed.

While this is indeed true, and I wouldn't be opposed to computing these on the fly, I could imagine some users getting tired of copying (or having to remember) boiler plate like the following to compute positions

fi = position["FACE"][0]
bcoords = positions["FACE"][1]
node_ids = field.grid.uxgrid.face_node_connectivity[fi, :].values
x = np.sum(field.grid.uxgrid.node_lon.values[node_ids] * bcoords, axis=-1)
y = np.sum(field.grid.uxgrid.node_lat.values[node_ids] * bcoords, axis=-1)

Having the partiticle positions as an argument is quite a nice convenience

[erikvansebille]

It changes things a little bit on our (the developer's) side in that the xgrid.search / uxgrid.search methods fill in the spatial component and time index searching is currently done withing field.eval; I suppose we could stitch the ti,tau information into the positions object within field.eval prior to calling interpolate.

Indeed, that would be the way to go. Not difficult, because the two function are called back-to-back in

Parcels/src/parcels/_core/field.py

Lines 219 to 220 in 4330de3

	tau, ti = _search_time_index(self, time)
	position = self.grid.search(z, y, x, ei=_ei)

Having the partiticle positions as an argument is quite a nice convenience

OK, I see your point about having access to the particle positions too, @fluidnumerics-joe. But in that case I propose we add these to the positions dictionary, as z, lat and lon (following the naming of the particle variables)

So positions is then a dictionary with 12 keys: time, ti, tau, z, zi, zeta, lat, yi, eta, lon, xi, xsi

Agree?

[fluidnumerics-joe]

So positions is then a dictionary with 12 keys: time, ti, tau, z, zi, zeta, lat, yi, eta, lon, xi, xsi
Agree?

Yes! I like the consolidation

[VeckoTheGecko]

Just catching up on all of this now - sorry for the delay.

I like the idea of encapsulating time indices and linear interpolation weights into the positions object.

+1

While this is indeed true, and I wouldn't be opposed to computing these on the fly, I could imagine some users getting tired of copying (or having to remember) boiler plate like the following to compute positions

One option would be for us to provide utilities that are able to compute physical positions from barycentric coordinates (at least hiding this boilerplate from the interpolator writers).

So positions is then a dictionary with 12 keys: time, ti, tau, z, zi, zeta, lat, yi, eta, lon, xi, xsi

Agree?

Should we make the distinction here between barycentric and physical positions by splitting these into two dictionaries? (otherwise its a mix of the two based on which keys are being accessed, which I think introduces some more complexity)

---

I think computing the physical coordinates in the interpolator (via boilerplace or helper functions) or pre-computing them are both viable options which we'll only know the impact of through benchmarking

[erikvansebille]

I just implemented the 'flat' dictionary approach in fe9844d. But I can see that this may be difficult (for example I can't easily get grid.localize() to work), so an alternative could be a nested dictionary

postions = {
    "T": {"time": time, "ti": ti, "tau": tau}, 
    "Z": {"z": z, "zi": zi, "zeta": zeta}, 
    "Y": {"lat": lat, "yi": yi, "eta": eta}, 
    "X": {"lon": lon, "xi": xi, "xsi", xsi"
}

Then for unstructured grids, the X and Y entries could be replaced by FACE

postions = {
    "T": {"time": time, "ti": ti, "tau": tau}, 
    "Z": {"z": z, "zi": zi, "zeta": zeta}, 
    "FACE": {"lat": lat, "lon": lon, "fi": fi, "bcoord": bcoord}
}

Would that make more sense (also for uxgrids, @fluidnumerics-joe)?

[erikvansebille]

As we discussed at the meeting, an alternative could be

positions = {
    "T": {"coord": time, "index": ti, "bcoord": tau}, 
    "Z": {"coord": z, "index": zi, "bcoord": zeta}, 
    "Y": {"coord": lat, "index": yi, "bcoord": eta}, 
    "X": {"coord": lon, "index": xi, "bcoord", xsi"},
}

[erikvansebille]

However, @VeckoTheGecko proposed to separate the particle position from the grid positions.

In my interpretation of that idea, this could then become

def Interpolator(
    particle_positions: dict,
    grid_positions: dict,
    field: Field,
) -> np.float32 | np.float64:

where then

particle_positions = {"time", time, "z", z, "lat", lat, "lon", lon}

and for XGrids

grid_positions = {
    "T": {"index": ti, "bcoord": tau}, 
    "Z": {"index": zi, "bcoord": zeta}, 
    "Y": {"index": yi, "bcoord": eta}, 
    "X": {"index": xi, "bcoord", xsi},
}

and for UXGrids

grid_positions = {
    "T": {"index": ti, "bcoord": tau}, 
    "Z": {"index": zi, "bcoord": zeta}, 
    "FACE": {"index": fi, "bcoord": bcoord}
}