Delayed starts#

In many applications, it is needed to ‘delay’ the start of particle advection. For example because particles need to be released at different times throughout an experiment. Or because particles need to be released at a constant rate from the same set of locations.

This tutorial will show how this can be done. We start with importing the relevant modules.

[1]:
from datetime import timedelta as delta

import matplotlib.pyplot as plt
import numpy as np
import xarray as xr
from IPython.display import HTML
from matplotlib.animation import FuncAnimation

from parcels import (
    AdvectionRK4,
    FieldSet,
    JITParticle,
    ParticleFile,
    ParticleSet,
    Variable,
    download_example_dataset,
)

First import a FieldSet (from the Peninsula example, in this case)

[2]:
example_dataset_folder = download_example_dataset("Peninsula_data")
fieldset = FieldSet.from_parcels(
    f"{example_dataset_folder}/peninsula", allow_time_extrapolation=True
)

Now, there are two ways to delay the start of particles. Either by defining the whole ParticleSet at initialisation and giving each particle its own time. Or by using the repeatdt argument. We will show both options here

Assigning each particle its own time#

The simplest way to delaye the start of a particle is to use the time argument for each particle

[3]:
npart = 10  # number of particles to be released
lon = 3e3 * np.ones(npart)
lat = np.linspace(3e3, 45e3, npart, dtype=np.float32)

# release every particle one hour later
time = np.arange(0, npart) * delta(hours=1).total_seconds()

pset = ParticleSet(fieldset=fieldset, pclass=JITParticle, lon=lon, lat=lat, time=time)

Then we can execute the pset as usual

[4]:
output_file = pset.ParticleFile(name="DelayParticle_time.zarr", outputdt=delta(hours=1))

pset.execute(
    AdvectionRK4, runtime=delta(hours=24), dt=delta(minutes=5), output_file=output_file
)
INFO: Output files are stored in DelayParticle_time.zarr.
100%|██████████| 86400.0/86400.0 [00:02<00:00, 35968.90it/s]

And then finally, we can show a movie of the particles. Note that the southern-most particles start to move first.

[5]:
%%capture

ds = xr.open_zarr("DelayParticle_time.zarr")

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

ax.set_ylabel("Meridional distance [m]")
ax.set_xlabel("Zonal distance [m]")
ax.set_xlim(0, 9e4)
ax.set_ylim(0, 5e4)

timerange = np.unique(ds["time"].values[np.isfinite(ds["time"])])

# Indices of the data where time = 0
time_id = np.where(ds["time"] == timerange[0])

sc = ax.scatter(ds["lon"].values[time_id], ds["lat"].values[time_id])

t = str(timerange[0].astype("timedelta64[h]"))
title = ax.set_title(f"Particles at t = {t}")


def animate(i):
    t = str(timerange[i].astype("timedelta64[h]"))
    title.set_text(f"Particles at t = {t}")

    time_id = np.where(ds["time"] == timerange[i])
    sc.set_offsets(np.c_[ds["lon"].values[time_id], ds["lat"].values[time_id]])


anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)
[6]:
HTML(anim.to_jshtml())
[6]:

Using the repeatdt argument#

The second method to delay the start of particle releases is to use the repeatdt argument when constructing a ParticleSet. This is especially useful if you want to repeatedly release particles from the same set of locations.

[7]:
npart = 10  # number of particles to be released
lon = 3e3 * np.ones(npart)
lat = np.linspace(3e3, 45e3, npart, dtype=np.float32)
repeatdt = delta(hours=3)  # release from the same set of locations every 3h

pset = ParticleSet(
    fieldset=fieldset, pclass=JITParticle, lon=lon, lat=lat, repeatdt=repeatdt
)

Now we again define an output file and execute the pset as usual.

[8]:
output_file = pset.ParticleFile(name="DelayParticle_releasedt", outputdt=delta(hours=1))

pset.execute(
    AdvectionRK4, runtime=delta(hours=24), dt=delta(minutes=5), output_file=output_file
)
INFO: Output files are stored in DelayParticle_releasedt.zarr.
  0%|          | 0/86400.0 [00:00<?, ?it/s]WARNING: ParticleFile chunks are set to (10, 1), but this may lead to a significant slowdown in Parcels when many calls to repeatdt. Consider setting a larger chunk size for your ParticleFile (e.g. chunks=(int(1e4), 1)).
100%|██████████| 86400.0/86400.0 [00:01<00:00, 49894.21it/s]

And we get an animation where a new particle is released every 3 hours from each start location

[9]:
%%capture

ds = xr.open_zarr("DelayParticle_releasedt.zarr")

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

ax.set_ylabel("Meridional distance [m]")
ax.set_xlabel("Zonal distance [m]")
ax.set_xlim(0, 9e4)
ax.set_ylim(0, 5e4)

timerange = np.unique(ds["time"].values[np.isfinite(ds["time"])])

# Indices of the data where time = 0
time_id = np.where(ds["time"] == timerange[0])

sc = ax.scatter(ds["lon"].values[time_id], ds["lat"].values[time_id])

t = str(timerange[0].astype("timedelta64[h]"))
title = ax.set_title(f"Particles at t = {t}")


def animate(i):
    t = str(timerange[i].astype("timedelta64[h]"))
    title.set_text(f"Particles at t = {t}")

    time_id = np.where(ds["time"] == timerange[i])
    sc.set_offsets(np.c_[ds["lon"].values[time_id], ds["lat"].values[time_id]])


anim = FuncAnimation(fig, animate, frames=len(timerange), interval=100)
[10]:
HTML(anim.to_jshtml())
[10]: