Treatment of particles started at different time

[philippemiron]

Hi everyone,

I have a question to understand what is done with particles initialized at different times.

For example, with a dt=30min and outputdt=60min and the following particles:

P1 = (lon1, lat1, t1=1:35)
P2 = (lon2, lat3, t2=2:50)
P3 = (lon3, lat3, t3=4:20)

My expectation was that the time in the output file would be independent for each particle, but instead, it is a function of the earliest release time. More precisely, P2's output time = [3:35, 4:35, 5:35, ...] and P3 starts with time = [4:35, 5:35, ...].

My question is: are particles independently evolved and then interpolated to a common time? If this is the case, why not simply avoid the interpolation (when the dt is a multiple of the outputdt) and output the positions obtained from the integration since the time variable in the output file is already a 2d matrix?

Thanks,

[erikvansebille]

HI @philippemiron, thanks for reporting. Can I ask which version of Parcels you're running? Because something related was recently addressed in #1095?

And then, could you please add a minimal working (breaking) example that shows this unexpected behaviour, so that we can use that to try bugfix it, if the issue isn't resolved already?

[philippemiron]

Hi @erikvansebille,

Thanks for quick answer. I'm currently using version 2.3.0.

Here is a minimal example (inspired by one of the tutorial):

from parcels import FieldSet, ParticleSet, JITParticle, plotTrajectoriesFile
from parcels import AdvectionRK4
from datetime import timedelta
import numpy as np
import xarray as xr

seconds_per_day = 60*60*24
a, b = 20, 10  # domain size
lon = np.linspace(0, a, 21, dtype=np.float32)
lat = np.linspace(0, b, 11, dtype=np.float32)
dx, dy = lon[2]-lon[1], lat[2]-lat[1]
times = np.linspace(0, 10, 11)*seconds_per_day

U = np.ones((times.size, lat.size, lon.size), dtype=np.float32)*1e-4
V = np.zeros((times.size, lat.size, lon.size), dtype=np.float32)

data = {'U': U, 'V': V}
dimensions = {'lon': lon, 'lat': lat, 'time': times}
fieldset = FieldSet.from_data(data, dimensions, mesh='flat')

X = np.array([1, 1, 1])
Y = np.array([1, 2, 3])
T = np.array([1, 2, 3])*60*60 # hours
pset = ParticleSet(fieldset, pclass=JITParticle, lon=X, lat=Y, time=T)

output_file = pset.ParticleFile(name='test.nc', outputdt=timedelta(hours=6))

kernel = pset.Kernel(AdvectionRK4)

pset.execute(kernel,
             dt=timedelta(hours=1),
             runtime=timedelta(days=1),
             output_file=output_file,
             verbose_progress=True)

output_file.export()
output_file.close()

with xr.open_dataset('test.nc') as df:
    for i in range(0, df.dims['traj']):
        print(df['time'][i].values/3600/1e9)

The velocity field is only a small positive value in x (so particles don't exit). The different parameters are:

• 3 particles: X=[1,1,1] Y=[1,2,3] T=[1,2,3] hours
• dt = 1 hours
• outputdt = 6 hours
• runtime = 1 day

The output times for each particles is:

• T1 = [ 1 7 13 19 25]
• T2 = [ 7 13 19 25 'NaT']
• T3 = [ 7 13 19 25 'NaT']

I would have expected: [2, 8, 14, 20] and [3, 9, 15, 21], respectively for T2 and T3. With the current output, I'm wondering how the particles are evolved.

[erikvansebille]

Thanks for this elaboration, @philippemiron. I have now looked into it; and in fact this is a design choice (not a bug) in Parcels. By synchronising output time, it is easier to create e.g. maps or animations of the state of all particles at a specific time.

Now, if you want to do so for all particles at a specific age, then you may need to adjust outputdt so that it is a divisor of all times, or do multiple simulations.

See also the extra section I added to the delaystart tutorial in #1103. Does this answer your question?

[philippemiron]

Hi @erikvansebille,

Ok, I understand the synchronization and why this is easier to later plot particles all together at specific times.

I now realize that my example didn't cover the second interrogation I had. See the following case, where the particles times are not multiple of the outputdt nor the dt:

from parcels import FieldSet, ParticleSet, JITParticle, plotTrajectoriesFile
from parcels import AdvectionRK4
from datetime import timedelta
import numpy as np
import xarray as xr

a, b = 20, 10  # domain size
lon = np.linspace(0, a, 21, dtype=np.float32)
lat = np.linspace(0, b, 11, dtype=np.float32)
dx, dy = lon[2]-lon[1], lat[2]-lat[1]
times = np.linspace(0, 10, 11)

U = np.ones((times.size, lat.size, lon.size), dtype=np.float32)
V = np.zeros((times.size, lat.size, lon.size), dtype=np.float32)

data = {'U': U, 'V': V}
dimensions = {'lon': lon, 'lat': lat, 'time': times}
fieldset = FieldSet.from_data(data, dimensions, mesh='flat')

X = np.array([1, 1, 1])
Y = np.array([1, 2, 3])
T = np.array([0, 1.1, 2.2])

pset = ParticleSet(fieldset, pclass=JITParticle, lon=X, lat=Y, time=T)

output_file = pset.ParticleFile(name='test.nc', outputdt=2)

kernel = pset.Kernel(AdvectionRK4)

pset.execute(kernel,
             dt=1,
             endtime=4,
             output_file=output_file,
             verbose_progress=True)

output_file.export()
output_file.close()

with xr.open_dataset('test.nc') as df:
    for i in range(0, df.dims['traj']):
        print(df['time'][i].values/1e9)

The output is:

[0 2 4]
[    2     4 'NaT']
[    2     4 'NaT']

1. For the second particle, is it advected with dt=1 so t2=[1.1, 2.1, 3.1, 4.1] and then interpolated in the output file at the same time as the first particle, i.e t2=[2, 4, NaT]? Or is it advected first with dt=0.9 to synchronize it with the first particle?
2. For the third particle, how can we get something in the output file at t=2 if it starts at t=2.2? If I move the initial time of the third particle to ]2.5, 4[, it is only output at [4, NaT, NaT].

[philippemiron]

Oh, I see! So in my case, if I'm forecasting particles from their last known positions, which are given at times during the day, I would have to integrate them all individually right (if I want to avoid rounding)? Or I guess if I decrease the dt, the error becomes negligible.

Thanks!

[erikvansebille]

Sorry, I now relaised that my answer above was not entirely correct. Parcels does deal non-multiples of dt in starttimes. See example below:

fieldset = FieldSet.from_data(data={'U':0., 'V':0.}, dimensions={'lon': 0, 'lat': 0})
times = [0, 1.1]
pset = pset_type[pset_mode]['pset'](fieldset, lon=[0]*len(times), lat=[0]*len(times), pclass=ptype[mode], time=times)

def Addlon(particle, fieldset, time):
    particle.lon += particle.dt

pset.execute(Addlon, dt=1, runtime=2)
assert np.allclose([p.lon for p in pset], [2 - t for t in times])
print(pset)

which returns

P[0](lon=2.000000, lat=0.000000, depth=0.000000, time=2.000000)
P[1](lon=0.900000, lat=0.000000, depth=0.000000, time=2.000000)

So the longitudes are now [2.0, 0.9], which means that the dt for the second particle is dt=0.9 on the first execute of that particle.

[philippemiron]

Oh, perfect then. It does a smaller dt to synchronize.

[erikvansebille]

Exactly!

[philippemiron]

Maybe this should be included in the delay start notebook as a future reference, otherwise, that resolves the issue. Thanks again.