Skip to content

Commit e9d136a

Browse files
Fixing backward time in RK45
1 parent 3db13dd commit e9d136a

1 file changed

Lines changed: 5 additions & 4 deletions

File tree

src/parcels/kernels/advection.py

Lines changed: 5 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -137,6 +137,7 @@ def AdvectionRK45(particles, fieldset): # pragma: no cover
137137
and doubled if error is smaller than 1/10th of tolerance.
138138
"""
139139
dt = particles.dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds
140+
sign_dt = np.sign(dt)
140141

141142
c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]
142143
A = [
@@ -189,17 +190,17 @@ def AdvectionRK45(particles, fieldset): # pragma: no cover
189190
)
190191
particles.dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt)
191192
particles.dt = np.where(
192-
particles.dt > fieldset.RK45_max_dt * np.timedelta64(1, "s"),
193-
fieldset.RK45_max_dt * np.timedelta64(1, "s"),
193+
np.abs(particles.next_dt) > np.abs(fieldset.RK45_max_dt * np.timedelta64(1, "s")),
194+
fieldset.RK45_max_dt * np.timedelta64(1, "s") * sign_dt,
194195
particles.dt,
195196
)
196197
particles.state = np.where(good_particles, StatusCode.Evaluate, particles.state)
197198

198199
repeat_particles = np.invert(good_particles)
199200
particles.dt = np.where(repeat_particles, particles.dt / 2, particles.dt)
200201
particles.dt = np.where(
201-
particles.dt < fieldset.RK45_min_dt * np.timedelta64(1, "s"),
202-
fieldset.RK45_min_dt * np.timedelta64(1, "s"),
202+
np.abs(particles.dt) < np.abs(fieldset.RK45_min_dt * np.timedelta64(1, "s")),
203+
fieldset.RK45_min_dt * np.timedelta64(1, "s") * sign_dt,
203204
particles.dt,
204205
)
205206
particles.state = np.where(repeat_particles, StatusCode.Repeat, particles.state)

0 commit comments

Comments
 (0)