diff --git a/src/parcels/kernels/advection.py b/src/parcels/kernels/advection.py index de0e78054..917402ea1 100644 --- a/src/parcels/kernels/advection.py +++ b/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) diff --git a/tests/test_advection.py b/tests/test_advection.py index 1f3e04df7..188e8343d 100644 --- a/tests/test_advection.py +++ b/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,12 +283,13 @@ 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") @@ -295,7 +297,7 @@ def truth_moving(x_0, y_0, 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)