Add next_dt in RK45

This is needed in case dt is signalled for increasing in the RK45 Kernel; otherwise the wrong dt gets added in the add_positionupdate_kernel

Author: erikvansebille

src/parcels/_core/kernel.py

@@ -123,6 +123,10 @@ def PositionUpdate(particles, fieldset):  # pragma: no cover
             particles.dlat = 0
             particles.dz = 0
 
+            if hasattr(self.fieldset, "RK45_tol"):
+                # Update dt in case it's increased in RK45 kernel
+                particles.dt = particles.next_dt
+
         self._pyfuncs = (PositionUpdate + self)._pyfuncs
 
     def check_fieldsets_in_kernels(self, pyfunc):  # TODO v4: this can go into another method? assert_is_compatible()?
@@ -138,6 +142,8 @@ def check_fieldsets_in_kernels(self, pyfunc):  # TODO v4: this can go into anoth
                 if self._fieldset.U.grid._gtype not in [GridType.CurvilinearZGrid, GridType.RectilinearZGrid]:
                     raise NotImplementedError("Analytical Advection only works with Z-grids in the vertical")
             elif pyfunc is AdvectionRK45:
+                if "next_dt" not in [v.name for v in self.ptype.variables]:
+                    raise ValueError('ParticleClass requires a "next_dt" for AdvectionRK45 Kernel.')
                 if not hasattr(self.fieldset, "RK45_tol"):
                     warnings.warn(
                         "Setting RK45 tolerance to 10 m. Use fieldset.add_constant('RK45_tol', [distance]) to change.",
@@ -240,6 +246,9 @@ def execute(self, pset, endtime, dt):
                 return StatusCode.Success
 
             # adapt dt to end exactly on endtime
+            if hasattr(self.fieldset, "RK45_tol"):
+                pset.next_dt = np.maximum(np.minimum(pset.next_dt, time_to_endtime), 0)
+
             if compute_time_direction == 1:
                 pset.dt = np.maximum(np.minimum(pset.dt, time_to_endtime), 0)
             else:

src/parcels/kernels/advection.py

@@ -188,11 +188,11 @@ def AdvectionRK45(particles, fieldset):  # pragma: no cover
     increase_dt_particles = (
         good_particles & (kappa <= fieldset.RK45_tol / 10) & (np.fabs(dt * 2) <= np.fabs(fieldset.RK45_max_dt))
     )
-    particles.dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt)
-    particles.dt = np.where(
+    particles.next_dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt)
+    particles.next_dt = np.where(
         np.abs(particles.next_dt) > np.abs(fieldset.RK45_max_dt * np.timedelta64(1, "s")),
         fieldset.RK45_max_dt * np.timedelta64(1, "s") * sign_dt,
-        particles.dt,
+        particles.next_dt,
     )
     particles.state = np.where(good_particles, StatusCode.Evaluate, particles.state)
 

tests/test_advection.py

@@ -285,8 +285,13 @@ def test_moving_eddy(kernel, rtol):
 
     if kernel == AdvectionRK45:
         fieldset.add_constant("RK45_tol", rtol)
+        MyParticle = Particle.add_variable(Variable("next_dt", dtype="timedelta64[s]", initial=dt))
+    else:
+        MyParticle = Particle
 
-    pset = ParticleSet(fieldset, lon=start_lon, lat=start_lat, z=start_z, time=np.timedelta64(0, "s"))
+    pset = ParticleSet(
+        fieldset, pclass=MyParticle, lon=start_lon, lat=start_lat, z=start_z, time=np.timedelta64(0, "s")
+    )
     pset.execute(kernel, dt=dt, endtime=np.timedelta64(1, "h"))
 
     def truth_moving(x_0, y_0, t):
@@ -325,8 +330,11 @@ def test_decaying_moving_eddy(kernel, rtol):
     if kernel == AdvectionRK45:
         fieldset.add_constant("RK45_tol", rtol)
         fieldset.add_constant("RK45_min_dt", 10 * 60)
+        MyParticle = Particle.add_variable(Variable("next_dt", dtype="timedelta64[s]", initial=dt))
+    else:
+        MyParticle = Particle
 
-    pset = ParticleSet(fieldset, lon=start_lon, lat=start_lat, time=np.timedelta64(0, "s"))
+    pset = ParticleSet(fieldset, pclass=MyParticle, lon=start_lon, lat=start_lat, time=np.timedelta64(0, "s"))
     pset.execute(kernel, dt=dt, endtime=np.timedelta64(23, "h"))
 
     def truth_moving(x_0, y_0, t):
@@ -373,13 +381,15 @@ def test_stommelgyre_fieldset(kernel, rtol, grid_type):
     start_lon = np.linspace(10e3, 100e3, npart)
     start_lat = np.ones_like(start_lon) * 5000e3
 
-    if kernel == AdvectionRK45:
-        fieldset.add_constant("RK45_tol", rtol)
-
     SampleParticle = Particle.add_variable(
         [Variable("p", initial=0.0, dtype=np.float32), Variable("p_start", initial=0.0, dtype=np.float32)]
     )
 
+    if kernel == AdvectionRK45:
+        fieldset.add_constant("RK45_tol", rtol)
+        fieldset.add_constant("RK45_min_dt", 5 * 60)
+        SampleParticle = SampleParticle.add_variable(Variable("next_dt", dtype="timedelta64[s]", initial=dt))
+
     def UpdateP(particles, fieldset):  # pragma: no cover
         particles.p = fieldset.P[particles.time, particles.z, particles.lat, particles.lon]
         particles.p_start = np.where(particles.time == 0, particles.p, particles.p_start)
@@ -394,7 +404,7 @@ def UpdateP(particles, fieldset):  # pragma: no cover
     [
         (AdvectionRK2, 2e-2),
         (AdvectionRK4, 5e-3),
-        (AdvectionRK45, 1e-4),
+        (AdvectionRK45, 1e-3),
     ],
 )
 @pytest.mark.parametrize("grid_type", ["A"])  # TODO also implement C-grid once available
@@ -413,13 +423,15 @@ def test_peninsula_fieldset(kernel, rtol, grid_type):
     start_lat = np.linspace(3e3, 47e3, npart)
     start_lon = 3e3 * np.ones_like(start_lat)
 
-    if kernel == AdvectionRK45:
-        fieldset.add_constant("RK45_tol", rtol)
-
     SampleParticle = Particle.add_variable(
         [Variable("p", initial=0.0, dtype=np.float32), Variable("p_start", initial=0.0, dtype=np.float32)]
     )
 
+    if kernel == AdvectionRK45:
+        fieldset.add_constant("RK45_tol", rtol)
+        fieldset.add_constant("RK45_min_dt", 5 * 60)
+        SampleParticle = SampleParticle.add_variable(Variable("next_dt", dtype="timedelta64[s]", initial=dt))
+
     def UpdateP(particles, fieldset):  # pragma: no cover
         particles.p = fieldset.P[particles.time, particles.z, particles.lat, particles.lon]
         particles.p_start = np.where(particles.time == 0, particles.p, particles.p_start)

tests/test_kernel.py

@@ -11,7 +11,7 @@
 )
 from parcels._datasets.structured.generic import datasets as datasets_structured
 from parcels.interpolators import XLinear
-from parcels.kernels import AdvectionRK4
+from parcels.kernels import AdvectionRK4, AdvectionRK45
 from tests.common_kernels import MoveEast, MoveNorth
 
 
@@ -88,6 +88,14 @@ def test_kernel_from_list_error_checking(fieldset):
         assert kernels_mixed.funcname == "AdvectionRK4MoveEastMoveNorth"
 
 
+def test_RK45Kernel_error_no_next_dt(fieldset):
+    """Tests that kernel throws error if Particle class does not have next_dt for RK45"""
+    pset = ParticleSet(fieldset, lon=[0.5], lat=[0.5])
+
+    with pytest.raises(ValueError, match='ParticleClass requires a "next_dt" for AdvectionRK45 Kernel.'):
+        pset.Kernel(AdvectionRK45)
+
+
 def test_kernel_signature(fieldset):
     pset = ParticleSet(fieldset, lon=[0.5], lat=[0.5])