Index search errors in NEMO-based run.

[christensen5]

Hi all. I am in the process of running a backwards-time simulation with 3D velocities obtained from a NEMO run. I am using the Parcels default RK4 3D kernel together with a sampling kernel (that simply samples the value of a temperature and a nutrient field at each timestep) and two recovery kernels (to resubmerge particles if they breach the surface and to delete particles if they run aground). During the run, I am (very infrequently) seeing the following error:

Correct cell not found for (-25588126015134105600.000000, -76.026820) after 1000000 iterations Debug info: old particle indices: (yi, xi) 45 1081 new particle indices: (yi, xi) 50 3552 Mesh 2d shape: 3059 4322 Relative particle position: (xsi, eta) -3.0705868351752969e+20 1.5449192168699499e-01

When this does occur, it seems to be for extremely high latitude particles (below 70 deg South, where I have a small number of particles in my run). I note the extremely negative xsi value (always on the order of -1e+20) in the debug info - maybe particles are being queried many cells to the "West" of where they were expected?

I am not inordinately worried, since it's not crashing the simulation, but I am curious about why this might be happening. As far as I can tell this error message originates in the reconnect_bound_indices method from index_search.h, but I am not familiar enough with Parcels' index search approach to know at a glance what is going on.

[christensen5]

Correction: it appears that this error is not exclusive to high-latitude particles. I have now observed it for equatorial particles too:

Correct cell not found for (-6291204100168089600.000000, -13.054622) after 1000000 iterations Debug info: old particle indices: (yi, xi) 1337 3930 new particle indices: (yi, xi) 1335 2080 Mesh 2d shape: 3059 4322 Relative particle position: (xsi, eta) -7.5492145364573094e+19 4.7120297077680956e-01

[reint-fischer]

Hi @christensen5,

Unfortunately, we have seen these errors on a few occansions before, especially with curvilinear grids and in high-latitude regions, but have not yet discovered what causes the index_search to fail. The error should originate at least from the search_indices_curvilinear method, and specifically appear here.

What is your idea for why it should originate in the reconnect_bnd_indices method? I believe that only affects xi when it is near 0 or xdim, which in your case seems to be 4322.

What caught my attention is that both large 'lon' values you get are near the order of magnitude of 2**64, which can be stored in a float64 or double. However, I am not sure yet how that would relate to the error at the moment.

[christensen5]

index_search_errors.csv
Hi @reint-fischer,

Thanks for this info. Looking again, you're quite right about the origin of the error message, I think I must have misread the call to reconnect_bnd_indices() in index_search.h line 357 as a definition.

I must say I don't really have a hypothesis for why this is happening in my case...in my current simulation involving about 17 million particles on a NEMO curvilinear grid, I have noticed the error about 20 times in the current run (about 10% complete). I attach a csv tabulating the debug info for each of these cases, in case it's helpful or interesting. The only thing I notice that all cases have in common is that the "new" xi index is always either 1283 or 1284. Once the run is finished I'll pull up the NEMO velocity fields and try to see if there are any strange values around there (I don't have enough memory spare to do this while the simulation is running unfortunately).

Regarding your observation about the 'lon' values and 2^64, all the cases I've seen so far have lon values in the range 2^61.1 - 2^64.5 . I'm afraid I don't have any more concrete idea than you at the moment about whether that is relevant or not.

I don't know how urgent this error is to the Parcels team - if it would be helpful I can try and create a MWE that reproduces the error without requiring 17 million particles or 100s of GB of NEMO velocity fields?

[erikvansebille]

Well, x=1283 is exactly the location of the dateline in the NEMO file (e.g. xr.load_dataset('means/ORCA0083-N06_20010105d05U.nc').nav_lon.plot(), see figure below
.

So this could very well be a dateline issue?

[christensen5]

Hi Erik, that looks like a rather probable explanation. My C isn't especially advanced but I can try and work through the logic of search_indices_curvilinear for a particle crossing the dateline to see if this would prevent the while loop therein from terminating, and thereby trigger the max iteration limit?