Merge branch 'v4-dev' into spatial_slip_interpolation

Author: erikvansebille

docs/v4/index.md

@@ -7,7 +7,7 @@ The key goals of this update are
 1. to support `Fields` on unstructured grids;
 2. to allow for user-defined interpolation methods (somewhat similar to user-defined kernels);
 3. to make the codebase more modular, easier to extend, and more maintainable;
-4. to align Parcels more with other tools in the [Pangeo ecosystemand](https://www.pangeo.io/#ecosystem), particularly by leveraging `xarray` more; and
+4. to align Parcels more with other tools in the [Pangeo ecosystem](https://www.pangeo.io/#ecosystem), particularly by leveraging `xarray` more; and
 5. to improve the performance of Parcels.
 
 The timeline for the release of Parcels v4 is not yet fixed, but we are aiming for a release of an 'alpha' version in September 2025. This v4-alpha will have support for unstructured grids and user-defined interpolation methods, but is not yet performance-optimised.

parcels/_index_search.py

@@ -5,9 +5,9 @@
 
 import numpy as np
 
+from parcels._typing import Mesh
 from parcels.tools.statuscodes import (
-    _raise_field_out_of_bound_error,
-    _raise_field_sampling_error,
+    _raise_grid_searching_error,
     _raise_time_extrapolation_error,
 )
 
@@ -64,12 +64,12 @@ def _search_indices_curvilinear_2d(
     # TODO: Re-enable in some capacity
     # if x < field.lonlat_minmax[0] or x > field.lonlat_minmax[1]:
     #     if grid.lon[0, 0] < grid.lon[0, -1]:
-    #         _raise_field_out_of_bound_error(y, x)
+    #         _raise_grid_searching_error(y, x)
     #     elif x < grid.lon[0, 0] and x > grid.lon[0, -1]:  # This prevents from crashing in [160, -160]
-    #         _raise_field_out_of_bound_error(z, y, x)
+    #         _raise_grid_searching_error(z, y, x)
 
     # if y < field.lonlat_minmax[2] or y > field.lonlat_minmax[3]:
-    #     _raise_field_out_of_bound_error(z, y, x)
+    #     _raise_grid_searching_error(z, y, x)
 
     while np.any(xsi < -tol) or np.any(xsi > 1 + tol) or np.any(eta < -tol) or np.any(eta > 1 + tol):
         px = np.array([grid.lon[yi, xi], grid.lon[yi, xi + 1], grid.lon[yi + 1, xi + 1], grid.lon[yi + 1, xi]])
@@ -99,30 +99,22 @@ def _search_indices_curvilinear_2d(
         it += 1
         if it > maxIterSearch:
             print(f"Correct cell not found after {maxIterSearch} iterations")
-            _raise_field_out_of_bound_error(0, y, x)
+            _raise_grid_searching_error(0, y, x)
     xsi = np.where(xsi < 0.0, 0.0, np.where(xsi > 1.0, 1.0, xsi))
     eta = np.where(eta < 0.0, 0.0, np.where(eta > 1.0, 1.0, eta))
 
     if np.any((xsi < 0) | (xsi > 1) | (eta < 0) | (eta > 1)):
-        _raise_field_sampling_error(y, x)
+        _raise_grid_searching_error(y, x)
     return (yi, eta, xi, xsi)
 
 
-def _reconnect_bnd_indices(yi: int, xi: int, ydim: int, xdim: int, sphere_mesh: bool):
-    if xi < 0:
-        if sphere_mesh:
-            xi = xdim - 2
-        else:
-            xi = 0
-    if xi > xdim - 2:
-        if sphere_mesh:
-            xi = 0
-        else:
-            xi = xdim - 2
-    if yi < 0:
-        yi = 0
-    if yi > ydim - 2:
-        yi = ydim - 2
-        if sphere_mesh:
-            xi = xdim - xi
+def _reconnect_bnd_indices(yi: int, xi: int, ydim: int, xdim: int, mesh: Mesh):
+    xi = np.where(xi < 0, (xdim - 2) if mesh == "spherical" else 0, xi)
+    xi = np.where(xi > xdim - 2, 0 if mesh == "spherical" else (xdim - 2), xi)
+
+    xi = np.where(yi > ydim - 2, xdim - xi if mesh == "spherical" else xi, xi)
+
+    yi = np.where(yi < 0, 0, yi)
+    yi = np.where(yi > ydim - 2, ydim - 2, yi)
+
     return yi, xi

parcels/field.py

@@ -238,10 +238,12 @@ def eval(self, time: datetime, z, y, x, particle=None, applyConversion=True):
 
         tau, ti = _search_time_index(self, time)
         position = self.grid.search(z, y, x, ei=_ei)
-        _update_particle_states(particle, position)
+        _update_particle_states_position(particle, position)
 
         value = self._interp_method(self, ti, position, tau, time, z, y, x)
 
+        _update_particle_states_interp_value(particle, value)
+
         if applyConversion:
             value = self.units.to_target(value, z, y, x)
         return value
@@ -318,16 +320,21 @@ def eval(self, time: datetime, z, y, x, particle=None, applyConversion=True):
 
         tau, ti = _search_time_index(self.U, time)
         position = self.grid.search(z, y, x, ei=_ei)
-        _update_particle_states(particle, position)
+        _update_particle_states_position(particle, position)
 
         if self._vector_interp_method is None:
             u = self.U._interp_method(self.U, ti, position, tau, time, z, y, x)
             v = self.V._interp_method(self.V, ti, position, tau, time, z, y, x)
             if "3D" in self.vector_type:
                 w = self.W._interp_method(self.W, ti, position, tau, time, z, y, x)
+            else:
+                w = 0.0
         else:
             (u, v, w) = self._vector_interp_method(self, ti, position, tau, time, z, y, x)
 
+        for vel in (u, v, w):
+            _update_particle_states_interp_value(particle, vel)
+
         if applyConversion:
             u = self.U.units.to_target(u, z, y, x)
             v = self.V.units.to_target(v, z, y, x)
@@ -349,14 +356,30 @@ def __getitem__(self, key):
             return _deal_with_errors(error, key, vector_type=self.vector_type)
 
 
-def _update_particle_states(particle, position):
+def _update_particle_states_position(particle, position):
     """Update the particle states based on the position dictionary."""
     if particle and "X" in position:  # TODO also support uxgrid search
-        particle.state = np.where(position["X"][0] < 0, StatusCode.ErrorOutOfBounds, particle.state)
-        particle.state = np.where(position["Y"][0] < 0, StatusCode.ErrorOutOfBounds, particle.state)
-        particle.state = np.where(position["Z"][0] == RIGHT_OUT_OF_BOUNDS, StatusCode.ErrorOutOfBounds, particle.state)
-        particle.state = np.where(
-            position["Z"][0] == LEFT_OUT_OF_BOUNDS, StatusCode.ErrorThroughSurface, particle.state
+        particle.state = np.maximum(
+            np.where(position["X"][0] < 0, StatusCode.ErrorOutOfBounds, particle.state), particle.state
+        )
+        particle.state = np.maximum(
+            np.where(position["Y"][0] < 0, StatusCode.ErrorOutOfBounds, particle.state), particle.state
+        )
+        particle.state = np.maximum(
+            np.where(position["Z"][0] == RIGHT_OUT_OF_BOUNDS, StatusCode.ErrorOutOfBounds, particle.state),
+            particle.state,
+        )
+        particle.state = np.maximum(
+            np.where(position["Z"][0] == LEFT_OUT_OF_BOUNDS, StatusCode.ErrorThroughSurface, particle.state),
+            particle.state,
+        )
+
+
+def _update_particle_states_interp_value(particle, value):
+    """Update the particle states based on the interpolated value, but only if state is not an Error already."""
+    if particle:
+        particle.state = np.maximum(
+            np.where(np.isnan(value), StatusCode.ErrorInterpolation, particle.state), particle.state
         )
 
 

parcels/kernel.py

@@ -15,9 +15,11 @@
 from parcels.basegrid import GridType
 from parcels.tools.statuscodes import (
     StatusCode,
+    _raise_field_interpolation_error,
     _raise_field_out_of_bound_error,
     _raise_field_out_of_bound_surface_error,
-    _raise_field_sampling_error,
+    _raise_general_error,
+    _raise_grid_searching_error,
     _raise_time_extrapolation_error,
 )
 from parcels.tools.warnings import KernelWarning
@@ -28,6 +30,16 @@
 __all__ = ["Kernel"]
 
 
+ErrorsToThrow = {
+    StatusCode.ErrorTimeExtrapolation: _raise_time_extrapolation_error,
+    StatusCode.ErrorOutOfBounds: _raise_field_out_of_bound_error,
+    StatusCode.ErrorThroughSurface: _raise_field_out_of_bound_surface_error,
+    StatusCode.ErrorInterpolation: _raise_field_interpolation_error,
+    StatusCode.ErrorGridSearching: _raise_grid_searching_error,
+    StatusCode.Error: _raise_general_error,
+}
+
+
 class Kernel:
     """Kernel object that encapsulates auto-generated code.
 
@@ -270,14 +282,7 @@ def execute(self, pset, endtime, dt):
             if np.any(pset.state == StatusCode.StopAllExecution):
                 return StatusCode.StopAllExecution
 
-            errors_to_throw = {
-                StatusCode.ErrorTimeExtrapolation: _raise_time_extrapolation_error,
-                StatusCode.ErrorOutOfBounds: _raise_field_out_of_bound_error,
-                StatusCode.ErrorThroughSurface: _raise_field_out_of_bound_surface_error,
-                StatusCode.Error: _raise_field_sampling_error,
-            }
-
-            for error_code, error_func in errors_to_throw.items():
+            for error_code, error_func in ErrorsToThrow.items():
                 if np.any(pset.state == error_code):
                     inds = pset.state == error_code
                     if error_code == StatusCode.ErrorTimeExtrapolation:

parcels/particleset.py

@@ -1,6 +1,7 @@
 import sys
 import warnings
 from collections.abc import Iterable
+from typing import Literal
 
 import numpy as np
 import xarray as xr
@@ -506,59 +507,17 @@ def execute(
         if output_file:
             output_file.metadata["parcels_kernels"] = self._kernel.name
 
-        if (dt is not None) and (not isinstance(dt, np.timedelta64)):
-            raise TypeError("dt must be a np.timedelta64 object")
-        if dt is None or np.isnat(dt):
+        if dt is None:
             dt = np.timedelta64(1, "s")
-        self._data["dt"][:] = dt
-        sign_dt = np.sign(dt).astype(int)
-        if sign_dt not in [-1, 1]:
-            raise ValueError("dt must be a positive or negative np.timedelta64 object")
 
-        if self.fieldset.time_interval is None:
-            start_time = np.timedelta64(0, "s")  # For the execution loop, we need a start time as a timedelta object
-            if runtime is None:
-                raise TypeError("The runtime must be provided when the time_interval is not defined for a fieldset.")
+        if not isinstance(dt, np.timedelta64) or np.isnat(dt) or (sign_dt := np.sign(dt).astype(int)) not in [-1, 1]:
+            raise ValueError(f"dt must be a positive or negative np.timedelta64 object, got {dt=!r}")
 
-            else:
-                if isinstance(runtime, np.timedelta64):
-                    end_time = runtime
-                else:
-                    raise TypeError("The runtime must be a np.timedelta64 object")
+        self._data["dt"][:] = dt
 
-        else:
-            if not np.isnat(self.time_nextloop).any():
-                if sign_dt > 0:
-                    start_time = self.time_nextloop.min()
-                else:
-                    start_time = self.time_nextloop.max()
-            else:
-                if sign_dt > 0:
-                    start_time = self.fieldset.time_interval.left
-                else:
-                    start_time = self.fieldset.time_interval.right
-
-            if runtime is None:
-                if endtime is None:
-                    raise ValueError(
-                        "Must provide either runtime or endtime when time_interval is defined for a fieldset."
-                    )
-                # Ensure that the endtime uses the same type as the start_time
-                if isinstance(endtime, self.fieldset.time_interval.left.__class__):
-                    if sign_dt > 0:
-                        if endtime < self.fieldset.time_interval.left:
-                            raise ValueError("The endtime must be after the start time of the fieldset.time_interval")
-                        end_time = min(endtime, self.fieldset.time_interval.right)
-                    else:
-                        if endtime > self.fieldset.time_interval.right:
-                            raise ValueError(
-                                "The endtime must be before the end time of the fieldset.time_interval when dt < 0"
-                            )
-                        end_time = max(endtime, self.fieldset.time_interval.left)
-                else:
-                    raise TypeError("The endtime must be of the same type as the fieldset.time_interval start time.")
-            else:
-                end_time = start_time + runtime * sign_dt
+        start_time, end_time = _get_simulation_start_and_end_times(
+            self.fieldset.time_interval, self._data["time_nextloop"], runtime, endtime, sign_dt
+        )
 
         # Set the time of the particles if it hadn't been set on initialisation
         if np.isnat(self._data["time"]).any():
@@ -619,15 +578,69 @@ def _warn_particle_times_outside_fieldset_time_bounds(release_times: np.ndarray,
 
     if isinstance(time.left, np.datetime64) and isinstance(release_times[0], np.timedelta64):
         release_times = np.array([t + time.left for t in release_times])
-    if np.any(release_times < time.left):
+    if np.any(release_times < time.left) or np.any(release_times > time.right):
         warnings.warn(
             "Some particles are set to be released outside the FieldSet's executable time domain.",
             ParticleSetWarning,
             stacklevel=2,
         )
-    if np.any(release_times > time.right):
-        warnings.warn(
-            "Some particles are set to be released after the fieldset's last time and the fields are not constant in time.",
-            ParticleSetWarning,
-            stacklevel=2,
+
+
+def _get_simulation_start_and_end_times(
+    time_interval: TimeInterval,
+    particle_release_times: np.ndarray,
+    runtime: np.timedelta64 | None,
+    endtime: np.datetime64 | None,
+    sign_dt: Literal[-1, 1],
+) -> tuple[np.datetime64, np.datetime64]:
+    if runtime is not None and endtime is not None:
+        raise ValueError(
+            f"runtime and endtime are mutually exclusive - provide one or the other. Got {runtime=!r}, {endtime=!r}"
         )
+
+    if runtime is None and time_interval is None:
+        raise ValueError("The runtime must be provided when the time_interval is not defined for a fieldset.")
+
+    if sign_dt == 1:
+        first_release_time = particle_release_times.min()
+    else:
+        first_release_time = particle_release_times.max()
+
+    start_time = _get_start_time(first_release_time, time_interval, sign_dt, runtime)
+
+    if endtime is None:
+        if not isinstance(runtime, np.timedelta64):
+            raise ValueError(f"The runtime must be a np.timedelta64 object. Got {type(runtime)}")
+
+        endtime = start_time + sign_dt * runtime
+
+    if time_interval is not None:
+        if type(endtime) != type(time_interval.left):  # noqa: E721
+            raise ValueError(
+                f"The endtime must be of the same type as the fieldset.time_interval start time. Got {endtime=!r} with {time_interval=!r}"
+            )
+        if endtime not in time_interval:
+            msg = (
+                f"Calculated/provided end time of {endtime!r} is not in fieldset time interval {time_interval!r}. Either reduce your runtime, modify your "
+                "provided endtime, or change your release timing."
+                "Important info:\n"
+                f"    First particle release: {first_release_time!r}\n"
+                f"    runtime: {runtime!r}\n"
+                f"    (calculated) endtime: {endtime!r}"
+            )
+            raise ValueError(msg)
+
+    return start_time, endtime
+
+
+def _get_start_time(first_release_time, time_interval, sign_dt, runtime):
+    if time_interval is None:
+        time_interval = TimeInterval(left=np.timedelta64(0, "s"), right=runtime)
+
+    if sign_dt == 1:
+        fieldset_start = time_interval.left
+    else:
+        fieldset_start = time_interval.right
+
+    start_time = first_release_time if not np.isnat(first_release_time) else fieldset_start
+    return start_time