Counting particles that exit the FieldSet domain

[martoconnh]

Hi!

Is there any way I can count the number of particles that go out of bounds in my simulation? I have used the following kernel in order to delete every particle that goes out of my domain (it's used in #1412 ):

def DeleteParticle(particle, fieldset, time):
    if particle.state == StatusCode.ErrorOutOfBounds:
        particle.delete()

I can't find a way to do this. I thought maybe I could create another kernel somehow but I can't come up with anything...

[erikvansebille]

The easiest solution is to add a Variable to each particle (e.g. particle.deleted) and then set that to 1 when the particle is deleted. Then, you can analyse which particles are deleted in the zarr file

So (note, not tested as I don’t have access to a laptop at the moment)

class CustomParticle(JITParticle):
    deleted = Variable(‘deleted’, initial=0)

def DeleteParticle(particle, fieldset, time):
    if particle.state == StatusCode.ErrorOutOfBounds:
        particle.deleted=1
        particle.delete()

[erikvansebille]

OK, so do you think that these 126 lost particles are deleted even before you start the pset.execute()? What is the RunTimeWarning you're getting?

[martoconnh]

I don't think so, I checked right before the pset.execute() and my ParticleSet has length = 3488. But when I check in the output file there are only 3362. The RuntimeWarning is the following (it was already on my last comment, sorry I didn't make it clear):

254: RuntimeWarning: invalid value encountered in cast
 flat_num_dates_ns_int = (flat_num_dates * _NS_PER_TIME_DELTA[delta]).astype(

[erikvansebille]

OK, that RunTimeWarning shouldn't be the problem.

So what if you do a pset.execute() with a dt=0 and a runtime=0? What is the len(pset) before and after you do the execute?

[martoconnh]

Then I just get this:

RuntimeError: dt and runtime are zero, or endtime is equal to Particle.time. ParticleSet.execute() will not do anything.

Besides, the number of particles depends entirely on the runtime I'm using, as I realease them in different times along my simulation. So setting dt=0 and runtime=0 would mean that none of my particles are realeased yet, as I only have 2 particles to be released in time=0.

In case this is not clear enough, I'm simulating particles released from 2 source points, released at different rates from each one. So for instance, if I simulate 1 day (what I have been doing), I will be releasing 3488 particles, but with 2 days it will be 6999 and so on if more days are simulated. Don't know if this is needed but I thought some context would be useful.

[erikvansebille]

Hmm, I didn’t realize you used such a complicated release scheme. Perhaps then invert the bug-finding: which of the particles you did expect to be in your output file did not turn up there? Is there a pattern to the 126 ‘lost’ particles?

[martoconnh]

Hi!

I have not tried yet what you said earlier, as right now I'm not really sure how can I know which particles are the ones missing, as they're supposed to be indistinguishable. But I have tried some other things that may be useful (I hope).

First of all, I have set 2 new source points so now there should not be any particle that goes out of bounds. I have done the pset.execute() , and the particleset lenght before and after is the same, 3488, so no out of bounds particles, as expected. But when I check the output file, I'm still missing these 126 particles. Also, I have tried modifying the dt argument in the pset.execute() and the outputdt argument in the output_file.

Previously I had (losing 126 particles):

dt = timedelta(minutes=15)
outputdt = timedelta(hours=1)

I have now tried setting both of them to 1 hour (losing 72 particles), and I also have tried 15 minutes for both (losing only 17 particles but increasing a LOT the execution time).

Finally, I have also checked how many particles I'm losing for 2 and 3 days simulation with my original dt and outputdt, obtaining:
· 2 days (6999 particles) --> 127 lost particles
· 3 days (11471 particles) --> 163 lost particles

I really hope this helps somehow...

Thanks once again!!

[erikvansebille]

Could it be that some particles are not in the output file because the simulation has not run that far yet? So they are the last few of the output?

[martoconnh]

Hi!

I have been trying by setting both dt = outputdt = 1 hour. Now I am missing 72 particles as I said before. I have checked the output file, and it turns out that in fact I'm missing the last particles released from each source point. So basically I'm just failing to save the data for the last time values in my simulation. Maybe I should just add more simulation time in order to get one more dt ?

[erikvansebille]

Ah yes, that could be the reason indeed. If you are using the latest Parcels v3.0, you could also use the solution here

pfile.write_latest_locations(pset, time=endtime)

[martoconnh]

Hi!

I think it worked, as now I am not losing any particles, but it adds 2 more values to 'trajectory' in my output file. Everything is normal, execept these 2 values, which do not make any sense. This is how I implemented what you told me, and below there is the last values of 'trajectory' in my output file:

endtime = pset.time[-1]
output_file.write_latest_locations(pset, time=endtime)

array([... , 3483, 3484, 3485, 3486, 3487, 9223372036854775807, 9223372036854775807], dtype=int64)

So I do have all 3488 particles, but maybe I'm getting these last 2 numbers (actually the same number) because I wrongly implemented the code you told me earlier?