What to do if particles are on land? about parcels HOT 5 CLOSED

I have been returning to this problem (and related out-of-domain issues) as it is becoming increasingly important for simulations involving diffusion or other, non-advection, forcing.

I'll attempt to summarise my understanding of our current functionality, to open discussion. Please correct me where if I am misunderstanding the code, particularly JIT mode where it's harder to debug the C-code converted kernels.
At present, the master branch:

• Interpolates sampled field values using interpolator2d
• If this interpolation involves reading values from outside the domain:
  • In Scipy mode this is caught, an error message printed, and a value of None returned for the interpolated value
  • In JIT mode, NaNs are returned, but no error or warning message is generated
• NaN values within the boundary demain are set to zero during initiation of individual Field objects

Is this an accurate summary?

This allows the user to check if sampled field values are None (for scipy) or nan (using math.isnan for JIT mode), and act accordingly for out-of-domain instances (e.g. calling particle.die() or perhaps moving the particle back to a previously stored position).

However, it seems that there is still no way of checking for presence on land masses within the model domain, aside from checking for zero field values within kernels (which in some instances may be a valid field value).
From the point of view of writing custom kernels, it would more useful to have a consistent NaN value returned, rather than a built in check-rejection test. It would allow users to, say, redraw the diffusion component of a time-step but keep the advection, for example.

I know NaNs cause problems with Scipy, but what options do we have here?

from parcels.

Ok, for this discussion I would like to point you to the ongoing work in the recovery-kernel-loop branch, which provides an explicit mechanism for dealing with "kernel errors", eg. "out of bounds" sampling or "NaN sampling". The idea here is to allow users to define "recovery kernels" that are run only once a specified error event occurs, and can range from throwing an error to "re-inserting" particles into the grid. Since the recovery mechanism for out-of-bounds and NaN cases can be specified separately, this should provide enough flexibility to do what you want, right?

The current implementation works for SciPy mode, but it does not yet expose the functionality through the execute() call. I will add this soon, including examples, so that we provide a choice for user how to deal with these special cases.

from parcels.

This sounds great @mlange05 ! I didn't realise this branch existed. I will take a proper look and try to begin implementing some recovery kernels for scipy in anticipation of the finished product. All my work arounds with the current master are becoming rather messy, but it's kind of critical for our applied tuna paper.

from parcels.

We've written some 'unbeaching kernels', see e.g. https://github.com/OceanParcels/Parcelsv2.0PaperNorthSeaScripts/blob/master/northsea_mp_kernels.py which is associated with the v2.0 paper in GMD

from parcels.

@erikvansebille and @Jacketless, are all NaNs converted to zero? I was hoping to recognize cells = NaN as beached (inside a kernel) and unbeach accordingly. Instead, I made all of my "land" cells = 0.00; but then there is the chance of actual 0 velocity occurring and being interpreted as beaching... Furthermore, I think I am running into spatial interpolation issues where the outermost land cells are not being recognized as land (since their interpolated velocity is > or < 0) and I end up with "longshore flow" on land. Is that possible?

from parcels.