Skip to content

Adding in the Runge-Kutta 2 advection schemes #2315

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
merged 7 commits into from
Nov 3, 2025
Merged

Adding in the Runge-Kutta 2 advection schemes #2315

Show file tree
Hide file tree
Changes from 5 commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
179d306
Implementing RK2 and RK2_3D advection schemes
michaeldenes Oct 8, 2025
202dc07
Adding the two new schemes to the list of schemes
michaeldenes Oct 8, 2025
3a95f67
Updating init to include the new kernels
michaeldenes Oct 8, 2025
4be9a2b
Adding RK2 and RK2_3D into (some of) the unit tests.
michaeldenes Oct 8, 2025
c107a5d
Updating tolerances for RK2 unit tests
erikvansebille Oct 8, 2025
a39a72c
Removing distionary from method to kernels
erikvansebille Nov 3, 2025
29b8d37
Merge branch 'v4-dev' into runge-kutta2
erikvansebille Nov 3, 2025
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
4 changes: 4 additions & 0 deletions src/parcels/kernels/__init__.py
View file
Open in desktop
Copy link
Member

@erikvansebille erikvansebille Oct 8, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@VeckoTheGecko, do we really need to add new Kernels at three locations? (two locations here in __init__.py and also at the top of advection.py). Can there not be less duplication?

Copy link
Contributor

@VeckoTheGecko VeckoTheGecko Oct 9, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We can make advection.py private (i.e., changing it to be _advection.py). Then we wouldn't need to define __all__ in advection.py. As it stands though, we need to have it in __all__ for advection.py and __init__.py

Original file line number Diff line number Diff line change
@@ -1,6 +1,8 @@
from .advection import (
AdvectionAnalytical,
AdvectionEE,
AdvectionRK2,
AdvectionRK2_3D,
AdvectionRK4,
AdvectionRK4_3D,
AdvectionRK4_3D_CROCO,
Expand All @@ -21,6 +23,8 @@
# advection
"AdvectionAnalytical",
"AdvectionEE",
"AdvectionRK2",
"AdvectionRK2_3D",
"AdvectionRK4_3D_CROCO",
"AdvectionRK4_3D",
"AdvectionRK4",
Expand Down
25 changes: 25 additions & 0 deletions src/parcels/kernels/advection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,13 +9,38 @@
__all__ = [
"AdvectionAnalytical",
"AdvectionEE",
"AdvectionRK2",
"AdvectionRK2_3D",
"AdvectionRK4",
"AdvectionRK4_3D",
"AdvectionRK4_3D_CROCO",
"AdvectionRK45",
]


def AdvectionRK2(particles, fieldset): # pragma: no cover
"""Advection of particles using second-order Runge-Kutta integration."""
dt = particles.dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds
(u1, v1) = fieldset.UV[particles]
lon1, lat1 = (particles.lon + u1 * 0.5 * dt, particles.lat + v1 * 0.5 * dt)
(u2, v2) = fieldset.UV[particles.time + 0.5 * particles.dt, particles.z, lat1, lon1, particles]
particles.dlon += u2 * dt
particles.dlat += v2 * dt


def AdvectionRK2_3D(particles, fieldset): # pragma: no cover
"""Advection of particles using second-order Runge-Kutta integration including vertical velocity."""
dt = particles.dt / np.timedelta64(1, "s")
(u1, v1, w1) = fieldset.UVW[particles]
lon1 = particles.lon + u1 * 0.5 * dt
lat1 = particles.lat + v1 * 0.5 * dt
z1 = particles.z + w1 * 0.5 * dt
(u2, v2, w2) = fieldset.UVW[particles.time + 0.5 * particles.dt, z1, lat1, lon1, particles]
particles.dlon += u2 * dt
particles.dlat += v2 * dt
particles.dz += w2 * dt


def AdvectionRK4(particles, fieldset): # pragma: no cover
"""Advection of particles using fourth-order Runge-Kutta integration."""
dt = particles.dt / np.timedelta64(1, "s") # TODO: improve API for converting dt to seconds
Expand Down
11 changes: 10 additions & 1 deletion tests/test_advection.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -18,6 +18,8 @@
AdvectionDiffusionEM,
AdvectionDiffusionM1,
AdvectionEE,
AdvectionRK2,
AdvectionRK2_3D,
AdvectionRK4,
AdvectionRK4_3D,
AdvectionRK45,
Expand All @@ -26,6 +28,8 @@

kernel = {
"EE": AdvectionEE,
"RK2": AdvectionRK2,
"RK2_3D": AdvectionRK2_3D,
"RK4": AdvectionRK4,
"RK4_3D": AdvectionRK4_3D,
"RK45": AdvectionRK45,
Expand Down Expand Up @@ -261,6 +265,8 @@ def test_radialrotation(npart=10):
("EE", 1e-2),
("AdvDiffEM", 1e-2),
("AdvDiffM1", 1e-2),
("RK2", 6e-5),
("RK2_3D", 6e-5),
("RK4", 1e-5),
("RK4_3D", 1e-5),
("RK45", 1e-4),
Expand All @@ -271,7 +277,7 @@ def test_moving_eddy(method, rtol):
grid = XGrid.from_dataset(ds)
U = Field("U", ds["U"], grid, interp_method=XLinear)
V = Field("V", ds["V"], grid, interp_method=XLinear)
if method == "RK4_3D":
if method in ["RK2_3D", "RK4_3D"]:
# Using W to test 3D advection (assuming same velocity as V)
W = Field("W", ds["V"], grid, interp_method=XLinear)
UVW = VectorField("UVW", U, V, W)
Expand Down Expand Up @@ -312,6 +318,7 @@ def truth_moving(x_0, y_0, t):
"method, rtol",
[
("EE", 1e-1),
("RK2", 3e-3),
("RK4", 1e-5),
("RK45", 1e-4),
],
Expand Down Expand Up @@ -356,6 +363,7 @@ def truth_moving(x_0, y_0, t):
@pytest.mark.parametrize(
"method, rtol",
[
("RK2", 0.1),
("RK4", 0.1),
("RK45", 0.1),
],
Expand Down Expand Up @@ -396,6 +404,7 @@ def UpdateP(particles, fieldset): # pragma: no cover
@pytest.mark.parametrize(
"method, rtol",
[
("RK2", 2e-2),
("RK4", 5e-3),
("RK45", 1e-4),
],
Expand Down