From ffddc51379290d4fdf63e09d6ca61a7a4efda9eb Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Thu, 25 Sep 2025 22:23:09 +0200 Subject: [PATCH 1/3] Adding unit test for subsecond dt This addresses #2241 --- tests/v4/test_particleset_execute.py | 15 +++++++++++++++ 1 file changed, 15 insertions(+) diff --git a/tests/v4/test_particleset_execute.py b/tests/v4/test_particleset_execute.py index 241868dd6..34f7be4ea 100644 --- a/tests/v4/test_particleset_execute.py +++ b/tests/v4/test_particleset_execute.py @@ -12,6 +12,7 @@ ParticleSet, StatusCode, UXPiecewiseConstantFace, + Variable, VectorField, ) from parcels._datasets.structured.generated import simple_UV_dataset @@ -150,6 +151,20 @@ def test_particleset_endtime_type(fieldset, endtime, expectation): pset.execute(endtime=endtime, dt=np.timedelta64(10, "m"), pyfunc=DoNothing) +@pytest.mark.parametrize( + "dt", [np.timedelta64(1, "s"), np.timedelta64(1, "ms"), np.timedelta64(10, "ms"), np.timedelta64(1, "ns")] +) +def test_pset_execute_subsecond_dt(fieldset, dt): + def AddDt(particles, fieldset): # pragma: no cover + dt = particles.dt / np.timedelta64(1, "s") + particles.added_dt += dt + + 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, 10.0 * dt / np.timedelta64(1, "s"), atol=1e-5) + + def test_pset_remove_particle_in_kernel(fieldset): npart = 100 pset = ParticleSet(fieldset, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart)) From 9c2e2d20153068f4b7125a3510d86eac010687f9 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 26 Sep 2025 08:36:03 +0200 Subject: [PATCH 2/3] Logic for checking dtype of particledata.dt against input dt --- parcels/particleset.py | 26 ++++++++++++++++++++++++-- tests/v4/test_particleset_execute.py | 7 +++++++ 2 files changed, 31 insertions(+), 2 deletions(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 03b2246c1..530079a71 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -425,6 +425,19 @@ def _num_error_particles(self): """ return np.sum(np.isin(self._data["state"], [StatusCode.Success, StatusCode.Evaluate], invert=True)) + def update_dt_dtype(self, dt_dtype: np.dtype): + """Update the dtype of dt + + Parameters + ---------- + dt_dtype : np.dtype + New dtype for dt. + """ + if dt_dtype not in [np.timedelta64, "timedelta64[ns]", "timedelta64[ms]", "timedelta64[s]"]: + raise ValueError(f"dt_dtype must be a numpy timedelta64 dtype. Got {dt_dtype=!r}") + + self._data["dt"] = self._data["dt"].astype(dt_dtype) + def set_variable_write_status(self, var, write_status): """Method to set the write status of a Variable. @@ -500,6 +513,17 @@ def execute( except (ValueError, AssertionError) as e: raise ValueError(f"dt must be a non-zero datetime.timedelta or np.timedelta64 object, got {dt=!r}") from e + # Check if particle dt has finer resolution than input dt + particle_resolution = np.timedelta64(1, np.datetime_data(self._data["dt"].dtype)) + input_resolution = np.timedelta64(1, np.datetime_data(dt.dtype)) + + if input_resolution >= particle_resolution: + self._data["dt"][:] = dt + else: + raise ValueError( + f"The dtype of dt ({dt.dtype}) is coarser than the dtype of the particle dt ({self._data['dt'].dtype}). Please use Particle.set_dt_dtype() to provide a dt with at least the same precision as the particle dt." + ) + if runtime is not None: try: runtime = maybe_convert_python_timedelta_to_numpy(runtime) @@ -508,8 +532,6 @@ def execute( f"The runtime must be a datetime.timedelta or np.timedelta64 object. Got {type(runtime)}" ) from e - self._data["dt"][:] = dt - start_time, end_time = _get_simulation_start_and_end_times( self.fieldset.time_interval, self._data["time_nextloop"], runtime, endtime, sign_dt ) diff --git a/tests/v4/test_particleset_execute.py b/tests/v4/test_particleset_execute.py index 34f7be4ea..fbad3d5ec 100644 --- a/tests/v4/test_particleset_execute.py +++ b/tests/v4/test_particleset_execute.py @@ -161,10 +161,17 @@ 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.update_dt_dtype(dt.dtype) pset.execute(AddDt, runtime=dt * 10, dt=dt) np.testing.assert_allclose(pset[0].added_dt, 10.0 * dt / np.timedelta64(1, "s"), atol=1e-5) +def test_pset_execute_subsecond_dt_error(fieldset): + pset = ParticleSet(fieldset, lon=0, lat=0) + with pytest.raises(ValueError, match="The dtype of dt"): + pset.execute(DoNothing, runtime=np.timedelta64(10, "ms"), dt=np.timedelta64(1, "ms")) + + def test_pset_remove_particle_in_kernel(fieldset): npart = 100 pset = ParticleSet(fieldset, lon=np.linspace(0, 1, npart), lat=np.linspace(1, 0, npart)) From 061f99642ba659dc48010a1c92928f79a50c7923 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 26 Sep 2025 08:59:46 +0200 Subject: [PATCH 3/3] Fixing docstring to use ParticleSet.set_dt_dtype() --- parcels/particleset.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/particleset.py b/parcels/particleset.py index 530079a71..0590dc067 100644 --- a/parcels/particleset.py +++ b/parcels/particleset.py @@ -521,7 +521,7 @@ def execute( self._data["dt"][:] = dt else: raise ValueError( - f"The dtype of dt ({dt.dtype}) is coarser than the dtype of the particle dt ({self._data['dt'].dtype}). Please use Particle.set_dt_dtype() to provide a dt with at least the same precision as the particle dt." + f"The dtype of dt ({dt.dtype}) is coarser than the dtype of the particle dt ({self._data['dt'].dtype}). Please use ParticleSet.set_dt_dtype() to provide a dt with at least the same precision as the particle dt." ) if runtime is not None: