diff --git a/docs/user_guide/examples/tutorial_fesom.ipynb b/docs/user_guide/examples/tutorial_fesom.ipynb new file mode 100644 index 000000000..2246ab8ed --- /dev/null +++ b/docs/user_guide/examples/tutorial_fesom.ipynb @@ -0,0 +1,258 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# 🖥️ FESOM tutorial\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Parcels v4 supports unstructured-grid model output via [uxarray](https://uxarray.readthedocs.io/). This tutorial walks through the minimum steps to advect particles in real [FESOM2](https://fesom.de/) output. The recipe is:\n", + "\n", + "1. Open the FESOM grid and data files with `uxarray`.\n", + "2. Rename FESOM-specific dimensions to Parcels' UGRID conventions with `parcels.convert.fesom_to_ugrid`.\n", + "3. Build a `FieldSet` with `parcels.FieldSet.from_ugrid_conventions`.\n", + "4. Run the simulation as on any structured grid.\n", + "\n", + "If you have not done so already, work through the [quickstart tutorial](../../getting_started/tutorial_quickstart.md) first to get familiar with `ParticleSet`, `Kernel`, and `ParticleFile`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from pathlib import Path\n", + "\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import uxarray as ux\n", + "\n", + "import parcels\n", + "import parcels.tutorial\n", + "from parcels.convert import fesom_to_ugrid" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Get the FESOM tutorial dataset\n", + "\n", + "We use a small periodic-channel snapshot from a FESOM2 simulation that ships with Parcels' tutorial data registry. As in the [quickstart](../../getting_started/tutorial_quickstart.md), `parcels.tutorial.open_dataset` downloads the files into a local cache on first use; subsequent calls just return the cached copy.\n", + "\n", + "`uxarray` expects file paths rather than an in-memory dataset, so we trigger the downloads and then point `ux.open_mfdataset` at the cached files:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for name in [\n", + " \"FESOM_periodic_channel/fesom_channel\", # grid description\n", + " \"FESOM_periodic_channel/u.fesom_channel\", # zonal velocity (face-registered)\n", + " \"FESOM_periodic_channel/v.fesom_channel\", # meridional velocity (face-registered)\n", + " \"FESOM_periodic_channel/w.fesom_channel\", # vertical velocity (node-registered)\n", + "]:\n", + " parcels.tutorial.open_dataset(name)\n", + "\n", + "from parcels._datasets.remote import _DATA_HOME\n", + "\n", + "data_dir = Path(_DATA_HOME) / \"data\" / \"FESOM_periodic_channel\"\n", + "\n", + "grid_path = str(data_dir / \"fesom_channel.nc\")\n", + "data_paths = [\n", + " str(data_dir / \"u.fesom_channel.nc\"),\n", + " str(data_dir / \"v.fesom_channel.nc\"),\n", + " str(data_dir / \"w.fesom_channel.nc\"),\n", + "]" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "```{note}\n", + "The `FESOM_periodic_channel` dataset is a single-time-step snapshot of an idealised channel (~4.4° × ~18° wide, 40 vertical layers). Working with multi-time FESOM output is identical, except you pass a glob or list of data files to `ux.open_mfdataset`.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Open the data with `uxarray`\n", + "\n", + "`ux.open_mfdataset(grid_path, data_paths)` reads the FESOM grid description and joins the velocity files on its grid. FESOM names its velocity variables `u`, `v`, `w` — we rename them to `U`, `V`, `W` so that `parcels.FieldSet.from_ugrid_conventions` recognises them as the velocity components:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = ux.open_mfdataset(grid_path, data_paths).rename_vars(\n", + " {\"u\": \"U\", \"v\": \"V\", \"w\": \"W\"}\n", + ")\n", + "ds" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Note the FESOM-specific dimension names: `elem` (number of triangular faces), `nod2` (number of nodes), `nz1` (vertical layer centres), and `nz` (layer interfaces). The horizontal velocities `U` and `V` live on face centres along `nz1`; the vertical velocity `W` lives on nodes along `nz`." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Convert to UGRID conventions\n", + "\n", + "Parcels works with a small UGRID-compliant dialect: nodes are `n_node`, faces are `n_face`, vertical centres are `zc`, and vertical interfaces are `zf`. The helper `parcels.convert.fesom_to_ugrid` does this rename in one call:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "ds = fesom_to_ugrid(ds)\n", + "print(\"dims:\", dict(ds.sizes))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Build the `FieldSet`\n", + "\n", + "With UGRID-compliant dimensions in place, `parcels.FieldSet.from_ugrid_conventions` builds the `FieldSet`. It detects `U`, `V`, `W`, assigns a `UxGrid` to each field, and picks an appropriate interpolator based on each variable's coordinate location (face- vs. node-registered, centre vs. interface). Use `mesh=\"spherical\"` so that velocities in m/s are correctly converted to deg/s along the lon/lat coordinates:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh=\"spherical\")\n", + "\n", + "for name, field in fieldset.fields.items():\n", + " interp = getattr(field, \"interp_method\", None)\n", + " interp_name = interp.__name__ if interp is not None else \"-\"\n", + " print(f\"{name:>4s} -> {type(field).__name__:<11s} interp={interp_name}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "`U` and `V` get face-registered interpolation, `W` gets node-registered linear interpolation. The combined vector fields `UV` and `UVW` are assembled automatically." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "## Release particles and advect\n\nWe seed particles on a grid of four latitudes spanning the channel and ten longitudes, and integrate for two days with RK4. Because this snapshot has only a single time level, `fieldset.time_interval` is `None` and we omit the `time=` argument so that Parcels treats the flow as constant in time:" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "lon_grid, lat_grid = np.meshgrid(\n", + " np.linspace(0.5, 4.0, 10),\n", + " np.linspace(3.0, 15.0, 4),\n", + ")\n", + "lon = lon_grid.ravel()\n", + "lat = lat_grid.ravel()\n", + "z = np.full(lon.size, 50.0) # release at 50 m depth\n", + "\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " pclass=parcels.Particle,\n", + " lon=lon,\n", + " lat=lat,\n", + " z=z,\n", + ")\n", + "\n", + "output_file = parcels.ParticleFile(\n", + " \"output-fesom.parquet\", outputdt=np.timedelta64(1, \"h\")\n", + ")\n", + "\n", + "pset.execute(\n", + " [parcels.kernels.AdvectionRK4],\n", + " runtime=np.timedelta64(2, \"D\"),\n", + " dt=np.timedelta64(5, \"m\"),\n", + " output_file=output_file,\n", + " verbose_progress=False,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": "## Plot the result\n\nEach particle trajectory is coloured by observation time so we can follow how the particles drift through the FESOM2 velocity field:" + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "df = parcels.read_particlefile(\"output-fesom.parquet\")\n", + "\n", + "fig, ax = plt.subplots(figsize=(10, 5))\n", + "ax.scatter(lon, lat, facecolors=\"none\", edgecolors=\"k\", s=60, label=\"release\")\n", + "sc = ax.scatter(\n", + " df[\"lon\"], df[\"lat\"], c=df[\"time\"].dt.total_seconds(), s=8, cmap=\"viridis\"\n", + ")\n", + "fig.colorbar(sc, ax=ax, label=\"time since release [s]\")\n", + "ax.set_xlabel(\"Longitude [deg E]\")\n", + "ax.set_ylabel(\"Latitude [deg N]\")\n", + "ax.legend(loc=\"upper right\")\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The particles drift through the channel following the FESOM2 velocity field. From here, the rest of Parcels — custom kernels, sampling fields onto particles, writing your own interpolators — works identically to structured grids. See the [interpolation tutorial](./tutorial_interpolation.ipynb) for the available `Ux*` interpolators and how to write a custom one." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "docs", + "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.14.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index 3303f1bc6..9e22f1236 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -28,6 +28,7 @@ examples/explanation_grids.md examples/tutorial_nemo.ipynb examples/tutorial_croco_3D.ipynb examples/tutorial_mitgcm.ipynb +examples/tutorial_fesom.ipynb examples/tutorial_velocityconversion.ipynb examples/tutorial_nestedgrids.ipynb examples/tutorial_manipulating_field_data.ipynb