What to do if particles are on land?

[erikvansebille]

We need to decide how to treat particles that are on land. @Abobie has worked hard on the boundaries branch where particles are 'pushed back' using ghost points, when they come very close to land. However, there can still be instances where particles end up on land. For example when the timestep is so large that particles step over multiple grid cells in one timestep. Or if random-walk diffusion is added

My thinking is that at the end of the kernel, we need a function that, for each timestep, assesses whether a particle is on land (and whether it has gone out of the domain, see #47). If that function return True, the user should have a few options:

1. Remove the particle (to simulate beaching)
2. Try computing the particle position again (essentially redrawing the random number in the diffusion; this should probably be defualt behaviour)
3. Move the particle to the closest ocean point (might be difficult/expensive to compute, though)

@Jacketless, you have been thinking about this too for your diffusion steps, right? Any comments?

And @mlange05 what do you think about a general test at the end of each particle timestep to accept/reject the move in that timestep? Would that be feasible?

[Jacketless]

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?

[mlange05]

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.

[Jacketless]

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.

[erikvansebille]

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

[hkhirsh]

@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?