Fix `tutorial_nested_grids.ipynb` - grid searching

[VeckoTheGecko]

OK - running into a separate (unrelated) problem surfaced by this change which I think runs deeper...

docs/user_guide/examples/tutorial_nestedgrids.ipynb combines different datasets together - but it was specifying "flat" for the gridIDs and "spherical" for the data itself (which although somehow worked - doesn't make sense to me, and doesn't align with the work we've been doing).

This is in lon/lat space, so I updated to use "spherical" and got grid searching errors.

---------------------------------------------------------------------------
GridSearchingError                        Traceback (most recent call last)
Cell In[10], [line 18](vscode-notebook-cell:?execution_count=10&line=18)
     14 def SampleGridID(particles, fieldset):
     15     particles.gridID = fieldset.GridID[particles]
     16 
     17 
---> [18](vscode-notebook-cell:?execution_count=10&line=18) pset.execute(
     19     SampleGridID,
     20     runtime=np.timedelta64(1, "s"),
     21     dt=np.timedelta64(1, "s"),

File ~/coding/repos/parcels/src/parcels/_core/particleset.py:453, in ParticleSet.execute(self, kernels, dt, endtime, runtime, output_file, verbose_progress)
    450 else:
    451     next_time = end_time
--> [453](https://file+.vscode-resource.vscode-cdn.net/Users/Hodgs004/coding/repos/parcels/docs/user_guide/examples/~/coding/repos/parcels/src/parcels/_core/particleset.py:453) self._kernel.execute(self, endtime=next_time, dt=dt)
    455 if next_output is not None:
    456     if np.abs(next_time - next_output) < 0.001:

File ~/coding/repos/parcels/src/parcels/_core/kernel.py:244, in Kernel.execute(self, pset, endtime, dt)
    242                 error_func(pset[inds].t)
    243             else:
--> [244](https://file+.vscode-resource.vscode-cdn.net/Users/Hodgs004/coding/repos/parcels/docs/user_guide/examples/~/coding/repos/parcels/src/parcels/_core/kernel.py:244)                 error_func(pset[inds].z, pset[inds].y, pset[inds].x)
    246 return pset

File ~/coding/repos/parcels/src/parcels/_core/statuscodes.py:85, in _raise_grid_searching_error(z, y, x)
     84 def _raise_grid_searching_error(z, y, x):
---> [85](https://file+.vscode-resource.vscode-cdn.net/Users/Hodgs004/coding/repos/parcels/docs/user_guide/examples/~/coding/repos/parcels/src/parcels/_core/statuscodes.py:85)     raise GridSearchingError(f"Grid searching failed at (z={z}, y={y}, x={x})")

GridSearchingError: Grid searching failed at (z=array([0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
       0.], dtype=float32), y=array([-18.        , -18.        , -18.        , -18.        ,
       -18.        , -15.052631  , -15.052631  , -15.052631  ,
       -15.052631  , -12.105263  , -12.105263  ,  -9.157895  ,
        -9.157895  ,  -9.157895  ,  -6.2105265 ,  -3.2631578 ,
        -3.2631578 ,  -3.2631578 ,  -3.2631578 ,  -3.2631578 ,
        -3.2631578 ,  -0.31578946,  -0.31578946,   2.631579  ,
         5.5789475 ,   8.526316  ,   8.526316  ,  11.473684  ,
        11.473684  ,  11.473684  ,  14.421053  ,  14.421053  ,
        14.421053  ,  17.368422  ,  17.368422  ,  17.368422  ,
        17.368422  ,  17.368422  ,  20.31579   ,  20.31579   ,
        23.263159  ,  26.210526  ,  29.157894  ,  29.157894  ,
        29.157894  ,  29.157894  ,  32.105263  ,  35.05263   ,
        35.05263   ,  38.        ,  38.        ,  38.        ],
      dtype=float32), x=array([ 0.25,  3.  ,  5.75,  8.5 , 11.25,  8.5 , 11.25, 14.  , 52.5 ,
       16.75, 49.75, -2.5 , 44.25, 47.  , 41.5 , -8.  ,  3.  , 16.75,
       19.5 , 30.5 , 38.75, -8.  , 33.25, -8.  , -8.  , -8.  , -5.25,
       -8.  , -5.25, 25.  , -8.  , -5.25, 25.  , -8.  , -5.25, -2.5 ,
       25.  , 44.25, 25.  , 47.  , 47.  , 49.75, 16.75, 22.25, 49.75,
       52.5 , 52.5 , 49.75, 55.25, 22.25, 55.25, 58.  ], dtype=float32))

Adapting the plot I was able to find the grid search failure locations

Using an LLM it was talking about there being a bug in the point-in-cell check, but to be honest I don't know how accurate that is.

LLM output

Spherical UxGrid search fails for large triangles

When using mesh="spherical" with UxGrid, grid searching fails for points inside large triangular faces. Two independent bugs:

Bug 1: SpatialHash Cartesian bounding box too tight

SpatialHash.__init__ (spatialhash.py:182-188) computes each face's bounding box from the Cartesian (x,y,z) coordinates of its vertices only. On a sphere, great-circle edges curve outward relative to the chord, so for large faces a query point on the sphere surface can lie inside the face but outside the vertex-only Cartesian bounding box. The hash never considers the face as a candidate.

Suggested fix: Include great-circle edge midpoints (average vertex Cartesian coords, renormalize to unit sphere) in the bounding box computation.

Bug 2: uxgrid_point_in_cell barycentric check too strict

uxgrid_point_in_cell (index_search.py:298-372) projects the query onto the face plane and computes barycentric coordinates using unsigned 3D triangle areas. For large spherical triangles, projection distortion causes sum(bcoords) to exceed 1 (e.g. 1.023 for a ~70°-wide face), failing the np.isclose(..., 1.0, rtol=1e-3) check at line 370.

Minimal reproducer

import numpy as np
import uxarray as ux
import xarray as xr
from shapely.geometry import Point, Polygon
from triangle import triangulate

from parcels._core.index_search import uxgrid_point_in_cell
from parcels._core.uxgrid import UxGrid

# Build a multi-face constrained Delaunay triangulation from three nested polygons
polygons = [
    np.array([(10, 15), (25, 10), (25, 25), (17, 36), (10, 32)]),
    np.array([(0, -5), (35, 0), (35, 25), (0, 20)]),
    np.array([(-10, -20), (60, -20), (60, 40), (-10, 40)]),
]
verts, segments = [], []
offset = 0
for poly in polygons:
    Ni = len(poly)
    verts.extend(poly.tolist())
    segments.extend([[offset + j, offset + ((j + 1) % Ni)] for j in range(Ni)])
    offset += Ni
B = triangulate(
    {"vertices": np.asarray(verts, dtype=float), "segments": np.asarray(segments, dtype=int)},
    "p",
)
points, face_tris = B["vertices"], B["triangles"].astype(int)

ds = xr.Dataset(
    {
        "node_lon": ("n_node", points[:, 0]),
        "node_lat": ("n_node", points[:, 1]),
        "face_node_connectivity": (("n_face", "n_max_face_nodes"), face_tris),
    },
    coords={
        "n_node": np.arange(points.shape[0]),
        "n_face": np.arange(face_tris.shape[0]),
    },
    attrs={"Conventions": "UGRID-1.0"},
)
uxgrid = ux.Grid(ds)
z = ux.UxDataArray(np.array([0.0, 1.0]), dims="z", uxgrid=uxgrid)
grid = UxGrid(uxgrid, z, mesh="spherical")

# --- Bug 1: SpatialHash misses points inside a large face ---
sh = grid.get_spatial_hash()
for lon, lat in [(0.25, -18.0), (3.0, -18.0), (5.75, -18.0)]:
    _, face_ids, _ = sh.query(
        np.array([lat], dtype=np.float32),
        np.array([lon], dtype=np.float32),
    )
    print(f"Bug 1 — ({lon}, {lat}): face_id={face_ids[0]} (expected >=0, got GRID_SEARCH_ERROR)")

# --- Bug 2: point-in-cell rejects a point near the edge of a large face ---
# Find which face contains (52.5, -15) by brute force — none pass the check
for fi in range(face_tris.shape[0]):
    is_in, coords = uxgrid_point_in_cell(
        grid,
        np.array([-15.0], dtype=np.float32),
        np.array([52.5], dtype=np.float32),
        np.array([fi]),
        np.array([fi]),
    )
    if abs(coords[0].sum() - 1.0) < 0.05:  # close but fails the rtol=1e-3 check
        print(
            f"Bug 2 — (52.5, -15) face {fi}: is_in={is_in[0]}, "
            f"bcoords sum={coords[0].sum():.4f} (fails isclose rtol=1e-3)"
        )

Output:

Bug 1 — (0.25, -18.0): face_id=-3 (expected >=0, got GRID_SEARCH_ERROR)
Bug 1 — (3.0, -18.0): face_id=-3 (expected >=0, got GRID_SEARCH_ERROR)
Bug 1 — (5.75, -18.0): face_id=-3 (expected >=0, got GRID_SEARCH_ERROR)
Bug 2 — (52.5, -15) face 12: is_in=0, bcoords sum=1.0234 (fails isclose rtol=1e-3)
Bug 2 — (52.5, -15) face 13: is_in=0, bcoords sum=1.0492 (fails isclose rtol=1e-3)

Bug 1 is the primary blocker — it prevents tutorial_nestedgrids.ipynb from running with mesh="spherical" on the triangulation FieldSet.

@erikvansebille did you run into something similar when authoring the notebook? cc @fluidnumericsJoe

Originally posted by @VeckoTheGecko in #2848 (comment)

---

cc @erikvansebille, @wyatt-fluidnumerics, @fluidnumericsJoe - nicer to continue here than in the original PR :)

[erikvansebille]

Thanks for opening this issue, @VeckoTheGecko. Did you make any progress yesterday, @fluidnumericsJoe and @wyatt-fluidnumerics? I think it's important to get this fixed before release, as nested grids is one of the key improvements of v4!

[fluidnumericsJoe]

@wyatt-fluidnumerics and I are fairly confident we've narrowed in on the issue here.. There's two bugs

First comes from the particle-in-cell search performing an orthogonal projection of the particle position onto the planar face of a triangular cell. When we have a spherical mesh type, we use the corner nodes of a cell to define a planar face to approximate the surface of the sphere. A particle's lat,lon position is converted to x,y,z using the spherical coordinate transformation. Then, we remove the component of the position that is in the face normal direction to "drag" the particle down onto the face. For particles near the edge of a face, this is actually quite catastrophic and can actually project the x,y,z position into a neighboring element erroneously. In the PIC search, this results in the sum of the barycentric coordinates being much greater than 1 (a condition we use to mark particles not in a cell).

What we want instead is a gnomonic projection, which drags the particle down into the planar face along a radial line that originates from the center of the unit sphere.

We had claude create a rendering that was slightly nicer than our chicken-scratch on a whiteboard :

For your mesh, with very large elements, this error is very noticeable; it was there for all of our test cases but apparently for the resolutions we were testing at, the error in the barycentric coordinate sum was less than our tolerance (1e-3). Using a gnomonic projection, we can drop that tolerance to machine precision now and shouldn't have to worry about this failure anymore.

The other failure mode comes from the hash table construction. When the mesh type is spherical, we construct the bounds of the hash grid from the bounding box formed by the corner nodes only. For regional domains with very large elements (like this case), this is problematic since, for example, converting a face center lat,lon to x,y,z could actually produce a position that is outside of the hash grid. If a particle is happens to be at these points that happen to lie just outside the hash grid, the morton code they end up associated with is not even a key in the hash table. To fix this, we're rethinking how we calculate the bounding box for the hash table.

We've created a few tests that confirm these two independent failure modes and have a fix for the first worked out. We're still sorting out the second one.

[erikvansebille]

Thanks for this update, @fluidnumericsJoe! Very interesting - and good then that we uncovered this because the nested grid tutorial uses such large cells. It will only make the code more accurate - and hopefully not much slower? How much more complex/involved is the gnomic projection?

Good luck with fixing the second failure mode!

[fluidnumericsJoe]

gnomonic projection actually looks like it will be fewer operations than the normal projection so I don't think there will be performance impacts here. For the other issue, @wyatt-fluidnumerics has got some clever idea for building the domain bounding box from face and edge centers rather than corner nodes. It should be fixable in initialization with a hardly noticeable cost.

[fluidnumericsJoe]

@wyatt-fluidnumerics - reproducing tests have been pushed to the bugfix/issue-2878 branch .