Repository: OceanParcels/Parcels (python)

Issue:
Backward in time outputs not written at desired frequency
### Parcels version

3.0.5

### Description

Looks like a regression from 3.0.4 to 3.0.5.

### Code sample

```python
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`:
```text
AssertionError: 
Arrays are not almost equal to 7 decimals
 ACTUAL: 3600000000000.0
 DESIRED: 86400000000000.0
```

With `v3.0.4` it works.


Pull request:
Fixing bug in backward-in-time-writing
This fixes #1722, by reverting a change in ed11e43

Thanks for the git-bisect hint, @VeckoTheGecko!

Patch:
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
diff --git a/pyproject.toml b/pyproject.toml
index 55ded6ff..5e5c2106 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -52,7 +52,7 @@ local_scheme = "no-local-version"
 
 [tool.pytest.ini_options]
 addopts = ["--strict-config", "--strict-markers"]
-testpaths = ["parcels/tests", "docs/examples"]
+testpaths = ["tests", "docs/examples"]
 python_files = ["test_*.py", "example_*.py", "*tutorial*"]
 minversion = "7"
 markers = [ # can be skipped by doing `pytest -m "not slow"` etc.


Test patch:
diff --git a/tests/test_particlefile.py b/tests/test_particlefile.py
index 050823d9..15ef80c6 100755
--- a/tests/test_particlefile.py
+++ b/tests/test_particlefile.py
@@ -1,4 +1,6 @@
 import os
+import tempfile
+from datetime import timedelta
 
 import cftime
 import numpy as np
@@ -6,9 +8,11 @@ import pytest
 import xarray as xr
 from zarr.storage import MemoryStore
 
+import parcels
 from parcels import (
     AdvectionRK4,
     Field,
+    FieldSet,
     JITParticle,
     ParticleSet,
     ScipyParticle,
@@ -344,3 +348,69 @@ def test_reset_dt(fieldset, mode, tmpdir):
     pset.execute(pset.Kernel(Update_lon), endtime=0.12, dt=0.02, output_file=ofile)
 
     assert np.allclose(pset.lon, 0.6)
+
+
+def setup_pset_execute(*, fieldset: FieldSet, outputdt: timedelta, execute_kwargs, particle_class=ScipyParticle):
+    npart = 10
+
+    if fieldset is None:
+        fieldset = create_fieldset_zeros_simple()
+
+    pset = ParticleSet(
+        fieldset,
+        pclass=particle_class,
+        lon=np.full(npart, fieldset.U.lon.mean()),
+        lat=np.full(npart, fieldset.U.lat.mean()),
+    )
+
+    with tempfile.TemporaryDirectory() as dir:
+        name = f"{dir}/test.zarr"
+        output_file = pset.ParticleFile(name=name, outputdt=outputdt)
+
+        pset.execute(DoNothing, output_file=output_file, **execute_kwargs)
+        ds = xr.open_zarr(name).load()
+
+    return ds
+
+
+def test_pset_execute_outputdt_forwards():
+    """Testing output data dt matches outputdt in forward time."""
+    outputdt = timedelta(hours=1)
+    runtime = timedelta(hours=5)
+    dt = timedelta(minutes=5)
+
+    ds = setup_pset_execute(
+        fieldset=create_fieldset_zeros_simple(), outputdt=outputdt, execute_kwargs=dict(runtime=runtime, dt=dt)
+    )
+
+    assert np.all(ds.isel(trajectory=0).time.diff(dim="obs").values == np.timedelta64(outputdt))
+
+
+def test_pset_execute_outputdt_backwards():
+    """Testing output data dt matches outputdt in backwards time."""
+    outputdt = timedelta(hours=1)
+    runtime = timedelta(days=2)
+    dt = -timedelta(minutes=5)
+
+    ds = setup_pset_execute(
+        fieldset=create_fieldset_zeros_simple(), outputdt=outputdt, execute_kwargs=dict(runtime=runtime, dt=dt)
+    )
+    file_outputdt = ds.isel(trajectory=0).time.diff(dim="obs").values
+    assert np.all(file_outputdt == np.timedelta64(-outputdt))
+
+
+def test_pset_execute_outputdt_backwards_fieldset_timevarying():
+    """test_pset_execute_outputdt_backwards() still passed despite #1722 as it doesn't account for time-varying fields,
+    which for some reason #1722
+    """
+    outputdt = timedelta(hours=1)
+    runtime = timedelta(days=2)
+    dt = -timedelta(minutes=5)
+
+    # TODO: Not ideal using the `download_example_dataset` here, but I'm struggling to recreate this error using the test suite fieldsets we have
+    example_dataset_folder = parcels.download_example_dataset("MovingEddies_data")
+    fieldset = parcels.FieldSet.from_parcels(f"{example_dataset_folder}/moving_eddies")
+
+    ds = setup_pset_execute(outputdt=outputdt, execute_kwargs=dict(runtime=runtime, dt=dt), fieldset=fieldset)
+    file_outputdt = ds.isel(trajectory=0).time.diff(dim="obs").values
+    assert np.all(file_outputdt == np.timedelta64(-outputdt)), (file_outputdt, np.timedelta64(-outputdt))
diff --git a/tests/test_particlesets.py b/tests/test_particlesets.py
index 5e0befd7..9973e1cd 100644
--- a/tests/test_particlesets.py
+++ b/tests/test_particlesets.py
@@ -22,6 +22,13 @@ def fieldset():
     return create_fieldset_zeros_simple()
 
 
+@pytest.fixture
+def pset(fieldset):
+    npart = 10
+    pset = ParticleSet(fieldset, pclass=JITParticle, lon=np.linspace(0, 1, npart), lat=np.zeros(npart))
+    return pset
+
+
 @pytest.mark.parametrize("mode", ["scipy", "jit"])
 def test_pset_create_lon_lat(fieldset, mode):
     npart = 100