OceanParcels / Parcels

Main code for Parcels (Probably A Really Computationally Efficient Lagrangian Simulator)
https://www.oceanparcels.org
MIT License
294 stars 135 forks source link

"Jumps" using different values of RK4_dt #1384

Closed luukrader closed 1 year ago

luukrader commented 1 year ago

Dear Erik or other reader,

In my notebook, I am faking an aquaplanet with a domain of 10S-89N and 90W-10W with a zonal velocity of 1 km/day only. Next, I am releasing particles in the domain of 10S-89N and 45W-40W with 1 degree in between them, as can be seen by the particleset. After the fieldset and particleset are written, I am executing the same experiment with the RK4 kernel, using 3 different timesteps of the kernel: 1 minute, 5 minutes and 30 minutes. Subsequently, the output data is called output_data_1, output_data_5 and output_data_30, representing the runs with RK4_dt = 1, 5 and 30 minutes respectively. Hereafter, the travelled distance after one day (which is zonal only) is computed in degrees and in kilometres and plotted.

As can be seen from the plots, the distance travelled by the particles as a function of latitude shows strange jumps when using smaller timesteps of the kernel, instead of showing a smooth profile.

What do you think might be the cause of these jumps? I will use OceanParcels to simulate plastics through (sub-)mesoscale dynamics. Therefore, I want to be able to choose a proper timestep for the RK4 kernel.

The Jupyter Notebook can be found via this link: https://github.com/luukrader/issue-OceanParcels/blob/main/Parcels_GH.ipynb

Thank you a lot in advance!

Cheers, Luuk

erikvansebille commented 1 year ago

Hi @luukrader, thanks for sharing this Issue and the notebook to reproduce it.

I checked it out this morning, and it seems an issue with floating-point accuracy of the particle.lon and particle.lat. The default accuracy for linear interpolation is np.float32 (see also the documentation here), but for very small time steps this is not precise enough.

The solution is quite simple, when you create your ParticleSet, just also specify

pset = ParticleSet.from_line(..., lonlatdepth_dtype=np.float64)

That will result in the figure below, where there is no more dependance on dt image

luukrader commented 1 year ago

Hi @erikvansebille,

That's very clear! Fun to see again how such a problem can be resolved with just one sentence of extra code :)

Thank you for the fast and clear answer, it helps me a lot.

Cheers, Luuk