Particle positions in output are lagged by one `dt` relative to time labels

[fluidnumerics-joe]

Parcels version

afb58da

Description

When writing particle trajectories via ParticleFile, the position recorded at each output time corresponds to the particle state one dt earlier than the time label. For a uniform velocity field u = U0 where the analytical trajectory is x(t) = x0 + U0*t, the final output at t = T contains the position at t = T - dt. The error is exactly proportional to dt.

Reproducer

import numpy as np
import parcels
import uxarray as ux
import xarray as xr

# --- Build a minimal Delaunay grid ---
NX = 10
lon, lat = np.meshgrid(
    np.linspace(0, 20, NX, dtype=np.float64),
    np.linspace(0, 20, NX, dtype=np.float64),
)
lon_flat, lat_flat = lon.ravel(), lat.ravel()
mask = (
    np.isclose(lon_flat, 0) | np.isclose(lon_flat, 20)
    | np.isclose(lat_flat, 0) | np.isclose(lat_flat, 20)
)
uxgrid = ux.Grid.from_points(
    (lon_flat, lat_flat),
    method="regional_delaunay",
    boundary_points=np.flatnonzero(mask),
)
uxgrid.attrs["Conventions"] = "UGRID-1.0"

# --- Uniform velocity field on face centers ---
U0 = 0.001  # degrees/s
V0 = 0.0
N_TIMESTEPS = 1
TIME = xr.date_range("2000-01-01", periods=N_TIMESTEPS, freq="1h")
zf = np.array([0.0, 1.0])
zc = np.array([0.5])

U = np.full((N_TIMESTEPS, 1, uxgrid.n_face), U0)
V = np.full((N_TIMESTEPS, 1, uxgrid.n_face), V0)
W = np.zeros((N_TIMESTEPS, 2, uxgrid.n_node))

ds = ux.UxDataset({
    "U": ux.UxDataArray(U, uxgrid=uxgrid, dims=["time", "zc", "n_face"],
             coords=dict(time=(["time"], TIME), zc=(["zc"], zc)),
             attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")),
    "V": ux.UxDataArray(V, uxgrid=uxgrid, dims=["time", "zc", "n_face"],
             coords=dict(time=(["time"], TIME), zc=(["zc"], zc)),
             attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")),
    "W": ux.UxDataArray(W, uxgrid=uxgrid, dims=["time", "zf", "n_node"],
             coords=dict(time=(["time"], TIME), nz=(["zf"], zf)),
             attrs=dict(location="node", mesh="delaunay", Conventions="UGRID-1.0")),
}, uxgrid=uxgrid)

# --- Run advection ---
T_TOTAL = np.timedelta64(3600, "s")
DT = np.timedelta64(300, "s")
dt_s = 300.0

fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat")
pset = parcels.ParticleSet(fieldset, lon=[5.0], lat=[5.0])
pfile = parcels.ParticleFile(store="test.zarr", outputdt=DT)

pset.execute(parcels.kernels.AdvectionRK4, runtime=T_TOTAL, dt=DT,
             output_file=pfile, verbose_progress=False)

# --- Check ---
lon_expected = 5.0 + U0 * 3600.0   # 8.6
lon_lagged   = 5.0 + U0 * 3300.0   # 8.3

print(f"pset.lon[0]        = {pset.lon[0]}")
print(f"analytical(T)      = {lon_expected}")
print(f"analytical(T - dt) = {lon_lagged}")
print(f"error vs T:        = {abs(pset.lon[0] - lon_expected):.6e}")
print(f"error vs T - dt:   = {abs(pset.lon[0] - lon_lagged):.6e}")
# Output:
#   error vs T      ≈ U0 * dt = 3e-1   (proportional to dt)
#   error vs T - dt ≈ 0                 (machine precision)

The issue reproduces with the AdvectionEE and AdvectionRK4 time integrators (I have not tested AdvectionRK45, though I don't believe the integrator choice matters). I suspect there is an issue in the particleset.execute() and file io logic.

[fluidnumerics-joe]

Interesting observation here.. When the particlefile is not included in the pset.execute() call, the particle position is correct.

[fluidnumerics-joe]

I ran a simple test that just forward steps twice with a time step size of 10 minutes. The particlefile output is configured to write particle positions at each time step. The output table below shows the time, time difference from initial time, longitude stored in the zarr output, the analytical solution for the longitude position, and the error. It looks like the initial condition is repeated after the first time step.

obs                        time    dt (s)    lon (stored)  lon (analytical)         error
     ------------------------------------------------------------------------------------------
       0  2000-01-01T00:00:00.000000000         0      5.00000000        5.00000000      0.00e+00
       1  2000-01-01T00:10:00.000000000       600      5.00000000        5.60000000     -6.00e-01
       2  2000-01-01T00:20:00.000000000      1200      5.59999990        6.20000000     -6.00e-01

Code to reproduce this output is below

"""

import shutil

import numpy as np
import parcels
from parcels.kernels import AdvectionEE
import uxarray as ux
import xarray as xr

# --- Build a minimal Delaunay grid ---
NX = 10
lon, lat = np.meshgrid(
    np.linspace(0, 20, NX, dtype=np.float64),
    np.linspace(0, 20, NX, dtype=np.float64),
)
lon_flat, lat_flat = lon.ravel(), lat.ravel()
mask = (
    np.isclose(lon_flat, 0) | np.isclose(lon_flat, 20)
    | np.isclose(lat_flat, 0) | np.isclose(lat_flat, 20)
)
uxgrid = ux.Grid.from_points(
    (lon_flat, lat_flat),
    method="regional_delaunay",
    boundary_points=np.flatnonzero(mask),
)
uxgrid.attrs["Conventions"] = "UGRID-1.0"

# --- Uniform velocity field on face centers ---
U0 = 0.001  # degrees/s
N_TIMES = 5
TIME = xr.date_range("2000-01-01", periods=N_TIMES, freq="1h")
zc = np.array([0.5])
zf = np.array([0.0, 1.0])

ds = ux.UxDataset({
    "U": ux.UxDataArray(
             np.full((N_TIMES, 1, uxgrid.n_face), U0),
             uxgrid=uxgrid, dims=["time", "zc", "n_face"],
             coords=dict(time=(["time"], TIME), zc=(["zc"], zc)),
             attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")),
    "V": ux.UxDataArray(
             np.full((N_TIMES, 1, uxgrid.n_face), 0.0),
             uxgrid=uxgrid, dims=["time", "zc", "n_face"],
             coords=dict(time=(["time"], TIME), zc=(["zc"], zc)),
             attrs=dict(location="face", mesh="delaunay", Conventions="UGRID-1.0")),
    "W": ux.UxDataArray(
             np.zeros((N_TIMES, 2, uxgrid.n_node)),
             uxgrid=uxgrid, dims=["time", "zf", "n_node"],
             coords=dict(time=(["time"], TIME), zf=(["zf"], zf)),
             attrs=dict(location="node", mesh="delaunay", Conventions="UGRID-1.0")),
}, uxgrid=uxgrid)

# --- Advect one particle for exactly 2 timesteps ---
DT = np.timedelta64(600, "s")
T_TOTAL = np.timedelta64(1200, "s")  # 2 * DT
x0 = 5.0

fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat")
pset = parcels.ParticleSet(fieldset, lon=[x0], lat=[10.0])
pfile = parcels.ParticleFile(store="lag_test.zarr", outputdt=DT)

pset.execute(AdvectionEE, runtime=T_TOTAL, dt=DT,
             output_file=pfile, verbose_progress=False)

# --- Read back and compare ---
traj = xr.open_zarr("lag_test.zarr")
t0 = traj["time"].values[0, 0]

print(f"Parcels version: {parcels.__version__}")
print()
print(f"{'obs':>3s}  {'time':>26s}  {'dt (s)':>8s}  {'lon (stored)':>14s}  {'lon (analytical)':>16s}  {'error':>12s}")
print("-" * 90)
for i in range(traj.sizes["obs"]):
    t = traj["time"].values[0, i]
    lon_stored = float(traj["lon"].values[0, i])
    elapsed_s = (t - t0) / np.timedelta64(1, "s")
    lon_exact = x0 + U0 * elapsed_s
    err = lon_stored - lon_exact
    print(f"{i:3d}  {str(t):>26s}  {elapsed_s:8.0f}  {lon_stored:14.8f}  {lon_exact:16.8f}  {err:12.2e}")

traj.close()
shutil.rmtree("lag_test.zarr", ignore_errors=True)

[fluidnumerics-joe]

In particleset.py lines 425-448. The outer loop starts with next_output = start_time (line 425) and time = start_time (line 427):

        next_output = start_time if output_file else None

        time = start_time
        while sign_dt * (time - end_time) < 0:
            if next_output is not None:
                f = min if sign_dt > 0 else max
                next_time = f(next_output, end_time)
            else:
                next_time = end_time

            self._kernel.execute(self, endtime=next_time, dt=dt)

For the reproducer shared in this issue, on the first iteration: next_output = start_time = 0, end_time = 3600, next_time = min(0, 3600) = 0. The kernel executes then with endtime = 0. Since particles are already at time = 0, time_to_endtime = 0, which clamps pset.dt to 0 inside kernel.execute() (line 210). Advection runs with dt=0, producing dlon=0, and PositionUpdate gets prepended at line 249.

In the next time iteration next_output = 600, next_time = min(600, 1200) = 600. Now, we're calling kernel.execute(endtime=600, dt=600) and the kernels are now [PositionUpdate, AdvectionEE]. Then,
Inner iteration 1 (kernel.py line 199):

• time_to_endtime = 600 - 0 = 600
• pset.dt clamped to min(600, 600) = 600
• PositionUpdate runs: lon += dlon (which is 0 from the initial step), time += 600 → time=600, dlon = 0
• AdvectionEE runs: computes dlon = U0 * 600 based on current position
• particles_endofloop: time == 600 == endtime → EndofLoop

After exiting the inner loop in kernel.py, we're back in the outer loop in pset.execute where we call output_file.write() at t=600. However, since PositionUpdate only applied the stale dlon=0 the wrong position is written.

@VeckoTheGecko - there's a few ways we can handle this. I'll suggest something in a branch...

[erikvansebille]

Thanks for raising this Issue, @fluidnumerics-joe. I probably introduced this bug during the kernel-loop refactoring, e.g. in #2333 or #2490 (or even before). The kernel loop is quite iffy and most of the testing so far was without output (which is probably why it does work as expected if you don't provide a ParticleFile)

If you want, I can also try working on a fix, but it will probably will have to wait until next week. Your reproducer is certainly very useful!

[fluidnumerics-joe]

I think I have a fix, just working through some broken tests at the moment. There will be a PR waiting for you to review at your leisure.