Adding support and test for RK45 advection

Author: erikvansebille

parcels/application_kernels/advection.py

@@ -113,7 +113,7 @@ def AdvectionRK45(particle, fieldset, time):  # 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 = min(particle.next_dt, fieldset.RK45_max_dt) / np.timedelta64(1, "s")  # noqa TODO improve API for converting dt to seconds
+    dt = min(particle.next_dt / np.timedelta64(1, "s"), fieldset.RK45_max_dt)  # noqa TODO improve API for converting dt to seconds
     c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]
     A = [
         [1.0 / 4.0, 0.0, 0.0, 0.0, 0.0],

parcels/kernel.py

@@ -384,13 +384,12 @@ def evaluate_particle(self, p, endtime):
                 return p
 
             pre_dt = p.dt
-            # TODO implement below later again
-            # try:  # Use next_dt from AdvectionRK45 if it is set
-            #     if abs(endtime - p.time_nextloop) < abs(p.next_dt) - 1e-6:
-            #         p.next_dt = abs(endtime - p.time_nextloop) * sign_dt
-            # except AttributeError:
-            if sign_dt * (endtime - p.time_nextloop) <= p.dt:
-                p.dt = sign_dt * (endtime - p.time_nextloop)
+            try:  # Use next_dt from AdvectionRK45 if it is set
+                if abs(endtime - p.time_nextloop) < abs(p.next_dt) - np.timedelta64(1000, "ns"):
+                    p.next_dt = sign_dt * (endtime - p.time_nextloop)
+            except KeyError:
+                if sign_dt * (endtime - p.time_nextloop) <= p.dt:
+                    p.dt = sign_dt * (endtime - p.time_nextloop)
             res = self._pyfunc(p, self._fieldset, p.time_nextloop)
 
             if res is None:

parcels/particleset.py

@@ -156,7 +156,7 @@ def __init__(
                 if isinstance(v.initial, attrgetter):
                     initial = v.initial(self)
                 else:
-                    initial = v.initial * np.ones(len(trajectory_ids), dtype=v.dtype)
+                    initial = [np.array(v.initial, dtype=v.dtype)] * len(trajectory_ids)
                 self._data[v.name] = initial
 
         # update initial values provided on ParticleSet creation

tests/v4/test_advection.py

@@ -2,7 +2,7 @@
 import pytest
 
 from parcels._datasets.structured.generic import simple_UV_dataset
-from parcels.application_kernels import AdvectionEE, AdvectionRK4, AdvectionRK4_3D
+from parcels.application_kernels import AdvectionEE, AdvectionRK4, AdvectionRK4_3D, AdvectionRK45
 from parcels.field import Field, VectorField
 from parcels.fieldset import FieldSet
 from parcels.particle import Particle, Variable
@@ -101,7 +101,7 @@ def TriLinear(  # TODO move to interpolation file
     "EE": AdvectionEE,
     "RK4": AdvectionRK4,
     "RK4_3D": AdvectionRK4_3D,
-    # "RK45": AdvectionRK45,
+    "RK45": AdvectionRK45,
     # "AA": AdvectionAnalytical,
     # "AdvDiffEM": AdvectionDiffusionEM,
     # "AdvDiffM1": AdvectionDiffusionM1,
@@ -222,7 +222,7 @@ def SubmergeParticle(particle, fieldset, time):  # pragma: no cover
         # ("AdvDiffM1", 1e-2),
         ("RK4", 1e-5),
         ("RK4_3D", 1e-5),
-        # ("RK45", 1e-5),
+        ("RK45", 1e-5),
     ],
 )
 def test_moving_eddy(method, rtol):
@@ -258,7 +258,10 @@ def truth_moving(x_0, y_0, t):
         fieldset = FieldSet([U, V, UV])
         start_depth = 0
 
-    RK45Particles = Particle.add_variable(Variable("next_dt", initial=dt))
+    if method == "RK45":
+        # Use RK45Particles to set next_dt
+        RK45Particles = Particle.add_variable(Variable("next_dt", initial=dt, dtype=np.timedelta64))
+        fieldset.add_constant("RK45_tol", 1e-6)
 
     pclass = RK45Particles if method == "RK45" else Particle
     pset = ParticleSet(