Merge pull request #2426 from Parcels-code/particle_particle_interaction_v4

Making particle-particle interaction tutorial v4-compatible

Author: erikvansebille

docs/user_guide/examples/tutorial_interaction.ipynb

@@ -54,6 +54,7 @@ examples/tutorial_statuscodes.ipynb
 examples/tutorial_gsw_density.ipynb
 examples/tutorial_Argofloats.ipynb
 examples/tutorial_diffusion.ipynb
+examples/tutorial_interaction.ipynb
 ```
 
 ```{toctree}
@@ -65,7 +66,6 @@ examples/tutorial_interpolation.ipynb
 ```
 
 <!-- examples/tutorial_particle_field_interaction.ipynb -->
-<!-- examples/tutorial_interaction.ipynb -->
 <!-- examples/tutorial_analyticaladvection.ipynb -->
 <!-- examples/tutorial_kernelloop.ipynb -->
 

docs/user_guide/examples_v3/tutorial_interaction.ipynb

@@ -14,6 +14,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat
 - The `particle` argument in the Kernel signature has been renamed to `particles`.
 - `math` functions should be replaced with array compatible equivalents (e.g., `math.sin` -> `np.sin`). Instead of `ParcelsRandom` you should use numpy's random functions.
 - `particle.depth` has been changed to `particles.z` to be consistent with the [CF conventions for trajectory data](https://cfconventions.org/cf-conventions/cf-conventions.html#trajectory-data), and to make Parcels also generalizable to atmospheric contexts.
+- The `InteractionKernel` class has been removed. Since normal Kernels now have access to _all_ particles, particle-particle interaction can be performed within normal Kernels.
 
 ## FieldSet
 

docs/user_guide/index.md

@@ -53,9 +53,6 @@ class ParticleSet:
         Optional list of "trajectory" values (integers) for the particle IDs
     partition_function :
         Function to use for partitioning particles over processors. Default is to use kMeans
-    periodic_domain_zonal :
-        Zonal domain size, used to apply zonally periodic boundaries for particle-particle
-        interaction. If None, no zonally periodic boundaries are applied
 
         Other Variables can be initialised using further arguments (e.g. v=... for a Variable named 'v')
     """
@@ -74,7 +71,6 @@ def __init__(
         self._data = None
         self._repeat_starttime = None
         self._kernel = None
-        self._interaction_kernel = None
 
         self.fieldset = fieldset
         lon = np.empty(shape=0) if lon is None else _convert_to_flat_array(lon)
@@ -228,8 +224,6 @@ def add(self, particles):
         for d in self._data:
             self._data[d] = np.concatenate((self._data[d], particles._data[d]))
 
-        # Adding particles invalidates the neighbor search structure.
-        self._dirty_neighbor = True
         return self
 
     def __iadd__(self, particles):
@@ -256,44 +250,6 @@ def remove_indices(self, indices):
         for d in self._data:
             self._data[d] = np.delete(self._data[d], indices, axis=0)
 
-    def _active_particles_mask(self, time, dt):
-        active_indices = (time - self._data["time"]) / dt >= 0
-        non_err_indices = np.isin(self._data["state"], [StatusCode.Success, StatusCode.Evaluate])
-        active_indices = np.logical_and(active_indices, non_err_indices)
-        self._active_particle_idx = np.where(active_indices)[0]
-        return active_indices
-
-    def _compute_neighbor_tree(self, time, dt):
-        active_mask = self._active_particles_mask(time, dt)
-
-        self._values = np.vstack(
-            (
-                self._data["z"],
-                self._data["lat"],
-                self._data["lon"],
-            )
-        )
-        if self._dirty_neighbor:
-            self._neighbor_tree.rebuild(self._values, active_mask=active_mask)
-            self._dirty_neighbor = False
-        else:
-            self._neighbor_tree.update_values(self._values, new_active_mask=active_mask)
-
-    def _neighbors_by_index(self, particle_idx):
-        neighbor_idx, distances = self._neighbor_tree.find_neighbors_by_idx(particle_idx)
-        neighbor_idx = self._active_particle_idx[neighbor_idx]
-        mask = neighbor_idx != particle_idx
-        neighbor_idx = neighbor_idx[mask]
-        if "horiz_dist" in self._data._ptype.variables:
-            self._data["vert_dist"][neighbor_idx] = distances[0, mask]
-            self._data["horiz_dist"][neighbor_idx] = distances[1, mask]
-        return True  # TODO fix for v4 ParticleDataIterator(self.particledata, subset=neighbor_idx)
-
-    def _neighbors_by_coor(self, coor):
-        neighbor_idx = self._neighbor_tree.find_neighbors_by_coor(coor)
-        neighbor_ids = self._data["trajectory"][neighbor_idx]
-        return neighbor_ids
-
     def populate_indices(self):
         """Pre-populate guesses of particle ei (element id) indices"""
         for i, grid in enumerate(self.fieldset.gridset):
@@ -359,13 +315,6 @@ def Kernel(self, pyfunc):
             pyfuncs=[pyfunc],
         )
 
-    def InteractionKernel(self, pyfunc_inter):
-        from parcels.interaction.interactionkernel import InteractionKernel
-
-        if pyfunc_inter is None:
-            return None
-        return InteractionKernel(self.fieldset, self._ptype, pyfunc=pyfunc_inter)
-
     def data_indices(self, variable_name, compare_values, invert=False):
         """Get the indices of all particles where the value of `variable_name` equals (one of) `compare_values`.