From 4d03c0e4f62e791a7074f479e840b13381f4b37f Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 24 Oct 2025 08:22:44 +0200 Subject: [PATCH 01/13] Adding unit tests that seemed to break on VirtualShip But work on v4-dev itself. Need to explore further --- tests/test_interpolation.py | 6 ++++++ tests/test_particleset_execute.py | 22 ++++++++++++++++++++++ 2 files changed, 28 insertions(+) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index 20ba10d1c5..d8f85d1948 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -91,6 +91,12 @@ def test_raw_2d_interpolation(field, func, t, z, y, x, expected): np.testing.assert_equal(value, expected) +def test_scalar_field_eval(field): + UV = VectorField("UV", field, field) + + UV.eval(np.timedelta64(2, "s"), 2, 1.5, 0.5) + + @pytest.mark.parametrize( "func, t, z, y, x, expected", [ diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index 163a2217eb..46515ece38 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -150,6 +150,28 @@ def test_particleset_endtime_type(fieldset, endtime, expectation): pset.execute(endtime=endtime, dt=np.timedelta64(10, "m"), pyfunc=DoNothing) +def test_particleset_run_to_endtime(fieldset): + starttime = fieldset.time_interval.left + endtime = fieldset.time_interval.right + + def SampleU(particles, fieldset): # pragma: no cover + _ = fieldset.U[particles] + + pset = ParticleSet(fieldset, lon=[0.2], lat=[5.0], time=[starttime]) + pset.execute(SampleU, endtime=endtime, dt=np.timedelta64(1, "D")) + assert pset[0].time + pset[0].dt == endtime + + +def test_particleset_interpolate_domainedge(fieldset): + def SampleU(particles, fieldset): # pragma: no cover + particles.dlon = fieldset.U[particles] + + print(fieldset.U.grid.lon) + pset = ParticleSet(fieldset, lon=fieldset.U.grid.lon[0], lat=fieldset.U.grid.lat[0]) + pset.execute(SampleU, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(1, "D")) + assert np.isfinite(pset[0].dlon) + + @pytest.mark.parametrize( "dt", [np.timedelta64(1, "s"), np.timedelta64(1, "ms"), np.timedelta64(10, "ms"), np.timedelta64(1, "ns")] ) From 833c6925237ab244efae682e0ec2412dd96bbb5c Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 24 Oct 2025 08:22:58 +0200 Subject: [PATCH 02/13] Adding note on gridset to migration guide --- docs/community/v4-migration.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/community/v4-migration.md b/docs/community/v4-migration.md index 81991c63fd..0fd553c6c3 100644 --- a/docs/community/v4-migration.md +++ b/docs/community/v4-migration.md @@ -34,3 +34,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat - Particlefiles should be created by `ParticleFile(...)` instead of `pset.ParticleFile(...)` - The `name` argument in `ParticleFile` has been replaced by `store` and can now be a string, a Path or a zarr store. + +## GridSet + +- `GridSet` is now a list, so change `fieldset.gridset.grids[0]` to `fieldset.gridset[0]`. From ced1413d52d357caf52aa7426d76052949443aeb Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 24 Oct 2025 11:32:08 +0200 Subject: [PATCH 03/13] Updating tests with breaking test at right boundaries --- tests/test_particleset_execute.py | 22 ++++++++++++++++++---- 1 file changed, 18 insertions(+), 4 deletions(-) diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index 46515ece38..c022fced05 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -162,14 +162,28 @@ def SampleU(particles, fieldset): # pragma: no cover assert pset[0].time + pset[0].dt == endtime -def test_particleset_interpolate_domainedge(fieldset): +def test_particleset_interpolate_on_domainedge(zonal_flow_fieldset): + fieldset = zonal_flow_fieldset + def SampleU(particles, fieldset): # pragma: no cover particles.dlon = fieldset.U[particles] - print(fieldset.U.grid.lon) - pset = ParticleSet(fieldset, lon=fieldset.U.grid.lon[0], lat=fieldset.U.grid.lat[0]) + pset = ParticleSet(fieldset, lon=fieldset.U.grid.lon[-1], lat=fieldset.U.grid.lat[-1]) pset.execute(SampleU, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(1, "D")) - assert np.isfinite(pset[0].dlon) + np.testing.assert_equal(pset[0].dlon, 1) + + +def test_particleset_interpolate_outside_domainedge(zonal_flow_fieldset): + fieldset = zonal_flow_fieldset + + def SampleU(particles, fieldset): # pragma: no cover + particles.dlon = fieldset.U[particles] + + dlat = 1e-3 + pset = ParticleSet(fieldset, lon=fieldset.U.grid.lon[-1], lat=fieldset.U.grid.lat[-1] + dlat) + + with pytest.raises(FieldOutOfBoundError): + pset.execute(SampleU, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(1, "D")) @pytest.mark.parametrize( From 2a893e0f9dd921df6f47728bdb7270100ba4b87f Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 24 Oct 2025 11:36:52 +0200 Subject: [PATCH 04/13] Fixing index_search to support particles exactly at right side --- src/parcels/_core/index_search.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/index_search.py b/src/parcels/_core/index_search.py index dcf4094518..6ff1cdb54f 100644 --- a/src/parcels/_core/index_search.py +++ b/src/parcels/_core/index_search.py @@ -44,7 +44,7 @@ def _search_1d_array( # TODO v4: We probably rework this to deal with 0D arrays before this point (as we already know field dimensionality) if len(arr) < 2: return np.zeros(shape=x.shape, dtype=np.int32), np.zeros_like(x) - index = np.searchsorted(arr, x, side="right") - 1 + index = np.clip(np.searchsorted(arr, x, side="right") - 1, 0, len(arr) - 2) # Use broadcasting to avoid repeated array access arr_index = arr[index] arr_next = arr[np.clip(index + 1, 1, len(arr) - 1)] # Ensure we don't go out of bounds @@ -57,7 +57,7 @@ def _search_1d_array( # bcoord = (x - arr[index]) / dx index = np.where(x < arr[0], LEFT_OUT_OF_BOUNDS, index) - index = np.where(x >= arr[-1], RIGHT_OUT_OF_BOUNDS, index) + index = np.where(x > arr[-1], RIGHT_OUT_OF_BOUNDS, index) return np.atleast_1d(index), np.atleast_1d(bcoord) From dc41fe76ae70cbc48e83eca96955ccef20024ebc Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 24 Oct 2025 11:40:43 +0200 Subject: [PATCH 05/13] Removing redundant print statements That were old debug statements --- tests/test_particleset_execute.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index c022fced05..c3dfb53349 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -365,7 +365,6 @@ def MoveRight(particles, fieldset): # pragma: no cover def MoveLeft(particles, fieldset): # pragma: no cover inds = np.where(particles.state == StatusCode.ErrorOutOfBounds) - print(inds, particles.state) particles[inds].dlon -= 1.0 particles[inds].state = StatusCode.Success @@ -400,7 +399,6 @@ def KernelCounter(particles, fieldset): # pragma: no cover pset = ParticleSet(fieldset, lon=np.zeros(1), lat=np.zeros(1)) pset.execute(KernelCounter, dt=np.timedelta64(2, "s"), runtime=np.timedelta64(5, "s")) assert pset.lon == 3 - print(pset.dt) assert pset.dt == np.timedelta64(2, "s") From b2da378878b447c4fbd5e3611dc52e6148d3c898 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 24 Oct 2025 12:49:14 +0200 Subject: [PATCH 06/13] Add field.eval info in migration guide --- docs/community/v4-migration.md | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/docs/community/v4-migration.md b/docs/community/v4-migration.md index 0fd553c6c3..3eccca5976 100644 --- a/docs/community/v4-migration.md +++ b/docs/community/v4-migration.md @@ -35,6 +35,10 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat - Particlefiles should be created by `ParticleFile(...)` instead of `pset.ParticleFile(...)` - The `name` argument in `ParticleFile` has been replaced by `store` and can now be a string, a Path or a zarr store. +## Field + +- `Field.eval()` returns an array of floats instead of a single float (related to the vectorization) + ## GridSet - `GridSet` is now a list, so change `fieldset.gridset.grids[0]` to `fieldset.gridset[0]`. From f3c59ab742caee8851f1221ccb2ef0fc671f6cfd Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 24 Oct 2025 17:03:05 +0200 Subject: [PATCH 07/13] Adding check on decreasing depth --- src/parcels/_core/xgrid.py | 10 +++++++++- tests/test_xgrid.py | 8 ++++++++ 2 files changed, 17 insertions(+), 1 deletion(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index d20d53f8b3..4e3d82daea 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -104,9 +104,12 @@ def __init__(self, grid: xgcm.Grid, mesh="flat"): if "lat" in ds: ds.set_coords("lat") - if len(set(grid.axes) & {"X", "Y", "Z"}) > 0: # Only if spatial grid is >0D (see #2054 for further development) + if len(set(grid.axes) & {"X", "Y"}) > 0: # Only if spatial grid is >0D (see #2054 for further development) assert_valid_lat_lon(ds["lat"], ds["lon"], grid.axes) + if "Z" in grid.axes: + assert_valid_depth(ds["depth"]) + assert_valid_mesh(mesh) self._ds = ds @@ -474,6 +477,11 @@ def assert_valid_lat_lon(da_lat, da_lon, axes: _XGCM_AXES): ) +def assert_valid_depth(da_depth): + if not np.all(np.diff(da_depth.values) > 0): + raise ValueError(f"Depth DataArray {da_depth.name!r} with dims {da_depth.dims} must be strictly increasing.") + + def _convert_center_pos_to_fpoint( *, index: int, bcoord: float, xgcm_position: _XGCM_AXIS_POSITION, f_points_xgcm_position: _XGCM_AXIS_POSITION ) -> tuple[int, float]: diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py index 40a203db54..2132d2af49 100644 --- a/tests/test_xgrid.py +++ b/tests/test_xgrid.py @@ -126,6 +126,14 @@ def test_invalid_lon_lat(): XGrid.from_dataset(ds) +def test_invalid_depth(): + ds = datasets["ds_2d_left"].copy() + ds = ds.reindex({"ZG": -ds.ZG}) + + with pytest.raises(ValueError, match="Depth DataArray .* must be strictly increasing*"): + XGrid.from_dataset(ds) + + @pytest.mark.parametrize( "ds", [ From 57d53bec5f52b88b09fb3e5c2b980ab5f1a1de0e Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 27 Oct 2025 15:29:26 +0100 Subject: [PATCH 08/13] Updating decreasing depth error message with HINT --- src/parcels/_core/xgrid.py | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/parcels/_core/xgrid.py b/src/parcels/_core/xgrid.py index 4e3d82daea..85b833fa55 100644 --- a/src/parcels/_core/xgrid.py +++ b/src/parcels/_core/xgrid.py @@ -479,7 +479,10 @@ def assert_valid_lat_lon(da_lat, da_lon, axes: _XGCM_AXES): def assert_valid_depth(da_depth): if not np.all(np.diff(da_depth.values) > 0): - raise ValueError(f"Depth DataArray {da_depth.name!r} with dims {da_depth.dims} must be strictly increasing.") + raise ValueError( + f"Depth DataArray {da_depth.name!r} with dims {da_depth.dims} must be strictly increasing. " + f'HINT: you may be able to use ds.reindex to flip depth - e.g., ds = ds.reindex({da_depth.name}=ds["{da_depth.name}"][::-1])' + ) def _convert_center_pos_to_fpoint( From 53286e62682d400e771eba05e1181551990c6420 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 27 Oct 2025 15:29:54 +0100 Subject: [PATCH 09/13] Adding info that field.eval does not throw an error message to migration guide --- docs/community/v4-migration.md | 1 + 1 file changed, 1 insertion(+) diff --git a/docs/community/v4-migration.md b/docs/community/v4-migration.md index 3eccca5976..2884049b10 100644 --- a/docs/community/v4-migration.md +++ b/docs/community/v4-migration.md @@ -38,6 +38,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat ## Field - `Field.eval()` returns an array of floats instead of a single float (related to the vectorization) +- `Field.eval()` does not throw OutOfBounds or other errors ## GridSet From eae899c5253559e5da6db916a0246f332a4f660d Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 27 Oct 2025 16:34:54 +0100 Subject: [PATCH 10/13] Update tests/test_xgrid.py Co-authored-by: Nick Hodgskin <36369090+VeckoTheGecko@users.noreply.github.com> --- tests/test_xgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py index 2132d2af49..7f3b85cdc2 100644 --- a/tests/test_xgrid.py +++ b/tests/test_xgrid.py @@ -128,7 +128,7 @@ def test_invalid_lon_lat(): def test_invalid_depth(): ds = datasets["ds_2d_left"].copy() - ds = ds.reindex({"ZG": -ds.ZG}) + ds = ds.reindex({"ZG": ds.ZG[::-1}) with pytest.raises(ValueError, match="Depth DataArray .* must be strictly increasing*"): XGrid.from_dataset(ds) From 9c2a7b5c4810b46b7819116f138d6a9a1cb4a13f Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Mon, 27 Oct 2025 16:39:02 +0100 Subject: [PATCH 11/13] Fixing parenthesis --- tests/test_xgrid.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_xgrid.py b/tests/test_xgrid.py index 7f3b85cdc2..361d027168 100644 --- a/tests/test_xgrid.py +++ b/tests/test_xgrid.py @@ -128,7 +128,7 @@ def test_invalid_lon_lat(): def test_invalid_depth(): ds = datasets["ds_2d_left"].copy() - ds = ds.reindex({"ZG": ds.ZG[::-1}) + ds = ds.reindex({"ZG": ds.ZG[::-1]}) with pytest.raises(ValueError, match="Depth DataArray .* must be strictly increasing*"): XGrid.from_dataset(ds) From 2faa3c0e219ae7c891cacbe139a718b44599ac08 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 30 Oct 2025 13:58:35 +0100 Subject: [PATCH 12/13] Removing redundant unit test --- tests/test_interpolation.py | 6 ------ 1 file changed, 6 deletions(-) diff --git a/tests/test_interpolation.py b/tests/test_interpolation.py index d8f85d1948..20ba10d1c5 100644 --- a/tests/test_interpolation.py +++ b/tests/test_interpolation.py @@ -91,12 +91,6 @@ def test_raw_2d_interpolation(field, func, t, z, y, x, expected): np.testing.assert_equal(value, expected) -def test_scalar_field_eval(field): - UV = VectorField("UV", field, field) - - UV.eval(np.timedelta64(2, "s"), 2, 1.5, 0.5) - - @pytest.mark.parametrize( "func, t, z, y, x, expected", [ From e01adfcd0cd92caa7d3f5fd2f51c623ec8165558 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 31 Oct 2025 15:22:08 +0100 Subject: [PATCH 13/13] fixing unit tests for new kernel loop --- src/parcels/_core/index_search.py | 4 ++-- tests/test_particleset_execute.py | 11 +++++++---- 2 files changed, 9 insertions(+), 6 deletions(-) diff --git a/src/parcels/_core/index_search.py b/src/parcels/_core/index_search.py index 91e9ef339e..656657e5bc 100644 --- a/src/parcels/_core/index_search.py +++ b/src/parcels/_core/index_search.py @@ -85,8 +85,8 @@ def _search_time_index(field: Field, time: datetime): if not field.time_interval.is_all_time_in_interval(time): _raise_outside_time_interval_error(time, field=None) - ti = np.searchsorted(field.data.time.data, time, side="right") - 1 - tau = (time - field.data.time.data[ti]) / (field.data.time.data[ti + 1] - field.data.time.data[ti]) + ti, tau = _search_1d_array(field.data.time.data, time) + return {"T": {"index": np.atleast_1d(ti), "bcoord": np.atleast_1d(tau)}} diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index b4111658f5..d33fb45963 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -159,18 +159,21 @@ def SampleU(particles, fieldset): # pragma: no cover pset = ParticleSet(fieldset, lon=[0.2], lat=[5.0], time=[starttime]) pset.execute(SampleU, endtime=endtime, dt=np.timedelta64(1, "D")) - assert pset[0].time + pset[0].dt == endtime + assert pset[0].time == endtime def test_particleset_interpolate_on_domainedge(zonal_flow_fieldset): fieldset = zonal_flow_fieldset + MyParticle = Particle.add_variable(Variable("var")) + def SampleU(particles, fieldset): # pragma: no cover - particles.dlon = fieldset.U[particles] + particles.var = fieldset.U[particles] - pset = ParticleSet(fieldset, lon=fieldset.U.grid.lon[-1], lat=fieldset.U.grid.lat[-1]) + print(fieldset.U.grid.lon) + pset = ParticleSet(fieldset, pclass=MyParticle, lon=fieldset.U.grid.lon[-1], lat=fieldset.U.grid.lat[-1]) pset.execute(SampleU, runtime=np.timedelta64(1, "D"), dt=np.timedelta64(1, "D")) - np.testing.assert_equal(pset[0].dlon, 1) + np.testing.assert_equal(pset[0].var, 1) def test_particleset_interpolate_outside_domainedge(zonal_flow_fieldset):