Fixed particle on boundary issue

CHANGELOG.md

@@ -16,10 +16,11 @@ and the project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0.
 - Multiplicity to the releaser module
 
 
-## [Unreleased] - 2024-09-19
+## [2.0.3] - 2024-09-19
 ### Fixed
 - Can import local IBM and gridforce modules
 - zROMS module now works with legacy config file
+- Simulation no longer breaks if particles reach domain boundary
 ### Changed
 - Logger now outputs current time
 

ladim/__init__.py

@@ -1,3 +1,3 @@
-__version__ = '2.0.2'
+__version__ = '2.0.3'
 
 from .main import main, run

ladim/gridforce/ROMS.py

@@ -182,6 +182,8 @@ def sample_depth(self, X, Y):
         """Return the depth of grid cells"""
         I = X.round().astype(int) - self.i0
         J = Y.round().astype(int) - self.j0
+        I = np.minimum(np.maximum(I, 0), self.H.shape[1] - 1)
+        J = np.minimum(np.maximum(J, 0), self.H.shape[0] - 1)
         return self.H[J, I]
 
     def lonlat(self, X, Y, method="bilinear"):
@@ -191,6 +193,8 @@ def lonlat(self, X, Y, method="bilinear"):
         # else: containing grid cell, less accurate
         I = X.round().astype("int") - self.i0
         J = Y.round().astype("int") - self.j0
+        I = np.minimum(np.maximum(I, 0), self.lon.shape[1] - 1)
+        J = np.minimum(np.maximum(J, 0), self.lon.shape[0] - 1)
         return self.lon[J, I], self.lat[J, I]
 
     def ingrid(self, X: np.ndarray, Y: np.ndarray) -> np.ndarray:
@@ -206,6 +210,11 @@ def onland(self, X, Y):
         """Returns True for points on land"""
         I = X.round().astype(int) - self.i0
         J = Y.round().astype(int) - self.j0
+
+        # Constrain to valid indices
+        I = np.minimum(np.maximum(I, 0), self.M.shape[-1] - 1)
+        J = np.minimum(np.maximum(J, 0), self.M.shape[-2] - 1)
+
         return self.M[J, I] < 1
 
     # Error if point outside
@@ -802,6 +811,11 @@ def sample3D(F, X, Y, K, A, method="bilinear"):
     # else:  method == 'nearest'
     I = X.round().astype("int")
     J = Y.round().astype("int")
+
+    # Constrain to valid indices
+    I = np.minimum(np.maximum(I, 0), F.shape[-1] - 1)
+    J = np.minimum(np.maximum(J, 0), F.shape[-2] - 1)
+
     return F[K, J, I]