error after pset.execute with runtime length

[keli559]

Running a model similar to the first two cases in:
https://nbviewer.org/github/OceanParcels/parcels/blob/master/parcels/examples/parcels_tutorial.ipynb
My ocean circulation model is different. It has been running OK, except when I increase the runtime in pset.execute. It shows pink with a progress bar. When it reaches 72%, It gives me the follow error (last few lines)
`
--> 492 return self._collection.toDictionary(pfile=pfile, time=time,
493 deleted_only=deleted_only)

File ~\Miniconda3\envs\py3_parcels\lib\site-packages\parcels\collection\collectionsoa.py:856, in ParticleCollectionSOA.toDictionary(self, pfile, time, deleted_only)
854 pset_errs = ((self._data['state'][indices_to_write] != OperationCode.Delete) & np.greater(np.abs(time - self._data['time'][indices_to_write]), 1e-3, where=np.isfinite(self._data['time'][indices_to_write])))
855 if np.count_nonzero(pset_errs) > 0:
--> 856 logger.warning_once('time argument in pfile.write() is {}, but particles have time {}'.format(time, self._data['time'][pset_errs]))
858 if len(pfile.var_names_once) > 0:
859 first_write = (_to_write_particles(self._data, time) & _is_particle_started_yet(self._data, time) & np.isin(self._data['id'], pfile.written_once, invert=True))

IndexError: boolean index did not match indexed array along dimension 0; dimension is 34 but corresponding boolean dimension is 12
72%|█████████████████████████████████████████████▉ | 3099600.0/4320000.0 [00:23<00:04, 284623.96it/s]
My code is:
output_file = pset.ParticleFile(name=savePartFn, outputdt=timedelta(hours=1))
pset.execute(AdvectionRK4,
#runtime=timedelta(days=dT/(24.0*3600.0)),
runtime=timedelta(days=50),
dt=timedelta(minutes=5),
output_file=output_file, recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})
`

when days = 8 or 20, all is good, when it is 50, it starts to report error.... Thanks

[erikvansebille]

Thanks for reporting this error, @keli559. I had not seen this bug before, and cannot reproduce it without having your code or data. However, note that we are about to release a new version of Parcels that has a complete overhaul of the way output is written (#1199); so if you wait a few weeks more this error might disappear?

[keli559]

Feel free to take a look at the following code, @erikvansebille. It is very similar to the first two cases in the second tutorial on the OceanParcels sight (I am still learning .... ) It has been working alright, I used the DeleteParticle function that you suggested before to take off the points that are advected outside of the domain.
This time, I set number of particles to be 50, released sequentially along the time axis of the eulerian data (ocean circulation output), The total time axis length of the eurlerian data is about 3 months. However, I only allow the pset.execute to run 50 days as runtime. That is where I got the error from. Thanks
`
import xarray
%matplotlib inline
from parcels import FieldSet, ParticleSet, Variable, JITParticle, AdvectionRK4, plotTrajectoriesFile, ErrorCode
import numpy as np
import math
from datetime import timedelta
import datetime
from operator import attrgetter
import matplotlib.pyplot as plt
folder = "C:\Users\kl591\parcels_examples\parcels_examples\cop_mer\"

def DeleteParticle(particle, fieldset, time):
particle.delete()
##################
ds = xarray.open_dataset(folder + "data_v2.nc")
T_start = np.datetime64('2022-05-01')
T_end = np.datetime64('2022-08-25')
Npart = 50
dT = np.timedelta64(T_end-T_start, 's')/np.timedelta64(1, 's')
ds1 = ds.sel(depth=0.494, method="nearest").drop("depth")
ds1 = ds1.sel(time=slice(T_start, T_end))
######################
saveFn = '{0}_{1}_v1.nc'.format(str(T_start)[4:].replace('-', ''), str(T_end)[4:].replace('-', ''))
savePartFn = 'part{0}.nc'.format(str(T_start)[0:7].replace('-', ''))
ds1.to_netcdf(folder+saveFn)
######################
dimensions = {"lat":"latitude", "lon":"longitude", "time": "time"}
variables = {"U":"utotal", "V":"vtotal"}
fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)
fieldset.V.show()
########################
pset = ParticleSet.from_line(fieldset=fieldset, pclass=JITParticle,
size=Npart, # releasing 5 particles
start=(14.1, -24.0), # releasing on a line: the start longitude and latitude
finish=(14.1, -24.0), # releasing on a line: the end longitude and latitude
time = np.linspace(0.0, dT, Npart))
print(pset)
pset.show(field=fieldset.U)
########################
output_file = pset.ParticleFile(name=savePartFn, outputdt=timedelta(hours=1))
pset.execute(AdvectionRK4,
#runtime=timedelta(days=dT/(24.0*3600.0)),
runtime=timedelta(days=50),
dt=timedelta(minutes=5),
output_file=output_file, recovery={ErrorCode.ErrorOutOfBounds: DeleteParticle})

...

...

...

...

`

[keli559]

Dear Dr. @erikvansebille, I am trying to track a series of particles released at a fixed (lat, lon) and each is released 6 hours after the other. I am using ParticleSet.from_line function. What would be a good function that you would recommend me using if I want to track 500 of them... many thanks.... so currently, as you see, I am using from_line... just trying to figure out a work-around. thanks

[VeckoTheGecko]

If you'd like to instantiate a particleset with arbitrary lon, lat and time, you can use the ParticleSet class directly as documented here in the documentation.

times = # Array of np.timedelta64 objects
lons = np.full(Npart, 14.1)
lats = np.full(Npart, -24.0)
ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=lons, lat=lats, time=times)

.from_line() should still work for your purposes, I'm not sure where you're encountering limitations with the method.

[keli559]

Hi all, I figured out the problem. this is what I did... 1) I switch from jupyter notebook to .py script (I am much more comfortable with .py) 2) I added allow_time_extrapolation..... I don't know how it is needed as the particle tracking time range is the same as the ocean circulation time axis length..... but it worked.... fieldset = FieldSet.from_netcdf(filenames, variables, dimensions, allow_time_extrapolation=True)

[keli559]

My last question here: how do I save the movies I made using plotTrajectoriesFile()? Many thanks!

[VeckoTheGecko]

My last question here: how do I save the movies I made using plotTrajectoriesFile()?

Taking a peak at the source code, you can provide mode="movie2d_notebook" which will return the corresponding matplotlib anim object. From here, its just a matter of using the same tooling that matplotlib provides in order to save animations.

eg.

anim = plotTrajectoriesFile(pfile, mode = "movie2d_notebook")

from matplotlib.animation import FFMpegWriter
writervideo = FFMpegWriter(fps=60)
anim.save("my_simulation.mp4", writer=writervideo)