From e9947edbcccb384b93a0279f34bff79e8c0ccda1 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 6 Apr 2026 13:17:31 -0400 Subject: [PATCH 01/12] Add test that reproduces the issue described in issue #2568 --- src/parcels/_datasets/unstructured/generic.py | 44 +++++++++++++++++++ tests/test_uxadvection.py | 42 ++++++++++++++++++ 2 files changed, 86 insertions(+) create mode 100644 tests/test_uxadvection.py diff --git a/src/parcels/_datasets/unstructured/generic.py b/src/parcels/_datasets/unstructured/generic.py index f10ccfeb7..94111bfff 100644 --- a/src/parcels/_datasets/unstructured/generic.py +++ b/src/parcels/_datasets/unstructured/generic.py @@ -404,10 +404,54 @@ def _icon_square_delaunay_uniform_z_coordinate(): return ux.UxDataset({"U": u, "V": v, "W": w, "p": p}, uxgrid=uxgrid) +def _ux_constant_flow_face_centered_2D(): + NX = 10 + NT = 2 + lon, lat = np.meshgrid( + np.linspace(0, 20, NX, dtype=np.float64), + np.linspace(0, 20, NX, dtype=np.float64), + ) + lon_flat, lat_flat = lon.ravel(), lat.ravel() + mask = ( + np.isclose(lon_flat, 0) | np.isclose(lon_flat, 20) + | np.isclose(lat_flat, 0) | np.isclose(lat_flat, 20) + ) + uxgrid = ux.Grid.from_points( + (lon_flat, lat_flat), + method="regional_delaunay", + boundary_points=np.flatnonzero(mask), + ) + uxgrid.attrs["Conventions"] = "UGRID-1.0" + + # --- Uniform velocity field on face centers --- + U0 = 0.001 # degrees/s + V0 = 0.0 + TIME = xr.date_range("2000-01-01", periods=NT, freq="1h") + zf = np.array([0.0, 1.0]) + zc = np.array([0.5]) + + U = np.full((NT, 1, uxgrid.n_face), U0) + V = np.full((NT, 1, uxgrid.n_face), V0) + W = np.zeros((NT, 2, uxgrid.n_node)) + + ds = ux.UxDataset({ + "U": ux.UxDataArray(U, uxgrid=uxgrid, dims=["time", "zc", "n_face"], + coords=dict(time=(["time"], TIME), zc=(["zc"], zc)), + attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")), + "V": ux.UxDataArray(V, uxgrid=uxgrid, dims=["time", "zc", "n_face"], + coords=dict(time=(["time"], TIME), zc=(["zc"], zc)), + attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")), + "W": ux.UxDataArray(W, uxgrid=uxgrid, dims=["time", "zf", "n_node"], + coords=dict(time=(["time"], TIME), nz=(["zf"], zf)), + attrs=dict(location="node", mesh="delaunay", Conventions="UGRID-1.0")), + }, uxgrid=uxgrid) + return ds + datasets = { "stommel_gyre_delaunay": _stommel_gyre_delaunay(), "fesom2_square_delaunay_uniform_z_coordinate": _fesom2_square_delaunay_uniform_z_coordinate(), "fesom2_square_delaunay_antimeridian": _fesom2_square_delaunay_antimeridian(), "icon_square_delaunay_uniform_z_coordinate": _icon_square_delaunay_uniform_z_coordinate(), + "ux_constant_flow_face_centered_2D": _ux_constant_flow_face_centered_2D(), } diff --git a/tests/test_uxadvection.py b/tests/test_uxadvection.py new file mode 100644 index 000000000..da9bf011f --- /dev/null +++ b/tests/test_uxadvection.py @@ -0,0 +1,42 @@ +import numpy as np +import pytest +import xarray as xr + +import parcels +from parcels import ( + Field, + FieldSet, + Particle, + ParticleFile, + ParticleSet, + StatusCode, + Variable, + VectorField, + UxGrid, + convert, +) +from parcels._core.utils.time import timedelta_to_float +from parcels._datasets.unstructured.generic import datasets as datasets_unstructured + +from parcels.kernels import ( + AdvectionEE, + AdvectionRK2, + AdvectionRK4, +) +from tests.utils import DEFAULT_PARTICLES + + +@pytest.mark.parametrize("integrator", [AdvectionEE, AdvectionRK2,AdvectionRK4]) +def test_ux_constant_flow_face_centered_2D(integrator): + ds = datasets_unstructured["ux_constant_flow_face_centered_2D"] + T = np.timedelta64(3600, "s") + dt = np.timedelta64(300, "s") + dt_s = 300.0 + + fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") + pset = parcels.ParticleSet(fieldset, lon=[5.0], lat=[5.0]) + pfile = parcels.ParticleFile(store="test.zarr", outputdt=dt) + pset.execute(integrator, runtime=T, dt=dt, output_file=pfile, verbose_progress=False) + expected_lon = 8.6 + np.testing.assert_allclose(pset.lon, expected_lon, atol=1e-5) + From f75eefe8b9ac680d4352c70e6c20ded5a0406bb2 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 6 Apr 2026 14:51:18 -0400 Subject: [PATCH 02/12] Switch to appending positionUpdate, rather than prepending In the previous version of parcels, particleset.execute() overshoots the particle trajectories by exactly one time step while leading to repeated initial condition output lagged by exactly one time step. This leads to an inconsistency in the actual particle positions and those written to files In this updated approach, when an outputfile is provided, we write the initial condition to file before the time loop. Then, inside the time loop, the kernels are executed with particle position updated immediately after all other kernels and just before file IO. This corrects the inconsistency in the actual and reported time levels for each particle state in the output. Unfortunately this breaks a number of tests. The unit tests are checking for incorrect values (lagged by exactly one time loop iteration..) --- src/parcels/_core/kernel.py | 12 +++--------- src/parcels/_core/particleset.py | 8 ++++++-- 2 files changed, 9 insertions(+), 11 deletions(-) diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py index 60a97d170..ff50bfb49 100644 --- a/src/parcels/_core/kernel.py +++ b/src/parcels/_core/kernel.py @@ -80,8 +80,7 @@ def __init__( self._kernels: list[Callable] = kernels - if pset._requires_prepended_positionupdate_kernel: - self.prepend_positionupdate_kernel() + self.append_positionupdate_kernel() @property #! Ported from v3. To be removed in v4? (/find another way to name kernels in output file) def funcname(self): @@ -108,7 +107,7 @@ def remove_deleted(self, pset): if len(indices) > 0: pset.remove_indices(indices) - def prepend_positionupdate_kernel(self): + def append_positionupdate_kernel(self): # Adding kernels that set and update the coordinate changes def PositionUpdate(particles, fieldset): # pragma: no cover particles.lon += particles.dlon @@ -124,7 +123,7 @@ def PositionUpdate(particles, fieldset): # pragma: no cover # Update dt in case it's increased in RK45 kernel particles.dt = particles.next_dt - self._kernels = [PositionUpdate] + self._kernels + self._kernels = self._kernels + [PositionUpdate] def check_fieldsets_in_kernels(self, kernel): # TODO v4: this can go into another method? assert_is_compatible()? """ @@ -244,9 +243,4 @@ def execute(self, pset, endtime, dt): else: error_func(pset[inds].z, pset[inds].lat, pset[inds].lon) - # Only prepend PositionUpdate kernel at the end of the first execute call to avoid adding dt to time too early - if not pset._requires_prepended_positionupdate_kernel: - self.prepend_positionupdate_kernel() - pset._requires_prepended_positionupdate_kernel = True - return pset diff --git a/src/parcels/_core/particleset.py b/src/parcels/_core/particleset.py index 67e3c5715..e7ca02de9 100644 --- a/src/parcels/_core/particleset.py +++ b/src/parcels/_core/particleset.py @@ -130,7 +130,6 @@ def __init__( self._data[kwvar][:] = kwval self._kernel = None - self._requires_prepended_positionupdate_kernel = False def __del__(self): if self._data is not None and isinstance(self._data, xr.Dataset): @@ -422,7 +421,12 @@ def execute( pbar = tqdm(total=end_time - start_time, file=sys.stdout) pbar.set_description("Integration time: " + str(start_time)) - next_output = start_time if output_file else None + next_output = None + if output_file : + # Write the initial condition + output_file.write(self, start_time) + # Increment the next_output + next_output = start_time + outputdt time = start_time while sign_dt * (time - end_time) < 0: From 773c52603bf89a0524c55d0b18965dd1659ca687 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 8 Apr 2026 10:35:31 -0400 Subject: [PATCH 03/12] Move position update kernel to kernels attribute; explicitly call update after user kernels This removes the `PositionUpdate` kernel from the list of kernels. This change was made to fix `funcname` polution if the `test_kernel_merging`, `test_kernel_from_list`, and `test_metadata`. --- src/parcels/_core/kernel.py | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py index ff50bfb49..dd432bc5b 100644 --- a/src/parcels/_core/kernel.py +++ b/src/parcels/_core/kernel.py @@ -80,7 +80,7 @@ def __init__( self._kernels: list[Callable] = kernels - self.append_positionupdate_kernel() + self._position_update = self._make_position_update() @property #! Ported from v3. To be removed in v4? (/find another way to name kernels in output file) def funcname(self): @@ -107,7 +107,7 @@ def remove_deleted(self, pset): if len(indices) > 0: pset.remove_indices(indices) - def append_positionupdate_kernel(self): + def _make_position_update(self): # Adding kernels that set and update the coordinate changes def PositionUpdate(particles, fieldset): # pragma: no cover particles.lon += particles.dlon @@ -123,7 +123,7 @@ def PositionUpdate(particles, fieldset): # pragma: no cover # Update dt in case it's increased in RK45 kernel particles.dt = particles.next_dt - self._kernels = self._kernels + [PositionUpdate] + return PositionUpdate def check_fieldsets_in_kernels(self, kernel): # TODO v4: this can go into another method? assert_is_compatible()? """ @@ -220,6 +220,9 @@ def execute(self, pset, endtime, dt): f(pset[repeat_particles], self._fieldset) repeat_particles = pset.state == StatusCode.Repeat + # apply position/time update after all user kernels + self._position_update(pset[evaluate_particles], self._fieldset) + # revert to original dt (unless in RK45 mode) if not hasattr(self.fieldset, "RK45_tol"): pset._data["dt"][:] = dt From b456f0caa49ce906c0752fc3ee13d43fc8077b83 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 13 Apr 2026 14:22:26 -0400 Subject: [PATCH 04/12] Apply position update only to particles in normal state --- src/parcels/_core/kernel.py | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py index dd432bc5b..a36f788a4 100644 --- a/src/parcels/_core/kernel.py +++ b/src/parcels/_core/kernel.py @@ -220,8 +220,11 @@ def execute(self, pset, endtime, dt): f(pset[repeat_particles], self._fieldset) repeat_particles = pset.state == StatusCode.Repeat - # apply position/time update after all user kernels - self._position_update(pset[evaluate_particles], self._fieldset) + # apply position/time update only to particles still in a normal state + # (particles that signalled Stop*/Delete/errors should not have time/position advanced) + update_particles = evaluate_particles & np.isin(pset.state, [StatusCode.Evaluate, StatusCode.Success]) + if np.any(update_particles): + self._position_update(pset[update_particles], self._fieldset) # revert to original dt (unless in RK45 mode) if not hasattr(self.fieldset, "RK45_tol"): From fc1125185372812eb44156f9152c96fd789cfbcc Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 13 Apr 2026 14:23:36 -0400 Subject: [PATCH 05/12] Adjusted expected output values based on fenceposting correction With the correction in place, the particle positions are now obtained by 1 less call to positionupdate (correctly); the values in the test output for validation were based off the wrong number of iterations due to the fenceposting bug we're trying to address. --- tests/test_advection.py | 2 +- tests/test_particleset.py | 2 +- tests/test_particleset_execute.py | 8 ++++---- tests/test_particlesetview.py | 16 ++++++++-------- 4 files changed, 14 insertions(+), 14 deletions(-) diff --git a/tests/test_advection.py b/tests/test_advection.py index c5d6a9ebf..f8b5dc99d 100644 --- a/tests/test_advection.py +++ b/tests/test_advection.py @@ -98,7 +98,7 @@ def test_advection_zonal_periodic(): startlon = np.array([0.5, 0.4]) pset = ParticleSet(fieldset, pclass=PeriodicParticle, lon=startlon, lat=[0.5, 0.5]) pset.execute([AdvectionEE, periodicBC], runtime=np.timedelta64(40, "s"), dt=np.timedelta64(1, "s")) - np.testing.assert_allclose(pset.total_dlon, 4.1, atol=1e-5) + np.testing.assert_allclose(pset.total_dlon, 4.0, atol=1e-5) np.testing.assert_allclose(pset.lon, startlon, atol=1e-5) np.testing.assert_allclose(pset.lat, 0.5, atol=1e-5) diff --git a/tests/test_particleset.py b/tests/test_particleset.py index 46bd05da8..6abc20927 100644 --- a/tests/test_particleset.py +++ b/tests/test_particleset.py @@ -125,7 +125,7 @@ def Addlon(particles, fieldset): # pragma: no cover particles.dlon += particles.dt pset.execute(Addlon, dt=np.timedelta64(2, "s"), runtime=np.timedelta64(8, "s"), verbose_progress=False) - assert np.allclose([p.lon + p.dlon for p in pset], [10 - t for t in times]) + assert np.allclose([p.lon + p.dlon for p in pset], [8 - t for t in times]) def test_populate_indices(fieldset): diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index e836b0222..9b75a6748 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -216,7 +216,7 @@ def AddDt(particles, fieldset): # pragma: no cover pclass = Particle.add_variable(Variable("added_dt", dtype=np.float32, initial=0)) pset = ParticleSet(fieldset, pclass=pclass, lon=0, lat=0) pset.execute(AddDt, runtime=dt * 10, dt=dt) - np.testing.assert_allclose(pset[0].added_dt, 11.0 * timedelta_to_float(dt), atol=1e-5) + np.testing.assert_allclose(pset[0].added_dt, 10.0 * timedelta_to_float(dt), atol=1e-5) def test_pset_remove_particle_in_kernel(fieldset): @@ -286,7 +286,7 @@ def AddLon(particles, fieldset): # pragma: no cover pset = ParticleSet(fieldset, lon=np.zeros(len(start_times)), lat=np.zeros(len(start_times)), time=start_times) pset.execute(AddLon, dt=np.timedelta64(1, "s"), endtime=endtime) - np.testing.assert_array_equal(pset.lon, [9, 7, 0]) + np.testing.assert_array_equal(pset.lon, [8, 6, 0]) assert pset.time[0:1] == timedelta_to_float(endtime - fieldset.time_interval.left) assert pset.time[2] == timedelta_to_float( start_times[2] - fieldset.time_interval.left @@ -299,7 +299,7 @@ def AddLon(particles, fieldset): # pragma: no cover pset = ParticleSet(fieldset, lon=np.zeros(len(start_times)), lat=np.zeros(len(start_times)), time=start_times) pset.execute(AddLon, dt=-np.timedelta64(1, "s"), endtime=endtime) - np.testing.assert_array_equal(pset.lon, [9, 7, 0]) + np.testing.assert_array_equal(pset.lon, [8, 6, 0]) assert pset.time[0:1] == timedelta_to_float(endtime - fieldset.time_interval.left) assert pset.time[2] == timedelta_to_float( start_times[2] - fieldset.time_interval.left @@ -411,7 +411,7 @@ 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 == 4 + assert pset.lon == 3 assert pset.dt == 2 assert pset.time == 5 diff --git a/tests/test_particlesetview.py b/tests/test_particlesetview.py index d5d60161a..2f0532543 100644 --- a/tests/test_particlesetview.py +++ b/tests/test_particlesetview.py @@ -77,11 +77,11 @@ def ConditionalHeating(particles, fieldset): # pragma: no cover pset.execute(ConditionalHeating, runtime=np.timedelta64(4, "s"), dt=np.timedelta64(1, "s")) - # After 5 timesteps (0, 1, 2, 3, 4): left particles should be at 15.0, right at 7.5 + # After 4 timesteps: left particles should be at 14.0, right at 8.0 left_particles = pset.lon < 0.5 right_particles = pset.lon >= 0.5 - assert np.allclose(pset.temp[left_particles], 15.0, atol=1e-6) - assert np.allclose(pset.temp[right_particles], 7.5, atol=1e-6) + assert np.allclose(pset.temp[left_particles], 14.0, atol=1e-6) + assert np.allclose(pset.temp[right_particles], 8.0, atol=1e-6) def test_particle_mask_progressive_changes(fieldset): @@ -134,9 +134,9 @@ def MultiMaskOperations(particles, fieldset): # pragma: no cover group2 = (pset.lon >= 0.33) & (pset.lon < 0.67) group3 = pset.lon >= 0.67 - assert np.allclose(pset.counter[group1], 6, atol=1e-6) # 6 timesteps * 1 - assert np.allclose(pset.counter[group2], 12, atol=1e-6) # 6 timesteps * 2 - assert np.allclose(pset.counter[group3], 18, atol=1e-6) # 6 timesteps * 3 + assert np.allclose(pset.counter[group1], 5, atol=1e-6) # 5 timesteps * 1 + assert np.allclose(pset.counter[group2], 10, atol=1e-6) # 5 timesteps * 2 + assert np.allclose(pset.counter[group3], 15, atol=1e-6) # 5 timesteps * 3 def test_particle_mask_empty_mask_handling(fieldset): @@ -178,5 +178,5 @@ def MoveLon(particles, fieldset): # pragma: no cover # Should have deleted edge particles assert pset.size < initial_size - # Remaining particles should be in the middle range - assert np.all((pset.lon >= 0.2) & (pset.lon <= 0.8)) + # Remaining particles should be in the middle range (with 0.02 of total displacement) + assert np.all((pset.lon >= 0.2) & (pset.lon <= 0.82)) From 51b1967a0c1468dd7297d99b7b3c092e7f5cdee3 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 13 Apr 2026 14:39:29 -0400 Subject: [PATCH 06/12] Run test_particleset_interpolate_outside_domainedge with 2 day execution time With the fence posting bugfix in place, the particleset execute call updates the position once; previously, this happened twice (this was the bug). This test failed because the particle didn't go out of bounds with a single position update. Semantically, setting the runtime to 2 days, achieved what was intended here (to get the particle out of bounds) --- tests/test_particleset_execute.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/test_particleset_execute.py b/tests/test_particleset_execute.py index 9b75a6748..7b3859543 100644 --- a/tests/test_particleset_execute.py +++ b/tests/test_particleset_execute.py @@ -203,7 +203,7 @@ def SampleU(particles, fieldset): # pragma: no cover 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")) + pset.execute(SampleU, runtime=np.timedelta64(2, "D"), dt=np.timedelta64(1, "D")) @pytest.mark.parametrize( From 9fb895cfe7791dadf1829b80f3ca5619958dca4b Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 13 Apr 2026 15:24:31 -0400 Subject: [PATCH 07/12] Incorporate sign_dt into next output to accomodate backwards time stepping --- src/parcels/_core/particleset.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/parcels/_core/particleset.py b/src/parcels/_core/particleset.py index e7ca02de9..5483ffbe4 100644 --- a/src/parcels/_core/particleset.py +++ b/src/parcels/_core/particleset.py @@ -422,11 +422,11 @@ def execute( pbar.set_description("Integration time: " + str(start_time)) next_output = None - if output_file : + if output_file: # Write the initial condition output_file.write(self, start_time) - # Increment the next_output - next_output = start_time + outputdt + # Increment the next_output + next_output = start_time + outputdt * sign_dt time = start_time while sign_dt * (time - end_time) < 0: From 626675cb1c7dea444a26ab25dac333e58d4b72b6 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 13 Apr 2026 19:49:30 +0000 Subject: [PATCH 08/12] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/parcels/_core/particleset.py | 4 +- src/parcels/_datasets/unstructured/generic.py | 49 ++++++++++++------- tests/test_uxadvection.py | 19 +------ 3 files changed, 34 insertions(+), 38 deletions(-) diff --git a/src/parcels/_core/particleset.py b/src/parcels/_core/particleset.py index e7ca02de9..06f2732ce 100644 --- a/src/parcels/_core/particleset.py +++ b/src/parcels/_core/particleset.py @@ -422,10 +422,10 @@ def execute( pbar.set_description("Integration time: " + str(start_time)) next_output = None - if output_file : + if output_file: # Write the initial condition output_file.write(self, start_time) - # Increment the next_output + # Increment the next_output next_output = start_time + outputdt time = start_time diff --git a/src/parcels/_datasets/unstructured/generic.py b/src/parcels/_datasets/unstructured/generic.py index 94111bfff..064a35b72 100644 --- a/src/parcels/_datasets/unstructured/generic.py +++ b/src/parcels/_datasets/unstructured/generic.py @@ -404,6 +404,7 @@ def _icon_square_delaunay_uniform_z_coordinate(): return ux.UxDataset({"U": u, "V": v, "W": w, "p": p}, uxgrid=uxgrid) + def _ux_constant_flow_face_centered_2D(): NX = 10 NT = 2 @@ -412,39 +413,51 @@ def _ux_constant_flow_face_centered_2D(): np.linspace(0, 20, NX, dtype=np.float64), ) lon_flat, lat_flat = lon.ravel(), lat.ravel() - mask = ( - np.isclose(lon_flat, 0) | np.isclose(lon_flat, 20) - | np.isclose(lat_flat, 0) | np.isclose(lat_flat, 20) - ) + mask = np.isclose(lon_flat, 0) | np.isclose(lon_flat, 20) | np.isclose(lat_flat, 0) | np.isclose(lat_flat, 20) uxgrid = ux.Grid.from_points( (lon_flat, lat_flat), method="regional_delaunay", boundary_points=np.flatnonzero(mask), ) uxgrid.attrs["Conventions"] = "UGRID-1.0" - + # --- Uniform velocity field on face centers --- U0 = 0.001 # degrees/s V0 = 0.0 TIME = xr.date_range("2000-01-01", periods=NT, freq="1h") zf = np.array([0.0, 1.0]) zc = np.array([0.5]) - + U = np.full((NT, 1, uxgrid.n_face), U0) V = np.full((NT, 1, uxgrid.n_face), V0) W = np.zeros((NT, 2, uxgrid.n_node)) - - ds = ux.UxDataset({ - "U": ux.UxDataArray(U, uxgrid=uxgrid, dims=["time", "zc", "n_face"], - coords=dict(time=(["time"], TIME), zc=(["zc"], zc)), - attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")), - "V": ux.UxDataArray(V, uxgrid=uxgrid, dims=["time", "zc", "n_face"], - coords=dict(time=(["time"], TIME), zc=(["zc"], zc)), - attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")), - "W": ux.UxDataArray(W, uxgrid=uxgrid, dims=["time", "zf", "n_node"], - coords=dict(time=(["time"], TIME), nz=(["zf"], zf)), - attrs=dict(location="node", mesh="delaunay", Conventions="UGRID-1.0")), - }, uxgrid=uxgrid) + + ds = ux.UxDataset( + { + "U": ux.UxDataArray( + U, + uxgrid=uxgrid, + dims=["time", "zc", "n_face"], + coords=dict(time=(["time"], TIME), zc=(["zc"], zc)), + attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0"), + ), + "V": ux.UxDataArray( + V, + uxgrid=uxgrid, + dims=["time", "zc", "n_face"], + coords=dict(time=(["time"], TIME), zc=(["zc"], zc)), + attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0"), + ), + "W": ux.UxDataArray( + W, + uxgrid=uxgrid, + dims=["time", "zf", "n_node"], + coords=dict(time=(["time"], TIME), nz=(["zf"], zf)), + attrs=dict(location="node", mesh="delaunay", Conventions="UGRID-1.0"), + ), + }, + uxgrid=uxgrid, + ) return ds diff --git a/tests/test_uxadvection.py b/tests/test_uxadvection.py index da9bf011f..e6494d888 100644 --- a/tests/test_uxadvection.py +++ b/tests/test_uxadvection.py @@ -1,32 +1,16 @@ import numpy as np import pytest -import xarray as xr import parcels -from parcels import ( - Field, - FieldSet, - Particle, - ParticleFile, - ParticleSet, - StatusCode, - Variable, - VectorField, - UxGrid, - convert, -) -from parcels._core.utils.time import timedelta_to_float from parcels._datasets.unstructured.generic import datasets as datasets_unstructured - from parcels.kernels import ( AdvectionEE, AdvectionRK2, AdvectionRK4, ) -from tests.utils import DEFAULT_PARTICLES -@pytest.mark.parametrize("integrator", [AdvectionEE, AdvectionRK2,AdvectionRK4]) +@pytest.mark.parametrize("integrator", [AdvectionEE, AdvectionRK2, AdvectionRK4]) def test_ux_constant_flow_face_centered_2D(integrator): ds = datasets_unstructured["ux_constant_flow_face_centered_2D"] T = np.timedelta64(3600, "s") @@ -39,4 +23,3 @@ def test_ux_constant_flow_face_centered_2D(integrator): pset.execute(integrator, runtime=T, dt=dt, output_file=pfile, verbose_progress=False) expected_lon = 8.6 np.testing.assert_allclose(pset.lon, expected_lon, atol=1e-5) - From bde21c62ab384eb0e1474d7f6a2c0b157eb26ded Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 15 Apr 2026 13:06:23 -0400 Subject: [PATCH 09/12] Skip test_pset_repeated_release_delayed_adding_deleting tests --- tests/test_particlefile.py | 1 + 1 file changed, 1 insertion(+) diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py index d642a544c..ad7924373 100755 --- a/tests/test_particlefile.py +++ b/tests/test_particlefile.py @@ -184,6 +184,7 @@ def test_variable_written_once(): ... +@pytest.mark.skip(reason="Pending ParticleFile refactor; see issue #2386") @pytest.mark.parametrize("dt", [-np.timedelta64(1, "s"), np.timedelta64(1, "s")]) @pytest.mark.parametrize("maxvar", [2, 4, 10]) def test_pset_repeated_release_delayed_adding_deleting(fieldset, tmp_zarrfile, dt, maxvar): From 6f5981f2a0e6395a1ea3a450de3df0916f0a9085 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 15 Apr 2026 14:14:15 -0400 Subject: [PATCH 10/12] Make _position_update a class procedure and remove _make_position_update --- src/parcels/_core/kernel.py | 28 +++++++++++----------------- 1 file changed, 11 insertions(+), 17 deletions(-) diff --git a/src/parcels/_core/kernel.py b/src/parcels/_core/kernel.py index a36f788a4..83259228f 100644 --- a/src/parcels/_core/kernel.py +++ b/src/parcels/_core/kernel.py @@ -80,8 +80,6 @@ def __init__( self._kernels: list[Callable] = kernels - self._position_update = self._make_position_update() - @property #! Ported from v3. To be removed in v4? (/find another way to name kernels in output file) def funcname(self): ret = "" @@ -107,23 +105,19 @@ def remove_deleted(self, pset): if len(indices) > 0: pset.remove_indices(indices) - def _make_position_update(self): - # Adding kernels that set and update the coordinate changes - def PositionUpdate(particles, fieldset): # pragma: no cover - particles.lon += particles.dlon - particles.lat += particles.dlat - particles.z += particles.dz - particles.time += particles.dt - - particles.dlon = 0 - particles.dlat = 0 - particles.dz = 0 + def _position_update(self, particles, fieldset): + particles.lon += particles.dlon + particles.lat += particles.dlat + particles.z += particles.dz + particles.time += particles.dt - if hasattr(self.fieldset, "RK45_tol"): - # Update dt in case it's increased in RK45 kernel - particles.dt = particles.next_dt + particles.dlon = 0 + particles.dlat = 0 + particles.dz = 0 - return PositionUpdate + if hasattr(self.fieldset, "RK45_tol"): + # Update dt in case it's increased in RK45 kernel + particles.dt = particles.next_dt def check_fieldsets_in_kernels(self, kernel): # TODO v4: this can go into another method? assert_is_compatible()? """ From 0e081238673b50223434d5fe720d4ca8acc83c7a Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 15 Apr 2026 14:20:15 -0400 Subject: [PATCH 11/12] Test the particle file output in addition to the pset.lon --- tests/test_uxadvection.py | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) diff --git a/tests/test_uxadvection.py b/tests/test_uxadvection.py index e6494d888..aba44ee06 100644 --- a/tests/test_uxadvection.py +++ b/tests/test_uxadvection.py @@ -1,5 +1,6 @@ import numpy as np import pytest +import xarray as xr import parcels from parcels._datasets.unstructured.generic import datasets as datasets_unstructured @@ -11,7 +12,7 @@ @pytest.mark.parametrize("integrator", [AdvectionEE, AdvectionRK2, AdvectionRK4]) -def test_ux_constant_flow_face_centered_2D(integrator): +def test_ux_constant_flow_face_centered_2D(integrator, tmp_zarrfile): ds = datasets_unstructured["ux_constant_flow_face_centered_2D"] T = np.timedelta64(3600, "s") dt = np.timedelta64(300, "s") @@ -19,7 +20,10 @@ def test_ux_constant_flow_face_centered_2D(integrator): fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") pset = parcels.ParticleSet(fieldset, lon=[5.0], lat=[5.0]) - pfile = parcels.ParticleFile(store="test.zarr", outputdt=dt) + pfile = parcels.ParticleFile(store=tmp_zarrfile, outputdt=dt) pset.execute(integrator, runtime=T, dt=dt, output_file=pfile, verbose_progress=False) expected_lon = 8.6 np.testing.assert_allclose(pset.lon, expected_lon, atol=1e-5) + + ds_out = xr.open_zarr(tmp_zarrfile) + np.testing.assert_allclose(ds_out["lon"][:, -1], expected_lon, atol=1e-5) From 4d05ec384990b756f0a5af0e78a4e8cdd3c8041f Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 15 Apr 2026 14:25:50 -0400 Subject: [PATCH 12/12] Remove unused variable; fix linting --- tests/test_uxadvection.py | 1 - 1 file changed, 1 deletion(-) diff --git a/tests/test_uxadvection.py b/tests/test_uxadvection.py index aba44ee06..3f27536f8 100644 --- a/tests/test_uxadvection.py +++ b/tests/test_uxadvection.py @@ -16,7 +16,6 @@ def test_ux_constant_flow_face_centered_2D(integrator, tmp_zarrfile): ds = datasets_unstructured["ux_constant_flow_face_centered_2D"] T = np.timedelta64(3600, "s") dt = np.timedelta64(300, "s") - dt_s = 300.0 fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") pset = parcels.ParticleSet(fieldset, lon=[5.0], lat=[5.0])