🖥️ Choosing a timestep and integration method#
In Parcels, we can simulate virtual particles through time. Kernels which compute a process as a function of time, such as advection, can do so by numerically integrating the process over each timestep dt. In this guide we explore how to control the numerical accuracy of a Parcels simulation by choosing a timestep dt and an integration method.
Note
This notebook is based on the workshop on writing custom Kernels from the 10-year Parcels event.
Timestep#
The choice of timestep affects the accuracy and duration of our simulation. The required accuracy of the integration ultimately depends on the desired application of the particle trajectories. In general, we want our timestep to be small enough to smoothly resolve the time-evolving process of interest, and not so small that our simulations take forever. Here, we will first look at how to estimate a priori what an appropriate timestep could be, and then at how to test a Parcels simulation to check the accuracy.
A priori estimation of dt#
In this example, we will estimate the appropriate timestep to compute advection of a virtual particle through an ocean current field.
Advection, the change in position \(\mathbf{x}(t) = (x(t), y(t))\) at time \(t\), can be described by the equation:
where \(\mathbf{v}(\mathbf{x},t) = (u(\mathbf{x},t), v(\mathbf{x},t))\) describes the ocean velocity field at position \(\mathbf{x}\) at time \(t\).
To estimate the timescale that we want to resolve, we can think about the scales at which advection varies. Here we use the daily velocity fields at 1/12th degree horizontal resolution from the Copernicus Marine Service. This means that the velocity will vary in time at scales >= 24 hours, and in space at scales >= 1/12th degree.
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import polars as pl
import parcels
import parcels.tutorial
# Load the CopernicusMarine data in the Agulhas region from the example_datasets
ds_fields = parcels.tutorial.open_dataset(
"CopernicusMarine_data_for_Argo_tutorial/data"
)
# Convert to SGRID-compliant dataset and create FieldSet
fields = {"U": ds_fields["uo"], "V": ds_fields["vo"]}
ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)
# Convert the FieldSet to windowed arrays for better performance
fieldset = fieldset.to_windowed_arrays()
# Check field resolution in time and space
print(
f"dt = {np.unique(np.diff(ds_fields.time.values).astype('timedelta64[h]'))} hours"
)
print(f"dy = {np.unique(np.diff(ds_fields.latitude.values))} degrees")
print(f"dx = {np.unique(np.diff(ds_fields.longitude.values))} degrees")
dt = [24] hours
dy = [0.08333206 0.08333588] degrees
dx = [0.08332825 0.08334351] degrees
In our Parcels simulation, we must ensure that:
dt< 24 hoursparticles.dx/particles.dy< 1/12th degree.
We can rewrite the second condition using the scale relation:
to express the condition as a function of dt:
We can use the maximum velocity (in degrees s-1) in our flow field as a conservative estimate for \(U\):
# Use the maximum velocity at the surface to scale dx/dy
U_max_surface = np.nanmax(np.hypot(ds_fields.uo, ds_fields.vo))
print(f"U_max = {U_max_surface:.2f} m s-1")
# convert to degrees s-1 (at lat = 30 deg S)
U_max_surface_deg = U_max_surface / (1852 * 60 * np.cos(np.deg2rad(-30)))
print(f" == {U_max_surface_deg * 1e5:.2f}e-5 degrees s-1")
U_max = 1.65 m s-1
== 1.71e-5 degrees s-1
🖥️ Spherical grids and unit converters
Our displacement occurs in units of longitude and latitde, but our velocity field is in m/s. That’s why we have to convert the units of velocity from m/s to degrees/s. In Parcels, this is done automatically for VectorFields when using a spherical grid. See the Spherical grids and velocity conversion guide guide for more information.
Using U_max_surface_deg, we find a second estimated limit of an appropriate dt:
Test timesteps#
Now that we have an estimate of an appropriate choice for dt, we can test a range of values in Parcels to understand the impact on the accuracy of our trajectories. To do so, we can run simulations using parcels.kernels.AdvectionRK2 for a limited number of particles and quantify the differences.
# Particle locations and initial time
npart = 5 # number of particles to be released
initial_release_lons = 32 * np.ones(npart)
initial_release_lats = np.linspace(-32.5, -32, npart, dtype=np.float32)
initial_release_times = np.repeat(
ds_fields.time.values[0], npart
) # release all particles at the start time of the fieldset
We simulate particles using a range of timesteps that differ by a factor 2-10, starting at (dt < 24 hours). We also keep track of the time it takes to run each simulation:
Warning
dt must be chosen such that \(k\) \(\text{d}t\) fits exactly in the FieldSet dt (24 hours), where \(k\) is an integer.
import time
runtime = np.timedelta64(7, "D") # Total simulation time
dt_choices = [
np.timedelta64(12, "h"),
np.timedelta64(6, "h"),
np.timedelta64(1, "h"),
np.timedelta64(20, "m"),
np.timedelta64(5, "m"),
]
sim_duration = np.zeros(len(dt_choices))
for i, dt in enumerate(dt_choices):
pset = parcels.ParticleSet(
fieldset=fieldset,
pclass=parcels.Particle,
t=initial_release_times,
y=initial_release_lats,
x=initial_release_lons,
)
outputdt = dt
pfile = parcels.ParticleFile(
path=f"AdvectionRK2_dt_{int(dt / np.timedelta64(1, 's'))}s.parquet",
outputdt=outputdt,
)
print(f"Begin simulation for dt = {int(dt / np.timedelta64(1, 's'))} s")
# time and run simulation
start_time = time.time()
pset.execute(
parcels.kernels.AdvectionRK2,
runtime=runtime,
dt=dt,
output_file=pfile,
verbose_progress=False,
)
sim_duration_i = time.time() - start_time
sim_duration[i] = sim_duration_i
print(
f"Simulation duration (dt = {int(dt / np.timedelta64(1, 's'))} s) = {np.round(sim_duration_i, 2)} seconds"
)
Begin simulation for dt = 43200 s
INFO: Output files are stored in AdvectionRK2_dt_43200s.parquet
Simulation duration (dt = 43200 s) = 0.2 seconds
Begin simulation for dt = 21600 s
INFO: Output files are stored in AdvectionRK2_dt_21600s.parquet
Simulation duration (dt = 21600 s) = 0.26 seconds
Begin simulation for dt = 3600 s
INFO: Output files are stored in AdvectionRK2_dt_3600s.parquet
Simulation duration (dt = 3600 s) = 1.0 seconds
Begin simulation for dt = 1200 s
INFO: Output files are stored in AdvectionRK2_dt_1200s.parquet
Simulation duration (dt = 1200 s) = 2.76 seconds
Begin simulation for dt = 300 s
INFO: Output files are stored in AdvectionRK2_dt_300s.parquet
Simulation duration (dt = 300 s) = 10.55 seconds
Let’s start by comparing the trajectories of the different simulations:
# To compare the results visually, we will need distinct colour maps for each timestep
dt_colours = np.linspace(0, 1, len(dt_choices), endpoint=True)
fig = plt.figure()
ax = plt.axes()
temperature = ds_fields.isel(time=0, depth=0).thetao.plot(cmap="Greys")
for j, dt in enumerate(dt_choices):
df = parcels.read_particlefile(
f"AdvectionRK2_dt_{int(dt / np.timedelta64(1, 's'))}s.parquet"
)
for i, traj in enumerate(df.partition_by("particle_id", maintain_order=True)):
ax.plot(
traj["x"],
traj["y"],
alpha=0.75,
color=plt.cm.viridis(dt_colours[j]),
label=f"dt = {dt}" if i == 0 else None,
)
df_start = df.filter(pl.col("t") == df["t"].min())
ax.scatter(df_start["x"], df_start["y"], c="r", marker="s", label="starting locations")
ax.legend()
ax.set_ylim(-32.7, -31.3)
ax.set_xlim(31, 32.4)
plt.show()
We see that as the timestep becomes smaller, the trajectories start to converge. As expected from our horizontal grid scaling, the 2 simulations with dt > 2 hours produce trajectories that end up in different grid cells. The other three simulations (dt = (1 h, 20 min, 5 min)) are more difficult to separate visually.
To quantify the precision of these simulations, we can compute the separation distance between the trajectories with the same starting location. Since we have no exact solution to these trajectories, we use our highest resolution simulation as a benchmark. The separation distance is a common metric to analyse particle trajectories, as described in the Lagrangian Diagnostics project.
def dist_km(lona, lonb, lata, latb):
"""
Function to calculate the distance between 2 points in km
Haversine formula used, which assumes the Earth is a sphere.
source: https://stackoverflow.com/questions/19412462/getting-distance-between-two-points-based-on-latitude-longitude
"""
R = 6371.0 # approximate radius of earth in km
lat1 = np.radians(lata)
lon1 = np.radians(lona)
lat2 = np.radians(latb)
lon2 = np.radians(lonb)
dlon = lon2 - lon1
dlat = lat2 - lat1
a = np.square(np.sin(dlat / 2)) + np.cos(lat1) * np.cos(lat2) * np.square(
np.sin(dlon / 2)
)
c = 2 * np.atan2(np.sqrt(a), np.sqrt(1 - a))
distance = R * c
return distance
dist_end = np.zeros((len(dt_choices) - 1, npart))
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(10, 5))
axs[0].set_title("Separation distance from dt = 5-min trajectory")
axs[0].set_xlabel("Time")
axs[0].tick_params("x", rotation=45)
axs[0].set_yscale("log")
axs[0].set_ylim(1e-2, 1e2)
axs[1].set_title("Linear y-scale")
axs[1].set_xlabel("Time")
axs[1].tick_params("x", rotation=45)
axs[1].set_ylim(0, 50)
# set 5 minute dt as benchmark
df_5min = parcels.read_particlefile("AdvectionRK2_dt_300s.parquet")
for i, dt in enumerate(dt_choices[:-1]):
df = parcels.read_particlefile(
f"AdvectionRK2_dt_{int(dt / np.timedelta64(1, 's'))}s.parquet"
)
# subset 5 minute data to match dt
x_5min_sub = df_5min.filter(pl.col("t").is_in(df["t"].implode()))["x"]
y_5min_sub = df_5min.filter(pl.col("t").is_in(df["t"].implode()))["y"]
# compute separation distance between each particle in km
dist = dist_km(df["x"], x_5min_sub, df["y"], y_5min_sub)
df = df.with_columns(pl.Series("dist", dist))
# plot
for j, traj in enumerate(df.partition_by("particle_id", maintain_order=True)):
axs[0].plot(
traj["t"],
traj["dist"],
alpha=0.75,
color=plt.cm.viridis(dt_colours[i]),
label=f"dt = {dt}" if j == 0 else None,
)
axs[1].plot(
traj["t"],
traj["dist"],
alpha=0.75,
color=plt.cm.viridis(dt_colours[i]),
label=f"dt = {dt}" if j == 0 else None,
)
dist_end[i] = df.filter(pl.col("t") == df["t"].max())["dist"]
axs[0].legend()
axs[0].set_ylabel("Distance (km)")
plt.show()
Warning
As you can see, the separation distance starts to increase rapidly with time. Small errors can grow, and the accuracy of the integration therefore depends strongly on the runtime and flow conditions (see the discussion on flow conditions at the end of this guide).
fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
for i, dt in enumerate(dt_choices[:-1]):
ax[0].scatter(
(dt / np.timedelta64(1, "m")).astype(int),
np.mean(dist_end[i]),
color=plt.cm.viridis(dt_colours[i]),
label=f"dt = {dt}",
)
ax[0].plot(
(dt_choices[:-1] / np.timedelta64(1, "m")).astype(int), np.mean(dist_end, axis=1)
)
ax[0].set_ylabel("Mean separation distance (km)")
ax[0].set_xlabel("dt (minutes)")
ax[0].grid()
ax[1].plot(dt_choices, sim_duration)
for i, dt in enumerate(dt_choices):
ax[1].scatter(
(dt / np.timedelta64(1, "m")).astype(int),
sim_duration[i],
color=plt.cm.viridis(dt_colours[i]),
label=f"dt = {dt}",
)
ax[1].set_ylabel("Simulation Duration (s)")
ax[1].set_xlabel("dt (minutes)")
ax[1].grid()
ax[1].legend()
plt.show()
We can see that in our simulation advecting particles for 7 days, the effect of dt on the precision of our simulation is approximately linear. The precision for a simulation with a timestep of 20 minutes is order of magnitude ~100 m. The effect on the time it takes to run a simulation is not linear in our case however; it increases sharply as we decrease our timestep. This may be optimised using more efficient chunking.
df = pd.DataFrame(
{
"dt": dt_choices,
"Mean separation distance (km)": np.append(
np.round(np.mean(dist_end, axis=1), 2), [None]
),
"Simulation duration (s)": np.round(sim_duration, 2),
}
)
df
| dt | Mean separation distance (km) | Simulation duration (s) | |
|---|---|---|---|
| 0 | 0 days 12:00:00 | 0.63 | 0.20 |
| 1 | 0 days 06:00:00 | 0.36 | 0.26 |
| 2 | 0 days 01:00:00 | 0.01 | 1.00 |
| 3 | 0 days 00:20:00 | 0.01 | 2.76 |
| 4 | 0 days 00:05:00 | None | 10.55 |
Integration schemes#
In addition to our choice of timestep, the method we use to integrate the process of interest can have an effect on the accuracy of our calculations. Parcels comes with a number of built-in advection Kernels. These kernels use different integration methods to compute the change in position from one timestep to the next. The higher-order Runge-Kutta (RK) methods are generally more accurate than the simple Explicit Euler (EE).
import inspect
# To get a sense of what the built-in Euler-forward scheme looks like, let's look at the source code!
print(inspect.getsource(parcels.kernels.AdvectionEE))
def AdvectionEE(particles, fieldset): # pragma: no cover
"""Advection of particles using Explicit Euler (aka Euler Forward) integration."""
(u1, v1) = fieldset.UV[particles]
particles.dx += u1 * particles.dt
particles.dy += v1 * particles.dt
# And the second-order Runge-Kutta (RK2) method
print(inspect.getsource(parcels.kernels._advection.rhs_uv))
print(inspect.getsource(parcels.kernels.AdvectionRK2))
print(inspect.getsource(parcels.kernels._integration.RK2))
def rhs_uv(fieldset, t, z, y, x, particles):
return fieldset.UV[t, z, y, x, particles]
def AdvectionRK2(particles, fieldset): # pragma: no cover
"""Advection of particles using second-order Runge-Kutta integration."""
u, v = RK2(particles, fieldset, rhs_uv)
particles.dx += u * particles.dt
particles.dy += v * particles.dt
def RK2(particles, fieldset, rhs):
z, y, x = particles.z, particles.y, particles.x
fields = rhs(fieldset, particles.t, z, y, x, particles)
_validate_rhs_output(fields, "RK2")
x = particles.x + fields[0] * 0.5 * particles.dt
y = particles.y + fields[1] * 0.5 * particles.dt
if len(fields) > 2:
z = particles.z + fields[2] * 0.5 * particles.dt
t = particles.t + 0.5 * particles.dt
return rhs(fieldset, t, z, y, x, particles)
# And RK4:
print(inspect.getsource(parcels.kernels._advection.rhs_uv))
print(inspect.getsource(parcels.kernels.AdvectionRK4))
print(inspect.getsource(parcels.kernels._integration.RK4))
def rhs_uv(fieldset, t, z, y, x, particles):
return fieldset.UV[t, z, y, x, particles]
def AdvectionRK4(particles, fieldset): # pragma: no cover
"""Advection of particles using fourth-order Runge-Kutta integration."""
u, v = RK4(particles, fieldset, rhs_uv)
particles.dx += u * particles.dt
particles.dy += v * particles.dt
def RK4(particles, fieldset, rhs):
z, y, x = particles.z, particles.y, particles.x
k1 = rhs(fieldset, particles.t, z, y, x, particles)
_validate_rhs_output(k1, "RK4")
x = particles.x + k1[0] * 0.5 * particles.dt
y = particles.y + k1[1] * 0.5 * particles.dt
if len(k1) > 2:
z = particles.z + k1[2] * 0.5 * particles.dt
t = particles.t + 0.5 * particles.dt
k2 = rhs(fieldset, t, z, y, x, particles)
x = particles.x + k2[0] * 0.5 * particles.dt
y = particles.y + k2[1] * 0.5 * particles.dt
if len(k2) > 2:
z = particles.z + k2[2] * 0.5 * particles.dt
k3 = rhs(fieldset, t, z, y, x, particles)
x = particles.x + k3[0] * particles.dt
y = particles.y + k3[1] * particles.dt
if len(k3) > 2:
z = particles.z + k3[2] * particles.dt
t = particles.t + particles.dt
k4 = rhs(fieldset, t, z, y, x, particles)
if len(k4) == 2:
return ((k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0, (k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0)
if len(k4) == 3:
return (
(k1[0] + 2 * k2[0] + 2 * k3[0] + k4[0]) / 6.0,
(k1[1] + 2 * k2[1] + 2 * k3[1] + k4[1]) / 6.0,
(k1[2] + 2 * k2[2] + 2 * k3[2] + k4[2]) / 6.0,
)
The higher-order methods use weighted intermediate steps in time and space to obtain a more accurate estimate of dy and dx for a given timestep.
Now let’s see how the different advection schemes affect our simulation. As before, there is no exact solution to the simulated trajectories, so we will use the highest order method as the benchmark.
advection_schemes = [
parcels.kernels.AdvectionEE,
parcels.kernels.AdvectionRK2,
parcels.kernels.AdvectionRK4,
]
sim_duration = np.zeros((len(advection_schemes), len(dt_choices)))
for i, advection_scheme in enumerate(advection_schemes):
for j, dt in enumerate(dt_choices):
pset = parcels.ParticleSet(
fieldset=fieldset,
pclass=parcels.Particle,
t=initial_release_times,
y=initial_release_lats,
x=initial_release_lons,
)
outputdt = dt
pfile = parcels.ParticleFile(
path=f"KernelCompare_{advection_scheme.__name__}_dt_{int(dt / np.timedelta64(1, 's'))}s.parquet",
outputdt=outputdt,
)
print(
f"Begin simulation for (scheme: {advection_scheme.__name__}, dt={int(dt / np.timedelta64(1, 's'))} s)"
)
start_time = time.time()
pset.execute(
advection_scheme,
runtime=runtime,
dt=dt,
output_file=pfile,
verbose_progress=False,
)
sim_duration_ij = time.time() - start_time
sim_duration[i, j] = sim_duration_ij
print(
f"Simulation duration (dt = {int(dt / np.timedelta64(1, 's'))} s) = {np.round(sim_duration_ij, 2)} seconds"
)
scheme_colours = np.linspace(0, 1, len(advection_schemes), endpoint=True)
# Now let's compare different advection schemes with the same timestep
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
for i, dt in enumerate(dt_choices):
m = i // 3
n = i % 3
axs[m, n].set_title(f"dt = {dt}")
axs[m, n].set_xlabel("Longitude")
for j, advection_scheme in enumerate(advection_schemes):
df = parcels.read_particlefile(
f"KernelCompare_{advection_scheme.__name__}_dt_{int(dt / np.timedelta64(1, 's'))}s.parquet"
)
for i, traj in enumerate(df.partition_by("particle_id", maintain_order=True)):
axs[m, n].plot(
traj["x"],
traj["y"],
alpha=0.75,
color=plt.cm.viridis(scheme_colours[j]),
label=f"{advection_scheme.__name__}" if i == 0 else None,
)
df_start = df.filter(pl.col("t") == df["t"].min())
axs[m, n].scatter(
df_start["x"], df_start["y"], c="r", marker="s", label="starting locations"
)
axs[m, n].grid()
axs[-1, -1].axis("off")
axs[0, 0].legend()
axs[0, 0].set_ylabel("Latitude")
axs[1, 0].set_ylabel("Latitude")
plt.show()
Clearly, for longer timesteps, the RK2 and RK4 schemes perform better. However, if the timestep is appropriate, as we have determined in the previous section, then the Explicit Euler scheme does not perform notably different.
dist_end = np.zeros((len(advection_schemes) - 1, len(dt_choices), npart))
fig, axs = plt.subplots(nrows=2, ncols=3, figsize=(15, 10))
for i, dt in enumerate(dt_choices):
m = i // 3
n = i % 3
axs[m, n].set_title(f"dt = {dt}")
axs[m, n].set_xlabel("Time")
axs[m, n].tick_params("x", rotation=45)
axs[m, n].set_yscale("log")
axs[m, n].set_ylim(1e-4, 1e1)
df_RK4 = parcels.read_particlefile(
f"KernelCompare_AdvectionRK4_dt_{int(dt / np.timedelta64(1, 's'))}s.parquet"
)
for j, advection_scheme in enumerate(advection_schemes[:-1]):
df = parcels.read_particlefile(
f"KernelCompare_{advection_scheme.__name__}_dt_{int(dt / np.timedelta64(1, 's'))}s.parquet"
)
dist = dist_km(df["x"], df_RK4["x"], df["y"], df_RK4["y"])
df = df.with_columns(pl.Series("dist", dist))
for k, traj in enumerate(df.partition_by("particle_id", maintain_order=True)):
axs[m, n].plot(
traj["t"],
traj["dist"],
alpha=0.75,
color=plt.cm.viridis(scheme_colours[j]),
label=f"|{advection_scheme.__name__} - AdvectionRK4|"
if k == 0
else None,
)
dist_end[j, i] = df.filter(pl.col("t") == df["t"].max())["dist"]
axs[m, n].grid()
axs[-1, -1].axis("off")
axs[0, 0].legend()
axs[0, 0].set_ylabel("Distance (km)")
axs[1, 0].set_ylabel("Distance (km)")
plt.tight_layout()
plt.show()
By quantifying the precision of the integration methods, we can see that for a given timestep the Runge-Kutta methods perform orders of magnitude better than the Explicit Euler method. In this example, the error associated with the selected integration methods is smaller than that of the range of timesteps.
fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 5))
axs[0].plot(
(dt_choices / np.timedelta64(1, "m")).astype(int),
np.nanmean(dist_end[0], axis=1),
"-o",
color=plt.cm.viridis(scheme_colours[0]),
label="AdvectionEE - AdvectionRK4",
)
axs[0].plot(
(dt_choices / np.timedelta64(1, "m")).astype(int),
np.nanmean(dist_end[1], axis=1),
"-o",
color=plt.cm.viridis(scheme_colours[1]),
label="AdvectionRK2 - AdvectionRK4",
)
axs[0].set_ylabel("Mean separation distance (km)")
axs[0].set_yscale("log")
axs[0].set_xlabel("dt (minutes)")
axs[0].legend()
axs[0].grid()
for i, advection_scheme in enumerate(advection_schemes):
axs[1].plot(
(dt_choices / np.timedelta64(1, "m")).astype(int),
sim_duration[i],
"-o",
color=plt.cm.viridis(scheme_colours[i]),
label=advection_scheme.__name__,
)
axs[1].set_ylabel("Simulation Duration (s)")
axs[1].set_yscale("log")
axs[1].set_xlabel("dt (minutes)")
axs[1].grid()
axs[1].legend()
plt.show()
In this last figure, we see that the improvement of RK2 by orders of magnitude with respect to EE comes at a small computational cost. Since the RK2 and RK4 methods are practically indistinguishable, using the RK2 method with a timestep of 1 hour or 20 minutes, depending on the application, would be an appropriate choice for this simulation.
Flow conditions#
The limitation to numerical accuracy of the integration introduces a positional difference of a certain order of magnitude (\(\delta\)), as seen in the figures above. After this initial deviation, the time evolution of the separation distance is a feature of the specific flow field. As with the Lorenz equations, two initially close particles often separate over time. The rate of separation from a finite perturbation (\(\delta\)) over a given timeperiod can be characterized by the Finite Size Lyapunov Exponent.
Finite Size Lyapunov Exponent references
Lagrangian Diagnostics example to calculate the related Finite Time Lyapunov Exponent
We can see that the separation rate of particles with more precise simulations does not depend mainly on the initial perturbation, but follows a rather specific time-evolution, as the particles enter different flow conditions. The final accuracy of your simulation may therefore strongly depend on your flow conditions and runtime.
🖥️ Lorenz equations in Parcels kernels
Check out this notebook to learn how to model the Lorenz equations with Parcels!
Writing custom Kernels with higher-order integrators#
Similarly to the built-in AdvectionRK2 Kernel, you can also write your own custom Kernels that use higher-order integrators. To do so, you can use the RK2 and RK4 integrators from parcels.integrators in your Kernel definition.
The ‘trick’ is to write a right-hand-side function (rhs) that computes the change in position at a given time and location, and then use this function in the RK2 or RK4 integrator to compute the change in position over a given timestep.
For example, let’s say we want to write a custom Kernel that computes advection and a simple windage effect. We can write the right-hand-side function as follows
def rhs(fieldset, t, z, y, x, particles):
uo, vo = fieldset.UV[t, z, y, x, particles]
uw, vw = fieldset.UVw[t, z, y, x, particles]
return (uo + fieldset.wind_coeff * (uw - uo), vo + fieldset.wind_coeff * (vw - vo))
def WindageRK2(particles, fieldset):
u, v = parcels.kernels._integration.RK2(particles, fieldset, rhs)
particles.dx += u * particles.dt
particles.dy += v * particles.dt
Now, if we also load in the wind data for the same region with the code below, we can execute this Kernel to compute the trajectories of particles that are advected by the ocean currents and affected by the wind. The solution will be second-order accurate in time, since we are using the RK2 integrator; but it is easy to swap in the RK4 integrator if you want to increase the accuracy even more.
ds_w = parcels.tutorial.open_dataset(
"CopernicusMarine_data_for_Argo_tutorial/data_wind"
)
ds_w = ds_w.expand_dims({"depth": [0]})
ds_w.load() # load the dataset into memory
fields = {"Uw": ds_w["eastward_wind"], "Vw": ds_w["northward_wind"]}
ds_wind = parcels.convert.copernicusmarine_to_sgrid(fields=fields)
fieldset_w = parcels.FieldSet.from_sgrid_conventions(
ds_wind, vector_fields={"UVw": ("Uw", "Vw")}
)
fieldset += fieldset_w
fieldset.describe()
| Name | Type | Grid number | Interp method / value | Parcels backend |
|:-------|:------------|--------------:|:------------------------|:------------------|
| U | Field | 0 | XLinear(...) | WindowedArray |
| V | Field | 0 | XLinear(...) | WindowedArray |
| UV | VectorField | 0 | XLinear_Velocity(...) | - |
| Uw | Field | 1 | XLinear(...) | NumPy |
| Vw | Field | 1 | XLinear(...) | NumPy |
| UVw | VectorField | 1 | XLinear_Velocity(...) | - |
mesh: SphericalMesh(radius=6366707.019493707)
time interval: (np.datetime64('2024-01-01T00:00:00.000000000'), np.datetime64('2024-02-01T00:00:00.000000000'))
fieldset.wind_coeff = 0.01 # windage coefficient (1% of the relative wind speed)
pset = parcels.ParticleSet(
fieldset=fieldset,
pclass=parcels.Particle,
t=initial_release_times,
z=np.ones_like(initial_release_times),
y=initial_release_lats,
x=initial_release_lons,
)
pfile = parcels.ParticleFile(
path="WindRK2.parquet",
outputdt=outputdt,
mode="w",
)
pset.execute(
WindageRK2,
runtime=runtime,
dt=np.timedelta64(1, "h"),
output_file=pfile,
verbose_progress=False,
)
INFO: Output files are stored in WindRK2.parquet
As you can see, the particles move very differently when we include the windage effect!
df = parcels.read_particlefile(pfile.path)
for traj in df.partition_by("particle_id"):
plt.plot(traj["x"], traj["y"], alpha=0.75)
plt.show()
Summary#
If you want to test the accuracy and efficiency of your own simulation, here is a brief list of things to consider:
Use the temporal and spatial resolution of the input to estimate the timescales you need to resolve.
dtmust fit into the FieldSet dt by a factor of an integer.Run the simulation for the full runtime. Many errors grow over time so to quantify these, we must take into account the full length of the simulation.
Consider which step introduces the largest error. If we compute advection and diffusion, increasing the accuracy of one process much more than the other will not improve the results.