From 38ed02619a37a06adc4a7ded1c853fdee6a43485 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 10 Jul 2026 11:49:00 +0200 Subject: [PATCH 1/3] Updating tutorials that combine fields The particle release was not very good; too close to the boundary. Also made a few other fixes, such as using RK2 for wind advection and improving wind field --- .../examples/explanation_kernelloop.md | 20 +++++++++------- .../tutorial_manipulating_field_data.ipynb | 24 ++++++++++++------- 2 files changed, 27 insertions(+), 17 deletions(-) diff --git a/docs/user_guide/examples/explanation_kernelloop.md b/docs/user_guide/examples/explanation_kernelloop.md index e6f1cccb0..e415e3250 100644 --- a/docs/user_guide/examples/explanation_kernelloop.md +++ b/docs/user_guide/examples/explanation_kernelloop.md @@ -62,7 +62,7 @@ ds_fields.load() # load the dataset into memory # Create an idealised wind field and add it to the dataset tdim, ydim, xdim = (len(ds_fields.time),len(ds_fields.latitude), len(ds_fields.longitude)) ds_fields["UWind"] = xr.DataArray( - data=np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values)[None, :, None], + data=0.5 * np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[None, :, None], coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude]) ds_fields["VWind"] = xr.DataArray( @@ -88,10 +88,12 @@ fieldset = parcels.FieldSet.from_sgrid_conventions( Now we define a wind Kernel that uses a forward Euler method to apply the wind forcing. Note that we update the `particles.dx` and `particles.dy` variables, rather than `particles.x` and `particles.y` directly. ```{code-cell} -def wind_kernel(particles, fieldset): - uwind, vwind = fieldset.Wind[particles] - particles.dx += uwind * particles.dt - particles.dy += vwind * particles.dt +def windRK2(particles, fieldset): + (u1, v1) = fieldset.Wind[particles] + x1, y1 = (particles.x + u1 * 0.5 * particles.dt, particles.y + v1 * 0.5 * particles.dt) + (u2, v2) = fieldset.Wind[particles.t + 0.5 * particles.dt, particles.z, y1, x1, particles] + particles.dx += u2 * particles.dt + particles.dy += v2 * particles.dt ``` First run a simulation where we apply Kernels as `[AdvectionRK2, wind_kernel]` @@ -100,7 +102,7 @@ First run a simulation where we apply Kernels as `[AdvectionRK2, wind_kernel]` :tags: [hide-output] npart = 10 z = np.repeat(ds_fields.depth[0].values, npart) -lons = np.repeat(31, npart) +lons = np.repeat(32.2, npart) lats = np.linspace(-32.5, -30.5, npart) pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, y=lats, x=lons) @@ -108,14 +110,14 @@ output_file = parcels.ParticleFile( path="advection_then_wind.parquet", outputdt=np.timedelta64(6,'h') ) pset.execute( - [parcels.kernels.AdvectionRK2, wind_kernel], + [parcels.kernels.AdvectionRK2, windRK2], runtime=np.timedelta64(5,'D'), dt=np.timedelta64(1,'h'), output_file=output_file, ) ``` -Then also run a simulation where we apply the Kernels in the reverse order as `[wind_kernel, AdvectionRK2]` +Then also run a simulation where we apply the Kernels in the reverse order as `[windRK2, AdvectionRK2]` ```{code-cell} :tags: [hide-output] @@ -126,7 +128,7 @@ output_file_reverse = parcels.ParticleFile( path="wind_then_advection.parquet", outputdt=np.timedelta64(6,"h") ) pset_reverse.execute( - [wind_kernel, parcels.kernels.AdvectionRK2], + [windRK2, parcels.kernels.AdvectionRK2], runtime=np.timedelta64(5,"D"), dt=np.timedelta64(1,"h"), output_file=output_file_reverse, diff --git a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb index 3d473ecd8..31e1230c0 100644 --- a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb +++ b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb @@ -52,19 +52,24 @@ "ds_fields.load() # load the dataset into memory\n", "\n", "# Create an idealised wind field and add it to the dataset\n", - "tdim, ydim, xdim = (\n", + "tdim, zdim, ydim, xdim = (\n", " len(ds_fields.time),\n", + " len(ds_fields.depth),\n", " len(ds_fields.latitude),\n", " len(ds_fields.longitude),\n", ")\n", "ds_fields[\"UWind\"] = xr.DataArray(\n", - " data=np.ones((tdim, ydim, xdim)) * np.sin(ds_fields.latitude.values)[None, :, None],\n", - " coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n", + " data=0.5\n", + " * np.ones((tdim, zdim, ydim, xdim))\n", + " * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[\n", + " None, None, :, None\n", + " ],\n", + " coords=[ds_fields.time, ds_fields.depth, ds_fields.latitude, ds_fields.longitude],\n", ")\n", "\n", "ds_fields[\"VWind\"] = xr.DataArray(\n", - " data=np.zeros((tdim, ydim, xdim)),\n", - " coords=[ds_fields.time, ds_fields.latitude, ds_fields.longitude],\n", + " data=np.zeros((tdim, zdim, ydim, xdim)),\n", + " coords=[ds_fields.time, ds_fields.depth, ds_fields.latitude, ds_fields.longitude],\n", ")" ] }, @@ -126,12 +131,14 @@ "\n", "npart = 10\n", "z = np.repeat(ds_fields.depth[0].values, npart)\n", - "lons = np.repeat(31, npart)\n", + "lons = np.repeat(32.2, npart)\n", "lats = np.linspace(-32.5, -30.5, npart)\n", "\n", "pset = parcels.ParticleSet(fieldset, pclass=parcels.Particle, z=z, y=lats, x=lons)\n", "output_file = parcels.ParticleFile(\n", - " path=\"summed_advection_wind.parquet\", outputdt=np.timedelta64(6, \"h\")\n", + " path=\"summed_advection_wind.parquet\",\n", + " outputdt=np.timedelta64(6, \"h\"),\n", + " mode=\"w\",\n", ")\n", "pset.execute(\n", " [parcels.kernels.AdvectionRK2],\n", @@ -158,8 +165,9 @@ "source": [ "# Plot the resulting particle trajectories overlapped for both cases\n", "summed_advection_wind = parcels.read_particlefile(\"summed_advection_wind.parquet\")\n", + "fig, ax = plt.subplots(figsize=(5, 3))\n", "for traj in summed_advection_wind.partition_by(\"particle_id\", maintain_order=True):\n", - " plt.plot(traj[\"x\"], traj[\"y\"], \"-\")\n", + " ax.plot(traj[\"x\"], traj[\"y\"], \"-\")\n", "plt.show()" ] } From b071b6def920834bea59741370eb09c56e0aece2 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 10 Jul 2026 11:50:16 +0200 Subject: [PATCH 2/3] Update tutorial_manipulating_field_data.ipynb --- docs/user_guide/examples/tutorial_manipulating_field_data.ipynb | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb index 31e1230c0..0be1a7106 100644 --- a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb +++ b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb @@ -138,7 +138,6 @@ "output_file = parcels.ParticleFile(\n", " path=\"summed_advection_wind.parquet\",\n", " outputdt=np.timedelta64(6, \"h\"),\n", - " mode=\"w\",\n", ")\n", "pset.execute(\n", " [parcels.kernels.AdvectionRK2],\n", From fe6dc512da88510868b2a56efddb9062a918c004 Mon Sep 17 00:00:00 2001 From: Erik van Sebille Date: Fri, 10 Jul 2026 11:57:39 +0200 Subject: [PATCH 3/3] Remove time and depth dimension from field --- .../tutorial_manipulating_field_data.ipynb | 19 ++++++------------- 1 file changed, 6 insertions(+), 13 deletions(-) diff --git a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb index 0be1a7106..a5f9dd104 100644 --- a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb +++ b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb @@ -52,24 +52,17 @@ "ds_fields.load() # load the dataset into memory\n", "\n", "# Create an idealised wind field and add it to the dataset\n", - "tdim, zdim, ydim, xdim = (\n", - " len(ds_fields.time),\n", - " len(ds_fields.depth),\n", - " len(ds_fields.latitude),\n", - " len(ds_fields.longitude),\n", - ")\n", + "ydim, xdim = len(ds_fields.latitude), len(ds_fields.longitude)\n", "ds_fields[\"UWind\"] = xr.DataArray(\n", " data=0.5\n", - " * np.ones((tdim, zdim, ydim, xdim))\n", - " * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[\n", - " None, None, :, None\n", - " ],\n", - " coords=[ds_fields.time, ds_fields.depth, ds_fields.latitude, ds_fields.longitude],\n", + " * np.ones((ydim, xdim))\n", + " * np.sin(ds_fields.latitude.values - ds_fields.latitude.values.mean())[:, None],\n", + " coords=[ds_fields.latitude, ds_fields.longitude],\n", ")\n", "\n", "ds_fields[\"VWind\"] = xr.DataArray(\n", - " data=np.zeros((tdim, zdim, ydim, xdim)),\n", - " coords=[ds_fields.time, ds_fields.depth, ds_fields.latitude, ds_fields.longitude],\n", + " data=np.zeros((ydim, xdim)),\n", + " coords=[ds_fields.latitude, ds_fields.longitude],\n", ")" ] },