Add fesom tutorial to documentation

docs/user_guide/examples/tutorial_fesom.ipynb

@@ -0,0 +1,308 @@
+{
+ "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 matplotlib.tri as mtri\n",
+    "import numpy as np\n",
+    "import uxarray as ux\n",
+    "\n",
+    "import parcels\n",
+    "import parcels.tutorial"
+   ]
+  },
+  {
+   "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 = parcels.convert.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 velocity field and trajectories\n\nWe plot the particle paths on top of the velocity field they advect through: triangle colour shows the speed at the release depth (≈50 m), black arrows show the velocity at face centres (drawn in lon/lat space, length proportional to speed), grey lines trace each particle's path, and the coloured dots mark the positions over time. Drawing the arrows with `angles=\"xy\"` keeps them aligned with the trajectories, so you can see the particles streak along the fast jets and barely move in the quiet bands between them:"
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "df = parcels.read_particlefile(\"output-fesom.parquet\")\n",
+    "\n",
+    "triang = mtri.Triangulation(\n",
+    "    ds.uxgrid.node_lon.values,\n",
+    "    ds.uxgrid.node_lat.values,\n",
+    "    triangles=ds.uxgrid.face_node_connectivity.values,\n",
+    ")\n",
+    "\n",
+    "depth_idx = int(np.argmin(np.abs(ds.zc.values - 50.0)))\n",
+    "U_face = np.asarray(ds[\"U\"].isel(zc=depth_idx)).squeeze()\n",
+    "V_face = np.asarray(ds[\"V\"].isel(zc=depth_idx)).squeeze()\n",
+    "speed = np.hypot(U_face, V_face)\n",
+    "\n",
+    "fig, ax = plt.subplots(figsize=(11, 5))\n",
+    "\n",
+    "# Background: speed at the release depth.\n",
+    "tpc = ax.tripcolor(triang, facecolors=speed, shading=\"flat\", cmap=\"Blues\")\n",
+    "fig.colorbar(tpc, ax=ax, label=\"speed [m/s]\", location=\"left\", shrink=0.85)\n",
+    "\n",
+    "# Velocity arrows at face centres. Drawing in lon/lat space (angles/scale_units\n",
+    "# \"xy\") keeps the arrows aligned with the trajectories; length tracks speed.\n",
+    "step = max(1, U_face.size // 400)\n",
+    "xf = ds.uxgrid.face_lon.values\n",
+    "yf = ds.uxgrid.face_lat.values\n",
+    "max_speed = float(np.nanmax(speed))\n",
+    "q = ax.quiver(\n",
+    "    xf[::step],\n",
+    "    yf[::step],\n",
+    "    U_face[::step],\n",
+    "    V_face[::step],\n",
+    "    angles=\"xy\",\n",
+    "    scale_units=\"xy\",\n",
+    "    scale=max_speed / 0.3,\n",
+    "    color=\"k\",\n",
+    "    width=0.0018,\n",
+    "    pivot=\"tail\",\n",
+    ")\n",
+    "ax.quiverkey(\n",
+    "    q, 0.86, 1.04, max_speed, f\"{max_speed:.2f} m/s\", labelpos=\"E\", coordinates=\"axes\"\n",
+    ")\n",
+    "\n",
+    "# Particle paths (grey lines) and positions coloured by time.\n",
+    "for traj in df.sort(\"time\").partition_by(\"particle_id\"):\n",
+    "    ax.plot(traj[\"lon\"], traj[\"lat\"], color=\"0.4\", linewidth=0.6, alpha=0.7, zorder=2)\n",
+    "ax.scatter(lon, lat, facecolors=\"none\", edgecolors=\"k\", s=60, label=\"release\", zorder=3)\n",
+    "sc = ax.scatter(\n",
+    "    df[\"lon\"],\n",
+    "    df[\"lat\"],\n",
+    "    c=df[\"time\"].dt.total_seconds(),\n",
+    "    s=6,\n",
+    "    cmap=\"viridis\",\n",
+    "    zorder=3,\n",
+    ")\n",
+    "fig.colorbar(sc, ax=ax, label=\"time since release [s]\", shrink=0.85)\n",
+    "\n",
+    "ax.set_xlabel(\"Longitude [deg E]\")\n",
+    "ax.set_ylabel(\"Latitude [deg N]\")\n",
+    "ax.set_title(\n",
+    "    f\"FESOM2 velocity at z ≈ {float(ds.zc.values[depth_idx]):.1f} m with particle trajectories\"\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
+}

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