Default ParticleFile output frequency results in infinite loop

[observingClouds]

Parcels version

3.1.2.dev4

Description

The default output frequency of ParticleFile (np.inf) causes the parcel simulation in the execution to end up in an infinite loop without progressing.

I would expect from the default value that the output at all timesteps is written to file.

Code sample

The following code snippet is identical to test_fieldset_from_xarray with the exception that the output_file attribute is set.

import numpy as np
import xarray as xr
from parcels import FieldSet, ParticleSet, JITParticle, AdvectionRK4, ParticleFile

def generate_dataset(xdim, ydim, zdim=1, tdim=1):
    lon = np.linspace(0.0, 12, xdim, dtype=np.float32)
    lat = np.linspace(0.0, 12, ydim, dtype=np.float32)
    depth = np.linspace(0.0, 20.0, zdim, dtype=np.float32)
    if tdim:
        time = np.linspace(0.0, 10, tdim, dtype=np.float64)
        Uxr = np.ones((tdim, zdim, ydim, xdim), dtype=np.float32)
        Vxr = np.ones((tdim, zdim, ydim, xdim), dtype=np.float32)
        for t in range(Uxr.shape[0]):
            Uxr[t, :, :, :] = t / 10.0
        coords = {"lat": lat, "lon": lon, "depth": depth, "time": time}
        dims = ("time", "depth", "lat", "lon")
    else:
        Uxr = np.ones((zdim, ydim, xdim), dtype=np.float32)
        Vxr = np.ones((zdim, ydim, xdim), dtype=np.float32)
        for z in range(Uxr.shape[0]):
            Uxr[z, :, :] = z / 2.0
        coords = {"lat": lat, "lon": lon, "depth": depth}
        dims = ("depth", "lat", "lon")
    return xr.Dataset(
        {"Uxr": xr.DataArray(Uxr, coords=coords, dims=dims), "Vxr": xr.DataArray(Vxr, coords=coords, dims=dims)}
    )


def test_fieldset_from_xarray(tdim):
    ds = generate_dataset(3, 3, 2, tdim)
    variables = {"U": "Uxr", "V": "Vxr"}
    if tdim:
        dimensions = {"lat": "lat", "lon": "lon", "depth": "depth", "time": "time"}
    else:
        dimensions = {"lat": "lat", "lon": "lon", "depth": "depth"}
    fieldset = FieldSet.from_xarray_dataset(ds, variables, dimensions, mesh="flat")
    assert fieldset.U._creation_log == "from_xarray_dataset"

    pset = ParticleSet(fieldset, JITParticle, 0, 0, depth=20)

    outputfile = ParticleFile("test_fieldset_from_xarray", pset)

    pset.execute(AdvectionRK4, dt=1, runtime=10, output_file=outputfile)
    if tdim == 10:
        assert np.allclose(pset.lon_nextloop[0], 4.5) and np.allclose(pset.lat_nextloop[0], 10)
    else:
        assert np.allclose(pset.lon_nextloop[0], 5.0) and np.allclose(pset.lat_nextloop[0], 10)

test_fieldset_from_xarray(10)

INFO: Output files are stored in test_fieldset_from_xarray.zarr.
 10%|█         | 1.0/10.0 [00:00<00:06,  1.49it/s]

Setting the output frequency to e.g. 1 (outputfile = ParticleFile("test_fieldset_from_xarray", pset, outputdt=1)) results in the expected results.

INFO: Output files are stored in test_fieldset_from_xarray.zarr.
100%|██████████| 10.0/10.0 [00:00<00:00, 18.33it/s]

[erikvansebille]

Thanks for reporting this bug, @observingClouds. I agree that that the np.inf default value for outputdt is confusing, so I suggest we remove it altogether. Have just created PR #1808 to fix this