[OUTDATED; DO NOT MERGE] Vectorial summing of particle displacement at end of kernel loop

[erikvansebille]

(note: the PR is now superseded by #1402, which was reset to caf7418 (before the merging of #1399), because Parquet and sqlite writing had in practice far worse performance to zarr)

Until now, the way that Kernel-concatenation worked in Parcels meant that the execution of multiple Kernels was generally not commutable: for example first advection a particle and then letting it sink would lead to a different result than first sinking and then advecting. That is because particle positions were updated within each kernel, instead of the vectorial summing of displacements that would be more appropriate (and would ensure commutability).

More problematically, while the particle positions would be updated within a kernel, the time was generally only updated at the end. This meant that Field-sampling could give unexpected results (e.g. when sampling was done at the new location but the old time).

To solve these issues, in this PR I propose that individual Kernels compute displacements (particle_dlon, particle_dlat and particle_ddepth) and that these displacements are only added to the particle locations when all kernels have been executed for a timestep dt.

This means that the workflow for each particle p and each time t becomes

1. Set particle_dlon = particle_dlat = particle_ddepth = 0
2. Execute all Kernels, including sampling Kernels, and let the Kernels update particle_dlon, particle_dlat and particle_ddepth when required
3. Write the location of the particle p at time t (if needed)
4. Update time to t += dt and particle locations to particle.lon += particle_dlon etc

Note that this workflow means that particles will be written at t=0 (and there is no need anymore for a separate initial values sampling step!), but not anymore at t=runtime (the last tilmestep), unless we also change that the particle loop is executed one time more at the end (tbd).

Issues/questions that need to be explored or implemented:

• a. Move the writing step to after the Kernel execution but before the location and time update
• b. Investigate whether it's wise to perform the loop standard one extra time to also write output at t=runtime (see comment above)
• c. Check if the particles_backups are stil needed now that particles are not moved in a Kernel
• d. Check whether RK45 still works
• e. Explore how to deal with Error Kernels
• f. Decide if the SummedFields class is still necessary
• g. Update this new workflow in documentation and tutorials

Since this is a major overhaul that can change results, I suggest to make this PR (when merged) part of v2.5.0

[erikvansebille]

I've done some exploration of how best to update the kernels. See this gist that explores this on a very minimal example ('nanoparcels')

Bottom line is that the cleanest working implementation (cell 6 in the gist) is something like (only for dlon, similar for other two directions)

def concatenated_kernels(particle, fieldset, time):
    # set dlon to zero at beginning of loop
    dlon = 0

    # sample fields at (time, location), and add to dlon when particle moves
    particle.p = fieldset.P(particle.time, particle.lon)
    dlon += fieldset.U(particle.time, particle.lon) * particle.dt
    
    # write the particle if needed (depending on outputdt)
    particle.write()

    # update the particle time and particle location
    particle.time += particle.dt
    particle.lon += dlon

So this requires a particle.write() function that can be called in a Kernel, something that may be feasible using the Parquet library. This is now being explored in #1399

[michaeldenes]

I've been playing around with his a bit today. The obvious major change is that custom kernels will need to compute displacements, rather than positions. That's easy enough (just subtract the current position particle.lon or particle.lat).

One minor issue is that an unbeaching kernel is still not commutable. You must first perform advection and custom kernels, and only then identify the updated position of the particle (particle.lon + particle_dlon for e.g.) and apply an unbeaching kernel.

For custom kernels that set custom variables (e.g. a variable to track if a particle goes below 500m), I think this would be even more prominent. Should you check an 'if statement' before any movement, and set the custom variable at the current timestep (based on the previous timestep position), or check the 'if statement' after all displacements and set at the current timestep? Any thoughts?

[erikvansebille]

Thanks @michaeldenes, for your feedback on this Draft PR. You're right that actions like unteaching will still not be commutable; kernels where the action depends on position will generally have that problem. The point is that with this PR some can commute. It also means we might be able to get rid of the SummedFields (which was essentially an implementation of vectorial summing for Advection only).

I think/hope that actions like unteaching could be treated in Error/RecoveryKernels, which need a major upgrade anyways. We still use the RecoveryKernel code from v0.9 essentially unchanged, and it is very slow. That's point e. in the list at the top of this page. But I'm first focussing on point a, since that's a sine qua non...

[erikvansebille]

Closing this PR as #1402 is now the working implementation of this idea