Merge pull request #1052 from OceanParcels/partial_slip_field_interpolation

Partialslip and freeslip interpolation

Author: erikvansebille

parcels/examples/documentation_unstuck_Agrid.ipynb

@@ -884,7 +884,7 @@ def interpolator2D(self, ti, z, y, x, particle=None):
             xii = xi if xsi <= .5 else xi+1
             yii = yi if eta <= .5 else yi+1
             return self.data[ti, yii, xii]
-        elif self.interp_method in ['linear', 'bgrid_velocity']:
+        elif self.interp_method in ['linear', 'bgrid_velocity', 'partialslip', 'freeslip']:
             val = (1-xsi)*(1-eta) * self.data[ti, yi, xi] + \
                 xsi*(1-eta) * self.data[ti, yi, xi+1] + \
                 xsi*eta * self.data[ti, yi+1, xi+1] + \
@@ -972,7 +972,7 @@ def interpolator3D(self, ti, z, y, x, time, particle=None):
                     xsi * eta * data[yi + 1, xi + 1] + \
                     (1 - xsi) * eta * data[yi + 1, xi]
                 return (1 - zeta) * f0 + zeta * f1
-        elif self.interp_method in ['linear', 'bgrid_velocity', 'bgrid_w_velocity']:
+        elif self.interp_method in ['linear', 'bgrid_velocity', 'bgrid_w_velocity', 'partialslip', 'freeslip']:
             if self.interp_method == 'bgrid_velocity':
                 if self.gridindexingtype == 'mom5':
                     zeta = 1.
@@ -1635,8 +1635,91 @@ def spatial_c_grid_interpolation3D(self, ti, z, y, x, time, particle=None):
             w = self.W.units.to_target(w, x, y, z)
         return (u, v, w)
 
+    def _is_land2D(self, di, yi, xi):
+        if self.U.data.ndim == 3:
+            if di < np.shape(self.U.data)[0]:
+                return np.isclose(self.U.data[di, yi, xi], 0.) and np.isclose(self.V.data[di, yi, xi], 0.)
+            else:
+                return True
+        else:
+            if di < self.U.grid.zdim and yi < np.shape(self.U.data)[-2] and xi < np.shape(self.U.data)[-1]:
+                return np.isclose(self.U.data[0, di, yi, xi], 0.) and np.isclose(self.V.data[0, di, yi, xi], 0.)
+            else:
+                return True
+
+    def spatial_slip_interpolation(self, ti, z, y, x, time, particle=None):
+        (xsi, eta, zeta, xi, yi, zi) = self.U.search_indices(x, y, z, ti, time, particle=particle)
+        di = ti if self.U.grid.zdim == 1 else zi  # general third dimension
+
+        f_u, f_v, f_w = 1, 1, 1
+        if self._is_land2D(di, yi, xi) and self._is_land2D(di, yi, xi+1) and self._is_land2D(di+1, yi, xi) \
+                and self._is_land2D(di+1, yi, xi+1) and eta > 0:
+            if self.U.interp_method == 'partialslip':
+                f_u = f_u * (.5 + .5 * eta) / eta
+                if self.vector_type == '3D':
+                    f_w = f_w * (.5 + .5 * eta) / eta
+            elif self.U.interp_method == 'freeslip':
+                f_u = f_u / eta
+                if self.vector_type == '3D':
+                    f_w = f_w / eta
+        if self._is_land2D(di, yi+1, xi) and self._is_land2D(di, yi+1, xi+1) and self._is_land2D(di+1, yi+1, xi) \
+                and self._is_land2D(di+1, yi+1, xi+1) and eta < 1:
+            if self.U.interp_method == 'partialslip':
+                f_u = f_u * (1 - .5 * eta) / (1 - eta)
+                if self.vector_type == '3D':
+                    f_w = f_w * (1 - .5 * eta) / (1 - eta)
+            elif self.U.interp_method == 'freeslip':
+                f_u = f_u / (1 - eta)
+                if self.vector_type == '3D':
+                    f_w = f_w / (1 - eta)
+        if self._is_land2D(di, yi, xi) and self._is_land2D(di, yi+1, xi) and self._is_land2D(di+1, yi, xi) \
+                and self._is_land2D(di+1, yi+1, xi) and xsi > 0:
+            if self.U.interp_method == 'partialslip':
+                f_v = f_v * (.5 + .5 * xsi) / xsi
+                if self.vector_type == '3D':
+                    f_w = f_w * (.5 + .5 * xsi) / xsi
+            elif self.U.interp_method == 'freeslip':
+                f_v = f_v / xsi
+                if self.vector_type == '3D':
+                    f_w = f_w / xsi
+        if self._is_land2D(di, yi, xi+1) and self._is_land2D(di, yi+1, xi+1) and self._is_land2D(di+1, yi, xi+1) \
+                and self._is_land2D(di+1, yi+1, xi+1) and xsi < 1:
+            if self.U.interp_method == 'partialslip':
+                f_v = f_v * (1 - .5 * xsi) / (1 - xsi)
+                if self.vector_type == '3D':
+                    f_w = f_w * (1 - .5 * xsi) / (1 - xsi)
+            elif self.U.interp_method == 'freeslip':
+                f_v = f_v / (1 - xsi)
+                if self.vector_type == '3D':
+                    f_w = f_w / (1 - xsi)
+        if self.U.grid.zdim > 1:
+            if self._is_land2D(di, yi, xi) and self._is_land2D(di, yi, xi+1) and self._is_land2D(di, yi+1, xi) \
+                    and self._is_land2D(di, yi+1, xi+1) and zeta > 0:
+                if self.U.interp_method == 'partialslip':
+                    f_u = f_u * (.5 + .5 * zeta) / zeta
+                    f_v = f_v * (.5 + .5 * zeta) / zeta
+                elif self.U.interp_method == 'freeslip':
+                    f_u = f_u / zeta
+                    f_v = f_v / zeta
+            if self._is_land2D(di+1, yi, xi) and self._is_land2D(di+1, yi, xi+1) and self._is_land2D(di+1, yi+1, xi) \
+                    and self._is_land2D(di+1, yi+1, xi+1) and zeta < 1:
+                if self.U.interp_method == 'partialslip':
+                    f_u = f_u * (1 - .5 * zeta) / (1 - zeta)
+                    f_v = f_v * (1 - .5 * zeta) / (1 - zeta)
+                elif self.U.interp_method == 'freeslip':
+                    f_u = f_u / (1 - zeta)
+                    f_v = f_v / (1 - zeta)
+
+        u = f_u * self.U.eval(time, z, y, x, particle)
+        v = f_v * self.V.eval(time, z, y, x, particle)
+        if self.vector_type == '3D':
+            w = f_w * self.W.eval(time, z, y, x, particle)
+            return u, v, w
+        else:
+            return u, v
+
     def eval(self, time, z, y, x, particle=None):
-        if self.U.interp_method != 'cgrid_velocity':
+        if self.U.interp_method not in ['cgrid_velocity', 'partialslip', 'freeslip']:
             u = self.U.eval(time, z, y, x, particle=particle, applyConversion=False)
             v = self.V.eval(time, z, y, x, particle=particle, applyConversion=False)
             u = self.U.units.to_target(u, x, y, z)
@@ -1648,19 +1731,22 @@ def eval(self, time, z, y, x, particle=None):
             else:
                 return (u, v)
         else:
+            interp = {'cgrid_velocity': {'2D': self.spatial_c_grid_interpolation2D, '3D': self.spatial_c_grid_interpolation3D},
+                      'partialslip': {'2D': self.spatial_slip_interpolation, '3D': self.spatial_slip_interpolation},
+                      'freeslip': {'2D': self.spatial_slip_interpolation, '3D': self.spatial_slip_interpolation}}
             grid = self.U.grid
             (ti, periods) = self.U.time_index(time)
             time -= periods*(grid.time_full[-1]-grid.time_full[0])
             if ti < grid.tdim-1 and time > grid.time[ti]:
                 t0 = grid.time[ti]
                 t1 = grid.time[ti + 1]
                 if self.vector_type == '3D':
-                    (u0, v0, w0) = self.spatial_c_grid_interpolation3D(ti, z, y, x, time, particle=particle)
-                    (u1, v1, w1) = self.spatial_c_grid_interpolation3D(ti + 1, z, y, x, time, particle=particle)
+                    (u0, v0, w0) = interp[self.U.interp_method]['3D'](ti, z, y, x, time, particle=particle)
+                    (u1, v1, w1) = interp[self.U.interp_method]['3D'](ti + 1, z, y, x, time, particle=particle)
                     w = w0 + (w1 - w0) * ((time - t0) / (t1 - t0))
                 else:
-                    (u0, v0) = self.spatial_c_grid_interpolation2D(ti, z, y, x, time, particle=particle)
-                    (u1, v1) = self.spatial_c_grid_interpolation2D(ti + 1, z, y, x, time, particle=particle)
+                    (u0, v0) = interp[self.U.interp_method]['2D'](ti, z, y, x, time, particle=particle)
+                    (u1, v1) = interp[self.U.interp_method]['2D'](ti + 1, z, y, x, time, particle=particle)
                 u = u0 + (u1 - u0) * ((time - t0) / (t1 - t0))
                 v = v0 + (v1 - v0) * ((time - t0) / (t1 - t0))
                 if self.vector_type == '3D':
@@ -1672,9 +1758,9 @@ def eval(self, time, z, y, x, particle=None):
                 # of the defined time range or if we have hit an
                 # exact value in the time array.
                 if self.vector_type == '3D':
-                    return self.spatial_c_grid_interpolation3D(ti, z, y, x, grid.time[ti], particle=particle)
+                    return interp[self.U.interp_method]['3D'](ti, z, y, x, grid.time[ti], particle=particle)
                 else:
-                    return self.spatial_c_grid_interpolation2D(ti, z, y, x, grid.time[ti], particle=particle)
+                    return interp[self.U.interp_method]['2D'](ti, z, y, x, grid.time[ti], particle=particle)
 
     def __getitem__(self, key):
         if _isParticle(key):

parcels/field.py

@@ -587,6 +587,8 @@ def from_c_grid_dataset(cls, filenames, variables, dimensions, indices=None, mes
         if 'U' in dimensions and 'W' in dimensions and dimensions['U'] != dimensions['W']:
             raise ValueError("On a C-grid, the dimensions of velocities should be the corners (f-points) of the cells, so the same for U, V and W. "
                              "See also https://nbviewer.jupyter.org/github/OceanParcels/parcels/blob/master/parcels/examples/documentation_indexing.ipynb")
+        if 'interp_method' in kwargs.keys():
+            raise TypeError("On a C-grid, the interpolation method for velocities should not be overridden")
 
         interp_method = {}
         for v in variables:

parcels/fieldset.py

@@ -18,7 +18,7 @@ typedef float type_coord;
 
 typedef enum
   {
-    LINEAR=0, NEAREST=1, CGRID_VELOCITY=2, CGRID_TRACER=3, BGRID_VELOCITY=4, BGRID_W_VELOCITY=5, BGRID_TRACER=6, LINEAR_INVDIST_LAND_TRACER=7
+    LINEAR=0, NEAREST=1, CGRID_VELOCITY=2, CGRID_TRACER=3, BGRID_VELOCITY=4, BGRID_W_VELOCITY=5, BGRID_TRACER=6, LINEAR_INVDIST_LAND_TRACER=7, PARTIALSLIP=8, FREESLIP=9
   } InterpCode;
 
 typedef enum