Implement support for preventing beaching (#160)

* Implement support for no_beaching in FVCOM

* Add entry on beaching to docs cookbook

doc/source/documentation/cookbook/beaching.ipynb

@@ -0,0 +1,55 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "a4abeb9c-4943-4b50-8b16-3a6a824b0d2f",
+   "metadata": {},
+   "source": [
+    "# Beaching\n",
+    "\n",
+    "In the ocean, objects that are transported by ocean currents may become beached or stranded on land.  If an object becomes beached in an intertidal zone, the beaching event may be short-lived, with the object being resuspended on the next flood tide. In contrast, if the object is deposited on the strand line during a spring tide, it may be several days or weeks before it has an opportunity to be resuspended. In addition, beaching will also be influenced by waves, the wind and storm surges. Objects may even become permanently or semi-permanently beached if deposited far inland during a storm surge. Beaching and resuspension dynamics are also dependent on the type of object under consideration.\n",
+    "\n",
+    "Particle tracking models handle beaching in different ways. The approach taken is strongly influenced by the type of driving data that is being used. Driving data from a global or broad scale ocean circulation model will not, in general, include wetting and drying with changes in the tide - the spatial resolution is too poor - meaning grid cells remain wet at all times. However, beaching may still be parameterised; for example, by imposing a probability that a particle becomes beached if it resides within a given distance of the model's land boundary, or if it crosses a land boundary during the simulation. Resuspension may be handled in a similar way. Alternatively, beaching may be ignored completely. In this scenario, if a simulated particle does cross a land boundary - which can happen for various reasons - it may be restored to its previous position or reflected back into the model domain. When using the output of fine scale models of coastal or estuarine systems which include wetting and drying, one is also forced to consider what happens to particles that become trapped in intertidal zones.\n",
+    "\n",
+    "## Beaching in PyLag\n",
+    "\n",
+    "Presently, PyLag does not include a paramaterisation of beaching in the absence of wetting and drying; if a land boundary is crossed, particles are either restored to their previous position or reflected back into the model domain depending on the type of boundary condition calculator used. However, it does include a parameterisation for wetting and drying. There are two options which are controlled through the binary flag `allow_beaching`, which should be set in the `[SIMULATION]` section of the PyLag run configuration file.\n",
+    "\n",
+    "If `allow_beaching=True`, then particles can become beached in intertidal zones. If this occurs, they are flagged as having beached using the particle attribute `is_beached`. Particles that are beached are passed over until the host cell they are in is flooded with water once more. If `allow_beaching=False`, particles are prevented from beaching by moving them to the centroid of the nearest wet element.\n",
+    "\n",
+    "Although in reality intertidal zones are free of water, this is not the case in models. Typically, models will impose a threshold water depth at which a cell is deemed \"dry\". If a particle is in an element when this threshold is met, it will become beached. The situation is complicated slightly when doing offline particle tracking simulations using discrete model outputs which are (necessarily!) saved infrequently relative to the model time step. In this event, the tidal cycle may be poorly resolved - or even aliased and the period for which the particle remains beached may become artificially long.\n",
+    "\n",
+    "When beaching is excluded, and the run is fully deterministic (i.e., diffusion is neglected), care should also be taken when particles beach in the same location at the same time. If this occurs, the selected particles will be moved to the same point in space - the centroid of the nearest wet element - and continue to move as one."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "id": "ecbed4b4-e69b-4f38-ba7e-06a3aa138fef",
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.9"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}

doc/source/documentation/cookbook/index.rst

@@ -11,6 +11,7 @@ In the cookbook, you can find hints and tips for working with PyLag; explanation
 
     initial_positions
     boundary_conditions
+    beaching
     fixed_vertical_positioning
     boundary_layer_dynamics
     settling

pylag/data_reader.pxd

@@ -33,6 +33,8 @@ cdef class DataReader:
                                    DTYPE_FLOAT_t particle_pathline[2],
                                    DTYPE_FLOAT_t intersection[2])
 
+    cdef find_nearest_wet_host(self, DTYPE_FLOAT_t time, Particle *particle)
+
     cdef set_default_location(self, Particle *particle)
 
     cdef set_local_coordinates(self, Particle *particle)

pylag/data_reader.pyx

@@ -178,6 +178,29 @@ cdef class DataReader:
                                    DTYPE_FLOAT_t intersection[2]):
         raise NotImplementedError
 
+    def find_nearest_wet_host_wrapper(self, DTYPE_FLOAT_t time,
+                                       ParticleSmartPtr particle):
+        """ Python wrapper for finding the nearest wet host element
+
+        This function is used to find the nearest wet host element to the
+        particle's current position.
+
+        Parameters
+        ----------
+        time : float
+            The current time.
+        
+        particle : pylag.particle_cpp_wrapper.ParticleSmartPtr
+            The particle at its current position.
+        
+        Returns
+        -------
+        """
+        return self.find_nearest_wet_host(time, particle.get_ptr())
+
+    cdef find_nearest_wet_host(self, DTYPE_FLOAT_t time, Particle *particle):
+        raise NotImplementedError
+
     def set_default_location_wrapper(self, ParticleSmartPtr particle):
         """ Python wrapper for setting the default location of a particle
 

pylag/fvcom_data_reader.pyx

@@ -28,7 +28,7 @@ from pylag.data_types_cython cimport DTYPE_INT_t, DTYPE_FLOAT_t
 
 from libcpp.string cimport string
 from libcpp.vector cimport vector
-from libc.math cimport log
+from libc.math cimport log, sqrt
 
 # PyLag cython imports
 from pylag.parameters cimport cartesian, geographic, deg_to_radians
@@ -463,6 +463,98 @@ cdef class FVCOMDataReader(DataReader):
 
         return
 
+    cdef find_nearest_wet_host(self, DTYPE_FLOAT_t time, Particle *particle):
+        """ Find the nearest wet host element
+
+        Find the nearest wet host element to the particle's position and relocate
+        the particle to the centroid of this element. Begin by searching the 
+        the particle's host element and then its immediate neighbours. If no wet
+        elements are found, the neighbours of neighbours are searched and so on.
+        If multiple wet elements are found, the one closest to the particle is
+        selected.
+        """
+        cdef vector[DTYPE_INT_t] current_nbe_checks
+        cdef vector[DTYPE_INT_t] next_nbe_checks
+        cdef vector[DTYPE_INT_t] elems_checked
+        cdef vector[DTYPE_INT_t] wet_elems
+
+        cdef DTYPE_FLOAT_t dist, min_dist
+        cdef DTYPE_INT_t current_host, nbe, n_wet_elems, new_host
+        cdef DTYPE_INT_t i, j 
+        cdef bint has_been_checked
+
+        # Confirm the particle currently resides in a dry element
+        if self.is_wet(time, particle) == 1:
+            return
+
+        # Save the current host
+        current_host = particle.get_host_horizontal_elem(self._name)
+
+        # Add the current host to the list of elements whose neighbours we will check next
+        next_nbe_checks.push_back(current_host)
+
+        # Flag that we have already checked the element itself
+        elems_checked.push_back(current_host)
+
+        # Perform an ever widening search to find the nearest wet element
+        while True:
+            current_nbe_checks = next_nbe_checks
+
+            # Clear the list of elements to check next - we will repopulate this as we go
+            next_nbe_checks.clear()
+
+            # Loop over the current list of elements whose neighbours we will check
+            for elem in current_nbe_checks:
+
+                # Search the element's immediate neighbours
+                for i in range(3):
+                    nbe = self._nbe[i, elem]
+
+                    # Check if sea point
+                    if nbe == LAND_BDY_CROSSED or nbe == OPEN_BDY_CROSSED:
+                        continue
+
+                    # Have we cheked this element before?
+                    has_been_checked = False
+                    for j in range(elems_checked.size()):
+                        if nbe == elems_checked[j]:
+                            has_been_checked = True
+                            break
+
+                    if has_been_checked:
+                        continue
+                    else:
+                        # Check if the neighbour is wet
+                        particle.set_host_horizontal_elem(self._name, nbe)
+
+                        if self.is_wet(time, particle) == 1:
+                            wet_elems.push_back(nbe)
+
+                        next_nbe_checks.push_back(nbe)
+                        elems_checked.push_back(nbe)
+
+                n_wet_elems = wet_elems.size()
+                if n_wet_elems == 0:
+                    continue
+                else:
+                    if n_wet_elems == 1:
+                        new_host = wet_elems[0]
+                    else:
+                        # If multiple wet elements are found, select the one closest to the particle
+                        min_dist = 1.e20
+                        for i in range(wet_elems.size()):
+                            # TODO - Will work with geographic coordinates?
+                            dist = sqrt((particle.get_x1() - self._xc[wet_elems[i]])**2 +
+                                        (particle.get_x2() - self._yc[wet_elems[i]])**2)
+                            if dist < min_dist:
+                                min_dist = dist
+                                new_host = wet_elems[i]
+
+                    particle.set_host_horizontal_elem(self._name, new_host)
+                    self.set_default_location(particle)
+
+                    return
+
     cdef set_local_coordinates(self, Particle *particle):
         """ Set local coordinates