Skip to content

Add example of limiting particle timestep by local cfl in schism tuto… #3016

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
fluidnumericsJoe wants to merge 4 commits into
base: main
Choose a base branch
from
+164 −29
Open

Add example of limiting particle timestep by local cfl in schism tuto… #3016

Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
276442d
Add example of limiting particle timestep by local cfl in schism tuto…
fluidnumericsJoe Oct 27, 2026
3d7d844
[pre-commit.ci] auto fixes from pre-commit.com hooks
pre-commit-ci[bot] Oct 27, 2026
681ca7a
Address review commments
fluidnumericsJoe Oct 29, 2026
d0db413
Merge branch 'main' into feature/cfl-limited-dt-tutorial-update
fluidnumericsJoe Oct 29, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
193 changes: 164 additions & 29 deletions docs/user_guide/examples/tutorial_schism.ipynb
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -174,8 +174,24 @@
"U = u.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(\n",
" rename\n",
") # reverse to surface-first\n",
"V = v.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename)\n",
"\n",
"V = v.isel(nSCHISM_vgrid_layers=slice(None, None, -1)).rename(rename)"
]
},
{
"cell_type": "markdown",
"id": "9",
"metadata": {},
"source": [
"Now, we can construct the `ux.UxDataset` object that we will use to create the `Fieldset` in Parcels with `Field` objects for `U` and `V`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "10",
"metadata": {},
"outputs": [],
"source": [
"uxds = ux.UxDataset(\n",
" xr.Dataset(\n",
" {\n",
Expand All @@ -191,7 +207,7 @@
},
{
"cell_type": "markdown",
"id": "9",
"id": "11",
"metadata": {},
"source": [
"## Build the `FieldSet`\n",
Expand All @@ -207,7 +223,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "10",
"id": "12",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -220,36 +236,38 @@
},
{
"cell_type": "markdown",
"id": "11",
"id": "13",
"metadata": {},
"source": [
"## Stop particles that end up below the bathymetry\n",
"\n",
"Because SCHISM's LSC2 grid keeps below-seabed levels as `NaN`, a particle whose (fixed) depth ends up\n",
"deeper than the local water column will sample `NaN`. Rather than feed it a fabricated velocity, we let\n",
"Parcels flag it: sampling `NaN` sets the particle state to {py:obj}`parcels.StatusCode.ErrorInterpolation`, and ending up outside\n",
"the mesh sets {py:obj}`parcels.StatusCode.ErrorOutOfBounds`. We add a small kernel that runs after advection, records that the\n",
"particle went out of bounds, and sets its state to {py:obj}`parcels.StatusCode.Delete` so it stops being advected (its trajectory\n",
"Parcels flag it: sampling `NaN` sets the particle state to `parcels.StatusCode.ErrorInterpolation`, and ending up outside\n",
"the mesh sets `parcels.StatusCode.ErrorOutOfBounds` or `parcels.StatusCode.ErrorGridSearching`. We add a small kernel that runs after advection, records that the\n",
"particle went out of bounds, and sets its state to `parcels.StatusCode.Delete` so it stops being advected (its trajectory\n",
"up to that point is kept).\n",
"\n",
"```{note}\n",
"Most of the Lake Ontario mesh is shallow nearshore water with only a couple of valid vertical levels;\n",
"the deep central basin has the full set. We therefore release particles in the deep basin (nodes with\n",
"many valid levels) so they start within the resolved water column.\n",
"the deep central basin has the full set. We therefore release most particles at 2 m depth in the deep\n",
"basin (nodes with many valid levels), and a smaller set at 1 m depth in the finely resolved nearshore\n",
"regions.\n",
"```"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "12",
"id": "14",
"metadata": {},
"outputs": [],
"source": [
"OUT_OF_BOUNDS_STATES = [\n",
" parcels.StatusCode.ErrorInterpolation, # sampled NaN (below the local seabed)\n",
" parcels.StatusCode.ErrorOutOfBounds, # left the horizontal mesh / below the grid\n",
" parcels.StatusCode.ErrorThroughSurface, # above the surface\n",
" parcels.StatusCode.ErrorGridSearching, # left the mesh during a step\n",
"]\n",
"\n",
"SchismParticle = parcels.Particle.add_variable(\n",
Expand All @@ -266,16 +284,104 @@
},
{
"cell_type": "markdown",
"id": "13",
"id": "15",
"metadata": {},
"source": [
"## Release particles and advect"
"## Limit the timestep by the local CFL number\n",
"\n",
"Unstructured meshes are often strongly refined near the coast, so a single `dt` that suits the large\n",
"cells in the basin can carry a particle across several small cells in one step. This is measured by\n",
"the CFL (Courant–Friedrichs–Lewy) number\n",
"\n",
"$$\n",
"C = \\frac{|u|\\,\\Delta t}{\\Delta x},\n",
"$$\n",
"\n",
"the number of cells of width $\\Delta x$ that a particle moving at speed $|u|$ crosses in one timestep\n",
"$\\Delta t$. Keeping $C \\lesssim 1$ ensures that particles do not skip over cells. See\n",
"[a priori estimation of `dt`](./tutorial_dt_integrators.ipynb#a-priori-estimation-of-dt) for the same\n",
"reasoning applied to a structured grid.\n",
"\n",
"To use the local cell width in a kernel, we add a `cell_width` field to the `FieldSet`. While there are a number of valid measures of the cell width of a triangle, here we use the shortest distance from a corner to the opposite edge. This is equivalent to twice the triangle's area divided by its longest edge. The cell width is defined on the faces and does not change in time or depth, so we broadcast it over the `time` and `zc` dimensions, build a `FieldSet` from it, and add that to the velocity `FieldSet`."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "14",
"id": "16",
"metadata": {},
"outputs": [],
"source": [
"# corner coordinates of each triangle, shape (n_face, 3, 2)\n",
"corners = np.stack([node_x[face_nodes], node_y[face_nodes]], axis=-1)\n",
"edges = np.roll(corners, -1, axis=1) - corners # edge vectors\n",
"edge_lengths = np.hypot(edges[..., 0], edges[..., 1])\n",
"cell_area = 0.5 * np.abs(\n",
" edges[:, 0, 0] * edges[:, 1, 1] - edges[:, 0, 1] * edges[:, 1, 0]\n",
")\n",
"cell_width = (\n",
" 2.0 * cell_area / edge_lengths.max(axis=1)\n",
") # smallest height of each triangle\n",
"\n",
"cell_width_da = xr.DataArray(\n",
" np.broadcast_to(cell_width, (U.sizes[\"time\"], zc.size, cell_width.size)),\n",
" dims=(\"time\", \"zc\", \"n_face\"),\n",
" coords={\"time\": U[\"time\"]},\n",
")\n",
"uxds_cell_width = ux.UxDataset(\n",
" xr.Dataset(\n",
" {\"cell_width\": cell_width_da}, coords={\"zf\": (\"zf\", zf), \"zc\": (\"zc\", zc)}\n",
" ),\n",
" uxgrid=uxgrid,\n",
")\n",
"fieldset = fieldset + parcels.FieldSet.from_ugrid_conventions(uxds_cell_width)\n",
"fieldset.describe()"
]
},
{
"cell_type": "markdown",
"id": "17",
"metadata": {},
"source": [
"The `CFLLimitedDt` kernel below shrinks each particle's `dt` so that it moves at most `CFL_target`\n",
"local cell widths per step: $\\Delta t \\le C\\,\\Delta x/|u|$, with $C$ = `CFL_target` and\n",
"$\\Delta x$ = `cell_width`. To expose the `CFL_target` as a value in the kernel, we use the `fieldset.add_context(name,value)` method to store a scalar value under the name `CFL_target`. The `dt` passed to `pset.execute()` then acts as the *maximum* timestep, and\n",
"Parcels sub-steps each particle until it reaches the next output time."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "18",
"metadata": {},
"outputs": [],
"source": [
"def CFLLimitedDt(particles, fieldset):\n",
" \"\"\"Shrink each particle's dt so it moves at most CFL_target cell widths per step.\"\"\"\n",
" cell_width = fieldset.cell_width[particles]\n",
" u, v = fieldset.UV[particles]\n",
" dt_cfl = (\n",
" fieldset.CFL_target * cell_width / np.maximum(np.sqrt(u * u + v * v), 1e-12)\n",
" )\n",
" particles.dt = np.sign(particles.dt) * np.minimum(np.abs(particles.dt), dt_cfl)\n",
"\n",
"\n",
"fieldset.add_context(\"CFL_target\", 1.0)"
]
},
{
"cell_type": "markdown",
"id": "19",
"metadata": {},
"source": [
"## Release particles and advect\n",
"Here, the initial particle positions are selected at the face centroids in the deep basin and in areas where the cells are expected to be small."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "20",
"metadata": {},
"outputs": [],
"source": [
Expand All @@ -290,19 +396,47 @@
"idx = np.where(deep)[0]\n",
"idx = idx[:: max(1, idx.size // n_max)][:n_max]\n",
"\n",
"x, y = cx[idx], cy[idx]\n",
"z = np.full(x.size, 2.0) # release at 2 m depth\n",
"print(f\"releasing {x.size} particles at z = 2 m in the deep basin\")\n",
"\n",
"# release 500 more on faces where the CFL limit at dt = 5 min is below 1 (small cells, fast flow), at 1 m depth\n",
"surf_speed = np.hypot(U.isel(time=0, zf=0).values, V.isel(time=0, zf=0).values)\n",
"dt_cfl = cell_width / np.maximum(\n",
" np.nan_to_num(surf_speed[face_nodes].mean(axis=1)), 1e-12\n",
")\n",
"fine = np.where(dt_cfl < 300)[0]\n",
"idx_fine = np.random.default_rng(0).choice(fine, 500, replace=False)\n",
"\n",
"x = np.concatenate([cx[idx], cx[idx_fine]])\n",
"y = np.concatenate([cy[idx], cy[idx_fine]])\n",
"z = np.concatenate([np.full(idx.size, 2.0), np.full(idx_fine.size, 1.0)])\n",
"print(\n",
" f\"releasing {idx.size} particles in the deep basin and {idx_fine.size} in fine nearshore cells\"\n",
")"
]
},
{
"cell_type": "markdown",
"id": "21",
"metadata": {},
"source": [
"When calling `pset.execute` we call the `CFLLimitedDt` kernel before the `parcels.kernels.AdvectionRK4`. With each execution loop, each particle's time step size is adjusted so that its local CFL estimate is at most `CFL_target` before being advected."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "22",
"metadata": {},
"outputs": [],
"source": [
"pset = parcels.ParticleSet(fieldset=fieldset, pclass=SchismParticle, x=x, y=y, z=z)\n",
"output_file = parcels.ParticleFile(\n",
" \"output-schism.parquet\", outputdt=np.timedelta64(30, \"m\")\n",
")\n",
"\n",
"pset.execute(\n",
" [parcels.kernels.AdvectionRK4, StopBelowBed],\n",
" [CFLLimitedDt, parcels.kernels.AdvectionRK4, StopBelowBed],\n",
" runtime=np.timedelta64(5, \"h\"), # the dataset spans 5 hours\n",
" dt=np.timedelta64(5, \"m\"),\n",
" dt=np.timedelta64(5, \"m\"), # maximum timestep\n",
" output_file=output_file,\n",
" verbose_progress=False,\n",
")\n",
Expand All @@ -311,7 +445,7 @@
},
{
"cell_type": "markdown",
"id": "15",
"id": "23",
"metadata": {},
"source": [
"## Plot the velocity field and trajectories\n",
Expand All @@ -325,7 +459,7 @@
{
"cell_type": "code",
"execution_count": null,
"id": "16",
"id": "24",
"metadata": {},
"outputs": [],
"source": [
Expand Down Expand Up @@ -388,7 +522,7 @@
},
{
"cell_type": "markdown",
"id": "17",
"id": "25",
"metadata": {},
"source": [
"The particles move with the SCHISM surface currents, and any that wander over shallow water where\n",
Expand All @@ -406,7 +540,7 @@
},
{
"cell_type": "markdown",
"id": "18",
"id": "26",
"metadata": {},
"source": [
"## Adapting this tutorial to your own SCHISM setup\n",
Expand All @@ -420,7 +554,8 @@
" $\\Delta x$ approximated by each triangle's shortest edge) is well above 1 in the narrow channel, so\n",
" particles can jump across several cells, or out of the mesh entirely, in a single step. Choose `dt`\n",
" from the *smallest* cells and *fastest* flows that your particles will visit, not from the domain\n",
" average. See [choosing a timestep and integration method](./tutorial_dt_integrators.ipynb).\n",
" average, or use the `CFLLimitedDt` kernel from this tutorial to adapt `dt` per particle. See\n",
" [choosing a timestep and integration method](./tutorial_dt_integrators.ipynb).\n",
"\n",
" ![CFL number for u and v near a SCHISM river inflow with dt = 5 minutes](../../_static/schism_inflow_cfl.png)\n",
"\n",
Expand All @@ -437,15 +572,15 @@
"\n",
"* **Particles leaving the mesh.** Because of the two points above, some particles near the coast may\n",
" still exit the mesh. Decide how to handle them: release particles away from the boundaries (as done\n",
" here with the `valid_levels` filter, whose threshold depends on your mesh), or catch the relevant\n",
" status codes (e.g. add `parcels.StatusCode.ErrorGridSearching` to `OUT_OF_BOUNDS_STATES`) and stop\n",
" those particles. See [working with status codes](./tutorial_statuscodes.md)."
" here for the deep-basin set with the `valid_levels` filter, whose threshold depends on your mesh),\n",
" or catch the relevant status codes (as done here with `OUT_OF_BOUNDS_STATES`) and stop those\n",
" particles. See [working with status codes](./tutorial_statuscodes.md)."
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Parcels:default (3.14.8)",
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
Expand All @@ -459,7 +594,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.14.8"
"version": "3.14.3"
},
"mystnb": {
"execution_in_temp": true
Expand Down
Loading