diff --git a/docs/user_guide/convert_schism.ipynb b/docs/user_guide/convert_schism.ipynb index 95fe1fc2c..5692fb20f 100644 --- a/docs/user_guide/convert_schism.ipynb +++ b/docs/user_guide/convert_schism.ipynb @@ -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", @@ -191,7 +207,7 @@ }, { "cell_type": "markdown", - "id": "9", + "id": "11", "metadata": {}, "source": [ "## Build the `FieldSet`\n", @@ -207,7 +223,7 @@ { "cell_type": "code", "execution_count": null, - "id": "10", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -220,29 +236,30 @@ }, { "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": [ @@ -250,6 +267,7 @@ " 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", @@ -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`](./simulation_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": [ @@ -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", @@ -311,7 +445,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "23", "metadata": {}, "source": [ "## Plot the velocity field and trajectories\n", @@ -325,7 +459,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -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", @@ -406,7 +540,7 @@ }, { "cell_type": "markdown", - "id": "18", + "id": "26", "metadata": {}, "source": [ "## Adapting this tutorial to your own SCHISM setup\n", @@ -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](./simulation_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](./simulation_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", @@ -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](./kernels_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](./kernels_statuscodes.md)." ] } ], "metadata": { "kernelspec": { - "display_name": "Parcels:default (3.14.8)", + "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, @@ -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 diff --git a/docs/user_guide/simulation_dt_integrators.ipynb b/docs/user_guide/simulation_dt_integrators.ipynb index d45e3464c..b1df99c24 100644 --- a/docs/user_guide/simulation_dt_integrators.ipynb +++ b/docs/user_guide/simulation_dt_integrators.ipynb @@ -174,6 +174,16 @@ "cell_type": "markdown", "id": "8", "metadata": {}, + "source": [ + "```{note}\n", + "Note that if your grid cells vary much in size (such as in an unstructured grid), the local CFL number can vary significantly across the domain. In such cases, using a fixed `dt` for all particles may lead to instability in regions with small cells and fast flows. Consider using the `CFLLimitedDt` kernel to adapt `dt` per particle, as documented in the [SCHISM tutorial](./convert_schism.ipynb#).\n", + "```" + ] + }, + { + "cell_type": "markdown", + "id": "9", + "metadata": {}, "source": [ "### Test timesteps\n", "\n", @@ -183,7 +193,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9", + "id": "10", "metadata": {}, "outputs": [], "source": [ @@ -198,7 +208,7 @@ }, { "cell_type": "markdown", - "id": "10", + "id": "11", "metadata": {}, "source": [ "We simulate particles using a range of timesteps that differ by a factor 2-10, starting at (dt < 24 hours). We also keep track of the time it takes to run each simulation:\n", @@ -211,7 +221,7 @@ { "cell_type": "code", "execution_count": null, - "id": "11", + "id": "12", "metadata": {}, "outputs": [], "source": [ @@ -262,7 +272,7 @@ }, { "cell_type": "markdown", - "id": "12", + "id": "13", "metadata": {}, "source": [ "Let's start by comparing the trajectories of the different simulations:" @@ -271,7 +281,7 @@ { "cell_type": "code", "execution_count": null, - "id": "13", + "id": "14", "metadata": {}, "outputs": [], "source": [ @@ -282,7 +292,7 @@ { "cell_type": "code", "execution_count": null, - "id": "14", + "id": "15", "metadata": {}, "outputs": [], "source": [ @@ -311,7 +321,7 @@ }, { "cell_type": "markdown", - "id": "15", + "id": "16", "metadata": {}, "source": [ "We see that as the timestep becomes smaller, the trajectories start to converge. As expected from our horizontal grid scaling, the 2 simulations with `dt` > 2 hours produce trajectories that end up in different grid cells. The other three simulations (`dt` = (1 h, 20 min, 5 min)) are more difficult to separate visually. \n", @@ -322,7 +332,7 @@ { "cell_type": "code", "execution_count": null, - "id": "16", + "id": "17", "metadata": {}, "outputs": [], "source": [ @@ -355,7 +365,7 @@ { "cell_type": "code", "execution_count": null, - "id": "17", + "id": "18", "metadata": {}, "outputs": [], "source": [ @@ -412,7 +422,7 @@ }, { "cell_type": "markdown", - "id": "18", + "id": "19", "metadata": {}, "source": [ "```{warning}\n", @@ -423,7 +433,7 @@ { "cell_type": "code", "execution_count": null, - "id": "19", + "id": "20", "metadata": {}, "outputs": [], "source": [ @@ -459,7 +469,7 @@ }, { "cell_type": "markdown", - "id": "20", + "id": "21", "metadata": {}, "source": [ "We can see that in our simulation advecting particles for 7 days, the effect of `dt` on the precision of our simulation is approximately linear. The precision for a simulation with a timestep of 20 minutes is order of magnitude ~100 m. The effect on the time it takes to run a simulation is not linear in our case however; it increases sharply as we decrease our timestep. This may be optimised using more efficient chunking." @@ -468,7 +478,7 @@ { "cell_type": "code", "execution_count": null, - "id": "21", + "id": "22", "metadata": {}, "outputs": [], "source": [ @@ -486,7 +496,7 @@ }, { "cell_type": "markdown", - "id": "22", + "id": "23", "metadata": {}, "source": [ "## Integration schemes\n", @@ -497,7 +507,7 @@ { "cell_type": "code", "execution_count": null, - "id": "23", + "id": "24", "metadata": {}, "outputs": [], "source": [ @@ -510,7 +520,7 @@ { "cell_type": "code", "execution_count": null, - "id": "24", + "id": "25", "metadata": {}, "outputs": [], "source": [ @@ -523,7 +533,7 @@ { "cell_type": "code", "execution_count": null, - "id": "25", + "id": "26", "metadata": {}, "outputs": [], "source": [ @@ -535,7 +545,7 @@ }, { "cell_type": "markdown", - "id": "26", + "id": "27", "metadata": {}, "source": [ "The higher-order methods use weighted intermediate steps in time and space to obtain a more accurate estimate of `dy` and `dx` for a given timestep.\n", @@ -546,7 +556,7 @@ { "cell_type": "code", "execution_count": null, - "id": "27", + "id": "28", "metadata": {}, "outputs": [], "source": [ @@ -560,7 +570,7 @@ { "cell_type": "code", "execution_count": null, - "id": "28", + "id": "29", "metadata": { "tags": [ "hide-output" @@ -607,7 +617,7 @@ { "cell_type": "code", "execution_count": null, - "id": "29", + "id": "30", "metadata": {}, "outputs": [], "source": [ @@ -645,7 +655,7 @@ }, { "cell_type": "markdown", - "id": "30", + "id": "31", "metadata": {}, "source": [ "Clearly, for longer timesteps, the RK2 and RK4 schemes perform better. However, if the timestep is appropriate, as we have determined in the previous section, then the Explicit Euler scheme does not perform notably different." @@ -654,7 +664,7 @@ { "cell_type": "code", "execution_count": null, - "id": "31", + "id": "32", "metadata": {}, "outputs": [], "source": [ @@ -702,7 +712,7 @@ }, { "cell_type": "markdown", - "id": "32", + "id": "33", "metadata": {}, "source": [ "By quantifying the precision of the integration methods, we can see that for a given timestep the Runge-Kutta methods perform orders of magnitude better than the Explicit Euler method. In this example, the error associated with the selected integration methods is smaller than that of the range of timesteps." @@ -711,7 +721,7 @@ { "cell_type": "code", "execution_count": null, - "id": "33", + "id": "34", "metadata": {}, "outputs": [], "source": [ @@ -753,7 +763,7 @@ }, { "cell_type": "markdown", - "id": "34", + "id": "35", "metadata": {}, "source": [ "In this last figure, we see that the improvement of RK2 by orders of magnitude with respect to EE comes at a small computational cost. Since the RK2 and RK4 methods are practically indistinguishable, using the RK2 method with a timestep of 1 hour or 20 minutes, depending on the application, would be an appropriate choice for this simulation." @@ -761,7 +771,7 @@ }, { "cell_type": "markdown", - "id": "35", + "id": "36", "metadata": {}, "source": [ "### Flow conditions\n", @@ -784,7 +794,7 @@ }, { "cell_type": "markdown", - "id": "36", + "id": "37", "metadata": {}, "source": [ "## Writing custom Kernels with higher-order integrators\n", @@ -799,7 +809,7 @@ { "cell_type": "code", "execution_count": null, - "id": "37", + "id": "38", "metadata": {}, "outputs": [], "source": [ @@ -819,7 +829,7 @@ }, { "cell_type": "markdown", - "id": "38", + "id": "39", "metadata": {}, "source": [ "Now, if we also load in the wind data for the same region with the code below, we can execute this Kernel to compute the trajectories of particles that are advected by the ocean currents and affected by the wind. The solution will be second-order accurate in time, since we are using the RK2 integrator; but it is easy to swap in the RK4 integrator if you want to increase the accuracy even more." @@ -828,7 +838,7 @@ { "cell_type": "code", "execution_count": null, - "id": "39", + "id": "40", "metadata": {}, "outputs": [], "source": [ @@ -851,7 +861,7 @@ { "cell_type": "code", "execution_count": null, - "id": "40", + "id": "41", "metadata": {}, "outputs": [], "source": [ @@ -883,7 +893,7 @@ }, { "cell_type": "markdown", - "id": "41", + "id": "42", "metadata": {}, "source": [ "As you can see, the particles move very differently when we include the windage effect!" @@ -892,7 +902,7 @@ { "cell_type": "code", "execution_count": null, - "id": "42", + "id": "43", "metadata": {}, "outputs": [], "source": [ @@ -904,7 +914,7 @@ }, { "cell_type": "markdown", - "id": "43", + "id": "44", "metadata": {}, "source": [ "## Summary\n",