{ "cells": [ { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "# Parcels structure overview\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The flexibility of Parcels allows for a wide range of applications and to build complex simulations. In order to help structuring your code, this tutorial describes the structure that a Parcels script uses.\n", "\n", "Code that uses Parcels is generally build up from four different components:\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "1. [**FieldSet**](#1.-FieldSet). Load and set up the fields. These can be velocity fields that are used to advect the particles, but it can also be e.g. temperature.\n", "2. [**ParticleSet**](#2.-ParticleSet). Define the type of particles. Also additional `Variables` can be added to the particles (e.g. temperature, to keep track of the temperature that particles experience).\n", "3. [**Kernels**](#3.-Kernels). Define and compile kernels. Kernels perform some specific operation on the particles every time step (e.g. interpolate the temperature from the temperature field to the particle location).\n", "4. [**Execution and output**](#4.-Execution-and-Output). Execute the simulation and write and store the output in a NetCDF file.\n", "\n", "We discuss each component in more detail below.\n", "\n", "![png](images/parcels_user_diagram.png)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 1. FieldSet\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Parcels provides a framework to simulate the movement of particles **within an existing flow field environment**. To start a parcels simulation we must define this environment with the `FieldSet` class. The minimal requirements for this Fieldset are that it must contain the `'U'` and `'V'` fields: the 2D hydrodynamic data that will move the particles in a horizontal direction. Additionally, it can contain e.g. a temperature or vertical flow field.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "A fieldset can be loaded with `FieldSet.from_netcdf`, if the model output with the fields is written in NetCDF files. This function requires `filenames`, `variables` and `dimensions`.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "In this example, we only load the `'U'` and `'V'` fields, which represent the zonal and meridional flow velocity. First, `fname` points to the location of the model output.\n" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "from parcels import download_example_dataset\n", "\n", "example_dataset_folder = download_example_dataset(\"GlobCurrent_example_data\")\n", "fname = f\"{example_dataset_folder}/*.nc\"\n", "filenames = {\"U\": fname, \"V\": fname}" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Second, we have to specify the `'U'` and `'V'` variable names, as given by the NetCDF files.\n" ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "variables = {\n", " \"U\": \"eastward_eulerian_current_velocity\",\n", " \"V\": \"northward_eulerian_current_velocity\",\n", "}" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Third, we specify the names of the variable dimensions, as given by the NetCDF files.\n" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "# In the GlobCurrent data the dimensions are also called 'lon', 'lat' and 'time\n", "dimensions = {\n", " \"U\": {\"lat\": \"lat\", \"lon\": \"lon\", \"time\": \"time\"},\n", " \"V\": {\"lat\": \"lat\", \"lon\": \"lon\", \"time\": \"time\"},\n", "}" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Finally, we load the fieldset.\n" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "from parcels import FieldSet\n", "\n", "fieldset = FieldSet.from_netcdf(filenames, variables, dimensions)" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### For more advanced tutorials on creating `FieldSets`:\n", "\n", "- [Implement **periodic boundaries**](https://docs.oceanparcels.org/en/latest/examples/tutorial_periodic_boundaries.html)\n", "- [How to **interpolate** field data for different fields](https://docs.oceanparcels.org/en/latest/examples/tutorial_interpolation.html)\n", "- [**Converting units** in the field data](https://docs.oceanparcels.org/en/latest/examples/tutorial_unitconverters.html)\n", "- [Working around incompatible **time coordinates**](https://docs.oceanparcels.org/en/latest/examples/tutorial_timestamps.html)\n", "\n", "### If you use interpolated output (such as from the [Copernicus Marine Service](https://marine.copernicus.eu)):\n", "\n", "- [Understanding **stuck particles** at land boundaries](https://docs.oceanparcels.org/en/latest/examples/documentation_stuck_particles.html)\n", "- [**Preventing stuck particles** at land boundaries](https://docs.oceanparcels.org/en/latest/examples/documentation_unstuck_Agrid.html)\n", "- Or download and use the [data on the **original C-grid**](https://www.mercator-ocean.eu/en/solutions-expertise/accessing-digital-data/) from [Mercator Ocean International](https://www.mercator-ocean.eu) directly\n", "\n", "### If you are working with field data on different grids:\n", "\n", "- [Grid **indexing** on different grids](https://docs.oceanparcels.org/en/latest/examples/documentation_indexing.html)\n", "- [Load field data from **Curvilinear grids**](https://docs.oceanparcels.org/en/latest/examples/tutorial_nemo_curvilinear.html)\n", "- [Load field data from **3D C-grids**](https://docs.oceanparcels.org/en/latest/examples/tutorial_nemo_3D.html)\n", "\n", "### If you want to combine different velocity fields:\n", "\n", "- [Nest velocity fields of different regions or resolutions in a **NestedField**](https://docs.oceanparcels.org/en/latest/examples/tutorial_NestedFields.html)\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 2. ParticleSet\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Once the environment has a `FieldSet` object, you can start defining your particles in a `ParticleSet` object. This object requires:\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "1. The `FieldSet` object in which the particles will be released.\n", "2. The type of `Particle`: Either a `JITParticle` or `ScipyParticle`.\n", "3. Initial conditions for each `Variable` defined in the `Particle`, most notably the release locations in `lon` and `lat`.\n" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "from parcels import JITParticle, ParticleSet, Variable\n", "\n", "# Define a new particleclass with Variable 'age' with initial value 0.\n", "AgeParticle = JITParticle.add_variable(Variable(\"age\", initial=0))\n", "\n", "pset = ParticleSet(\n", " fieldset=fieldset, # the fields that the particleset uses\n", " pclass=AgeParticle, # define the type of particle\n", " lon=29, # release longitude\n", " lat=-33, # release latitude\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### For more advanced tutorials on how to setup your `ParticleSet`:\n", "\n", "- [**Releasing particles** at different times](https://docs.oceanparcels.org/en/latest/examples/tutorial_delaystart.html)\n", "- [The difference between **JITParticles and ScipyParticles**](https://docs.oceanparcels.org/en/latest/examples/tutorial_jit_vs_scipy.html)\n", "\n", "For more information on how to implement `Particle` types with specific behaviour, see the [section on writing your own kernels](#For-more-advanced-tutorials-on-writing-custom-kernels-that-work-on-custom-particles:).\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 3. Kernels\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "Kernels are little snippets of code, which are applied to every particle in the `ParticleSet`, for every time step during a simulation.\n", "Basic kernels are [included in Parcels](https://parcels.readthedocs.io/en/latest/reference/predefined_kernels.html#predefined-kernels), among which several types of advection kernels. `AdvectionRK4` is the main advection kernel for two-dimensional advection, which is also used in this example.\n", "\n", "One can also write custom kernels, to add certain types of 'behaviour' to the particles. Kernels can then be combined with the `+` operator (where at least one of the kernels needs to be cast to a `pset.Kernel()` object), or by creating a `list` of the kernels. Note that the kernels are executed in order.\n" ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "# load basic advection kernels\n", "from parcels import AdvectionRK4\n", "\n", "# Create a custom kernel which displaces each particle southward\n", "\n", "\n", "def NorthVel(particle, fieldset, time):\n", " if time > 10 * 86400 and time < 10.2 * 86400:\n", " vvel = -1e-4\n", " particle_dlat += vvel * particle.dt\n", "\n", "\n", "# Create a custom kernel which keeps track of the particle age (minutes)\n", "\n", "\n", "def Age(particle, fieldset, time):\n", " particle.age += particle.dt / 3600\n", "\n", "\n", "# define all kernels to be executed on particles using an (ordered) list\n", "kernels = [Age, NorthVel, AdvectionRK4]" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "\n", "Some key limitations exist to the Kernels that everyone who wants to write their own should be aware of:\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "- Every Kernel must be a function with the following (and only those) arguments: `(particle, fieldset, time)`\n", "\n", "- In order to run successfully in JIT mode, Kernel definitions can only contain the following types of commands:\n", " - Basic arithmetical operators (`+`, `-`, `*`, `/`, `**`) and assignments (`=`).\n", "\n", " - Basic logical operators (`<`, `==`, `!=`, `>`, `&`, `|`). Note that you can use a statement like `particle.lon != particle.lon` to check if `particle.lon` is NaN (since `math.nan != math.nan`).\n", "\n", " - `if` and `while` loops, as well as `break` statements. Note that `for`-loops are not supported in JIT mode.\n", " \n", " - Interpolation of a `Field` from the `FieldSet` at a `[time, depth, lat, lon]` point, using square brackets notation.\n", " For example, to interpolate the zonal velocity (U) field at the particle location, use the following statement:\n", " ```python\n", " value = fieldset.U[time, particle.depth, particle.lat, particle.lon]\n", " ```\n", " or simple\n", " ```python\n", " value = fieldset.U[particle]\n", " ```\n", " \n", " - Functions from the maths standard library.\n", " \n", " - Functions from the custom ParcelsRandom library at parcels.rng. Note that these have to be used as ParcelsRandom.random(), ParcelsRandom.uniform() etc for the code to compile.\n", " \n", " - Simple `print` statements, such as:\n", " - `print(\"Some print\")`\n", " - `print(particle.lon)`\n", " - `print(f\"particle id: {particle.id}\")`\n", " - `print(f\"lon: {particle.lon}, lat: {particle.lat}\")`\n", " \n", " - Local variables can be used in Kernels, and these variables will be accessible in all concatenated Kernels. Note that these local variables are not shared between particles, and also not between time steps.\n", " \n", " - It is advised _not_ to update the particle location (`particle.lon`, `particle.lat`, `particle.depth`, and/or `particle.time`) directly, as that can negatively interfere with the way that particle movements by different kernels are vectorially added. Use `particle_dlon`, `particle_dlat`, `particle_ddepth`, and/or `particle_dtime` instead. See also the [kernel loop tutorial](https://docs.oceanparcels.org/en/latest/examples/tutorial_kernelloop.html).\n", " \n", " - Note that one has 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 for example wind data. In that case, a custom rotation function will have to be written.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### For more advanced tutorials on writing custom kernels that work on custom particles:\n", "\n", "- [Sample other fields like temperature](https://docs.oceanparcels.org/en/latest/examples/tutorial_sampling.html).\n", "- [Mimic the behaviour of ARGO floats](https://docs.oceanparcels.org/en/latest/examples/tutorial_Argofloats.html).\n", "- [Adding diffusion to approximate subgrid-scale processes and unresolved physics](https://docs.oceanparcels.org/en/latest/examples/tutorial_diffusion.html).\n", "- [Converting between units in m/s and degree/s](https://docs.oceanparcels.org/en/latest/examples/tutorial_unitconverters.html).\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "## 4. Execution and output\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "The final part executes the simulation, given the `ParticleSet`, `FieldSet` and `Kernels`, that have been defined in the previous steps. If you like to store the particle data generated in the simulation, you define the `ParticleFile` to which the output of the kernel execution will be written. Then, on the `ParticleSet` you have defined, you can use the method `ParticleSet.execute()` which requires the following arguments:\n", "\n", "1. The kernels to be executed.\n", "2. The `runtime` defining how long the execution loop runs. Alternatively, you may define the `endtime` at which the execution loop stops.\n", "3. The timestep `dt` at which to execute the kernels.\n", "4. (Optional) The `ParticleFile` object to write the output to.\n" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "INFO: Output files are stored in GCParticles.zarr.\n", "100%|██████████| 2073600.0/2073600.0 [00:04<00:00, 444884.89it/s]\n" ] } ], "source": [ "output_file = pset.ParticleFile(\n", " name=\"GCParticles.zarr\", # the name of the output file\n", " outputdt=3600, # the time period between consecutive out output steps\n", " chunks=(1, 10), # the chunking of the output file (number of particles, timesteps)\n", ")\n", "\n", "pset.execute(\n", " kernels, # the kernel (which defines how particles move)\n", " runtime=86400 * 24, # the total length of the run in seconds\n", " dt=300, # the timestep of the kernel in seconds\n", " output_file=output_file,\n", ")" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "
\n", "\n", "### A note on output chunking\n", "\n", "Note the use of the `chunks` argument in the `pset.ParticleFile()` above. This controls the 'chunking' of the output file, which is a way to optimize the writing of the output file. The default chunking for the output in Parcels is `(number of particles in initial particleset, 1)`. \n", "Note that this default may not be very efficient if \n", "1. you use `repeatdt` to release a relatively small number of particles _many_ times during the simulation and/or\n", "2. you expect to output _a lot of timesteps_ (e.g. more than 1000).\n", "\n", "In the first case, it is best to increase the first argument of `chunks` to 10 to 100 times the size of your initial particleset. In the second case, it is best to increase the second argument of `chunks` to 10 to 1000, depending a bit on the size of your initial particleset.\n", "\n", "In either case, it will generally be much more efficient if `chunks[0]*chunks[1]` is (much) greater than several thousand.\n", "\n", "See also [the advanced output in zarr format tutorial](https://docs.oceanparcels.org/en/latest/examples/documentation_advanced_zarr.html) for more information on this. The details will depend on the nature of the filesystem the data is being written to, so it is worth to optimise this parameter in your runs, as it can significantly speed up the writing of the output file and thus the runtime of `pset.execution()`.\n", "\n", "
" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "After running your simulation, you probably want to analyze the output. You can use the [trajan](https://opendrift.github.io/trajan/index.html) package for some simple plotting. However, we **recommend you write your own code** to analyze your specific output and you can probably separate the analysis from the simulation.\n" ] }, { "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ "### For more tutorials on the parcels output:\n", "\n", "- [How the output is structured and how to start your own analysis](https://docs.oceanparcels.org/en/latest/examples/tutorial_output.html)\n" ] } ], "metadata": { "celltoolbar": "Metagegevens bewerken", "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": 4 }