Set the runtime based on each particle lifetime

[iuryt]

Hi all,

I am struggling to setup an experiment that seems reasonable, but I could not find a good example from the documentation.

Let's say I have thousands of particles that will be released at different times, but I just want to track them until they reach a certain age.

In my case, I want to track particles for 30 days of lifetime with hourly samples, i.e. I expect to have only 30*24 observations of each particle.

If I set the runtime for 30 days, we will only track those particles that are relased in the first 30 days, but my simulation has 20 years of particles being released every 5 days and I just want to integrate each of them for 30 days. If I set the runtime for the whole 20 years, I will track most of the particles unnecessarily for longer than I need.

I could subset the velocity field and run many simulations (one per month), but I was just wondering if there is a more straightforward way to do that. If not, that could be an interesting feature to implement to the tool. What do you think?

[michaeldenes]

Hi @iuryt, one possible solution is that you could define an 'Age' variable for the particle, that only increments once the particle is released, via a kernel, perhaps. Within that kernel, once the particle reaches an age of 30 days, you can delete the particle (with particle.delete()). This way, after the 30 days, the particle will no longer be tracked, and the lon/lat/depth in the data should show up as NaNs.

[michaeldenes]

If you delete particles after 30 days, and the total number of particles that are iterated over decreases, then you should expect a speed up. However, in your case, as a set of particles are deleted, a new set will start (as they are initialised in the future), and so you're not decreasing the number of particles you are iterating over, it remains roughly the same. So I'm not entirely sure if you should expect a speed up or not. Please do let us know how it goes though, I'm interested to see if it does in fact speed things up!

[iuryt]

I ran for 2 months. If I dont delete the particles
100%|██████████| 5184000.0/5184000.0 [01:19<00:00, 65323.84it/s]
Now if I delete the particles, it is a little bit faster:
100%|██████████| 5184000.0/5184000.0 [01:00<00:00, 86339.17it/s]
However, I keep getting WARNING: ParticleSet is empty on writing as array at time .... warnings.
This the part of the code I am using to run this

def Age(particle, fieldset, time):
    particle.age += particle.dt / 86400

def DeleteParticle(particle, fieldset, time):
    if particle.age > 30:
        particle.delete()


AgeParticle = parcels.JITParticle.add_variable(parcels.Variable("age", initial=0))

ploni, plati, timei = [], [], []
for di in np.arange(0, 30., 5):
    ind = np.random.randint(0, plon.size-1, nparticles)
    ploni.append(plon[ind])
    plati.append(plat[ind])
    timei.append([timedelta(days = di).total_seconds()] * nparticles)
ploni, plati, timei = np.concatenate(ploni), np.concatenate(plati), np.concatenate(timei)


pset = parcels.ParticleSet(fieldset=fieldset, pclass=AgeParticle, lon=ploni, lat=plati, time=timei)

output_fname = f'parcels_{time.strftime("%Y_%m_%d")}.zarr'
output_file = pset.ParticleFile(name=os.path.join(path, output_fname), outputdt=timedelta(hours=3))

pset.execute(
    [Age, parcels.AdvectionRK4, smagdiff, DeleteParticle], runtime=timedelta(days=60), dt = timedelta(hours=1), output_file = output_file, verbose_progress = True,
)

[michaeldenes]

That's odd, if there are some particles, that warning shouldn't pop up. Do you. get the desired output even with the warnings?

Perhaps you could keep a second variable to store the "particle start time", and then in the Age kernel, use an if statement for (if particle.dt >+ particle.start_time ...).

[iuryt]

Yes, I get the desired output. I see what you are thinking. Probably the error is occurring before the particle start time. I will try that.

[iuryt]

I think this is something else @michaeldenes

I tried:

def Age(particle, fieldset, time):
    if time>=particle.time:
        particle.age += particle.dt / 86400

def DeleteParticle(particle, fieldset, time):
    if particle.age > 30:
        particle.delete()

AgeParticle = parcels.JITParticle.add_variable(parcels.Variable("age", initial=0))

and I am still getting:

WARNING: ParticleSet is empty on writing as array at time 4.752e+06
WARNING: ParticleSet is empty on writing as array at time 4.7628e+06
WARNING: ParticleSet is empty on writing as array at time 4.7736e+06
WARNING: ParticleSet is empty on writing as array at time 4.7844e+06
WARNING: ParticleSet is empty on writing as array at time 4.7952e+06
WARNING: ParticleSet is empty on writing as array at time 4.806e+06
WARNING: ParticleSet is empty on writing as array at time 4.8168e+06
WARNING: ParticleSet is empty on writing as array at time 4.8276e+06
WARNING: ParticleSet is empty on writing as array at time 4.8384e+06
WARNING: ParticleSet is empty on writing as array at time 4.8492e+06
WARNING: ParticleSet is empty on writing as array at time 4.86e+06
WARNING: ParticleSet is empty on writing as array at time 4.8708e+06
WARNING: ParticleSet is empty on writing as array at time 4.8816e+06
WARNING: ParticleSet is empty on writing as array at time 4.8924e+06
WARNING: ParticleSet is empty on writing as array at time 4.9032e+06
WARNING: ParticleSet is empty on writing as array at time 4.914e+06
WARNING: ParticleSet is empty on writing as array at time 4.9248e+06
WARNING: ParticleSet is empty on writing as array at time 4.9356e+06
WARNING: ParticleSet is empty on writing as array at time 4.9464e+06
WARNING: ParticleSet is empty on writing as array at time 4.9572e+06
WARNING: ParticleSet is empty on writing as array at time 4.968e+06
WARNING: ParticleSet is empty on writing as array at time 4.9788e+06
WARNING: ParticleSet is empty on writing as array at time 4.9896e+06
WARNING: ParticleSet is empty on writing as array at time 5.0004e+06
WARNING: ParticleSet is empty on writing as array at time 5.0112e+06
WARNING: ParticleSet is empty on writing as array at time 5.022e+06
WARNING: ParticleSet is empty on writing as array at time 5.0328e+06
WARNING: ParticleSet is empty on writing as array at time 5.0436e+06
WARNING: ParticleSet is empty on writing as array at time 5.0544e+06
WARNING: ParticleSet is empty on writing as array at time 5.0652e+06
WARNING: ParticleSet is empty on writing as array at time 5.076e+06
WARNING: ParticleSet is empty on writing as array at time 5.0868e+06
WARNING: ParticleSet is empty on writing as array at time 5.0976e+06
WARNING: ParticleSet is empty on writing as array at time 5.1084e+06
WARNING: ParticleSet is empty on writing as array at time 5.1192e+06
WARNING: ParticleSet is empty on writing as array at time 5.13e+06
WARNING: ParticleSet is empty on writing as array at time 5.1408e+06
WARNING: ParticleSet is empty on writing as array at time 5.1516e+06
WARNING: ParticleSet is empty on writing as array at time 5.1624e+06
WARNING: ParticleSet is empty on writing as array at time 5.1732e+06

I believe this is likely because my runtime exceeds the time it takes to run all particles for 30 days; that is, there are timesteps when no particles are present in pset, which explains the warnings.
I conducted a test with a runtime of 50 days, and there were no warnings. How can I stop pset.execute() when no particles are left?