Backward in time outputs not written at desired frequency

[willirath]

Parcels version

3.0.5

Description

Looks like a regression from 3.0.4 to 3.0.5.

Code sample

from datetime import timedelta
import xarray as xr
import numpy as np
import parcels

example_dataset_folder = parcels.download_example_dataset("MovingEddies_data")
fieldset = parcels.FieldSet.from_parcels(f"{example_dataset_folder}/moving_eddies")

pset = parcels.ParticleSet.from_list(
    fieldset=fieldset,
    pclass=parcels.JITParticle,
    lon=[3.3e5, 3.3e5],
    lat=[1e5, 2.8e5],
)

# run forward

output_file = pset.ParticleFile(
    name="EddyParticles.zarr",
    outputdt=timedelta(hours=1),
)
pset.execute(
    parcels.AdvectionRK4,
    runtime=timedelta(days=2),
    dt=timedelta(minutes=5),
    output_file=output_file,
)

# run backward

output_file = pset.ParticleFile(
    name="EddyParticles_Bwd.zarr",
    outputdt=timedelta(hours=1),
)
pset.execute(
    parcels.AdvectionRK4,
    runtime=timedelta(days=2),
    dt=-timedelta(minutes=5),
    output_file=output_file,
)

# check that timesteps match (up to a sign)

ds = xr.open_zarr("EddyParticles.zarr/")
ds_bwd = xr.open_zarr("EddyParticles_Bwd.zarr/")

np.testing.assert_almost_equal(
    ds.time.isel(trajectory=0).diff("obs").mean().compute().astype(float).data[()],
    - ds_bwd.time.isel(trajectory=0).diff("obs").mean().compute().astype(float).data[()],
)

With v3.0.5:

AssertionError: 
Arrays are not almost equal to 7 decimals
 ACTUAL: 3600000000000.0
 DESIRED: 86400000000000.0

With v3.0.4 it works.

[willirath]

Further hints: The output times appear to match the time steps of the fieldset.

[erikvansebille]

Thanks for reporting @willirath, I'll take a look asap. This could be the same issue that @sruehs previously encountered (personal communication); this minimal example is super-useful!

[VeckoTheGecko]

Git bisect points to ed11e43173fd246c2d1043df9afc8be45e742d12

[VeckoTheGecko]

Ah, this was the change

- next_callback = starttime + callbackdt if dt > 0 else starttime - callbackdt
+ next_callback = starttime * np.sign(dt)

but we need

diff --git a/parcels/particleset.py b/parcels/particleset.py
index d7179cf3..72696eb6 100644
--- a/parcels/particleset.py
+++ b/parcels/particleset.py
@@ -1137,7 +1137,7 @@ class ParticleSet:
             next_output = starttime + dt
         else:
             next_output = np.inf * np.sign(dt)
-        next_callback = starttime * np.sign(dt)
+        next_callback = starttime + callbackdt * np.sign(dt)

         tol = 1e-12
         time = starttime