{
"cells": [
{
"attachments": {},
"cell_type": "markdown",
"id": "6ad92baa",
"metadata": {},
"source": [
"# Particle-Particle interaction\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "fe5d7b1c",
"metadata": {},
"source": [
"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:\n",
"\n",
"- Particle-particle interaction only works in Scipy mode\n",
"- The type of interactions that are supported is still limited\n",
"\n",
"Interactions are implemented through `InteractionKernels`, which are similar to normal `Kernels`. The `InteractionKernels` are applied between particles that are located closer to each other than a specified `interaction_distance`. In general, the code structure needs three adaptations to apply particle-particle interaction:\n",
"\n",
"1. The `ParticleSet` requires an `interaction_distance` argument upon creation, to define the `interaction_distance`.\n",
"2. `ParticleSet.execute()` requires the `pyfunc_inter` argument, which contains the `InteractionKernels` that will be executed, similarly to the `pyfunc` argument for normal `Kernels`.\n",
"3. `InteractionKernels` have two additional arguments compared to normal `Kernels`:\n",
"\n",
"```python\n",
"def InteractionKernel(particle, fieldset, time, neighbors, mutator)\n",
"```\n",
"\n",
"The `neighbors` argument provides a list of the particles that are within a neighborhood (i.e. closer than the `interaction_distance` argument in `ParticleSet` creation).\n",
"\n",
"The `mutator` argument is an initially empty list with all the mutations that need to be performed on particles at the end of running all `InteractionKernels` on all particles.\n",
"This `mutator` argument is required, because otherwise the order at which interactions are applied has implications for the simulation. As a consequence, the simulation will likely be dependent on the order of the particle list if no mutator list is used.\n"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "66b91c29",
"metadata": {},
"source": [
"## Pulling particles\n",
"\n",
"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.\n"
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "da408a23",
"metadata": {},
"outputs": [],
"source": [
"%matplotlib notebook\n",
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"import xarray as xr\n",
"from IPython.display import HTML\n",
"from matplotlib.animation import FuncAnimation\n",
"\n",
"from parcels import (\n",
" DiffusionUniformKh,\n",
" FieldSet,\n",
" MergeWithNearestNeighbor,\n",
" NearestNeighborWithinRange,\n",
" ParticleSet,\n",
" ScipyInteractionParticle,\n",
" ScipyParticle,\n",
" Variable,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "9485e2e9",
"metadata": {},
"outputs": [],
"source": [
"def Pull(particle, fieldset, time, neighbors, mutator):\n",
" \"\"\"InterActionKernel that \"pulls\" all neighbor particles\n",
" toward the attracting particle with a constant velocity\"\"\"\n",
" distances = []\n",
" na_neighbors = []\n",
" # only execute kernel for particles that are 'attractor'\n",
" if not particle.attractor:\n",
" return StateCode.Success\n",
" for n in neighbors:\n",
" if n.attractor:\n",
" continue\n",
" x_p = np.array([particle.depth, particle.lat, particle.lon])\n",
" x_n = np.array([n.depth, n.lat, n.lon])\n",
"\n",
" # compute distance between attracted and attracting particle\n",
" distances.append(np.linalg.norm(x_p - x_n))\n",
" na_neighbors.append(n)\n",
"\n",
" velocity = 0.04 # predefined attracting velocity\n",
" for n in na_neighbors:\n",
" assert n.dt == particle.dt\n",
" dx = np.array(\n",
" [particle.lat - n.lat, particle.lon - n.lon, particle.depth - n.depth]\n",
" )\n",
" dx_norm = np.linalg.norm(dx)\n",
"\n",
" # calculate vector of position change\n",
" distance = velocity * n.dt\n",
" d_vec = distance * dx / dx_norm\n",
"\n",
" # define mutation function for mutator\n",
" def f(n, dlat, dlon, ddepth):\n",
" n.lat_nextloop += (\n",
" dlat # note that we need to change the locations for the next loop\n",
" )\n",
" n.lon_nextloop += dlon\n",
" n.depth_nextloop += ddepth\n",
"\n",
" # add mutation to the mutator\n",
" mutator[n.id].append((f, d_vec))"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "451638ae",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"INFO: Output files are stored in InteractingParticles.zarr.\n",
" 0%| | 0/60.0 [00:00 /var/folders/1n/500ln6w97859_nqq86vwpl000000gr/T/parcels-504/parcels_random_fec11e5e-713b-43cc-b561-8d45de86d6f2.c\n",
"WARNING: Some InteractionKernel was not completed succesfully, likely because a Particle threw an error that was not captured.\n",
"100%|██████████| 60.0/60.0 [00:02<00:00, 28.27it/s]\n"
]
}
],
"source": [
"npart = 11\n",
"\n",
"X, Y = np.meshgrid(np.linspace(-1, 1, npart), np.linspace(-1, 1, npart))\n",
"\n",
"# Define a fieldset without flow\n",
"fieldset = FieldSet.from_data({\"U\": 0, \"V\": 0}, {\"lon\": 0, \"lat\": 0}, mesh=\"flat\")\n",
"fieldset.add_constant_field(\"Kh_zonal\", 0.0005, mesh=\"flat\")\n",
"fieldset.add_constant_field(\"Kh_meridional\", 0.0005, mesh=\"flat\")\n",
"\n",
"\n",
"# Create custom particle class with extra variable that indicates\n",
"# whether the interaction kernel should be executed on this particle.\n",
"InteractingParticle = ScipyParticle.add_variable(\n",
" \"attractor\", dtype=np.bool_, to_write=\"once\"\n",
")\n",
"\n",
"\n",
"attractor = [\n",
" True if i in [int(npart * npart / 2 - 3), int(npart * npart / 2 + 3)] else False\n",
" for i in range(npart * npart)\n",
"]\n",
"pset = ParticleSet(\n",
" fieldset=fieldset,\n",
" pclass=InteractingParticle,\n",
" lon=X,\n",
" lat=Y,\n",
" interaction_distance=0.5, # note the interaction_distance argument here\n",
" attractor=attractor,\n",
")\n",
"\n",
"output_file = pset.ParticleFile(name=\"InteractingParticles.zarr\", outputdt=1)\n",
"\n",
"pset.execute(\n",
" pyfunc=DiffusionUniformKh,\n",
" pyfunc_inter=Pull, # note the pyfunc_inter here\n",
" runtime=60,\n",
" dt=1,\n",
" output_file=output_file,\n",
")"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "36c8f905",
"metadata": {},
"outputs": [],
"source": [
"%%capture\n",
"data_xarray = xr.open_zarr(\"InteractingParticles.zarr\")\n",
"data_attr = data_xarray.where(data_xarray[\"attractor\"].compute() == 1, drop=True)\n",
"data_other = data_xarray.where(data_xarray[\"attractor\"].compute() == 0, drop=True)\n",
"\n",
"timerange = np.arange(\n",
" np.nanmin(data_xarray[\"time\"].values),\n",
" np.nanmax(data_xarray[\"time\"].values),\n",
" np.timedelta64(1, \"s\"), # timerange in nanoseconds\n",
")\n",
"\n",
"fig = plt.figure(figsize=(4, 4), constrained_layout=True)\n",
"ax = fig.add_subplot()\n",
"\n",
"ax.set_ylabel(\"Meridional distance [m]\")\n",
"ax.set_xlabel(\"Zonal distance [m]\")\n",
"ax.set_xlim(-1.1, 1.1)\n",
"ax.set_ylim(-1.1, 1.1)\n",
"\n",
"time_id = np.where(\n",
" data_other[\"time\"] == timerange[0]\n",
") # Indices of the data where time = 0\n",
"time_id_attr = np.where(\n",
" data_attr[\"time\"] == timerange[0]\n",
") # Indices of the data where time = 0\n",
"\n",
"scatter = ax.scatter(\n",
" data_other[\"lon\"].values[time_id],\n",
" data_other[\"lat\"].values[time_id],\n",
" c=\"b\",\n",
" s=5,\n",
" zorder=1,\n",
")\n",
"scatter_attr = ax.scatter(\n",
" data_attr[\"lon\"].values[time_id_attr],\n",
" data_attr[\"lat\"].values[time_id_attr],\n",
" c=\"r\",\n",
" s=40,\n",
" zorder=2,\n",
")\n",
"\n",
"circs = []\n",
"for lon_a, lat_a in zip(\n",
" data_attr[\"lon\"].values[time_id_attr], data_attr[\"lat\"].values[time_id_attr]\n",
"):\n",
" circs.append(\n",
" ax.add_patch(\n",
" plt.Circle(\n",
" (lon_a, lat_a), 0.25, facecolor=\"None\", edgecolor=\"r\", linestyle=\"--\"\n",
" )\n",
" )\n",
" )\n",
"\n",
"t = str(timerange[0].astype(\"timedelta64[s]\"))\n",
"title = ax.set_title(\"Particles at t = \" + t + \" (Red particles are attractors)\")\n",
"\n",
"\n",
"def animate(i):\n",
" t = str(timerange[i].astype(\"timedelta64[s]\"))\n",
" title.set_text(\"Particles at t = \" + t + \"\\n (Red particles are attractors)\")\n",
"\n",
" time_id = np.where(data_other[\"time\"] == timerange[i])\n",
" time_id_attr = np.where(data_attr[\"time\"] == timerange[i])\n",
" scatter.set_offsets(\n",
" np.c_[data_other[\"lon\"].values[time_id], data_other[\"lat\"].values[time_id]]\n",
" )\n",
" scatter_attr.set_offsets(\n",
" np.c_[\n",
" data_attr[\"lon\"].values[time_id_attr], data_attr[\"lat\"].values[time_id_attr]\n",
" ]\n",
" )\n",
" for c, lon_a, lat_a in zip(\n",
" circs,\n",
" data_attr[\"lon\"].values[time_id_attr],\n",
" data_attr[\"lat\"].values[time_id_attr],\n",
" ):\n",
" c.center = (lon_a, lat_a)\n",
" return (\n",
" scatter,\n",
" scatter_attr,\n",
" circs,\n",
" )\n",
"\n",
"\n",
"anim = FuncAnimation(fig, animate, frames=len(timerange), interval=200, blit=True)\n",
"data_xarray.close()"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "cc7ced7f",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"\n",
"\n",
"\n",
"\n",
"\n",
"\n",
"
\n",
"\n",
"\n",
"\n"
],
"text/plain": [
""
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"HTML(anim.to_jshtml())"
]
},
{
"attachments": {},
"cell_type": "markdown",
"id": "cbaf5190",
"metadata": {},
"source": [
"## Interacting with the FieldSet\n",
"\n",
"An important feature of Parcels is to evaluate a `Field` at the `Particle` location using the square bracket notation: `particle.Temperature = fieldset.T[time, depth, lat, lon]`. These types of particle-field interactions are recommended to be implemented in standard `Kernels`, since the `InteractionKernels` do not report the `StateCodes` that are used to flag particles that encounter an error in the particle-field interaction, e.g. `OutOfBoundsError`. Any variable that is needed in the particle-particle interaction can be stored in a `Variable` by sampling the field in a `Kernel` before executing the `InteractionKernel`.\n"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.11.6"
}
},
"nbformat": 4,
"nbformat_minor": 5
}