Particle-particle interaction example

In this notebook, we show an example of the new 'particle-particle-interaction' functionality in Parcels. Note that this functionality is still in development, and the implementation is fairly rudimentary and slow for now. Importantly:

  • Particle-particle interaction only works in Scipy mode
  • The type of interactions that are supported is still limited

Interactions are implemented through InteractionKernels, which are similar to normal Kernels but with a slightly extended signature:

def InteractionKernel(particle, fieldset, time, neighbors, mutator)

where the last two arguments are different from normal Kernel signature.

The neighbors argument provided a list of the particles that are within a neighborhood (defined by the new interaction_distance argument in ParticleSet creation).

The mutator argument is an initially empty list with all the mutations that need to be performed on particles at the end of running the InteractionKernel on all particles. Such mutations can typically be moving a particle, but we can't do that 'on-the-fly' as the neighburs-list could then change and the result of the InteractionKernel might depend on the order of the particles list.

The two other changes for users are the above mentioned interaction_distance argument in ParticleSet creation and the pyfunc_inter argument in ParticleSet.execute().

Pulling particles

Below is an example of what can be done with particle-particle interaction. We create a square grid of $N\times N$ particles, which are all subject to Brownian Motion (via the built-in DiffusionUniformKh Kernel). Furthermore, some of the particles also 'attract' other particles that are within the interaction distance: these attracted particles move with a constant velocity to the attracting particles.

In [1]:
%matplotlib notebook
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt

from parcels import FieldSet
from parcels import ParticleSet
from parcels import Variable
from parcels import ScipyParticle, ScipyInteractionParticle
from parcels import DiffusionUniformKh
from parcels import NearestNeighborWithinRange, MergeWithNearestNeighbor

from matplotlib.animation import FuncAnimation
from IPython.display import HTML
INFO: Compiled ParcelsRandom ==> C:\Users\GEBRUI~1\AppData\Local\Temp\parcels-tmp\
In [4]:
def Pull(particle, fieldset, time, neighbors, mutator):
    '''InterActionKernel that "pulls" all neighbor particles 
    toward the attracting particle with a constant velocity '''
    distances = []
    na_neighbors = []
    if not particle.attractor:  # only execute kernel for particles that are 'attractor'
        return StateCode.Success
    for n in neighbors:
        if n.attractor:
        x_p = np.array([particle.depth,, particle.lon])
        x_n = np.array([ n.depth,, n.lon])
        distances.append(np.linalg.norm(x_p-x_n))  # compute distance between attracted and attracting particle

    velocity = 0.04  # predefined attracting velocity
    for n in na_neighbors:
        assert n.dt == particle.dt
        dx = np.array([, particle.lon-n.lon,
        dx_norm = np.linalg.norm(dx)

        distance = velocity*n.dt
        d_vec = distance*dx/dx_norm  # calculate vector of position change

        def f(n, dlat, dlon, ddepth):  # define mutation function for mutator
   += dlat
            n.lon += dlon
            n.depth += ddepth

        mutator[].append((f, d_vec))  # add mutation to the mutator
In [5]:
npart = 11

X, Y = np.meshgrid(np.linspace(-1, 1, npart), np.linspace(-1, 1, npart))

# Define a fieldset without flow
fieldset = FieldSet.from_data({'U': 0, 'V': 0}, {'lon': 0, 'lat': 0}, mesh='flat')
fieldset.add_constant_field("Kh_zonal", 0.0005, mesh="flat")
fieldset.add_constant_field("Kh_meridional", 0.0005, mesh="flat")

# Create custom particle class with extra variable that indicates
# whether the interaction kernel should be executed on this particle.
class InteractingParticle(ScipyParticle):
    attractor = Variable('attractor', dtype=np.bool_, to_write='once')

attractor = [True if i in [int(npart*npart/2-3), int(npart*npart/2+3)] else False for i in range(npart*npart)]
pset = ParticleSet(fieldset=fieldset, pclass=InteractingParticle,
                   lon=X, lat=Y,
                   interaction_distance=0.5,  # note the interaction_distance argument here

output_file = pset.ParticleFile(name="", outputdt=1)

             pyfunc_inter=Pull,  # note the pyfunc_inter here
             runtime=60, dt=1, output_file=output_file)

In [6]:
data_xarray = xr.open_dataset('')
data_attr = data_xarray.where(data_xarray['attractor']==1, drop=True)
data_other = data_xarray.where(data_xarray['attractor']==0, drop=True)

timerange = np.arange(np.nanmin(data_xarray['time'].values),
                      np.nanmax(data_xarray['time'].values), np.timedelta64(1, 's'))  # timerange in nanoseconds

fig = plt.figure(figsize=(4, 4),constrained_layout=True)
ax = fig.add_subplot()

ax.set_ylabel('Meridional distance [m]')
ax.set_xlabel('Zonal distance [m]')
ax.set_xlim(-1.1, 1.1)
ax.set_ylim(-1.1, 1.1)

time_id = np.where(data_other['time'] == timerange[0]) # Indices of the data where time = 0
time_id_attr = np.where(data_attr['time'] == timerange[0]) # Indices of the data where time = 0

scatter = ax.scatter(data_other['lon'].values[time_id], data_other['lat'].values[time_id], c='b', s=5,zorder=1)
scatter_attr = ax.scatter(data_attr['lon'].values[time_id_attr], data_attr['lat'].values[time_id_attr], c='r', s=40,zorder=2)

circs = []
for lon_a, lat_a in zip(data_attr['lon'].values[time_id_attr], data_attr['lat'].values[time_id_attr]):
    circs.append(ax.add_patch(plt.Circle((lon_a, lat_a), 0.25, facecolor='None', edgecolor='r', linestyle='--')))

t = str(timerange[0].astype('timedelta64[s]'))
title = ax.set_title('Particles at t = '+t+' (Red particles are attractors)')

def animate(i):
    t = str(timerange[i].astype('timedelta64[s]'))
    title.set_text('Particles at t = '+t+'\n (Red particles are attractors)')
    time_id = np.where(data_other['time'] == timerange[i])
    time_id_attr = np.where(data_attr['time'] == timerange[i])
    scatter.set_offsets(np.c_[data_other['lon'].values[time_id], data_other['lat'].values[time_id]])
    scatter_attr.set_offsets(np.c_[data_attr['lon'].values[time_id_attr], data_attr['lat'].values[time_id_attr]])
    for c, lon_a, lat_a in zip(circs, data_attr['lon'].values[time_id_attr], data_attr['lat'].values[time_id_attr]): = (lon_a, lat_a)
    return scatter, scatter_attr, circs,
anim = FuncAnimation(fig, animate, frames = len(timerange), interval=500, blit=True)
In [7]: