Merge pull request #2397 from Parcels-code/fix_RK45

Fixes for the RK45 Kernel

Author: erikvansebille

src/parcels/kernels/advection.py

@@ -129,14 +129,14 @@ def AdvectionEE(particles, fieldset):  # pragma: no cover
 def AdvectionRK45(particles, fieldset):  # pragma: no cover
     """Advection of particles using adaptive Runge-Kutta 4/5 integration.
 
-    Note that this kernel requires a Particle Class that has an extra Variable 'next_dt'
-    and a FieldSet with constants 'RK45_tol' (in meters), 'RK45_min_dt' (in seconds)
-    and 'RK45_max_dt' (in seconds).
+    Note that this kernel requires a FieldSet with constants 'RK45_tol' (in meters),
+    'RK45_min_dt' (in seconds) and 'RK45_max_dt' (in seconds).
 
     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 = particles.dt / np.timedelta64(1, "s")  # TODO: improve API for converting dt to seconds
+    sign_dt = np.sign(dt)
 
     c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]
     A = [
@@ -189,17 +189,17 @@ def AdvectionRK45(particles, fieldset):  # pragma: no cover
     )
     particles.dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt)
     particles.dt = np.where(
-        particles.dt > fieldset.RK45_max_dt * np.timedelta64(1, "s"),
-        fieldset.RK45_max_dt * np.timedelta64(1, "s"),
+        np.abs(particles.dt) > np.abs(fieldset.RK45_max_dt * np.timedelta64(1, "s")),
+        fieldset.RK45_max_dt * np.timedelta64(1, "s") * sign_dt,
         particles.dt,
     )
-    particles.state = np.where(good_particles, StatusCode.Success, particles.state)
+    particles.state = np.where(good_particles, StatusCode.Evaluate, particles.state)
 
     repeat_particles = np.invert(good_particles)
     particles.dt = np.where(repeat_particles, particles.dt / 2, particles.dt)
     particles.dt = np.where(
-        particles.dt < fieldset.RK45_min_dt * np.timedelta64(1, "s"),
-        fieldset.RK45_min_dt * np.timedelta64(1, "s"),
+        np.abs(particles.dt) < np.abs(fieldset.RK45_min_dt * np.timedelta64(1, "s")),
+        fieldset.RK45_min_dt * np.timedelta64(1, "s") * sign_dt,
         particles.dt,
     )
     particles.state = np.where(repeat_particles, StatusCode.Repeat, particles.state)

tests/test_advection.py

@@ -235,11 +235,12 @@ def test_radialrotation(npart=10):
     lon = np.linspace(32, 50, npart)
     lat = np.ones(npart) * 30
     starttime = np.arange(np.timedelta64(0, "s"), npart * dt, dt)
+    endtime = np.timedelta64(10, "m")
 
     pset = parcels.ParticleSet(fieldset, lon=lon, lat=lat, time=starttime)
-    pset.execute(parcels.kernels.AdvectionRK4, endtime=np.timedelta64(10, "m"), dt=dt)
+    pset.execute(parcels.kernels.AdvectionRK4, endtime=endtime, dt=dt)
 
-    theta = 2 * np.pi * (pset.time - starttime) / np.timedelta64(24 * 3600, "s")
+    theta = 2 * np.pi * (endtime - starttime) / np.timedelta64(24 * 3600, "s")
     true_lon = (lon - 30.0) * np.cos(theta) + 30.0
     true_lat = -(lon - 30.0) * np.sin(theta) + 30.0
 
@@ -282,20 +283,21 @@ def test_moving_eddy(kernel, rtol):
 
     start_lon, start_lat, start_z = 12000, 12500, 12500
     dt = np.timedelta64(30, "m")
+    endtime = np.timedelta64(1, "h")
 
     if kernel == AdvectionRK45:
         fieldset.add_constant("RK45_tol", rtol)
 
     pset = ParticleSet(fieldset, lon=start_lon, lat=start_lat, z=start_z, time=np.timedelta64(0, "s"))
-    pset.execute(kernel, dt=dt, endtime=np.timedelta64(1, "h"))
+    pset.execute(kernel, dt=dt, endtime=endtime)
 
     def truth_moving(x_0, y_0, t):
         t /= np.timedelta64(1, "s")
         lat = y_0 - (ds.u_0 - ds.u_g) / ds.f * (1 - np.cos(ds.f * t))
         lon = x_0 + ds.u_g * t + (ds.u_0 - ds.u_g) / ds.f * np.sin(ds.f * t)
         return lon, lat
 
-    exp_lon, exp_lat = truth_moving(start_lon, start_lat, pset.time[0])
+    exp_lon, exp_lat = truth_moving(start_lon, start_lat, endtime)
     np.testing.assert_allclose(pset.lon, exp_lon, rtol=rtol)
     np.testing.assert_allclose(pset.lat, exp_lat, rtol=rtol)
     if kernel == AdvectionRK4_3D:
@@ -321,13 +323,14 @@ def test_decaying_moving_eddy(kernel, rtol):
 
     start_lon, start_lat = 10000, 10000
     dt = np.timedelta64(60, "m")
+    endtime = np.timedelta64(23, "h")
 
     if kernel == AdvectionRK45:
         fieldset.add_constant("RK45_tol", rtol)
         fieldset.add_constant("RK45_min_dt", 10 * 60)
 
     pset = ParticleSet(fieldset, lon=start_lon, lat=start_lat, time=np.timedelta64(0, "s"))
-    pset.execute(kernel, dt=dt, endtime=np.timedelta64(23, "h"))
+    pset.execute(kernel, dt=dt, endtime=endtime)
 
     def truth_moving(x_0, y_0, t):
         t /= np.timedelta64(1, "s")
@@ -343,7 +346,7 @@ def truth_moving(x_0, y_0, t):
         )
         return lon, lat
 
-    exp_lon, exp_lat = truth_moving(start_lon, start_lat, pset.time[0])
+    exp_lon, exp_lat = truth_moving(start_lon, start_lat, endtime)
     np.testing.assert_allclose(pset.lon, exp_lon, rtol=rtol)
     np.testing.assert_allclose(pset.lat, exp_lat, rtol=rtol)