Implementing View for KernelParticle/ParticleSet

docs/user_guide/examples/tutorial_Argofloats.ipynb

@@ -26,70 +26,62 @@
    "source": [
     "import numpy as np\n",
     "\n",
-    "# Define the new Kernels that mimic Argo vertical movement\n",
+    "# Define the new Kernel that mimics Argo vertical movement\n",
     "driftdepth = 1000  # maximum depth in m\n",
     "maxdepth = 2000  # maximum depth in m\n",
     "vertical_speed = 0.10  # sink and rise speed in m/s\n",
     "cycletime = 10 * 86400  # total time of cycle in seconds\n",
     "drifttime = 9 * 86400  # time of deep drift in seconds\n",
     "\n",
     "\n",
-    "def ArgoPhase1(particles, fieldset):\n",
-    "    def SinkingPhase(p):\n",
-    "        \"\"\"Phase 0: Sinking with vertical_speed until depth is driftdepth\"\"\"\n",
-    "        p.dz += vertical_speed * particles.dt\n",
-    "        p.cycle_phase = np.where(p.z + p.dz >= driftdepth, 1, p.cycle_phase)\n",
-    "        p.dz = np.where(p.z + p.dz >= driftdepth, driftdepth - p.z, p.dz)\n",
+    "def ArgoVerticalMovement(particles, fieldset):\n",
+    "    # Split particles based on their current cycle_phase\n",
+    "    ptcls0 = particles[particles.cycle_phase == 0]\n",
+    "    ptcls1 = particles[particles.cycle_phase == 1]\n",
+    "    ptcls2 = particles[particles.cycle_phase == 2]\n",
+    "    ptcls3 = particles[particles.cycle_phase == 3]\n",
+    "    ptcls4 = particles[particles.cycle_phase == 4]\n",
+    "\n",
+    "    # Phase 0: Sinking with vertical_speed until depth is driftdepth\n",
+    "    ptcls0.dz += vertical_speed * ptcls0.dt\n",
+    "    ptcls0.cycle_phase = np.where(\n",
+    "        ptcls0.z + ptcls0.dz >= driftdepth, 1, ptcls0.cycle_phase\n",
+    "    )\n",
+    "    ptcls0.dz = np.where(\n",
+    "        ptcls0.z + ptcls0.dz >= driftdepth, driftdepth - ptcls0.z, ptcls0.dz\n",
+    "    )\n",
+    "\n",
+    "    # Phase 1: Drifting at depth for drifttime seconds\n",
+    "    ptcls1.drift_age += ptcls1.dt\n",
+    "    ptcls1.cycle_phase = np.where(ptcls1.drift_age >= drifttime, 2, ptcls1.cycle_phase)\n",
+    "    ptcls1.drift_age = np.where(ptcls1.drift_age >= drifttime, 0, ptcls1.drift_age)\n",
+    "\n",
+    "    # Phase 2: Sinking further to maxdepth\n",
+    "    ptcls2.dz += vertical_speed * ptcls2.dt\n",
+    "    ptcls2.cycle_phase = np.where(\n",
+    "        ptcls2.z + ptcls2.dz >= maxdepth, 3, ptcls2.cycle_phase\n",
+    "    )\n",
+    "    ptcls2.dz = np.where(\n",
+    "        ptcls2.z + ptcls2.dz >= maxdepth, maxdepth - ptcls2.z, ptcls2.dz\n",
+    "    )\n",
+    "\n",
+    "    # Phase 3: Rising with vertical_speed until at surface\n",
+    "    ptcls3.dz -= vertical_speed * ptcls3.dt\n",
+    "    ptcls3.temp = fieldset.thetao[ptcls3.time, ptcls3.z, ptcls3.lat, ptcls3.lon]\n",
+    "    ptcls3.cycle_phase = np.where(\n",
+    "        ptcls3.z + ptcls3.dz <= fieldset.mindepth, 4, ptcls3.cycle_phase\n",
+    "    )\n",
+    "    ptcls3.dz = np.where(\n",
+    "        ptcls3.z + ptcls3.dz <= fieldset.mindepth,\n",
+    "        fieldset.mindepth - ptcls3.z,\n",
+    "        ptcls3.dz,\n",
+    "    )\n",
+    "\n",
+    "    # Phase 4: Transmitting at surface until cycletime is reached\n",
+    "    ptcls4.cycle_phase = np.where(ptcls4.cycle_age >= cycletime, 0, ptcls4.cycle_phase)\n",
+    "    ptcls4.cycle_age = np.where(ptcls4.cycle_age >= cycletime, 0, ptcls4.cycle_age)\n",
+    "    ptcls4.temp = np.nan  # no temperature measurement when at surface\n",
     "\n",
-    "    SinkingPhase(particles[particles.cycle_phase == 0])\n",
-    "\n",
-    "\n",
-    "def ArgoPhase2(particles, fieldset):\n",
-    "    def DriftingPhase(p):\n",
-    "        \"\"\"Phase 1: Drifting at depth for drifttime seconds\"\"\"\n",
-    "        p.drift_age += particles.dt\n",
-    "        p.cycle_phase = np.where(p.drift_age >= drifttime, 2, p.cycle_phase)\n",
-    "        p.drift_age = np.where(p.drift_age >= drifttime, 0, p.drift_age)\n",
-    "\n",
-    "    DriftingPhase(particles[particles.cycle_phase == 1])\n",
-    "\n",
-    "\n",
-    "def ArgoPhase3(particles, fieldset):\n",
-    "    def SecondSinkingPhase(p):\n",
-    "        \"\"\"Phase 2: Sinking further to maxdepth\"\"\"\n",
-    "        p.dz += vertical_speed * particles.dt\n",
-    "        p.cycle_phase = np.where(p.z + p.dz >= maxdepth, 3, p.cycle_phase)\n",
-    "        p.dz = np.where(p.z + p.dz >= maxdepth, maxdepth - p.z, p.dz)\n",
-    "\n",
-    "    SecondSinkingPhase(particles[particles.cycle_phase == 2])\n",
-    "\n",
-    "\n",
-    "def ArgoPhase4(particles, fieldset):\n",
-    "    def RisingPhase(p):\n",
-    "        \"\"\"Phase 3: Rising with vertical_speed until at surface\"\"\"\n",
-    "        p.dz -= vertical_speed * particles.dt\n",
-    "        p.temp = fieldset.thetao[p.time, p.z, p.lat, p.lon]\n",
-    "        p.cycle_phase = np.where(p.z + p.dz <= fieldset.mindepth, 4, p.cycle_phase)\n",
-    "        p.dz = np.where(\n",
-    "            p.z + p.dz <= fieldset.mindepth,\n",
-    "            fieldset.mindepth - p.z,\n",
-    "            p.dz,\n",
-    "        )\n",
-    "\n",
-    "    RisingPhase(particles[particles.cycle_phase == 3])\n",
-    "\n",
-    "\n",
-    "def ArgoPhase5(particles, fieldset):\n",
-    "    def TransmittingPhase(p):\n",
-    "        \"\"\"Phase 4: Transmitting at surface until cycletime is reached\"\"\"\n",
-    "        p.cycle_phase = np.where(p.cycle_age >= cycletime, 0, p.cycle_phase)\n",
-    "        p.cycle_age = np.where(p.cycle_age >= cycletime, 0, p.cycle_age)\n",
-    "        p.temp = np.nan  # no temperature measurement when at surface\n",
-    "\n",
-    "    TransmittingPhase(particles[particles.cycle_phase == 4])\n",
-    "\n",
-    "\n",
-    "def ArgoPhase6(particles, fieldset):\n",
     "    particles.cycle_age += particles.dt  # update cycle_age"
    ]
   },
@@ -136,9 +128,7 @@
     "ArgoParticle = parcels.Particle.add_variable(\n",
     "    [\n",
     "        parcels.Variable(\"cycle_phase\", dtype=np.int32, initial=0.0),\n",
-    "        parcels.Variable(\n",
-    "            \"cycle_age\", dtype=np.float32, initial=0.0\n",
-    "        ),  # TODO update to \"timedelta64[s]\"\n",
+    "        parcels.Variable(\"cycle_age\", dtype=np.float32, initial=0.0),\n",
     "        parcels.Variable(\"drift_age\", dtype=np.float32, initial=0.0),\n",
     "        parcels.Variable(\"temp\", dtype=np.float32, initial=np.nan),\n",
     "    ]\n",
@@ -155,12 +145,7 @@
     "\n",
     "# combine Argo vertical movement kernel with built-in Advection kernel\n",
     "kernels = [\n",
-    "    ArgoPhase1,\n",
-    "    ArgoPhase2,\n",
-    "    ArgoPhase3,\n",
-    "    ArgoPhase4,\n",
-    "    ArgoPhase5,\n",
-    "    ArgoPhase6,\n",
+    "    ArgoVerticalMovement,\n",
     "    parcels.kernels.AdvectionRK4,\n",
     "]\n",
     "\n",

docs/user_guide/examples/tutorial_interaction.ipynb

@@ -293,18 +293,9 @@
     "    larger_idx = np.where(mass_j > mass_i, pair_j, pair_i)\n",
     "    smaller_idx = np.where(mass_j > mass_i, pair_i, pair_j)\n",
     "\n",
-    "    # perform transfer and mark deletions\n",
-    "    # TODO note that we use temporary arrays for indexing because of KernelParticle bug (GH #2143)\n",
-    "    masses = particles.mass\n",
-    "    states = particles.state\n",
-    "\n",
     "    # transfer mass from smaller to larger and mark smaller for deletion\n",
-    "    masses[larger_idx] += particles.mass[smaller_idx]\n",
-    "    states[smaller_idx] = parcels.StatusCode.Delete\n",
-    "\n",
-    "    # TODO use particle variables directly after KernelParticle bug (GH #2143) is fixed\n",
-    "    particles.mass = masses\n",
-    "    particles.state = states"
+    "    particles.mass[larger_idx] += particles.mass[smaller_idx]\n",
+    "    particles.state[smaller_idx] = parcels.StatusCode.Delete"
    ]
   },
   {

src/parcels/__init__.py

@@ -17,7 +17,7 @@
     Variable,
     Particle,
     ParticleClass,
-    KernelParticle,  # ? remove?
+    ParticleSetView,
 )
 from parcels._core.field import Field, VectorField
 from parcels._core.basegrid import BaseGrid
@@ -87,8 +87,8 @@
     "logger",
     "download_example_dataset",
     "list_example_datasets",
-    # (marked for potential removal)
-    "KernelParticle",
+    # Helpers for particle handling
+    "ParticleSetView",
 ]
 
 _stdlib_warnings.warn(

src/parcels/_core/field.py

@@ -14,7 +14,7 @@
     _unitconverters_map,
 )
 from parcels._core.index_search import GRID_SEARCH_ERROR, LEFT_OUT_OF_BOUNDS, RIGHT_OUT_OF_BOUNDS, _search_time_index
-from parcels._core.particle import KernelParticle
+from parcels._core.particle import ParticleSetView
 from parcels._core.statuscodes import (
     AllParcelsErrorCodes,
     StatusCode,
@@ -35,9 +35,9 @@
 
 
 def _deal_with_errors(error, key, vector_type: VectorType):
-    if isinstance(key, KernelParticle):
+    if isinstance(key, ParticleSetView):
         key.state = AllParcelsErrorCodes[type(error)]
-    elif isinstance(key[-1], KernelParticle):
+    elif isinstance(key[-1], ParticleSetView):
         key[-1].state = AllParcelsErrorCodes[type(error)]
     else:
         raise RuntimeError(f"{error}. Error could not be handled because particles was not part of the Field Sampling.")
@@ -229,7 +229,7 @@ def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True):
     def __getitem__(self, key):
         self._check_velocitysampling()
         try:
-            if isinstance(key, KernelParticle):
+            if isinstance(key, ParticleSetView):
                 return self.eval(key.time, key.z, key.lat, key.lon, key)
             else:
                 return self.eval(*key)
@@ -330,7 +330,7 @@ def eval(self, time: datetime, z, y, x, particles=None, applyConversion=True):
 
     def __getitem__(self, key):
         try:
-            if isinstance(key, KernelParticle):
+            if isinstance(key, ParticleSetView):
                 return self.eval(key.time, key.z, key.lat, key.lon, key)
             else:
                 return self.eval(*key)

src/parcels/_core/index_search.py

@@ -239,18 +239,20 @@ def uxgrid_point_in_cell(grid, y: np.ndarray, x: np.ndarray, yi: np.ndarray, xi:
             axis=-1,
         )
 
-        # Get projection points onto element plane
-        # for the projection, all points are computed relative to v0
-        r1 = np.squeeze(face_vertices[:, 1, :] - face_vertices[:, 0, :])  # (M,3)
-        r2 = np.squeeze(face_vertices[:, 2, :] - face_vertices[:, 0, :])  # (M,3)
+        # Get projection points onto element plane. Keep the leading
+        # dimension even for single-face queries so shapes remain (M,3).
+        r1 = face_vertices[:, 1, :] - face_vertices[:, 0, :]
+        r2 = face_vertices[:, 2, :] - face_vertices[:, 0, :]
         nhat = np.cross(r1, r2)
         norm = np.linalg.norm(nhat, axis=-1)
+        # Avoid division by zero for degenerate faces
+        norm = np.where(norm == 0.0, 1.0, norm)
         nhat = nhat / norm[:, None]
         # Calculate the component of the points in the direction of nhat
-        ptilde = points - np.squeeze(face_vertices[:, 0, :])
+        ptilde = points - face_vertices[:, 0, :]
         pdotnhat = np.sum(ptilde * nhat, axis=-1)
         # Reconstruct points with normal component removed.
-        points = ptilde - pdotnhat[:, None] * nhat + np.squeeze(face_vertices[:, 0, :])
+        points = ptilde - pdotnhat[:, None] * nhat + face_vertices[:, 0, :]
 
     else:
         nids = grid.uxgrid.face_node_connectivity[xi].values