Merge branch 'v4-dev' into implement_ei_in_search

Author: erikvansebille

docs/examples/tutorial_kernelloop.ipynb

@@ -26,17 +26,7 @@
     "\n",
     "When you run a Parcels simulation (i.e. a call to `pset.execute()`), the Kernel loop is the main part of the code that is executed. This part of the code loops through all particles and executes the Kernels that are defined for each particle.\n",
     "\n",
-    "In order to make sure that the displacements of a particle in the different Kernels can be summed, all Kernels add to a _change_ in position (`particle_dlon`, `particle_dlat`, and `particle_ddepth`). This is important, because there are situations where movement kernels would otherwise not commute. Take the example of advecting particles by currents _and_ winds. If the particle would first be moved by the currents and then by the winds, the result could be different from first moving by the winds and then by the currents. Instead, by adding the changes in position, the ordering of the Kernels has no consequence on the particle displacement."
-   ]
-  },
-  {
-   "cell_type": "markdown",
-   "metadata": {},
-   "source": [
-    "<div class=\"alert alert-info\">\n",
-    "\n",
-    "__Note__ that the variables `particle_dlon`, `particle_dlat`, and `particle_ddepth` are defined by the wrapper function that Parcels generates for each Kernel. This is why you don't have to define these variables yourself when writing a Kernel. See [here](https://github.com/OceanParcels/parcels/blob/daa4b062ed8ae0b2be3d87367d6b45599d6f95db/parcels/kernel.py#L277-L294) for the implementation of the wrapper functions.\n",
-    "</div>"
+    "In order to make sure that the displacements of a particle in the different Kernels can be summed, all Kernels add to a _change_ in position (`particles.dlon`, `particles.dlat`, and `particles.ddepth`). This is important, because there are situations where movement kernels would otherwise not commute. Take the example of advecting particles by currents _and_ winds. If the particle would first be moved by the currents and then by the winds, the result could be different from first moving by the winds and then by the currents. Instead, by adding the changes in position, the ordering of the Kernels has no consequence on the particle displacement."
    ]
   },
   {
@@ -52,29 +42,25 @@
    "source": [
     "Below is a structured overview of the Kernel loop is implemented. Note that this is for longitude only, but the same process is applied for latitude and depth.\n",
     "\n",
-    "1. Define an extra variable `particle.lon_nextloop` for each particle, which is the longitude at the end of the Kernel loop. Inititalise it to `particle.lon`.\n",
-    "\n",
-    "2. Also define an extra variable `particle.time_nextloop` for each particle, which is the time at the end of the Kernel loop. Inititalise it to `particle.time`.\n",
+    "1. Initialise an extra Variable `particles.lon=0` and `particles.time_nextloop = particles.time`\n",
     "\n",
-    "3. Within the Kernel loop, for each particle:<br>\n",
+    "2. Within the Kernel loop, for each particle:<br>\n",
     "\n",
-    "    1. Update `particle.lon` with `particle.lon_nextloop`<br>\n",
+    "    1. Update `particles.lon += particles.dlon`<br>\n",
     "\n",
-    "    2. Update `particle.time` with `particle.time_nextloop`<br>\n",
+    "    2. Set variable `particles.dlon = 0`<br>\n",
     "\n",
-    "    3. Set local variable `particle_dlon = 0`<br>\n",
+    "    3. Update `particles.time = particles.time_nextloop`\n",
     "\n",
     "    4. For each Kernel in the list of Kernels:\n",
     "        \n",
     "        1. Execute the Kernel\n",
     "        \n",
-    "        2. Update `particle_dlon` by adding the change in longitude, if needed<br>\n",
+    "        2. Update `particles.dlon` by adding the change in longitude, if needed<br>\n",
     "\n",
-    "    5. Update `particle.lon_nextloop` with `particle.lon + particle_dlon`<br>\n",
-    "    \n",
-    "    6. Update `particle.time_nextloop` with `particle.time + particle.dt`<br>\n",
+    "    5. Update `particles.time_nextloop += particles.dt`<br>\n",
     "\n",
-    "    7. If `outputdt` is a multiple of `particle.time`, write `particle.lon` and `particle.time` to zarr output file<br>\n",
+    "    6. If `outputdt` is a multiple of `particle.time`, write `particle.lon` and `particle.time` to zarr output file<br>\n",
     "\n",
     "Besides having commutable Kernels, the main advantage of this implementation is that, when using Field Sampling with e.g. `particle.temp = fieldset.Temp[particle.time, particle.depth, particle.lat, particle.lon]`, the particle location stays the same throughout the entire Kernel loop. Additionally, this implementation ensures that the particle location is the same as the location of the sampled field in the output file."
    ]
@@ -268,12 +254,12 @@
     "### 1. Avoid updating particle locations directly in Kernels\n",
     "It is better not to update `particle.lon` directly in a Kernel, as it can interfere with the loop above. Assigning a value to `particle.lon` in a Kernel will throw a warning. \n",
     "\n",
-    "Instead, update the local variable `particle_dlon`.\n",
+    "Instead, update the local variable `particle.dlon`.\n",
     "\n",
     "### 2. Be careful with updating particle variables that do not depend on Fields.\n",
     "While assigning the interpolated value of a `Field` to a Particle goes well in the loop above, this is not necessarily so for assigning other attributes. For example, a line like `particle.age += particle.dt` is executed directly so may result in the age being `dt` at `time = 0` in the output file. \n",
     "\n",
-    "A workaround is to either initialise the age to `-dt`, or to increase the `age` only when `particle.time > 0` (using an `if` statement).\n",
+    "A workaround is to either initialise the age to `-dt`, or to increase the `age` only when `particle.time > 0` (using an `np.where` statement).\n",
     "\n",
     "\n",
     "### 3. The last time is not written to file\n",

parcels/basegrid.py

@@ -54,9 +54,9 @@ def search(self, z: float, y: float, x: float, ei=None) -> dict[str, tuple[int,
             - Unstructured grid: {"Z": (zi, zeta), "FACE": (fi, bcoords)}
 
             Where:
-            - index (int): The cell position of a particle along the given axis
+            - index (int): The cell position of the particles along the given axis
             - barycentric_coordinates (float or np.ndarray): The coordinates defining
-              a particle's position within the grid cell. For structured grids, this
+              the particles positions within the grid cell. For structured grids, this
               is a single coordinate per axis; for unstructured grids, this can be
               an array of coordinates for the face polygon.
 

parcels/field.py

@@ -41,7 +41,7 @@ def _deal_with_errors(error, key, vector_type: VectorType):
     elif isinstance(key[-1], KernelParticle):
         key[-1].state = AllParcelsErrorCodes[type(error)]
     else:
-        raise RuntimeError(f"{error}. Error could not be handled because particle was not part of the Field Sampling.")
+        raise RuntimeError(f"{error}. Error could not be handled because particles was not part of the Field Sampling.")
 
     if vector_type and "3D" in vector_type:
         return (0, 0, 0)
@@ -205,26 +205,26 @@ def _check_velocitysampling(self):
                 stacklevel=2,
             )
 
-    def eval(self, time: datetime, z, y, x, particle=None, applyConversion=True):
+    def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True):
         """Interpolate field values in space and time.
 
         We interpolate linearly in time and apply implicit unit
         conversion to the result. Note that we defer to
         scipy.interpolate to perform spatial interpolation.
         """
-        if particle is None:
+        if particles is None:
             _ei = None
         else:
-            _ei = particle.ei[:, self.igrid]
+            _ei = particles.ei[:, self.igrid]
 
         tau, ti = _search_time_index(self, time)
         position = self.grid.search(z, y, x, ei=_ei)
-        _update_particles_ei(particle, position, self)
-        _update_particle_states_position(particle, position)
+        _update_particles_ei(particles, position, self)
+        _update_particle_states_position(particles, position)
 
         value = self._interp_method(self, ti, position, tau, time, z, y, x)
 
-        _update_particle_states_interp_value(particle, value)
+        _update_particle_states_interp_value(particles, value)
 
         if applyConversion:
             value = self.units.to_target(value, z, y, x)
@@ -289,22 +289,22 @@ def vector_interp_method(self, method: Callable):
         _assert_same_function_signature(method, ref=ZeroInterpolator_Vector, context="Interpolation")
         self._vector_interp_method = method
 
-    def eval(self, time: datetime, z, y, x, particle=None, applyConversion=True):
+    def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True):
         """Interpolate field values in space and time.
 
         We interpolate linearly in time and apply implicit unit
         conversion to the result. Note that we defer to
         scipy.interpolate to perform spatial interpolation.
         """
-        if particle is None:
+        if particles is None:
             _ei = None
         else:
-            _ei = particle.ei[:, self.igrid]
+            _ei = particles.ei[:, self.igrid]
 
         tau, ti = _search_time_index(self.U, time)
         position = self.grid.search(z, y, x, ei=_ei)
-        _update_particles_ei(particle, position, self)
-        _update_particle_states_position(particle, position)
+        _update_particles_ei(particles, position, self)
+        _update_particle_states_position(particles, position)
 
         if self._vector_interp_method is None:
             u = self.U._interp_method(self.U, ti, position, tau, time, z, y, x)
@@ -322,7 +322,7 @@ def eval(self, time: datetime, z, y, x, particle=None, applyConversion=True):
             (u, v, w) = self._vector_interp_method(self, ti, position, tau, time, z, y, x, applyConversion)
 
         for vel in (u, v, w):
-            _update_particle_states_interp_value(particle, vel)
+            _update_particle_states_interp_value(particles, vel)
 
         if applyConversion and ("3D" in self.vector_type):
             w = self.W.units.to_target(w, z, y, x) if self.W else 0.0
@@ -362,34 +362,34 @@ def _update_particles_ei(particles, position, field):
             )
 
 
-def _update_particle_states_position(particle, position):
+def _update_particle_states_position(particles, position):
     """Update the particle states based on the position dictionary."""
-    if particle:  # TODO also support uxgrid search
+    if particles:  # TODO also support uxgrid search
         for dim in ["X", "Y"]:
             if dim in position:
-                particle.state = np.maximum(
-                    np.where(position[dim][0] == -1, StatusCode.ErrorOutOfBounds, particle.state), particle.state
+                particles.state = np.maximum(
+                    np.where(position[dim][0] == -1, StatusCode.ErrorOutOfBounds, particles.state), particles.state
                 )
-                particle.state = np.maximum(
-                    np.where(position[dim][0] == GRID_SEARCH_ERROR, StatusCode.ErrorGridSearching, particle.state),
-                    particle.state,
+                particles.state = np.maximum(
+                    np.where(position[dim][0] == GRID_SEARCH_ERROR, StatusCode.ErrorGridSearching, particles.state),
+                    particles.state,
                 )
         if "Z" in position:
-            particle.state = np.maximum(
-                np.where(position["Z"][0] == RIGHT_OUT_OF_BOUNDS, StatusCode.ErrorOutOfBounds, particle.state),
-                particle.state,
+            particles.state = np.maximum(
+                np.where(position["Z"][0] == RIGHT_OUT_OF_BOUNDS, StatusCode.ErrorOutOfBounds, particles.state),
+                particles.state,
             )
-            particle.state = np.maximum(
-                np.where(position["Z"][0] == LEFT_OUT_OF_BOUNDS, StatusCode.ErrorThroughSurface, particle.state),
-                particle.state,
+            particles.state = np.maximum(
+                np.where(position["Z"][0] == LEFT_OUT_OF_BOUNDS, StatusCode.ErrorThroughSurface, particles.state),
+                particles.state,
             )
 
 
-def _update_particle_states_interp_value(particle, value):
+def _update_particle_states_interp_value(particles, 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
+    if particles:
+        particles.state = np.maximum(
+            np.where(np.isnan(value), StatusCode.ErrorInterpolation, particles.state), particles.state
         )
 
 

parcels/kernel.py

@@ -10,6 +10,7 @@
 
 from parcels.application_kernels.advection import (
     AdvectionAnalytical,
+    AdvectionRK4,
     AdvectionRK45,
 )
 from parcels.basegrid import GridType
@@ -24,7 +25,6 @@
     _raise_time_extrapolation_error,
 )
 from parcels.tools.warnings import KernelWarning
-from tests.common_kernels import DoNothing
 
 if TYPE_CHECKING:
     from collections.abc import Callable
@@ -69,7 +69,7 @@ def __init__(
         for f in pyfuncs:
             if not isinstance(f, types.FunctionType):
                 raise TypeError(f"Argument pyfunc should be a function or list of functions. Got {type(f)}")
-            _assert_same_function_signature(f, ref=DoNothing, context="Kernel")
+            _assert_same_function_signature(f, ref=AdvectionRK4, context="Kernel")
 
         if len(pyfuncs) == 0:
             raise ValueError("List of `pyfuncs` should have at least one function.")
@@ -118,21 +118,20 @@ def add_positionupdate_kernels(self):
         def Setcoords(particles, fieldset):  # pragma: no cover
             import numpy as np  # noqa
 
+            particles.lon += particles.dlon
+            particles.lat += particles.dlat
+            particles.depth += particles.ddepth
+
             particles.dlon = 0
             particles.dlat = 0
             particles.ddepth = 0
-            particles.lon = particles.lon_nextloop
-            particles.lat = particles.lat_nextloop
-            particles.depth = particles.depth_nextloop
+
             particles.time = particles.time_nextloop
 
-        def Updatecoords(particles, fieldset):  # pragma: no cover
-            particles.lon_nextloop = particles.lon + particles.dlon
-            particles.lat_nextloop = particles.lat + particles.dlat
-            particles.depth_nextloop = particles.depth + particles.ddepth
+        def UpdateTime(particles, fieldset):  # pragma: no cover
             particles.time_nextloop = particles.time + particles.dt
 
-        self._pyfuncs = (Setcoords + self + Updatecoords)._pyfuncs
+        self._pyfuncs = (Setcoords + self + UpdateTime)._pyfuncs
 
     def check_fieldsets_in_kernels(self, pyfunc):  # TODO v4: this can go into another method? assert_is_compatible()?
         """

parcels/particle.py

@@ -170,13 +170,11 @@ def get_default_particle(spatial_dtype: np.float32 | np.float64) -> ParticleClas
                 dtype=spatial_dtype,
                 attrs={"standard_name": "longitude", "units": "degrees_east", "axis": "X"},
             ),
-            Variable("lon_nextloop", dtype=spatial_dtype, to_write=False),
             Variable(
                 "lat",
                 dtype=spatial_dtype,
                 attrs={"standard_name": "latitude", "units": "degrees_north", "axis": "Y"},
             ),
-            Variable("lat_nextloop", dtype=spatial_dtype, to_write=False),
             Variable(
                 "depth",
                 dtype=spatial_dtype,
@@ -185,7 +183,6 @@ def get_default_particle(spatial_dtype: np.float32 | np.float64) -> ParticleClas
             Variable("dlon", dtype=spatial_dtype, to_write=False),
             Variable("dlat", dtype=spatial_dtype, to_write=False),
             Variable("ddepth", dtype=spatial_dtype, to_write=False),
-            Variable("depth_nextloop", dtype=spatial_dtype, to_write=False),
             Variable(
                 "time",
                 dtype=_SAME_AS_FIELDSET_TIME_INTERVAL.VALUE,

parcels/particlefile.py

@@ -262,9 +262,9 @@ def write_latest_locations(self, pset, time):
         time :
             Time at which to write ParticleSet. Note that typically this would be pset.time_nextloop
         """
-        for var in ["lon", "lat", "depth", "time"]:
-            pset._data[f"{var}"] = pset._data[f"{var}_nextloop"]
-
+        for var in ["lon", "lat", "depth"]:
+            pset._data[f"{var}"] += pset._data[f"d{var}"]
+        pset._data["time"] = pset._data["time_nextloop"]
         self.write(pset, time)
 
 

parcels/particleset.py

@@ -22,12 +22,12 @@
 
 
 class ParticleSet:
-    """Class for storing particle and executing kernel over them.
+    """Class for storing particles and executing kernel over them.
 
     Please note that this currently only supports fixed size particle sets, meaning that the particle set only
     holds the particles defined on construction. Individual particles can neither be added nor deleted individually,
-    and individual particles can only be deleted as a set procedurally (i.e. by 'particle.delete()'-call during
-    kernel execution).
+    and individual particles can only be deleted as a set procedurally (i.e. by changing their state to 'StatusCode.Delete'
+    during kernel execution).
 
     Parameters
     ----------
@@ -128,9 +128,6 @@ def __init__(
                 lat=lat,
                 depth=depth,
                 time=time,
-                lon_nextloop=lon,
-                lat_nextloop=lat,
-                depth_nextloop=depth,
                 time_nextloop=time,
                 trajectory=trajectory_ids,
             ),

parcels/tools/statuscodes.py

@@ -110,10 +110,10 @@ def _raise_time_extrapolation_error(time: float, field=None):
 
 
 class KernelError(RuntimeError):
-    """General particle kernel error with optional custom message."""
+    """General particles kernel error with optional custom message."""
 
-    def __init__(self, particle, fieldset=None, msg=None):
-        message = f"{particle.state}\nParticle {particle}\nTime: {particle.time}\ntimestep dt: {particle.dt}\n"
+    def __init__(self, particles, fieldset=None, msg=None):
+        message = f"{particles.state}\nParticle {particles}\nTime: {particles.time}\ntimestep dt: {particles.dt}\n"
         if msg:
             message += msg
         super().__init__(message)

parcels/xgrid.py

@@ -508,7 +508,7 @@ def _search_1d_array(
     x: float,
 ) -> tuple[int, int]:
     """
-    Searches for the particle location in a 1D array and returns barycentric coordinate along dimension.
+    Searches for particle locations in a 1D array and returns barycentric coordinate along dimension.
 
     Assumptions:
     - array is strictly monotonically increasing.