`Field.eval` silently returns wrong values for out-of-domain queries

[fluidnumericsJoe]

Parcels version

sha: bc20889

Description

Summary

When Field.eval is called directly (without a ParticleSet), the search-error sentinels returned by the index-search routines are passed straight through to the interpolator's corner-lookup, where they are used as raw indices. numpy negative-indexing then wraps them to the opposite edge of the grid, and the interpolator silently returns a weighted blend of corners pulled from unrelated parts of the field — no exception, no NaN, no warning.

This was discovered while investigating #2521. The PIC bug there has been fixed; this sentinel-leak path is an independent issue.

Affected sentinels

All three search-error sentinels share the same defect:

Sentinel	Source	Triggered by
LEFT_OUT_OF_BOUNDS = -2	_search_1d_array (rectilinear)	query < arr[0]
RIGHT_OUT_OF_BOUNDS = -1	_search_1d_array (rectilinear)	query > arr[-1]
GRID_SEARCH_ERROR = -3	_search_indices_curvilinear_2d (curvilinear)	spatial hash + PIC find no containing cell

Suggested acceptance criteria for a fix

• Field.eval(...) on an out-of-domain query without a ParticleSet raises an explicit error or returns NaN deterministically (not via index wrap).
• Behavior is uniform across all three sentinels (LEFT_OUT_OF_BOUNDS, RIGHT_OUT_OF_BOUNDS, GRID_SEARCH_ERROR) and across rectilinear and curvilinear paths.
• A regression test covers the direct-call path on a query that lies outside the mesh.

Code sample

import warnings

import numpy as np

from parcels import FieldSet
from parcels._datasets.structured.generated import simple_UV_dataset

warnings.filterwarnings("ignore")

# Tiny flat-mesh dataset so the corner indices are easy to read.
# Defaults: dims=(time=360, depth=2, ydim=30, xdim=4); we'll keep ydim, xdim small
# but informative.
ydim, xdim = 6, 5
ds = simple_UV_dataset(dims=(2, 1, ydim, xdim), mesh="flat")

# Functional U so that U[t, z, j, i] = j + 0.001 * i.
# The integer part identifies the j index, the fractional part the i index,
# so the value returned by Field.eval reveals which cells were sampled.
U = np.zeros(ds["U"].shape, dtype=float)
for j in range(ydim):
    for i in range(xdim):
        U[:, :, j, i] = j + 0.001 * i
ds["U"].data = U

print(f"Domain: lon in [{ds['lon'].values.min():.2e}, {ds['lon'].values.max():.2e}]"
      f"  lat in [{ds['lat'].values.min():.2e}, {ds['lat'].values.max():.2e}]")
print(f"Grid:   ydim={ydim}, xdim={xdim}, U[t,z,j,i] = j + 0.001*i\n")

fieldset = FieldSet.from_sgrid_conventions(ds, mesh="flat")

# --- Case 1: query east of the domain (xi sentinel = -1) ---
t = np.array([0.0])
z = np.array([0.0])
y = np.array([0.0])  # mid-domain lat
x = np.array([5e6])  # 5x past the east edge
val = fieldset.U.eval(t, z, y, x)
expected_j = ydim // 2  # mid-y bilinear, ~2.5 if exactly between rows
print(f"Case 1: eval at lon=5e6 (5x past east edge), lat=0  -> U = {val[0]:.6f}")

# --- Case 2: query west of the domain (xi sentinel = -2) ---
x = np.array([-5e6])
val = fieldset.U.eval(t, z, y, x)
print(f"Case 2: eval at lon=-5e6 (5x past west edge), lat=0 -> U = {val[0]:.6f}")

# --- Case 3: query simultaneously out of bounds in both lat and lon ---
y = np.array([5e6])
x = np.array([5e6])
val = fieldset.U.eval(t, z, y, x)
print(f"Case 3: eval at lon=5e6, lat=5e6 (both past edges)  -> U = {val[0]:.6f}")

Output shown below :

Domain: lon in [-1.00e+06, 1.00e+06]  lat in [-1.00e+06, 1.00e+06]
Grid:   ydim=6, xdim=5, U[t,z,j,i] = j + 0.001*i

Case 1: eval at lon=5e6 (5x past east edge), lat=0  -> U = 2.468000
Case 2: eval at lon=-5e6 (5x past west edge), lat=0 -> U = 2.527000
Case 3: eval at lon=5e6, lat=5e6 (both past edges)  -> U = -50.032000