Skip to content

Commit 0fefbac

Browse files
Fixing AdvectionRK4 to only run for particles within time_interval
1 parent 12108cf commit 0fefbac

File tree

1 file changed

+19
-10
lines changed

1 file changed

+19
-10
lines changed

src/parcels/kernels/advection.py

Lines changed: 19 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -43,16 +43,25 @@ def AdvectionRK2_3D(particles, fieldset): # pragma: no cover
4343

4444
def AdvectionRK4(particles, fieldset): # pragma: no cover
4545
"""Advection of particles using fourth-order Runge-Kutta integration."""
46-
dt = particles.dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds
47-
(u1, v1) = fieldset.UV[particles]
48-
lon1, lat1 = (particles.lon + u1 * 0.5 * dt, particles.lat + v1 * 0.5 * dt)
49-
(u2, v2) = fieldset.UV[particles.time + 0.5 * particles.dt, particles.z, lat1, lon1, particles]
50-
lon2, lat2 = (particles.lon + u2 * 0.5 * dt, particles.lat + v2 * 0.5 * dt)
51-
(u3, v3) = fieldset.UV[particles.time + 0.5 * particles.dt, particles.z, lat2, lon2, particles]
52-
lon3, lat3 = (particles.lon + u3 * dt, particles.lat + v3 * dt)
53-
(u4, v4) = fieldset.UV[particles.time + particles.dt, particles.z, lat3, lon3, particles]
54-
particles.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt
55-
particles.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt
46+
47+
def Update_locations(p):
48+
dt = p.dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds
49+
(u1, v1) = fieldset.UV[p]
50+
lon1, lat1 = (p.lon + u1 * 0.5 * dt, p.lat + v1 * 0.5 * dt)
51+
(u2, v2) = fieldset.UV[p.time + 0.5 * p.dt, p.z, lat1, lon1, p]
52+
lon2, lat2 = (p.lon + u2 * 0.5 * dt, p.lat + v2 * 0.5 * dt)
53+
(u3, v3) = fieldset.UV[p.time + 0.5 * p.dt, p.z, lat2, lon2, p]
54+
lon3, lat3 = (p.lon + u3 * dt, p.lat + v3 * dt)
55+
(u4, v4) = fieldset.UV[p.time + p.dt, p.z, lat3, lon3, p]
56+
p.dlon += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt
57+
p.dlat += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt
58+
59+
Update_locations(
60+
particles[
61+
(particles.time + particles.dt <= fieldset.time_interval.right)
62+
& (particles.time + particles.dt >= fieldset.time_interval.left)
63+
]
64+
)
5665

5766

5867
def AdvectionRK4_3D(particles, fieldset): # pragma: no cover

0 commit comments

Comments
 (0)