FEAT: Add mesh properties, fail on mismatched meshes

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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

[VeckoTheGecko]

Description

Changes:

• Added mesh property to ModelData base class
• Added mesh property to FieldSet
• Updated FieldSet init to call assert_models_have_same_mesh (which raises IncompatibleMeshesException)
• Add test

Checklist

• Contributes to Error out when combining fieldsets of different mesh types #2844 (followup from Adds a warning whenever a _mesh is different for a field when added to a FieldSet, or during Fieldset creation #2761 )
• Tests added
• This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

None used

[VeckoTheGecko]

if you have capacity @wyatt-fluidnumerics , it would be great to have your review

[VeckoTheGecko]

With this change we can't have fieldsets with only constant fields (which is very reasonable), since that would result in an ambigious fieldset mesh.

Hence I also removed test_fieldset_time_interval_constant_fields

[VeckoTheGecko]

I realise now that our docs relied on "empty, constant only" fieldsets.

docs/user_guide/examples/tutorial_interaction.ipynb

I assume that you wrote this @erikvansebille . I'll postpone merging, and leave this open until you're back

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[VeckoTheGecko]

If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

I'm not quite sure I understand this - I feel that this is scope creep for the project. Why would users want to use Parcels as an engine for agent based modelling? (instead of other libraries in Python for this)

[erikvansebille]

OK, you're right that this may be project scope creep. But what do we do then about the particle-particle interaction tutorial? That just needs a flat Field with a constant value for diffusivity everywhere

[VeckoTheGecko]

Perhaps we can recontextualize it in the context of FADs or a usecase with flow fields? Or just add an grid with no attached data to show off the interactivity, and mention that you would want to specify the flow fields yourself as well.

Would you like to work on a diff or would you like me to?

[erikvansebille]

I just pushed 2bd7bf7 - which now uses a very simple UV dataset. Hope this fixes the issue

[VeckoTheGecko]

That seems to have worked

[wyatt-fluidnumerics]

This all looks great to me. I think this would also be a good time to remove the mesh argument from add_constant_field like you mentioned in the original issue.

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

[VeckoTheGecko]

remove the mesh argument from add_constant_field

Yeah sure, I can wrap that into here :)

Actually, I remember why I didn't - I wanted to do it after #2845 as that would make things easier (I also wanted to expand testing for add_constant_field since it doesn't work at the moment for SphericalMesh(some_custom_radius)

Separate PR would be better I think

[erikvansebille]

Looks good - except that I think we should support empty FieldSets. but that could perhaps be annother PR? See below

[erikvansebille]

Actually, I think there are use cases where we have empty fieldsets - like indeed in the particle-particle interaction tutorial. If users want to use parcels as an Agent-Based-Model engine, they typically don't have/use Fields.

How problematic is it to support empty fieldsets? Perhaps we should change the interaction-tutorial to not have any fields at all (instead of only constant fields)? Or will that be a major endeavour (for another PR)?

[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

importnumpyasnpimportuxarrayasuximportxarrayasxrfromshapely.geometryimportPoint, Polygonfromtriangleimporttriangulatefromparcels._core.index_searchimportuxgrid_point_in_cellfromparcels._core.uxgridimportUxGrid# Build a multi-face constrained Delaunay triangulation from three nested polygonspolygons= [
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=0forpolyinpolygons:
Ni=len(poly)
verts.extend(poly.tolist())
segments.extend([[offset+j, offset+ ((j+1) %Ni)] forjinrange(Ni)])
offset+=NiB=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()
forlon, latin [(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 checkforfiinrange(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]),
)
ifabs(coords[0].sum() -1.0) <0.05: # close but fails the rtol=1e-3 checkprint(
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