Best pratice to modify p.dt inside a kernel ?

[quai20]

Hello all !

Within my custom kernel, I would like to change particle.dt because in some part of the advection, I want to refine or enlarge the sampling rate. I see that I can do particle.update_next_dt(new_dt) but :

• does it apply to all the kernels ? (I have kernels = custom_kernel + AdvectionRK4 + periodicBC)
• How do I set the file writing ? because my dt is not constant, so what about the ParticleFile outputdt ? Is it possible to write all the data, when irregularly sampled ?

I tried combining particle.update_next_dt(new_dt) in my kernel, and leaving outputdt empty but I only get the first and last values of the simulation in my output file, butI think it's normal with the default outputdt=np.infty.

I you have any recommandation of example of varying, thanks very much.

Kevin

[erikvansebille]

Dear @quai20. Yes, I think that what you are doing is fine, and that Parcels should be able to indeed adapt dt on the fly. outputs will be respected no matter how you change dt, I think. But there's only one way to find out I guess, by just trying and see what happens?

[quai20]

Yes I'm in the process of doing that 😄 . I managed to bypass the outputdt of particleFile and outputing every dt in my file but I didnt succeed (yet) with the particle.update_next_dt(new_dt).
When printing dt within the pset.execute(), it does not change.

[erikvansebille]

Hmm, maybe @CKehl can help and explain what might be going on?

[CKehl]

Hi,

1. you have here 3 different variables: dt set globally and applicable to all particles within the pset.execute(...) main loop. This is set external and cannot be modified inside the kernel. You have outputdt, which is not connected to dt and cannot be influenced from inside the kernel. What you can influence in the kernel is the temporary particle.dt of each individual particle (as the kernel applies to each particle individually).
2. This particle.dt is updated with the method you describe, but is not written back to pset.execute(...) - obviously, because all particles may have different dt's then, and thus setting the global dt from the set of particle.dt is undefined.
3. The essential reason for even having that update_next_dt(...) is that when particles cross a cell boundaries during advection (i.e. during one dt), the kernel is stopped at t * particle.dt, the new field cell is loaded, and the advection continues then with the remaining (1-t) * particle.dt. Hence, adapting that particle.dt is only there to implicitly split a larger dt into smaller ones. Once the kernel is executed for all particles for 1 whole dt, it returns to pset.execute(...), and in the next timestep begins.
4. dt inside pset.execute(...) is consistent for all particles, but not constant. The dt is the difference between time and the next timestep, which is the moment of the next action (e.g. compute, input, plot, output writing, etc.).
5. in short, what you want isn't easily feasible unless sampling here refers to internal advection steps (i.e. advection steps that are not written to file). There is of course the possibility to delete(...) a particle at each sub-step in your advection, or cause a TimeInterpolation error in your kernel, for which can you then have a recovery procedure where you save the particle state, delete the particle and re-add the particle again. Note that this will dramatically slow down the simulation. Another option: do your actual time advance and advection outside pset.execute(...) in an external loop.
6. word of caution: auto-adapting particle.dt to very small or very large time steps (i.e. your particle.dt drops below 10^-6 or rising beyond 10^8) will cause numerical instability because the unit for time steps is seconds stored in float64 and there is a limit on how many floating point numbers you can represent with it. example: a simulation over one year has 31536000 (3.1 x 10^7) seconds, if your particle.dt is 1 x 10^-6, your numbers need to cover a range 10^14 - that is a range of 2^46. The hard max. one can represent is a range of 10^19, and that also only if your dt is a rational number in the base-2 system. Exceeding that numerical precision limit can result in seemingly-infinite loops cause the particles get stuck updating their own dt because particle.dt becomes non-monotonous.