Skip to content

Remove _constrain_dt_to_within_time_interval function #2754

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
erikvansebille merged 3 commits into from
Jul 16, 2026
+76 −99
Merged

Remove _constrain_dt_to_within_time_interval function #2754

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
5db1fbf
Removing the _constrain_dt_to_within_time_interval helper function
erikvansebille Jul 15, 2026
31dd2bd
Using particles.dt directly in sigmagrids
erikvansebille Jul 15, 2026
14547ae
Merge branch 'main' into remove_constrain_dt_to_within_time_interval
erikvansebille Jul 16, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
146 changes: 63 additions & 83 deletions src/parcels/kernels/_advection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -17,76 +17,62 @@
]


def _constrain_dt_to_within_time_interval(time_interval, time, dt):
"""Helper function to make sure dt does not go outside time_interval.

This is especially relevant for higher-order RK methods (RK2, RK4, RK45),
which require interpolations at time + dt. If time is at the edges of the
time_interval (typically the last integration step), such an operation would
lead to an OutofTimeError.
"""
if time_interval:
dt = np.where(time + dt <= time_interval.time_length_as_flt, dt, time_interval.time_length_as_flt - time)
dt = np.where(time + dt >= 0, dt, time)
return dt


def AdvectionRK2(particles, fieldset): # pragma: no cover
"""Advection of particles using second-order Runge-Kutta integration."""
dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
(u1, v1) = fieldset.UV[particles]
x1, y1 = (particles.x + u1 * 0.5 * dt, particles.y + v1 * 0.5 * dt)
(u2, v2) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y1, x1, particles]
particles.dx += u2 * dt
particles.dy += v2 * dt
x1 = particles.x + u1 * 0.5 * particles.dt
y1 = particles.y + v1 * 0.5 * particles.dt
(u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles]
particles.dx += u2 * particles.dt
particles.dy += v2 * particles.dt


def AdvectionRK2_3D(particles, fieldset): # pragma: no cover
"""Advection of particles using second-order Runge-Kutta integration including vertical velocity."""
dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
(u1, v1, w1) = fieldset.UVW[particles]
x1 = particles.x + u1 * 0.5 * dt
y1 = particles.y + v1 * 0.5 * dt
z1 = particles.z + w1 * 0.5 * dt
(u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * dt, z1, y1, x1, particles]
particles.dx += u2 * dt
particles.dy += v2 * dt
particles.dz += w2 * dt
x1 = particles.x + u1 * 0.5 * particles.dt
y1 = particles.y + v1 * 0.5 * particles.dt
z1 = particles.z + w1 * 0.5 * particles.dt
(u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * particles.dt, z1, y1, x1, particles]
particles.dx += u2 * particles.dt
particles.dy += v2 * particles.dt
particles.dz += w2 * particles.dt


def AdvectionRK4(particles, fieldset): # pragma: no cover
"""Advection of particles using fourth-order Runge-Kutta integration."""
dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
(u1, v1) = fieldset.UV[particles]
x1, y1 = (particles.x + u1 * 0.5 * dt, particles.y + v1 * 0.5 * dt)
(u2, v2) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y1, x1, particles]
x2, y2 = (particles.x + u2 * 0.5 * dt, particles.y + v2 * 0.5 * dt)
(u3, v3) = fieldset.UV[particles.t + 0.5 * dt, particles.z, y2, x2, particles]
x3, y3 = (particles.x + u3 * dt, particles.y + v3 * dt)
(u4, v4) = fieldset.UV[particles.t + dt, particles.z, y3, x3, particles]
particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * dt
particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * dt
x1 = particles.x + u1 * 0.5 * particles.dt
y1 = particles.y + v1 * 0.5 * particles.dt
(u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles]
x2 = particles.x + u2 * 0.5 * particles.dt
y2 = particles.y + v2 * 0.5 * particles.dt
(u3, v3) = fieldset.UV[particles.t + 0.5 * particles.dt, particles.z, y2, x2, particles]
x3 = particles.x + u3 * particles.dt
y3 = particles.y + v3 * particles.dt
(u4, v4) = fieldset.UV[particles.t + particles.dt, particles.z, y3, x3, particles]
particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6.0 * particles.dt
particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6.0 * particles.dt


def AdvectionRK4_3D(particles, fieldset): # pragma: no cover
"""Advection of particles using fourth-order Runge-Kutta integration including vertical velocity."""
dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
(u1, v1, w1) = fieldset.UVW[particles]
x1 = particles.x + u1 * 0.5 * dt
y1 = particles.y + v1 * 0.5 * dt
z1 = particles.z + w1 * 0.5 * dt
(u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * dt, z1, y1, x1, particles]
x2 = particles.x + u2 * 0.5 * dt
y2 = particles.y + v2 * 0.5 * dt
z2 = particles.z + w2 * 0.5 * dt
(u3, v3, w3) = fieldset.UVW[particles.t + 0.5 * dt, z2, y2, x2, particles]
x3 = particles.x + u3 * dt
y3 = particles.y + v3 * dt
z3 = particles.z + w3 * dt
(u4, v4, w4) = fieldset.UVW[particles.t + dt, z3, y3, x3, particles]
particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * dt
particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * dt
particles.dz += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * dt
x1 = particles.x + u1 * 0.5 * particles.dt
y1 = particles.y + v1 * 0.5 * particles.dt
z1 = particles.z + w1 * 0.5 * particles.dt
(u2, v2, w2) = fieldset.UVW[particles.t + 0.5 * particles.dt, z1, y1, x1, particles]
x2 = particles.x + u2 * 0.5 * particles.dt
y2 = particles.y + v2 * 0.5 * particles.dt
z2 = particles.z + w2 * 0.5 * particles.dt
(u3, v3, w3) = fieldset.UVW[particles.t + 0.5 * particles.dt, z2, y2, x2, particles]
x3 = particles.x + u3 * particles.dt
y3 = particles.y + v3 * particles.dt
z3 = particles.z + w3 * particles.dt
(u4, v4, w4) = fieldset.UVW[particles.t + particles.dt, z3, y3, x3, particles]
particles.dx += (u1 + 2 * u2 + 2 * u3 + u4) / 6 * particles.dt
particles.dy += (v1 + 2 * v2 + 2 * v3 + v4) / 6 * particles.dt
particles.dz += (w1 + 2 * w2 + 2 * w3 + w4) / 6 * particles.dt


def AdvectionEE(particles, fieldset): # pragma: no cover
Expand All @@ -105,8 +91,7 @@ def AdvectionRK45(particles, fieldset): # 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 = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
sign_dt = np.sign(dt)
sign_dt = np.sign(particles.dt)

c = [1.0 / 4.0, 3.0 / 8.0, 12.0 / 13.0, 1.0, 1.0 / 2.0]
A = [
Expand All @@ -120,42 +105,37 @@ def AdvectionRK45(particles, fieldset): # pragma: no cover
b5 = [16.0 / 135.0, 0.0, 6656.0 / 12825.0, 28561.0 / 56430.0, -9.0 / 50.0, 2.0 / 55.0]

(u1, v1) = fieldset.UV[particles]
x1, y1 = (particles.x + u1 * A[0][0] * dt, particles.y + v1 * A[0][0] * dt)
(u2, v2) = fieldset.UV[particles.t + c[0] * dt, particles.z, y1, x1, particles]
x2, y2 = (
particles.x + (u1 * A[1][0] + u2 * A[1][1]) * dt,
particles.y + (v1 * A[1][0] + v2 * A[1][1]) * dt,
)
(u3, v3) = fieldset.UV[particles.t + c[1] * dt, particles.z, y2, x2, particles]
x3, y3 = (
particles.x + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * dt,
particles.y + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * dt,
)
(u4, v4) = fieldset.UV[particles.t + c[2] * dt, particles.z, y3, x3, particles]
x4, y4 = (
particles.x + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * dt,
particles.y + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * dt,
)
(u5, v5) = fieldset.UV[particles.t + c[3] * dt, particles.z, y4, x4, particles]
x5, y5 = (
particles.x + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * dt,
particles.y + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * dt,
)
(u6, v6) = fieldset.UV[particles.t + c[4] * dt, particles.z, y5, x5, particles]

x_4th = (u1 * b4[0] + u2 * b4[1] + u3 * b4[2] + u4 * b4[3] + u5 * b4[4]) * dt
y_4th = (v1 * b4[0] + v2 * b4[1] + v3 * b4[2] + v4 * b4[3] + v5 * b4[4]) * dt
x_5th = (u1 * b5[0] + u2 * b5[1] + u3 * b5[2] + u4 * b5[3] + u5 * b5[4] + u6 * b5[5]) * dt
y_5th = (v1 * b5[0] + v2 * b5[1] + v3 * b5[2] + v4 * b5[3] + v5 * b5[4] + v6 * b5[5]) * dt
x1 = particles.x + u1 * A[0][0] * particles.dt
y1 = particles.y + v1 * A[0][0] * particles.dt
(u2, v2) = fieldset.UV[particles.t + c[0] * particles.dt, particles.z, y1, x1, particles]
x2 = particles.x + (u1 * A[1][0] + u2 * A[1][1]) * particles.dt
y2 = particles.y + (v1 * A[1][0] + v2 * A[1][1]) * particles.dt
(u3, v3) = fieldset.UV[particles.t + c[1] * particles.dt, particles.z, y2, x2, particles]
x3 = particles.x + (u1 * A[2][0] + u2 * A[2][1] + u3 * A[2][2]) * particles.dt
y3 = particles.y + (v1 * A[2][0] + v2 * A[2][1] + v3 * A[2][2]) * particles.dt
(u4, v4) = fieldset.UV[particles.t + c[2] * particles.dt, particles.z, y3, x3, particles]
x4 = particles.x + (u1 * A[3][0] + u2 * A[3][1] + u3 * A[3][2] + u4 * A[3][3]) * particles.dt
y4 = particles.y + (v1 * A[3][0] + v2 * A[3][1] + v3 * A[3][2] + v4 * A[3][3]) * particles.dt
(u5, v5) = fieldset.UV[particles.t + c[3] * particles.dt, particles.z, y4, x4, particles]
x5 = particles.x + (u1 * A[4][0] + u2 * A[4][1] + u3 * A[4][2] + u4 * A[4][3] + u5 * A[4][4]) * particles.dt
y5 = particles.y + (v1 * A[4][0] + v2 * A[4][1] + v3 * A[4][2] + v4 * A[4][3] + v5 * A[4][4]) * particles.dt
(u6, v6) = fieldset.UV[particles.t + c[4] * particles.dt, particles.z, y5, x5, particles]

x_4th = (u1 * b4[0] + u2 * b4[1] + u3 * b4[2] + u4 * b4[3] + u5 * b4[4]) * particles.dt
y_4th = (v1 * b4[0] + v2 * b4[1] + v3 * b4[2] + v4 * b4[3] + v5 * b4[4]) * particles.dt
x_5th = (u1 * b5[0] + u2 * b5[1] + u3 * b5[2] + u4 * b5[3] + u5 * b5[4] + u6 * b5[5]) * particles.dt
y_5th = (v1 * b5[0] + v2 * b5[1] + v3 * b5[2] + v4 * b5[3] + v5 * b5[4] + v6 * b5[5]) * particles.dt

kappa = np.sqrt(np.pow(x_5th - x_4th, 2) + np.pow(y_5th - y_4th, 2))

good_particles = (kappa <= fieldset.RK45_tol) | (np.fabs(dt) <= np.fabs(fieldset.RK45_min_dt))
good_particles = (kappa <= fieldset.RK45_tol) | (np.fabs(particles.dt) <= np.fabs(fieldset.RK45_min_dt))
particles.dx += np.where(good_particles, x_5th, 0)
particles.dy += np.where(good_particles, y_5th, 0)

increase_dt_particles = (
good_particles & (kappa <= fieldset.RK45_tol / 10) & (np.fabs(dt * 2) <= np.fabs(fieldset.RK45_max_dt))
good_particles
& (kappa <= fieldset.RK45_tol / 10)
& (np.fabs(particles.dt * 2) <= np.fabs(fieldset.RK45_max_dt))
)
particles.next_dt = np.where(increase_dt_particles, particles.dt * 2, particles.dt)
particles.next_dt = np.where(
Expand Down
29 changes: 13 additions & 16 deletions src/parcels/kernels/_sigmagrids.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,8 +2,6 @@

import numpy as np

from parcels.kernels._advection import _constrain_dt_to_within_time_interval


def convert_z_to_sigma_croco(fieldset, t, z, y, x, particle):
"""Calculate local sigma level of the particles, by linearly interpolating the
Expand Down Expand Up @@ -49,27 +47,26 @@ def AdvectionRK2_3D_CROCO(particles, fieldset): # pragma: no cover
category=RuntimeWarning,
) # Needed because of linear sampling of W with sigma conversion

dt = _constrain_dt_to_within_time_interval(fieldset.time_interval, particles.t, particles.dt)
sigma = particles.z / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x]

sig = convert_z_to_sigma_croco(fieldset, particles.t, particles.z, particles.y, particles.x, particles)
(u1, v1) = fieldset.UV[particles.t, sig, particles.y, particles.x, particles]
w1 = fieldset.W[particles.t, sig, particles.y, particles.x, particles]
w1 *= sigma / fieldset.h[particles.t, np.zeros_like(particles.z), particles.y, particles.x]
x1 = particles.x + u1 * 0.5 * dt
y1 = particles.y + v1 * 0.5 * dt
sig_dep1 = sigma + w1 * 0.5 * dt
x1 = particles.x + u1 * 0.5 * particles.dt
y1 = particles.y + v1 * 0.5 * particles.dt
sig_dep1 = sigma + w1 * 0.5 * particles.dt
dep1 = sig_dep1 * fieldset.h[particles.t, np.zeros_like(particles.z), y1, x1]

sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * dt, dep1, y1, x1, particles)
(u2, v2) = fieldset.UV[particles.t + 0.5 * dt, sig1, y1, x1, particles]
w2 = fieldset.W[particles.t + 0.5 * dt, sig1, y1, x1, particles]
w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y1, x1]
x2 = particles.x + u2 * 0.5 * dt
y2 = particles.y + v2 * 0.5 * dt
sig_dep2 = sigma + w2 * 0.5 * dt
dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * dt, np.zeros_like(particles.z), y2, x2]
sig1 = convert_z_to_sigma_croco(fieldset, particles.t + 0.5 * particles.dt, dep1, y1, x1, particles)
(u2, v2) = fieldset.UV[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles]
w2 = fieldset.W[particles.t + 0.5 * particles.dt, sig1, y1, x1, particles]
w2 *= sig_dep1 / fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y1, x1]
x2 = particles.x + u2 * 0.5 * particles.dt
y2 = particles.y + v2 * 0.5 * particles.dt
sig_dep2 = sigma + w2 * 0.5 * particles.dt
dep2 = sig_dep2 * fieldset.h[particles.t + 0.5 * particles.dt, np.zeros_like(particles.z), y2, x2]

particles.dx += u2 * dt
particles.dy += v2 * dt
particles.dx += u2 * particles.dt
particles.dy += v2 * particles.dt
particles.dz += (dep1 - particles.z) + (dep2 - particles.z)