|
| 1 | +--- |
| 2 | +file_format: mystnb |
| 3 | +kernelspec: |
| 4 | + name: python3 |
| 5 | +--- |
| 6 | + |
| 7 | +# 📖 Parcels conceptual workflow |
| 8 | + |
| 9 | +Parcels is a set of Python classes and methods to create customisable particle tracking simulations using gridded output from (ocean) circulation models. |
| 10 | + |
| 11 | +Here, we will explain the most important classes and functions. This overview can be useful to start understanding the different components we use in Parcels, and to structure the code in a simulation script. |
| 12 | + |
| 13 | +A Parcels simulation is generally built up from four different components: |
| 14 | + |
| 15 | +1. [**FieldSet**](#1-fieldset). The input dataset of gridded fields (e.g. ocean current velocity, temperature) in which virtual particles are defined. |
| 16 | +2. [**ParticleSet**](#2-particleset). The dataset of virtual particles. These always contain time, z, lat, and lon, for which initial values must be defined. The ParticleSet may also contain other, custom variables. |
| 17 | +3. [**Kernels**](#3-kernels). Kernels perform some specific operation on the particles every time step (e.g. advect the particles with the three-dimensional flow; or interpolate the temperature field to the particle location). |
| 18 | +4. [**Execute**](#4-execute). Execute the simulation. The core method which integrates the operations defined in Kernels for a given runtime and timestep, and writes output to a ParticleFile. |
| 19 | + |
| 20 | +We discuss each component in more detail below. The subsections titled **"Learn how to"** link to more detailed [how-to guide notebooks](../user_guide/index.md) and more detailed _explanations_ of Parcels functionality are included under **"Read more about"** subsections. The full list of classes and methods is in the [API reference](../reference.md). If you want to learn by doing, check out the [quickstart tutorial](./tutorial_quickstart.md) to start creating your first Parcels simulation. |
| 21 | + |
| 22 | +```{figure} ../_static/concepts_diagram.png |
| 23 | +:alt: Parcels concepts diagram |
| 24 | +:width: 100% |
| 25 | +
|
| 26 | +Parcels concepts diagram with key classes in blue boxes |
| 27 | +``` |
| 28 | + |
| 29 | +## 1. FieldSet |
| 30 | + |
| 31 | +Parcels provides a framework to simulate particles **within a set of fields**, such as flow velocities and temperature. To start a parcels simulation we must define this dataset with the **`parcels.FieldSet`** class. |
| 32 | + |
| 33 | +The input dataset from which to create a `parcels.FieldSet` can be an [`xarray.Dataset`](https://docs.xarray.dev/en/stable/user-guide/data-structures.html#dataset) with output from a hydrodynamic model or reanalysis. Such a dataset usually contains a number of gridded variables (e.g. `"U"`), which in Parcels become `parcels.Field` objects. A list of `parcels.Field` objects is stored in a `parcels.FieldSet` in an analoguous way to how `xarray.DataArray` objects combine to make an `xarray.Dataset`. |
| 34 | + |
| 35 | +For several common input datasets, such as the Copernicus Marine Service analysis products, Parcels has a specific method to read and parse the data correctly: |
| 36 | + |
| 37 | +```python |
| 38 | +dataset = xr.open_mfdataset("insert_copernicus_data_files.nc") |
| 39 | +fieldset = parcels.FieldSet.from_copernicusmarine(dataset) |
| 40 | +``` |
| 41 | + |
| 42 | +In some cases, we might want to combine `parcels.Field`s from different sources in the same `parcels.FieldSet`, such as ocean currents from one dataset and Stokes drift from another. This is possible in Parcels by adding each `parcels.Field` separately: |
| 43 | + |
| 44 | +```python |
| 45 | +dataset1 = xr.dataset("insert_current_data_files.nc") |
| 46 | +dataset2 = xr.dataset("insert_stokes_data_files.nc") |
| 47 | + |
| 48 | +Ucurrent = parcels.Field(name="Ucurrent", data=dataset1["Ucurrent"], grid=parcels.XGrid.from_dataset(dataset1), interp_method=parcels.interpolators.XLinear) |
| 49 | +Ustokes = parcels.Field(name="Ustokes", data=dataset2["Ustokes"], grid=parcels.XGrid.from_dataset(dataset2), interp_method=parcels.interpolators.XLinear) |
| 50 | + |
| 51 | +fieldset = parcels.FieldSet([Ucurrent, Ustokes]) |
| 52 | +``` |
| 53 | + |
| 54 | +### Grid |
| 55 | + |
| 56 | +Each `parcels.Field` is defined on a grid. With Parcels, we can simulate particles in fields on both structured (**`parcels.XGrid`**) and unstructured (**`parcels.UxGrid`**) grids. The grid is defined by the coordinates of grid cell nodes, edges, and faces. `parcels.XGrid` objects are based on [`xgcm.Grid`](https://xgcm.readthedocs.io/en/latest/grids.html), while `parcels.UxGrid` objects are based on [`uxarray.Grid`](https://uxarray.readthedocs.io/en/stable/generated/uxarray.Grid.html#uxarray.Grid) objects. |
| 57 | + |
| 58 | +```{admonition} 📖 Read more about grids |
| 59 | +:class: seealso |
| 60 | +- [Grids explanation](../user_guide/examples/explanation_grids.md) |
| 61 | +``` |
| 62 | + |
| 63 | +### Interpolation |
| 64 | + |
| 65 | +To find the value of a `parcels.Field` at any particle location, Parcels interpolates the gridded field. Depending on the variable, grid, and required accuracy, different interpolation methods may be appropriate. Parcels comes with a number of built-in **`parcels.interpolators`**. |
| 66 | + |
| 67 | +```{admonition} 📖 Read more about interpolation |
| 68 | +:class: seealso |
| 69 | +- [Interpolation explanation](../user_guide/examples/explanation_interpolation.md) |
| 70 | +``` |
| 71 | + |
| 72 | +```{admonition} 🖥️ Learn how to use Parcels interpolators |
| 73 | +:class: seealso |
| 74 | +- [Interpolators guide](../user_guide/examples/tutorial_interpolation.ipynb) |
| 75 | +``` |
| 76 | + |
| 77 | +## 2. ParticleSet |
| 78 | + |
| 79 | +Once the environment has a `parcels.FieldSet` object, you can start defining your particles in a **`parcels.ParticleSet`** object. This object requires: |
| 80 | + |
| 81 | +1. The `parcels.FieldSet` object in which the particles will be released. |
| 82 | +2. The type of `parcels.Particle`: A default `Particle` or a custom `Particle`-type with additional `Variable`s (see the [custom kernel example](custom-kernel)). |
| 83 | +3. Initial conditions for each `Variable` defined in the `Particle`, most notably the release coordinates of `time`, `z`, `lat` and `lon`. |
| 84 | + |
| 85 | +```python |
| 86 | +time = np.array([0]) |
| 87 | +z = np.array([0]) |
| 88 | +lat = np.array([0]) |
| 89 | +lon = np.array([0]) |
| 90 | + |
| 91 | +# Create a ParticleSet |
| 92 | +pset = parcels.ParticleSet(fieldset=fieldset, pclass=parcels.Particle, time=time, z=z, lat=lat, lon=lon) |
| 93 | +``` |
| 94 | + |
| 95 | +```{admonition} 🖥️ Learn more about how to create ParticleSets |
| 96 | +:class: seealso |
| 97 | +- [Release particles at different times](../user_guide/examples/tutorial_delaystart.ipynb) |
| 98 | +``` |
| 99 | + |
| 100 | +## 3. Kernels |
| 101 | + |
| 102 | +A **`parcels.Kernel`** object is a little snippet of code, which is applied to the particles in the `ParticleSet`, for every time step during a simulation. Kernels define the computation or numerical integration done by Parcels, and can represent many processes such as advection, ageing, growth, or simply the sampling of a field. |
| 103 | + |
| 104 | +Advection of a particle by the flow, the change in position $\mathbf{x}(t) = (lon(t), lat(t))$ at time $t$, can be described by the equation: |
| 105 | + |
| 106 | +$$ |
| 107 | +\begin{aligned} |
| 108 | +\frac{\text{d}\mathbf{x}(t)}{\text{d}t} = \mathbf{v}(\mathbf{x}(t),t), |
| 109 | +\end{aligned} |
| 110 | +$$ |
| 111 | + |
| 112 | +where $\mathbf{v}(\mathbf{x},t) = (u(\mathbf{x},t), v(\mathbf{x},t))$ describes the ocean velocity field at position $\mathbf{x}$ at time $t$. |
| 113 | + |
| 114 | +In Parcels, we can write a kernel function which integrates this equation at each timestep `particles.dt`. To do so, we need the ocean velocity field `fieldset.UV` at the `particles` location, and compute the change in position, `particles.dlon` and `particles.dlat`. |
| 115 | + |
| 116 | +```python |
| 117 | +def AdvectionEE(particles, fieldset): |
| 118 | + """Advection of particles using Explicit Euler (aka Euler Forward) integration.""" |
| 119 | + (u1, v1) = fieldset.UV[particles] |
| 120 | + particles.dlon += u1 * particles.dt |
| 121 | + particles.dlat += v1 * particles.dt |
| 122 | +``` |
| 123 | + |
| 124 | +Basic kernels are included in Parcels to compute advection and diffusion. The standard advection kernel is `parcels.kernels.AdvectionRK2`, a [second-order Runge-Kutta integrator](https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#The_Runge%E2%80%93Kutta_method) of the advection function. |
| 125 | + |
| 126 | +```{warning} |
| 127 | +It is advised _not_ to update the particle coordinates (`particles.time`, `particles.z`, `particles.lat`, or `particles.lon`) directly within a Kernel, as that can negatively interfere with the way that particle movements by different kernels are vectorially added. Use a change in the coordinates: `particles.dlon`, `particles.dlat` and/or `particles.dz`. Read the [kernel loop tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_kernelloop.html) to understand why. |
| 128 | +``` |
| 129 | + |
| 130 | +(custom-kernel)= |
| 131 | +We can write custom kernels to add many different types of 'behaviour' to the particles. To do so, we write a function with two arguments: `particles` and `fieldset`. We can then write any computation as a function of any variables defined in the `Particle` and any `Field` defined in the `FieldSet`. Kernels can then be combined by creating a `list` of the kernels. The kernels are executed in order: |
| 132 | + |
| 133 | +```python |
| 134 | +# Create a custom Particle object with an "age" variable |
| 135 | +CustomParticle = parcels.Particle.add_variable( |
| 136 | + parcels.Variable("age", initial=0) |
| 137 | +) |
| 138 | + |
| 139 | +# Create a custom kernel which keeps track of the particle age |
| 140 | +def Age(particles, fieldset): |
| 141 | + particles.age += particles.dt |
| 142 | + |
| 143 | +# define all kernels to be executed on particles using an (ordered) list |
| 144 | +kernels = [Age, parcels.kernels.AdvectionRK4] |
| 145 | +``` |
| 146 | + |
| 147 | +```{note} |
| 148 | +Every Kernel must be a function with the following (and only those) arguments: `(particles, fieldset)` |
| 149 | +``` |
| 150 | + |
| 151 | +```{warning} |
| 152 | +We have to be careful with writing kernels for vector fields on Curvilinear grids. While Parcels automatically rotates the "U" and "V" field when necessary, this is not the case for other fields such as Stokes drift. [This guide](../user_guide/examples/tutorial_nemo_curvilinear.ipynb) describes how to use a curvilinear grid in Parcels. |
| 153 | +``` |
| 154 | + |
| 155 | +```{admonition} 📖 Read more about the Kernel loop |
| 156 | +:class: seealso |
| 157 | +- [The Kernel loop](../user_guide/examples/explanation_kernelloop.md) |
| 158 | +``` |
| 159 | + |
| 160 | +```{admonition} 🖥️ Learn how to write kernels |
| 161 | +:class: seealso |
| 162 | +- [Sample fields like temperature](../user_guide/examples/tutorial_sampling.ipynb). |
| 163 | +- [Mimic the behaviour of ARGO floats](../user_guide/examples/tutorial_Argofloats.ipynb). |
| 164 | +- [Add diffusion to approximate subgrid-scale processes and unresolved physics](../user_guide/examples/tutorial_diffusion.ipynb). |
| 165 | +- [Convert between units in m/s and degrees](../user_guide/examples/tutorial_unitconverters.ipynb). |
| 166 | +``` |
| 167 | + |
| 168 | +## 4. Execute |
| 169 | + |
| 170 | +The execution of the simulation is done using the method **`parcels.ParticleSet.execute()`**, given the `FieldSet`, `ParticleSet`, and `Kernels` defined in the previous steps. This method requires the following arguments: |
| 171 | + |
| 172 | +1. The kernels to be executed. |
| 173 | +2. The `runtime` defining how long the execution loop runs. Alternatively, you may define the `endtime` at which the execution loop stops. |
| 174 | +3. The timestep `dt` at which to execute the kernels. |
| 175 | +4. (Optional) The `ParticleFile` object to write the output to. |
| 176 | + |
| 177 | +```python |
| 178 | +dt = np.timedelta64(5, "m") |
| 179 | +runtime = np.timedelta64(1, "D") |
| 180 | + |
| 181 | +# Run the simulation |
| 182 | +pset.execute(pyfunc=kernels, dt=dt, runtime=runtime) |
| 183 | +``` |
| 184 | + |
| 185 | +### Output |
| 186 | + |
| 187 | +To analyse the particle data generated in the simulation, we need to define a `parcels.ParticleFile` and add it as an argument to `parcels.ParticleSet.execute()`. The output will be written in a [zarr format](https://zarr.readthedocs.io/en/stable/), which can be opened as an `xarray.Dataset`. The dataset will contain the particle data with at least `time`, `z`, `lat` and `lon`, for each particle at timesteps defined by the `outputdt` argument. |
| 188 | + |
| 189 | +There are many ways to analyze particle output, and although we provide [a short tutorial to get started](./tutorial_output.ipynb), we recommend writing your own analysis code and checking out other projects such as [trajan](https://opendrift.github.io/trajan/index.html) and [Lagrangian Diagnostics](https://lagrangian-diags.readthedocs.io/en/latest/). |
| 190 | + |
| 191 | +```{admonition} 🖥️ Learn how to run a simulation |
| 192 | +:class: seealso |
| 193 | +- [Choose an appropriate timestep and integrator](../user_guide/examples/tutorial_numerical_accuracy.ipynb) |
| 194 | +- [Work with Parcels output](./tutorial_output.ipynb) |
| 195 | +``` |
0 commit comments