📖 The sensitivity of grid and velocity interpolation#
The notebook below shows how sensitive trajectories are to the combination of grid and velocity interpolation schemes. We will use the idealised Peninsula flowfield and carefully compare how trajectories on the A-grid and C-grid version of this flow differ. To enhance the differences, we deliberately use a very coarse-resolution grid.
from dataclasses import dataclass
from typing import Literal
import matplotlib.pyplot as plt
import numpy as np
import scipy
from matplotlib.lines import Line2D
import parcels
We start with defining the FieldSet for the Peninsula flow and the grid. Note that on the A-grid we can directly calculate the velocities on the grid points from the streamfunction P. On the C-grid, we need to differentiate the streamfunction P to get the velocities on the staggered U and V grid points.
from parcels._datasets.structured.generated import peninsula_dataset
Next, we define some helper Kernels and the ParticleType.
def SampleP(particles, fieldset):
particles.p = fieldset.P[particles]
def DeleteParticle(particles, fieldset):
any_error = particles.state >= 50 # This captures all Errors
particles[any_error].state = parcels.StatusCode.Delete
pclass = parcels.Particle.add_variable(parcels.Variable("p", np.float32))
Now we run six different experiments, each with a different combination of grid type (“A” or “C”) and velocity interpolation scheme (parcels.interpolators.XLinear_Velocity, parcels.interpolators.XFreeslip, and parcels.interpolators.CGrid_Velocity).
@dataclass
class Experiment:
"""Utility class to consolidate experiment parameters."""
grid_type: Literal["A", "C"]
interp_method: Literal[
"linear", "freeslip", "cgrid_velocity"
] # , "analytical"] TODO add analtical when implemented
@property
def fieldset_interp_method(self):
"""Which FieldSet interpolation method to use for each experiment."""
if self.interp_method == "analytical":
return "cgrid_velocity"
return self.interp_method
@property
def advection_kernel(self):
if self.interp_method == "analytical":
return parcels.kernels.AdvectionAnalytical
return parcels.kernels.AdvectionRK4
@property
def file_name(self):
return f"Trajs_{self.grid_type}_{self.interp_method}.parquet"
@property
def plot_title(self):
if str(self.interp_method) == "analytical":
interp_name = "analytical advection"
else:
interp_name = f"{str(self.interp_method).replace('_Velocity', '').replace('(...)', '')} interpolation"
return f"{self.grid_type}grid & {interp_name}"
exps = [
Experiment("A", parcels.interpolators.XLinear_Velocity()),
Experiment("A", parcels.interpolators.XFreeslip()),
Experiment("C", parcels.interpolators.CGrid_Velocity()),
# Experiment("C", "analytical"),
Experiment("A", parcels.interpolators.CGrid_Velocity()),
Experiment("C", parcels.interpolators.XLinear_Velocity()),
]
dt = np.timedelta64(1000, "s") # setting output and execution timestep to same value
for exp in exps:
ds = peninsula_dataset(xdim=14, ydim=16, grid_type=exp.grid_type)
fieldset = parcels.FieldSet.from_sgrid_conventions(ds)
fieldset.UV.interp_method = exp.fieldset_interp_method
pset = parcels.ParticleSet(
fieldset, pclass=pclass, x=1e3 * np.ones(7), y=np.linspace(1e3, 10e3, 7)
)
outfile = parcels.ParticleFile(exp.file_name, outputdt=dt, mode="w")
pset.execute(
[exp.advection_kernel, SampleP, DeleteParticle],
runtime=np.timedelta64(100_000, "s"),
dt=dt,
output_file=outfile,
verbose_progress=False,
)
INFO: Output files are stored in Trajs_A_XLinear_Velocity(...).parquet
INFO: Output files are stored in Trajs_A_XFreeslip(...).parquet
INFO: Output files are stored in Trajs_C_CGrid_Velocity(...).parquet
INFO: Output files are stored in Trajs_A_CGrid_Velocity(...).parquet
INFO: Output files are stored in Trajs_C_XLinear_Velocity(...).parquet
We carefully compute the landmask (i.e. where the peninsula field has land), which is important for the interpretation of land that is consistent with Parcels (see also the Preventing stuck particles tutorial).
ds = peninsula_dataset(xdim=14, ydim=16, grid_type="C")
landmask = np.ma.masked_values(ds.U.data[:, :], 0)
landmask = landmask.mask.astype("int")
lonplot, latplot = np.meshgrid(ds.lon, ds.lat)
fl = scipy.interpolate.RectBivariateSpline(ds.lat, ds.lon, landmask, kx=1, ky=1)
x = ds.lon[:-1] + np.diff(ds.lon) / 2
y = ds.lat[:-1] + np.diff(ds.lat) / 2
lon_centers, lat_centers = np.meshgrid(x, y)
l_centers = fl(lat_centers[:, 0], lon_centers[0, :])
lmask = np.ma.masked_values(l_centers, 1)
Finally, we plot the trajectories for all six experiments. We see that the trajectories are very sensitive to the combination of grid and velocity interpolation scheme. Note for example that using cgrid_velocity interpolation on an A-grid (fifth panel) leads to particles getting grounded on the land, which is not physical.
fig, axs = plt.subplots(1, len(exps), figsize=(15, 4), sharey="row")
def plot_experiment(ax, exp: Experiment) -> None:
df = parcels.read_particlefile(exp.file_name)
ax.pcolormesh(lonplot, latplot, lmask.mask, cmap="Blues")
ax.scatter(
lonplot,
latplot,
c=landmask,
s=20,
cmap="Blues",
vmin=-0.05,
vmax=0.05,
edgecolors="k",
)
for traj in df.partition_by("particle_id", maintain_order=True):
ax.plot(traj["x"], traj["y"], linewidth=2)
ax.set_title(exp.plot_title)
# Set the same limits for all subplots
ax.set_xlim([lonplot.min(), lonplot.max()])
ax.set_ylim([0, 23e3])
m2km = lambda x, _: f"{x / 1000:.1f}"
ax.xaxis.set_major_formatter(m2km)
ax.yaxis.set_major_formatter(m2km)
ax.set_xlabel("x [km]")
for i, (ax, exp) in enumerate(zip(axs, exps, strict=True)):
if i == 0:
ax.set_ylabel("y [km]")
plot_experiment(ax, exp)
plt.tight_layout()
plt.show()
Because we know that particles must follow lines of constant streamfunction, we can also plot the streamlines of the flow field. When we overlay the “Agrid & linear interpolation” and “Cgrid & cgrid interpolation” trajectories onto these streamlines, we see that they are subtly different.
fig, ax = plt.subplots(1, 1, figsize=(5, 4))
ax.pcolormesh(lonplot, latplot, lmask.mask, cmap="Blues")
ax.scatter(
lonplot,
latplot,
c=landmask,
s=20,
cmap="Blues",
vmin=-0.05,
vmax=0.05,
edgecolors="k",
)
ax.contour(lonplot, latplot, ds.P.data, colors="k", levels=50)
exps_of_interest = filter(
lambda exp: (
exp.plot_title
in ["Agrid & XLinear interpolation", "Cgrid & CGrid interpolation"]
),
exps,
)
handles = []
for i, exp in enumerate(exps_of_interest):
df = parcels.read_particlefile(exp.file_name)
for traj in df.partition_by("particle_id", maintain_order=True):
ax.plot(traj["x"], traj["y"], linewidth=2, color=f"C{i:02d}")
# Manually save label for legend
handles.append(Line2D([0], [0], label=exp.plot_title, color=f"C{i:02d}"))
ax.legend(handles=handles)
ax.set_xlim([lonplot.min(), lonplot.max()])
ax.set_ylim([0, 23e3])
ax.set_ylabel("y [km]")
ax.set_xlabel("x [km]")
m2km = lambda x, _: f"{x / 1000:.1f}"
ax.xaxis.set_major_formatter(m2km)
ax.yaxis.set_major_formatter(m2km)
plt.tight_layout()
plt.show()
To further investigate the differences, we can also plot the change in the streamfunction value sampled by the particle (particle.p) between the final and start location for each experiment. We see that the difference is larger for particles that start farther south and thus get closer to the peninsula. The difference in streamfunction change for the particle that started farthest north is very close to zero for all of the first four experiments; and much larger for the fifth and sixth experiments (where the grid and interpolation scheme do not match up).
fig, ax = plt.subplots(1, 1, figsize=(6, 4))
for i, exp in enumerate(exps):
df = parcels.read_particlefile(exp.file_name)
for j, traj in enumerate(df.partition_by("particle_id", maintain_order=True)):
label = exp.plot_title if j == 0 else None
ax.plot(
traj["y"][0] / 1e3 + i / 10,
traj["p"][-1] - traj["p"][1],
"o",
alpha=0.7,
label=label,
markeredgecolor="k",
color=f"C{i:02d}",
)
ax.legend()
ax.hlines(0, 3, 11, linestyles="--", color="k")
ax.set_xlabel("y-location at start [km]")
ax.set_xlim([3, 11])
ax.set_ylabel("$p_{end} - p_{start}$")
plt.tight_layout()
plt.show()
This clearly shows why it is important to consider the combination of grid and interpolation scheme when running Parcels simulations.