Skip to content

Add FACE to dims to check in _update_particle_states_position #2688

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
fluidnumericsJoe merged 1 commit into from
Jun 19, 2026
+26 −2
Merged

Add FACE to dims to check in _update_particle_states_position #2688

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/parcels/_core/field.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -337,8 +337,8 @@ def _update_particles_ei(particles, grid_positions: dict, field: Field):

def _update_particle_states_position(particles, grid_positions: dict):
"""Update the particle states based on the position dictionary."""
if particles: # TODO also support uxgrid search
for dim in ["X", "Y"]:
if particles:
for dim in ["X", "Y", "FACE"]:
if dim in grid_positions:
particles.state = np.maximum(
np.where(grid_positions[dim]["index"] == -1, StatusCode.ErrorOutOfBounds, particles.state),
Expand Down
24 changes: 24 additions & 0 deletions tests/test_particleset_execute.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,7 @@
Variable,
VectorField,
)
from parcels._core.statuscodes import GridSearchingError
from parcels._core.utils.time import timedelta_to_float
from parcels._datasets.structured.generated import simple_UV_dataset
from parcels._datasets.structured.generic import datasets as datasets_structured
Expand Down Expand Up @@ -592,3 +593,26 @@ def test_uxstommelgyre_pset_execute_output():
pset.execute(
runtime=np.timedelta64(10, "m"), dt=np.timedelta64(60, "s"), kernels=AdvectionEE, output_file=output_file
)


def test_uxgrid_particle_leaving_domain_raises():
"""A particle advecting out of an unstructured (UxGrid) domain must be flagged.

Once the particle crosses the outflow boundary the FACE search returns
``GRID_SEARCH_ERROR`` and ``pset.execute`` should raise a GridSearchingError.
"""
ds = datasets_unstructured["ux_constant_flow_face_centered_2D"]
grid = UxGrid(grid=ds.uxgrid, z=ds.coords["zc"], mesh="flat")
U = Field(name="U", data=ds.U, grid=grid, interp_method=UxConstantFaceConstantZC)
V = Field(name="V", data=ds.V, grid=grid, interp_method=UxConstantFaceConstantZC)
UV = VectorField(name="UV", U=U, V=V, interp_method=Ux_Velocity)
fieldset = FieldSet([UV, UV.U, UV.V])

# Uniform eastward flow (U0 = 0.001 deg/s); release 0.1 deg inside the eastern
# outflow boundary (domain spans lon, lat in [0, 20]). The particle crosses
# lon=20 after 100 s of travel.
lon_max = float(ds.uxgrid.node_lon.max())
pset = ParticleSet(fieldset, lon=[lon_max - 0.1], lat=[10.0], z=[0.5], pclass=Particle)

with pytest.raises(GridSearchingError):

@erikvansebille erikvansebille Jun 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does this not throw an ErrorOutOfBounds? Would that not be more similar to how it's done in Sgrid?

@fluidnumericsJoe fluidnumericsJoe Jun 19, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The grid search returns a value of GRID_SEARCH_ERROR on the face index when it's not found in the grid; this is done in spatialhash.query()

@fluidnumericsJoe fluidnumericsJoe Jun 19, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

ErrorOutOfBounds is thrown when the search index is -1

GRID_SEARCH_ERROR=-3 which triggers this line

@fluidnumericsJoe fluidnumericsJoe Jun 19, 2026

Copy link
Copy Markdown
Contributor Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why does the Sgrid return a value of grid index -1, corresponding to RIGHT_OUT_OF_BOUNDS ?

GRID_SEARCH_ERROR = -3
LEFT_OUT_OF_BOUNDS = -2
RIGHT_OUT_OF_BOUNDS = -1

@erikvansebille erikvansebille Jun 19, 2026

Copy link
Copy Markdown
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ah I see; I was confused

We use these LEFT and RIGHT only in the 1d_search:

Parcels/src/parcels/_core/index_search.py

Lines 59 to 60 in 481decc

index = np.where(x < arr[0], LEFT_OUT_OF_BOUNDS, index)
index = np.where(x > arr[-1], RIGHT_OUT_OF_BOUNDS, index)

This is especially useful for z/depth, where users may want to distinguish between a particle moving into the air or into the ground. And it just so happens we can also use it on rectilinear grids in the horizontal, because these use 1d_search

On CurviLinear grids, Xgrid would also not be able to disinguish between LEFT and RIGHT. So your test is fine as is

pset.execute(AdvectionEE, runtime=np.timedelta64(120, "s"), dt=np.timedelta64(10, "s"))
Loading