Remove _constrain_dt_to_within_time_interval function (#2754)

* Removing the _constrain_dt_to_within_time_interval helper function

As not needed anymore with recent kernel loop change

This simplifies the Kernels and reverts #2378

* Using particles.dt directly in sigmagrids

Author: erikvansebille

src/parcels/kernels/_advection.py

@@ -17,76 +17,62 @@
 ]
 
 
-def _constrain_dt_to_within_time_interval(time_interval, time, dt):
-    """Helper function to make sure dt does not go outside time_interval.
-
-    This is especially relevant for higher-order RK methods (RK2, RK4, RK45),
-    which require interpolations at time + dt. If time is at the edges of the
-    time_interval (typically the last integration step), such an operation would
-    lead to an OutofTimeError.
-    """
-    if time_interval:
-        dt = np.where(time + dt <= time_interval.time_length_as_flt, dt, time_interval.time_length_as_flt - time)
-        dt = np.where(time + dt >= 0, dt, time)
-    return dt
-
-
 def AdvectionRK2(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration."""
-    dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
     (u1, v1) = fieldset.UV[particles]
-    x1, y1 = (particles.x + u1 * 0.5 * dt, particles.y + v1 * 0.5 * dt)
-    (u2, v2) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y1, x1, particles]
-    particles.dx += u2 * dt
-    particles.dy += v2 * dt
+    x1 = particles.x + u1 * 0.5 * particles.dt
+    y1 = particles.y + v1 * 0.5 * particles.dt
+    (u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles]
+    particles.dx += u2 * particles.dt
+    particles.dy += v2 * particles.dt
 
 
 def AdvectionRK2_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using second-order Runge-Kutta integration including vertical velocity."""
-    dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
     (u1, v1, w1) = fieldset.UVW[particles]
-    x1 = particles.x + u1 * 0.5 * dt
-    y1 = particles.y + v1 * 0.5 * dt
-    z1 = particles.z + w1 * 0.5 * dt
-    (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * dt, z1, y1, x1, particles]
-    particles.dx += u2 * dt
-    particles.dy += v2 * dt
-    particles.dz += w2 * dt
+    x1 = particles.x + u1 * 0.5 * particles.dt
+    y1 = particles.y + v1 * 0.5 * particles.dt
+    z1 = particles.z + w1 * 0.5 * particles.dt
+    (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * particles.dt, z1, y1, x1, particles]
+    particles.dx += u2 * particles.dt
+    particles.dy += v2 * particles.dt
+    particles.dz += w2 * particles.dt
 
 
 def AdvectionRK4(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration."""
-    dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
     (u1, v1) = fieldset.UV[particles]
-    x1, y1 = (particles.x + u1 * 0.5 * dt, particles.y + v1 * 0.5 * dt)
-    (u2, v2) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y1, x1, particles]
-    x2, y2 = (particles.x + u2 * 0.5 * dt, particles.y + v2 * 0.5 * dt)
-    (u3, v3) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y2, x2, particles]
-    x3, y3 = (particles.x + u3 * dt, particles.y + v3 * dt)
-    (u4, v4) = fieldset.UV[particles.t + dt, particles.z, y3, x3, particles]
-    particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt
-    particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt
+    x1 = particles.x + u1 * 0.5 * particles.dt
+    y1 = particles.y + v1 * 0.5 * particles.dt
+    (u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles]
+    x2 = particles.x + u2 * 0.5 * particles.dt
+    y2 = particles.y + v2 * 0.5 * particles.dt
+    (u3, v3) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y2, x2, particles]
+    x3 = particles.x + u3 * particles.dt
+    y3 = particles.y + v3 * particles.dt
+    (u4, v4) = fieldset.UV[particles.t + particles.dt, particles.z, y3, x3, particles]
+    particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * particles.dt
+    particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particles.dt
 
 
 def AdvectionRK4_3D(particles, fieldset):  # pragma: no cover
     """Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
-    dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
     (u1, v1, w1) = fieldset.UVW[particles]
-    x1 = particles.x + u1 * 0.5 * dt
-    y1 = particles.y + v1 * 0.5 * dt
-    z1 = particles.z + w1 * 0.5 * dt
-    (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * dt, z1, y1, x1, particles]
-    x2 = particles.x + u2 * 0.5 * dt
-    y2 = particles.y + v2 * 0.5 * dt
-    z2 = particles.z + w2 * 0.5 * dt
-    (u3, v3, w3) = fieldset.UVW[particles.t + 0.5 * dt, z2, y2, x2, particles]
-    x3 = particles.x + u3 * dt
-    y3 = particles.y + v3 * dt
-    z3 = particles.z + w3 * dt
-    (u4, v4, w4) = fieldset.UVW[particles.t + dt, z3, y3, x3, particles]
-    particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
-    particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
-    particles.dz += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * dt
+    x1 = particles.x + u1 * 0.5 * particles.dt
+    y1 = particles.y + v1 * 0.5 * particles.dt
+    z1 = particles.z + w1 * 0.5 * particles.dt
+    (u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * particles.dt, z1, y1, x1, particles]
+    x2 = particles.x + u2 * 0.5 * particles.dt
+    y2 = particles.y + v2 * 0.5 * particles.dt
+    z2 = particles.z + w2 * 0.5 * particles.dt
+    (u3, v3, w3) = fieldset.UVW[particles.t + 0.5 * particles.dt, z2, y2, x2, particles]
+    x3 = particles.x + u3 * particles.dt
+    y3 = particles.y + v3 * particles.dt
+    z3 = particles.z + w3 * particles.dt
+    (u4, v4, w4) = fieldset.UVW[particles.t + particles.dt, z3, y3, x3, particles]
+    particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * particles.dt
+    particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * particles.dt
+    particles.dz += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * particles.dt
 
 
 def AdvectionEE(particles, fieldset):  # pragma: no cover
@@ -105,8 +91,7 @@ def AdvectionRK45(particles, fieldset):  # pragma: no cover
     Time-step dt is halved if error is larger than fieldset.RK45_tol,
     and doubled if error is smaller than 1/10th of tolerance.
     """
-    dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
-    sign_dt = np.sign(dt)
+    sign_dt = np.sign(particles.dt)
 
     c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]
     A = [
@@ -120,42 +105,37 @@ def AdvectionRK45(particles, fieldset):  # pragma: no cover
     b5 = [16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0]
 
     (u1, v1) = fieldset.UV[particles]
-    x1, y1 = (particles.x + u1 * A[0][0] * dt, particles.y + v1 * A[0][0] * dt)
-    (u2, v2) = fieldset.UV[particles.t + c[0] * dt, particles.z, y1, x1, particles]
-    x2, y2 = (
-        particles.x + (u1 * A[1][0] + u2 * A[1][1]) * dt,
-        particles.y + (v1 * A[1][0] + v2 * A[1][1]) * dt,
-    )
-    (u3, v3) = fieldset.UV[particles.t + c[1] * dt, particles.z, y2, x2, particles]
-    x3, y3 = (
-        particles.x + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * dt,
-        particles.y + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * dt,
-    )
-    (u4, v4) = fieldset.UV[particles.t + c[2] * dt, particles.z, y3, x3, particles]
-    x4, y4 = (
-        particles.x + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * dt,
-        particles.y + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * dt,
-    )
-    (u5, v5) = fieldset.UV[particles.t + c[3] * dt, particles.z, y4, x4, particles]
-    x5, y5 = (
-        particles.x + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * dt,
-        particles.y + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * dt,
-    )
-    (u6, v6) = fieldset.UV[particles.t + c[4] * dt, particles.z, y5, x5, particles]
-
-    x_4th = (u1 * b4[0] + u2 * b4[1] + u3 * b4[2] + u4 * b4[3] + u5 * b4[4]) * dt
-    y_4th = (v1 * b4[0] + v2 * b4[1] + v3 * b4[2] + v4 * b4[3] + v5 * b4[4]) * dt
-    x_5th = (u1 * b5[0] + u2 * b5[1] + u3 * b5[2] + u4 * b5[3] + u5 * b5[4] + u6 * b5[5]) * dt
-    y_5th = (v1 * b5[0] + v2 * b5[1] + v3 * b5[2] + v4 * b5[3] + v5 * b5[4] + v6 * b5[5]) * dt
+    x1 = particles.x + u1 * A[0][0] * particles.dt
+    y1 = particles.y + v1 * A[0][0] * particles.dt
+    (u2, v2) = fieldset.UV[particles.t + c[0] * particles.dt, particles.z, y1, x1, particles]
+    x2 = particles.x + (u1 * A[1][0] + u2 * A[1][1]) * particles.dt
+    y2 = particles.y + (v1 * A[1][0] + v2 * A[1][1]) * particles.dt
+    (u3, v3) = fieldset.UV[particles.t + c[1] * particles.dt, particles.z, y2, x2, particles]
+    x3 = particles.x + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * particles.dt
+    y3 = particles.y + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * particles.dt
+    (u4, v4) = fieldset.UV[particles.t + c[2] * particles.dt, particles.z, y3, x3, particles]
+    x4 = particles.x + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * particles.dt
+    y4 = particles.y + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * particles.dt
+    (u5, v5) = fieldset.UV[particles.t + c[3] * particles.dt, particles.z, y4, x4, particles]
+    x5 = particles.x + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * particles.dt
+    y5 = particles.y + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * particles.dt
+    (u6, v6) = fieldset.UV[particles.t + c[4] * particles.dt, particles.z, y5, x5, particles]
+
+    x_4th = (u1 * b4[0] + u2 * b4[1] + u3 * b4[2] + u4 * b4[3] + u5 * b4[4]) * particles.dt
+    y_4th = (v1 * b4[0] + v2 * b4[1] + v3 * b4[2] + v4 * b4[3] + v5 * b4[4]) * particles.dt
+    x_5th = (u1 * b5[0] + u2 * b5[1] + u3 * b5[2] + u4 * b5[3] + u5 * b5[4] + u6 * b5[5]) * particles.dt
+    y_5th = (v1 * b5[0] + v2 * b5[1] + v3 * b5[2] + v4 * b5[3] + v5 * b5[4] + v6 * b5[5]) * particles.dt
 
     kappa = np.sqrt(np.pow(x_5th - x_4th, 2) + np.pow(y_5th - y_4th, 2))
 
-    good_particles = (kappa <= fieldset.RK45_tol) | (np.fabs(dt) <= np.fabs(fieldset.RK45_min_dt))
+    good_particles = (kappa <= fieldset.RK45_tol) | (np.fabs(particles.dt) <= np.fabs(fieldset.RK45_min_dt))
     particles.dx += np.where(good_particles, x_5th, 0)
     particles.dy += np.where(good_particles, y_5th, 0)
 
     increase_dt_particles = (
-        good_particles & (kappa <= fieldset.RK45_tol / 10) & (np.fabs(dt * 2) <= np.fabs(fieldset.RK45_max_dt))
+        good_particles
+        & (kappa <= fieldset.RK45_tol / 10)
+        & (np.fabs(particles.dt * 2) <= np.fabs(fieldset.RK45_max_dt))
     )
     particles.next_dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt)
     particles.next_dt = np.where(

src/parcels/kernels/_sigmagrids.py

@@ -2,8 +2,6 @@
 
 import numpy as np
 
-from parcels.kernels._advection import _constrain_dt_to_within_time_interval
-
 
 def convert_z_to_sigma_croco(fieldset, t, z, y, x, particle):
     """Calculate local sigma level of the particles, by linearly interpolating the
@@ -49,27 +47,26 @@ def AdvectionRK2_3D_CROCO(particles, fieldset):  # pragma: no cover
             category=RuntimeWarning,
         )  # Needed because of linear sampling of W with sigma conversion
 
-        dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
         sigma = particles.z / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x]
 
         sig = convert_z_to_sigma_croco(fieldset, particles.t, particles.z, particles.y, particles.x, particles)
         (u1, v1) = fieldset.UV[particles.t, sig, particles.y, particles.x, particles]
         w1 = fieldset.W[particles.t, sig, particles.y, particles.x, particles]
         w1 *= sigma / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x]
-        x1 = particles.x + u1 * 0.5 * dt
-        y1 = particles.y + v1 * 0.5 * dt
-        sig_dep1 = sigma + w1 * 0.5 * dt
+        x1 = particles.x + u1 * 0.5 * particles.dt
+        y1 = particles.y + v1 * 0.5 * particles.dt
+        sig_dep1 = sigma + w1 * 0.5 * particles.dt
         dep1 = sig_dep1 * fieldset.h[particles.t, np.zeros_like(particles.z), y1, x1]
 
-        sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * dt, dep1, y1, x1, particles)
-        (u2, v2) = fieldset.UV[particles.t + 0.5 * dt, sig1, y1, x1, particles]
-        w2 = fieldset.W[particles.t + 0.5 * dt, sig1, y1, x1, particles]
-        w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y1, x1]
-        x2 = particles.x + u2 * 0.5 * dt
-        y2 = particles.y + v2 * 0.5 * dt
-        sig_dep2 = sigma + w2 * 0.5 * dt
-        dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y2, x2]
+        sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * particles.dt, dep1, y1, x1, particles)
+        (u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles]
+        w2 = fieldset.W[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles]
+        w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y1, x1]
+        x2 = particles.x + u2 * 0.5 * particles.dt
+        y2 = particles.y + v2 * 0.5 * particles.dt
+        sig_dep2 = sigma + w2 * 0.5 * particles.dt
+        dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y2, x2]
 
-        particles.dx += u2 * dt
-        particles.dy += v2 * dt
+        particles.dx += u2 * particles.dt
+        particles.dy += v2 * particles.dt
         particles.dz += (dep1 - particles.z) + (dep2 - particles.z)