Fixing XPartialslip

Author: erikvansebille

parcels/application_kernels/interpolation.py

@@ -183,10 +183,10 @@ def is_land(ti: int, zi: int, yi: int, xi: int):
     f_v = np.ones_like(eta)
 
     if lenZ == 1:
-        f_u = np.where(is_land(0, 0, 0, 0) & is_land(0, 0, 0, 1) & (eta > 0), f_u / eta, f_u)
-        f_u = np.where(is_land(0, 0, 1, 0) & is_land(0, 0, 1, 1) & (eta < 1), f_u / (1 - eta), f_u)
-        f_v = np.where(is_land(0, 0, 0, 0) & is_land(0, 0, 1, 0) & (xsi > 0), f_v / xsi, f_v)
-        f_v = np.where(is_land(0, 0, 0, 1) & is_land(0, 0, 1, 1) & (xsi < 1), f_v / (1 - xsi), f_v)
+        f_u = np.where(is_land(0, 0, 0, 0) & is_land(0, 0, 0, 1) & (eta > 0), f_u * (a + b * eta) / eta, f_u)
+        f_u = np.where(is_land(0, 0, 1, 0) & is_land(0, 0, 1, 1) & (eta < 1), f_u * (1 - b * eta) / (1 - eta), f_u)
+        f_v = np.where(is_land(0, 0, 0, 0) & is_land(0, 0, 1, 0) & (xsi > 0), f_v * (a + b * xsi) / xsi, f_v)
+        f_v = np.where(is_land(0, 0, 0, 1) & is_land(0, 0, 1, 1) & (xsi < 1), f_v * (1 - b * xsi) / (1 - xsi), f_v)
     else:
         f_u = np.where(
             is_land(0, 0, 0, 0) & is_land(0, 0, 0, 1) & is_land(0, 1, 0, 0) & is_land(0, 1, 0, 1) & (eta > 0),
@@ -195,7 +195,7 @@ def is_land(ti: int, zi: int, yi: int, xi: int):
         )
         f_u = np.where(
             is_land(0, 0, 1, 0) & is_land(0, 0, 1, 1) & is_land(0, 1, 1, 0) & is_land(0, 1, 1, 1) & (eta < 1),
-            f_u * (a - b * eta) / (1 - eta),
+            f_u * (1 - b * eta) / (1 - eta),
             f_u,
         )
         f_v = np.where(
@@ -205,7 +205,7 @@ def is_land(ti: int, zi: int, yi: int, xi: int):
         )
         f_v = np.where(
             is_land(0, 0, 0, 1) & is_land(0, 0, 1, 1) & is_land(0, 1, 0, 1) & is_land(0, 1, 1, 1) & (xsi < 1),
-            f_v * (a - b * xsi) / (1 - xsi),
+            f_v * (1 - b * xsi) / (1 - xsi),
             f_v,
         )
         f_u = np.where(
@@ -215,7 +215,7 @@ def is_land(ti: int, zi: int, yi: int, xi: int):
         )
         f_u = np.where(
             is_land(0, 1, 0, 0) & is_land(0, 1, 0, 1) & is_land(0, 1, 1, 0 & is_land(0, 1, 1, 1) & (zeta < 1)),
-            f_u * (a - b * zeta) / (1 - zeta),
+            f_u * (1 - b * zeta) / (1 - zeta),
             f_u,
         )
         f_v = np.where(
@@ -225,7 +225,7 @@ def is_land(ti: int, zi: int, yi: int, xi: int):
         )
         f_v = np.where(
             is_land(0, 1, 0, 0) & is_land(0, 1, 0, 1) & is_land(0, 1, 1, 0 & is_land(0, 1, 1, 1) & (zeta < 1)),
-            f_v * (a - b * zeta) / (1 - zeta),
+            f_v * (1 - b * zeta) / (1 - zeta),
             f_v,
         )
 

tests/v4/test_interpolation.py

@@ -11,6 +11,7 @@
     XFreeslip,
     XLinear,
     XNearest,
+    XPartialslip,
     ZeroInterpolator,
 )
 from parcels.field import Field, VectorField
@@ -88,7 +89,11 @@ def test_raw_2d_interpolation(field, func, t, z, y, x, expected):
 @pytest.mark.parametrize(
     "func, t, z, y, x, expected",
     [
+        (XPartialslip, np.timedelta64(1, "s"), 0, 0, 0.0, [[1], [1]]),
         (XFreeslip, np.timedelta64(1, "s"), 0, 0.5, 1.5, [[1], [0.5]]),
+        (XPartialslip, np.timedelta64(1, "s"), 0, 2.5, 1.5, [[0.75], [0.5]]),
+        (XFreeslip, np.timedelta64(1, "s"), 0, 2.5, 1.5, [[1], [0.5]]),
+        (XPartialslip, np.timedelta64(1, "s"), 0, 1.5, 0.5, [[0.5], [0.75]]),
         (XFreeslip, np.timedelta64(1, "s"), 0, 1.5, 0.5, [[0.5], [1]]),
         (
             XFreeslip,
@@ -103,7 +108,6 @@ def test_raw_2d_interpolation(field, func, t, z, y, x, expected):
 def test_spatial_slip_interpolation(field, func, t, z, y, x, expected):
     field.data[:] = 1.0
     field.data[:, :, 1:3, 1:3] = 0.0  # Set zero land value to test spatial slip
-    print(field.data[0, 0, :, :])
     U = field
     V = field
     UV = VectorField("UV", U, V, vector_interp_method=func)