Supporting particle.delete in Kernel

Author: erikvansebille

parcels/kernel.py

@@ -72,8 +72,9 @@ def remove_deleted(self, pset):
         """Utility to remove all particles that signalled deletion."""
         bool_indices = pset._data["state"] == StatusCode.Delete
         indices = np.where(bool_indices)[0]
-        if len(indices) > 0 and self.fieldset.particlefile is not None:
-            self.fieldset.particlefile.write(pset, None, indices=indices)
+        # TODO v4: need to implement ParticleFile writing of deleted particles
+        # if len(indices) > 0 and self.fieldset.particlefile is not None:
+        #     self.fieldset.particlefile.write(pset, None, indices=indices)
         pset.remove_indices(indices)
 
 

parcels/particle.py

@@ -2,6 +2,8 @@
 
 import numpy as np
 
+from parcels.tools.statuscodes import StatusCode
+
 __all__ = ["InteractionParticle", "Particle", "Variable"]
 
 
@@ -136,6 +138,10 @@ def __repr__(self):
                 p_string += f"{var}={getattr(self, var):f}, "
         return p_string + f"time={time_string})"
 
+    def delete(self):
+        """Signal the particle for deletion."""
+        self.state = StatusCode.Delete
+
     @classmethod
     def add_variable(cls, var, *args, **kwargs):
         """Add a new variable to the Particle class

tests/v4/test_particleset.py

@@ -218,6 +218,27 @@ def test_pset_merge_inplace(fieldset, npart=100):
     assert pset1.size == 2 * npart
 
 
+def test_pset_remove_index(fieldset, npart=100):
+    lon = np.linspace(0, 1, npart)
+    lat = np.linspace(1, 0, npart)
+    pset = ParticleSet(fieldset, lon=lon, lat=lat, lonlatdepth_dtype=np.float64)
+    indices_to_remove = [0, 10, 20]
+    pset.remove_indices(indices_to_remove)
+    assert pset.size == 97
+    assert not np.any(np.in1d(pset.trajectory, indices_to_remove))
+
+
+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))
+
+    def DeleteKernel(particle, fieldset, time):  # pragma: no cover
+        if particle.lon >= 0.4:
+            particle.delete()
+
+    pset.execute(pset.Kernel(DeleteKernel), runtime=np.timedelta64(1, "s"), dt=np.timedelta64(1, "s"))
+    assert pset.size == 40
+
+
 def test_pset_iterator(fieldset):
     npart = 10
     pset = ParticleSet(fieldset, lon=np.zeros(npart), lat=np.ones(npart))