Repository: OceanParcels/parcels (python)

Issue:
Error in Scipy mode when OutOfBounds in depth dimension.
It seems that in 3D Scipy mode, a particle moving outside of the grid in the depth dimension is not handled correctly, throwing a generic KernelError rather than passing to the recovery kernel as expected.

In the minimal working example below, a 10cell x 10cell x 10cell grid is created, followed by a velocity field with constant flow in the lon and depth directions (and zero flow in the lat direction). I also define a recovery kernel called `DeleteParticle` which simply deletes a particle from its particleset when called.

When a particle `p` crosses the lon limit (i.e. when `p.lon` > 10), it is correctly passed to the recovery kernel and deleted. But when a particle crosses depth limit (`p.depth` > 10), Parcels throws a generic KernelError, rather than calling the recovery kernel.

This issue is not present in JIT mode, where the second particle is deleted as expected.

```
#!/usr/bin/env
from parcels import *

# define recovery kernel
def DeleteParticle(particle, fieldset, time):  # delete particles who run out of bounds.
    print("Particle %d deleted at (%f, %f, %f)" % (particle.id, particle.lon, particle.lat, particle.depth))
    particle.delete()

# create grid
lon = np.linspace(0, 10, 11, dtype=np.float32)
lat = np.linspace(0, 10, 11, dtype=np.float32)
depth = np.linspace(0, 10, 11, dtype=np.float32)
grid = RectilinearZGrid(lon=lon, lat=lat, depth=depth, time=np.zeros(1), mesh='flat')

# create fields
U = Field(name='U', data=np.ones([len(depth), len(lat), len(lon)]), grid=grid,
          allow_time_extrapolation=True)
V = Field(name='V', data=np.zeros([len(depth), len(lat), len(lon)]), grid=grid,
          allow_time_extrapolation=True)
W = Field(name='W', data=np.ones([len(depth), len(lat), len(lon)]), grid=grid,
          allow_time_extrapolation=True)
fieldset = FieldSet(U=U, V=V, fields={'W': W})

# show expected OutOfBounds behaviour when particle crosses lon limit.
particleset_lon = ParticleSet.from_list(fieldset=fieldset, pclass=ScipyParticle,
                                    lon=[9.9],
		                    lat=[5],
		                    depth=[5])
particleset_lon.execute(AdvectionRK4_3D, runtime=1., dt=.5,
                        recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

# show erroneous OutOfBounds behaviour when particle crosses depth limit.
particleset_depth = ParticleSet.from_list(fieldset=fieldset, pclass=ScipyParticle,
                                          lon=[5],
		                          lat=[5],
		                          depth=[9.9])
particleset_depth.execute(AdvectionRK4_3D, runtime=1., dt=.25,
                          recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})
```



Pull request description:
Changing interpolating sanity check to throw FieldSamplingError
Also adding unit test (and removing old ‘main’ in test_advection.py to bring in line with all other unit tests)

This fixes #572

Patch:
diff --git a/parcels/field.py b/parcels/field.py
index e77485a9..d2a82481 100644
--- a/parcels/field.py
+++ b/parcels/field.py
@@ -1,7 +1,7 @@
 from parcels.tools.loggers import logger
 from parcels.tools.converters import unitconverters_map, UnitConverter, Geographic, GeographicPolar
 from parcels.tools.converters import TimeConverter
-from parcels.tools.error import FieldSamplingError, TimeExtrapolationError
+from parcels.tools.error import FieldSamplingError, FieldOutOfBoundError, TimeExtrapolationError
 import parcels.tools.interpolation_utils as i_u
 import collections
 from py import path
@@ -428,6 +428,8 @@ class Field(object):
     def search_indices_vertical_z(self, z):
         grid = self.grid
         z = np.float32(z)
+        if z < grid.depth[0] or z > grid.depth[-1]:
+            raise FieldOutOfBoundError(0, 0, z, field=self)
         depth_index = grid.depth <= z
         if z >= grid.depth[-1]:
             zi = len(grid.depth) - 2
@@ -469,7 +471,7 @@ class Field(object):
         else:
             zi = depth_index.argmin() - 1 if z >= depth_vector[0] else 0
         if z < depth_vector[zi] or z > depth_vector[zi+1]:
-            raise FieldSamplingError(x, y, z, field=self)
+            raise FieldOutOfBoundError(x, y, z, field=self)
         zeta = (z - depth_vector[zi]) / (depth_vector[zi+1]-depth_vector[zi])
         return (zi, zeta)
 
@@ -498,9 +500,9 @@ class Field(object):
 
         if not grid.zonal_periodic:
             if x < grid.lonlat_minmax[0] or x > grid.lonlat_minmax[1]:
-                raise FieldSamplingError(x, y, z, field=self)
+                raise FieldOutOfBoundError(x, y, z, field=self)
         if y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]:
-            raise FieldSamplingError(x, y, z, field=self)
+            raise FieldOutOfBoundError(x, y, z, field=self)
 
         if grid.mesh != 'spherical':
             lon_index = grid.lon < x
@@ -553,16 +555,18 @@ class Field(object):
         if grid.zdim > 1 and not search2D:
             if grid.gtype == GridCode.RectilinearZGrid:
                 # Never passes here, because in this case, we work with scipy
-                (zi, zeta) = self.search_indices_vertical_z(z)
+                try:
+                    (zi, zeta) = self.search_indices_vertical_z(z)
+                except FieldOutOfBoundError:
+                    raise FieldOutOfBoundError(x, y, z, field=self)
             elif grid.gtype == GridCode.RectilinearSGrid:
                 (zi, zeta) = self.search_indices_vertical_s(x, y, z, xi, yi, xsi, eta, ti, time)
         else:
             zi = -1
             zeta = 0
 
-        assert(xsi >= 0 and xsi <= 1)
-        assert(eta >= 0 and eta <= 1)
-        assert(zeta >= 0 and zeta <= 1)
+        if not ((0 <= xsi <= 1) and (0 <= eta <= 1) and (0 <= zeta <= 1)):
+            raise FieldSamplingError(x, y, z, field=self)
 
         return (xsi, eta, zeta, xi, yi, zi)
 
@@ -578,11 +582,11 @@ class Field(object):
         if not grid.zonal_periodic:
             if x < grid.lonlat_minmax[0] or x > grid.lonlat_minmax[1]:
                 if grid.lon[0, 0] < grid.lon[0, -1]:
-                    raise FieldSamplingError(x, y, z, field=self)
+                    raise FieldOutOfBoundError(x, y, z, field=self)
                 elif x < grid.lon[0, 0] and x > grid.lon[0, -1]:  # This prevents from crashing in [160, -160]
-                    raise FieldSamplingError(x, y, z, field=self)
+                    raise FieldOutOfBoundError(x, y, z, field=self)
         if y < grid.lonlat_minmax[2] or y > grid.lonlat_minmax[3]:
-            raise FieldSamplingError(x, y, z, field=self)
+            raise FieldOutOfBoundError(x, y, z, field=self)
 
         while xsi < 0 or xsi > 1 or eta < 0 or eta > 1:
             px = np.array([grid.lon[yi, xi], grid.lon[yi, xi+1], grid.lon[yi+1, xi+1], grid.lon[yi+1, xi]])
@@ -610,9 +614,9 @@ class Field(object):
             else:
                 xsi = (x-a[0]-a[2]*eta) / (a[1]+a[3]*eta)
             if xsi < 0 and eta < 0 and xi == 0 and yi == 0:
-                raise FieldSamplingError(x, y, 0, field=self)
+                raise FieldOutOfBoundError(x, y, 0, field=self)
             if xsi > 1 and eta > 1 and xi == grid.xdim-1 and yi == grid.ydim-1:
-                raise FieldSamplingError(x, y, 0, field=self)
+                raise FieldOutOfBoundError(x, y, 0, field=self)
             if xsi < 0:
                 xi -= 1
             elif xsi > 1:
@@ -625,20 +629,22 @@ class Field(object):
             it += 1
             if it > maxIterSearch:
                 print('Correct cell not found after %d iterations' % maxIterSearch)
-                raise FieldSamplingError(x, y, 0, field=self)
+                raise FieldOutOfBoundError(x, y, 0, field=self)
 
         if grid.zdim > 1 and not search2D:
             if grid.gtype == GridCode.CurvilinearZGrid:
-                (zi, zeta) = self.search_indices_vertical_z(z)
+                try:
+                    (zi, zeta) = self.search_indices_vertical_z(z)
+                except FieldOutOfBoundError:
+                    raise FieldOutOfBoundError(x, y, z, field=self)
             elif grid.gtype == GridCode.CurvilinearSGrid:
                 (zi, zeta) = self.search_indices_vertical_s(x, y, z, xi, yi, xsi, eta, ti, time)
         else:
             zi = -1
             zeta = 0
 
-        assert(xsi >= 0 and xsi <= 1)
-        assert(eta >= 0 and eta <= 1)
-        assert(zeta >= 0 and zeta <= 1)
+        if not ((0 <= xsi <= 1) and (0 <= eta <= 1) and (0 <= zeta <= 1)):
+            raise FieldSamplingError(x, y, z, field=self)
 
         return (xsi, eta, zeta, xi, yi, zi)
 
@@ -733,7 +739,7 @@ class Field(object):
             val = self.interpolator3D(ti, z, y, x, time)
         if np.isnan(val):
             # Detect Out-of-bounds sampling and raise exception
-            raise FieldSamplingError(x, y, z, field=self)
+            raise FieldOutOfBoundError(x, y, z, field=self)
         else:
             return val
 
@@ -768,7 +774,7 @@ class Field(object):
     def depth_index(self, depth, lat, lon):
         """Find the index in the depth array associated with a given depth"""
         if depth > self.grid.depth[-1]:
-            raise FieldSamplingError(lon, lat, depth, field=self)
+            raise FieldOutOfBoundError(lon, lat, depth, field=self)
         depth_index = self.grid.depth <= depth
         if depth_index.all():
             # If given depth == largest field depth, use the second-last
@@ -1346,7 +1352,7 @@ class NestedField(list):
                 try:
                     val = list.__getitem__(self, iField).eval(*key)
                     break
-                except FieldSamplingError:
+                except (FieldOutOfBoundError, FieldSamplingError):
                     if iField == len(self)-1:
                         raise
                     else:
diff --git a/parcels/include/index_search.h b/parcels/include/index_search.h
index 948bba04..ebc53634 100644
--- a/parcels/include/index_search.h
+++ b/parcels/include/index_search.h
@@ -221,9 +221,9 @@ static inline ErrorCode search_indices_rectilinear(type_coord x, type_coord y, t
   else
     *zeta = 0;
 
-  if ( (*xsi < 0) || (*xsi > 1) ) return ERROR_OUT_OF_BOUNDS;
-  if ( (*eta < 0) || (*eta > 1) ) return ERROR_OUT_OF_BOUNDS;
-  if ( (*zeta < 0) || (*zeta > 1) ) return ERROR_OUT_OF_BOUNDS;
+  if ( (*xsi < 0) || (*xsi > 1) ) return ERROR;
+  if ( (*eta < 0) || (*eta > 1) ) return ERROR;
+  if ( (*zeta < 0) || (*zeta > 1) ) return ERROR;
 
   return SUCCESS;
 }
@@ -357,9 +357,9 @@ static inline ErrorCode search_indices_curvilinear(type_coord x, type_coord y, t
   else
     *zeta = 0;
 
-  if ( (*xsi < 0) || (*xsi > 1) ) return ERROR_OUT_OF_BOUNDS;
-  if ( (*eta < 0) || (*eta > 1) ) return ERROR_OUT_OF_BOUNDS;
-  if ( (*zeta < 0) || (*zeta > 1) ) return ERROR_OUT_OF_BOUNDS;
+  if ( (*xsi < 0) || (*xsi > 1) ) return ERROR;
+  if ( (*eta < 0) || (*eta > 1) ) return ERROR;
+  if ( (*zeta < 0) || (*zeta > 1) ) return ERROR;
 
   return SUCCESS;
 }
diff --git a/parcels/kernel.py b/parcels/kernel.py
index 1f4729e8..8a3131d9 100644
--- a/parcels/kernel.py
+++ b/parcels/kernel.py
@@ -1,7 +1,7 @@
 from parcels.codegenerator import KernelGenerator, LoopGenerator
 from parcels.compiler import get_cache_dir
 from parcels.tools.error import ErrorCode, recovery_map as recovery_base_map
-from parcels.field import FieldSamplingError, Field, SummedField, NestedField
+from parcels.field import FieldOutOfBoundError, Field, SummedField, NestedField
 from parcels.tools.loggers import logger
 from parcels.kernels.advection import AdvectionRK4_3D
 from os import path, remove
@@ -226,7 +226,7 @@ class Kernel(object):
                 try:
                     p.dt = sign_dt * dt_pos
                     res = self.pyfunc(p, pset.fieldset, p.time)
-                except FieldSamplingError as fse:
+                except FieldOutOfBoundError as fse:
                     res = ErrorCode.ErrorOutOfBounds
                     p.exception = fse
                 except Exception as e:
diff --git a/parcels/tools/error.py b/parcels/tools/error.py
index 12544ab4..8542b1ee 100644
--- a/parcels/tools/error.py
+++ b/parcels/tools/error.py
@@ -2,7 +2,7 @@
 from enum import IntEnum
 
 
-__all__ = ['ErrorCode', 'FieldSamplingError', 'TimeExtrapolationError',
+__all__ = ['ErrorCode', 'FieldSamplingError', 'FieldOutOfBoundError', 'TimeExtrapolationError',
            'KernelError', 'OutOfBoundsError', 'recovery_map']
 
 
@@ -29,6 +29,20 @@ class FieldSamplingError(RuntimeError):
         super(FieldSamplingError, self).__init__(message)
 
 
+class FieldOutOfBoundError(RuntimeError):
+    """Utility error class to propagate out-of-bound field sampling in Scipy mode"""
+
+    def __init__(self, x, y, z, field=None):
+        self.field = field
+        self.x = x
+        self.y = y
+        self.z = z
+        message = "%s sampled out-of-bound, at (%f, %f, %f)" % (
+            field.name if field else "Field", self.x, self.y, self.z
+        )
+        super(FieldOutOfBoundError, self).__init__(message)
+
+
 class TimeExtrapolationError(RuntimeError):
     """Utility error class to propagate erroneous time extrapolation sampling in Scipy mode"""
 


Test patch:
diff --git a/tests/test_advection.py b/tests/test_advection.py
index 216350f3..0eaf884d 100644
--- a/tests/test_advection.py
+++ b/tests/test_advection.py
@@ -1,10 +1,9 @@
-from parcels import FieldSet, ParticleSet, ScipyParticle, JITParticle
+from parcels import FieldSet, ParticleSet, ScipyParticle, JITParticle, ErrorCode
 from parcels import AdvectionEE, AdvectionRK4, AdvectionRK45, AdvectionRK4_3D
 import numpy as np
 import pytest
 import math
 from datetime import timedelta as delta
-from argparse import ArgumentParser
 
 
 ptype = {'scipy': ScipyParticle, 'jit': JITParticle}
@@ -102,6 +101,25 @@ def test_advection_3D(mode, npart=11):
     assert np.allclose([p.depth*time for p in pset], [p.lon for p in pset], atol=1.e-1)
 
 
+@pytest.mark.parametrize('mode', ['jit', 'scipy'])
+def test_advection_3D_outofbounds(mode):
+    xdim = ydim = zdim = 2
+    dimensions = {'lon': np.linspace(0., 1, xdim, dtype=np.float32),
+                  'lat': np.linspace(0., 1, ydim, dtype=np.float32),
+                  'depth': np.linspace(0., 1, zdim, dtype=np.float32)}
+    data = {'U': np.zeros((xdim, ydim, zdim), dtype=np.float32),
+            'V': np.zeros((xdim, ydim, zdim), dtype=np.float32),
+            'W': np.ones((xdim, ydim, zdim), dtype=np.float32)}
+    fieldset = FieldSet.from_data(data, dimensions, mesh='flat')
+
+    def DeleteParticle(particle, fieldset, time):
+        particle.delete()
+
+    pset = ParticleSet(fieldset=fieldset, pclass=ptype[mode], lon=0.5, lat=0.5, depth=0.9)
+    pset.execute(AdvectionRK4_3D, runtime=1., dt=1,
+                 recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})
+
+
 def periodicfields(xdim, ydim, uvel, vvel):
     dimensions = {'lon': np.linspace(0., 1., xdim+1, dtype=np.float32)[1:],  # don't include both 0 and 1, for periodic b.c.
                   'lat': np.linspace(0., 1., ydim+1, dtype=np.float32)[1:]}
@@ -316,38 +334,3 @@ def test_decaying_eddy(fieldset_decaying, mode, method, rtol, npart=1):
     exp_lat = [truth_decaying(x, y, endtime)[1] for x, y, in zip(lon, lat)]
     assert np.allclose(np.array([p.lon for p in pset]), exp_lon, rtol=rtol)
     assert np.allclose(np.array([p.lat for p in pset]), exp_lat, rtol=rtol)
-
-
-if __name__ == "__main__":
-    p = ArgumentParser(description="""
-Example of particle advection around an idealised peninsula""")
-    p.add_argument('mode', choices=('scipy', 'jit'), nargs='?', default='jit',
-                   help='Execution mode for performing computation')
-    p.add_argument('-p', '--particles', type=int, default=1,
-                   help='Number of particles to advect')
-    p.add_argument('-v', '--verbose', action='store_true', default=False,
-                   help='Print particle information before and after execution')
-    p.add_argument('--fieldset', choices=('stationary', 'moving', 'decaying'),
-                   default='stationary', help='Generate fieldset file with given dimensions')
-    p.add_argument('-m', '--method', choices=('RK4', 'EE', 'RK45'), default='RK4',
-                   help='Numerical method used for advection')
-    args = p.parse_args()
-    filename = 'analytical_eddies'
-
-    # Generate fieldset files according to chosen test setup
-    if args.fieldset == 'stationary':
-        fieldset = fieldset_stationary()
-    elif args.fieldset == 'moving':
-        fieldset = fieldset_moving()
-    elif args.fieldset == 'decaying':
-        fieldset = fieldset_decaying()
-
-    npart = args.particles
-    pset = ParticleSet(fieldset, pclass=ptype[args.mode],
-                       lon=np.linspace(4000, 21000, npart),
-                       lat=np.linspace(12500, 12500, npart))
-    if args.verbose:
-        print("Initial particle positions:\n%s" % pset)
-    pset.execute(kernel[args.method], dt=delta(minutes=3), runtime=delta(hours=6))
-    if args.verbose:
-        print("Final particle positions:\n%s" % pset)