Repository: OceanParcels/Parcels (python)

Issue:
Warning when particles are initialised outside time domain
@sruehs was encountering a bug with her simulation where the number of the particles in the output weren't the same as in the particle set. Turns out not all particles in the particleset were initialised within the time domain of `pset.execute(runtime=...)`. Having an informative warning here would be helpful. Note that the warning would be raised if `pset.execute()` is used in a loop, but that can now be silenced with a filter from `warnings`.

cc @michaeldenes 

Pull request description:
Implement warning for particles initialised outside time domain
This PR adds a warning when some particles are set to be released outside the `FieldSet`'s time bounds

It also adds a new `ParticleSet` warning to the list of Parcels custom warnings

This fixes #1694

Patch:
diff --git a/parcels/particleset.py b/parcels/particleset.py
index 6f7b5299..a5e058f1 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -14,7 +14,7 @@ from tqdm import tqdm
 from parcels._compat import MPI
 from parcels.application_kernels.advection import AdvectionRK4
 from parcels.compilation.codecompiler import GNUCompiler
-from parcels.field import NestedField
+from parcels.field import Field, NestedField
 from parcels.grid import CurvilinearGrid, GridType
 from parcels.interaction.interactionkernel import InteractionKernel
 from parcels.interaction.neighborsearch import (
@@ -32,7 +32,7 @@ from parcels.tools.converters import _get_cftime_calendars, convert_to_flat_arra
 from parcels.tools.global_statics import get_package_dir
 from parcels.tools.loggers import logger
 from parcels.tools.statuscodes import StatusCode
-from parcels.tools.warnings import FileWarning
+from parcels.tools.warnings import ParticleSetWarning
 
 __all__ = ["ParticleSet"]
 
@@ -174,6 +174,8 @@ class ParticleSet:
             raise NotImplementedError("If fieldset.time_origin is not a date, time of a particle must be a double")
         time = np.array([self.time_origin.reltime(t) if _convert_to_reltime(t) else t for t in time])
         assert lon.size == time.size, "time and positions (lon, lat, depth) do not have the same lengths."
+        if isinstance(fieldset.U, Field) and (not fieldset.U.allow_time_extrapolation):
+            _warn_particle_times_outside_fieldset_time_bounds(time, fieldset.U.grid.time_full)
 
         if lonlatdepth_dtype is None:
             lonlatdepth_dtype = self.lonlatdepth_dtype_from_field_interp_method(fieldset.U)
@@ -792,7 +794,7 @@ class ParticleSet:
                 f"Note that the `repeatdt` argument is not retained from {filename}, and that "
                 "setting a new repeatdt will start particles from the _new_ particle "
                 "locations.",
-                FileWarning,
+                ParticleSetWarning,
                 stacklevel=2,
             )
 
@@ -1247,6 +1249,22 @@ def _warn_outputdt_release_desync(outputdt: float, starttime: float, release_tim
             "Some of the particles have a start time difference that is not a multiple of outputdt. "
             "This could cause the first output of some of the particles that start later "
             "in the simulation to be at a different time than expected.",
-            FileWarning,
+            ParticleSetWarning,
             stacklevel=2,
         )
+
+
+def _warn_particle_times_outside_fieldset_time_bounds(release_times: np.ndarray, time_full: np.ndarray):
+    if np.any(release_times):
+        if np.any(release_times < time_full[0]):
+            warnings.warn(
+                "Some particles are set to be released before the fieldset's first time and allow_time_extrapolation is set to False.",
+                ParticleSetWarning,
+                stacklevel=2,
+            )
+        if np.any(release_times > time_full[-1]):
+            warnings.warn(
+                "Some particles are set to be released after the fieldset's last time and allow_time_extrapolation is set to False.",
+                ParticleSetWarning,
+                stacklevel=2,
+            )
diff --git a/parcels/tools/warnings.py b/parcels/tools/warnings.py
index 2f3cfae7..7dbcbbbc 100644
--- a/parcels/tools/warnings.py
+++ b/parcels/tools/warnings.py
@@ -1,6 +1,6 @@
 import warnings
 
-__all__ = ["FieldSetWarning", "FileWarning", "KernelWarning"]
+__all__ = ["FieldSetWarning", "FileWarning", "KernelWarning", "ParticleSetWarning"]
 
 
 class FieldSetWarning(UserWarning):
@@ -13,6 +13,12 @@ class FieldSetWarning(UserWarning):
     pass
 
 
+class ParticleSetWarning(UserWarning):
+    """Warning that is raised when there are issues in the construction of the ParticleSet."""
+
+    pass
+
+
 class FileWarning(UserWarning):
     """Warning that is raised when there are issues with input or output files.
 


Test patch:
diff --git a/tests/test_particlesets.py b/tests/test_particlesets.py
index 21906118..0678df36 100644
--- a/tests/test_particlesets.py
+++ b/tests/test_particlesets.py
@@ -7,6 +7,7 @@ from parcels import (
     FieldSet,
     JITParticle,
     ParticleSet,
+    ParticleSetWarning,
     ScipyParticle,
     StatusCode,
     Variable,
@@ -175,6 +176,14 @@ def test_pset_create_with_time(fieldset, mode):
     assert np.allclose([p.time for p in pset], time, rtol=1e-12)
 
 
+@pytest.mark.parametrize("mode", ["scipy", "jit"])
+def test_pset_create_outside_time(mode):
+    fieldset = create_fieldset_zeros_simple(withtime=True)
+    time = [-1, 0, 1, 20 * 86400]
+    with pytest.warns(ParticleSetWarning, match="Some particles are set to be released*"):
+        ParticleSet(fieldset, pclass=ptype[mode], lon=[0] * len(time), lat=[0] * len(time), time=time)
+
+
 @pytest.mark.parametrize("mode", ["scipy", "jit"])
 def test_pset_not_multipldt_time(fieldset, mode):
     times = [0, 1.1]
diff --git a/tests/tools/test_warnings.py b/tests/tools/test_warnings.py
index dc3588fd..1979f9c8 100644
--- a/tests/tools/test_warnings.py
+++ b/tests/tools/test_warnings.py
@@ -9,9 +9,9 @@ from parcels import (
     AdvectionRK45,
     FieldSet,
     FieldSetWarning,
-    FileWarning,
     KernelWarning,
     ParticleSet,
+    ParticleSetWarning,
     ScipyParticle,
 )
 from tests.utils import TEST_DATA
@@ -62,7 +62,7 @@ def test_file_warnings(tmp_zarrfile):
     )
     pset = ParticleSet(fieldset=fieldset, pclass=ScipyParticle, lon=[0, 0], lat=[0, 0], time=[0, 1])
     pfile = pset.ParticleFile(name=tmp_zarrfile, outputdt=2)
-    with pytest.warns(FileWarning, match="Some of the particles have a start time difference.*"):
+    with pytest.warns(ParticleSetWarning, match="Some of the particles have a start time difference.*"):
         pset.execute(AdvectionRK4, runtime=3, dt=1, output_file=pfile)
 
 
diff --git a/tests/utils.py b/tests/utils.py
index a1f55e99..c20d95cc 100644
--- a/tests/utils.py
+++ b/tests/utils.py
@@ -72,15 +72,24 @@ def create_flat_positions(n_particle):
     return np.random.rand(n_particle * 3).reshape(3, n_particle)
 
 
-def create_fieldset_zeros_simple(xdim=40, ydim=100):
-    U = np.zeros((ydim, xdim), dtype=np.float32)
-    V = np.zeros((ydim, xdim), dtype=np.float32)
+def create_fieldset_zeros_simple(xdim=40, ydim=100, withtime=False):
     lon = np.linspace(0, 1, xdim, dtype=np.float32)
     lat = np.linspace(-60, 60, ydim, dtype=np.float32)
     depth = np.zeros(1, dtype=np.float32)
-    data = {"U": np.array(U, dtype=np.float32), "V": np.array(V, dtype=np.float32)}
     dimensions = {"lat": lat, "lon": lon, "depth": depth}
-    return FieldSet.from_data(data, dimensions)
+    if withtime:
+        tdim = 10
+        time = np.linspace(0, 86400, tdim, dtype=np.float64)
+        dimensions["time"] = time
+        datadims = (tdim, ydim, xdim)
+        allow_time_extrapolation = False
+    else:
+        datadims = (ydim, xdim)
+        allow_time_extrapolation = True
+    U = np.zeros(datadims, dtype=np.float32)
+    V = np.zeros(datadims, dtype=np.float32)
+    data = {"U": np.array(U, dtype=np.float32), "V": np.array(V, dtype=np.float32)}
+    return FieldSet.from_data(data, dimensions, allow_time_extrapolation=allow_time_extrapolation)
 
 
 def assert_empty_folder(path: Path):