diff --git a/docs/getting_started/explanation_concepts.md b/docs/getting_started/explanation_concepts.md index c3117ed776..8153149725 100644 --- a/docs/getting_started/explanation_concepts.md +++ b/docs/getting_started/explanation_concepts.md @@ -143,7 +143,7 @@ def Age(particles, fieldset): particles.age += particles.dt # define all kernels to be executed on particles using an (ordered) list -kernels = [Age, parcels.kernels.AdvectionRK4] +kernels = [Age, parcels.kernels.AdvectionRK2] ``` ```{note} @@ -151,7 +151,7 @@ Every Kernel must be a function with the following (and only those) arguments: ` ``` ```{warning} -We have to be careful with writing kernels for vector fields on Curvilinear grids. While Parcels automatically rotates the "U" and "V" field when necessary, this is not the case for other fields such as Stokes drift. [This guide](../user_guide/examples/tutorial_nemo_curvilinear.ipynb) describes how to use a curvilinear grid in Parcels. +We have to be careful with kernels that sample velocities on "spherical" grids (so with longitude and latitude in degrees). Parcels can automatically convert velocities from m s-1 to degrees s-1, but only when using `VectorFields`. [This guide](../user_guide/examples/tutorial_velocityconversion.ipynb) describes how to use velocities on a "spherical" grid in Parcels. ``` ```{admonition} 📖 Read more about the Kernel loop @@ -159,7 +159,7 @@ We have to be careful with writing kernels for vector fields on Curvilinear grid - [The Kernel loop](../user_guide/examples/explanation_kernelloop.md) ``` -```{admonition} 🖥️ Learn how to write kernels +```{admonition} 🖥️ Learn how to write Kernels :class: seealso - [Sample fields like temperature](../user_guide/examples/tutorial_sampling.ipynb). - [Mimic the behaviour of ARGO floats](../user_guide/examples/tutorial_Argofloats.ipynb). diff --git a/docs/getting_started/tutorial_output.ipynb b/docs/getting_started/tutorial_output.ipynb index acf54cafac..0b7def37c1 100644 --- a/docs/getting_started/tutorial_output.ipynb +++ b/docs/getting_started/tutorial_output.ipynb @@ -120,7 +120,7 @@ "outputs": [], "source": [ "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", + " parcels.kernels.AdvectionRK2,\n", " runtime=np.timedelta64(48, \"h\"),\n", " dt=np.timedelta64(5, \"m\"),\n", " output_file=output_file,\n", diff --git a/docs/getting_started/tutorial_quickstart.md b/docs/getting_started/tutorial_quickstart.md index 348309eaf0..60d5ced099 100644 --- a/docs/getting_started/tutorial_quickstart.md +++ b/docs/getting_started/tutorial_quickstart.md @@ -54,6 +54,13 @@ ds_fset = parcels.convert.copernicusmarine_to_sgrid(fields=fields) fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset) ``` +You can inspect the `fieldset` by simply printing it: + +```{code-cell} +:tags: [hide-output] +print(fieldset) +``` + The subset contains a region of the Agulhas current along the southeastern coast of Africa: ```{code-cell} @@ -85,6 +92,15 @@ pset = parcels.ParticleSet( ) ``` +Again, you can inspect the `pset` by printing it: + +```{code-cell} +:tags: [hide-output] +print(pset) +``` + +And you can plot the particles on top of the temperature and velocity field: + ```{code-cell} temperature = ds_fields.isel(time=0, depth=0).thetao.plot(cmap="magma") velocity = ds_fields.isel(time=0, depth=0).plot.quiver(x="longitude", y="latitude", u="uo", v="vo") diff --git a/docs/user_guide/examples/explanation_kernelloop.md b/docs/user_guide/examples/explanation_kernelloop.md index 607566c4c8..bb2d20743b 100644 --- a/docs/user_guide/examples/explanation_kernelloop.md +++ b/docs/user_guide/examples/explanation_kernelloop.md @@ -4,17 +4,17 @@ kernelspec: name: python3 --- -# 📖 The Parcels Kernel loop +# 📖 Kernel loop -On this page we discuss Parcels' execution loop, and what happens under the hood when you combine multiple Kernels. +On this page we discuss how Parcels executes the Kernel loop, and what happens under the hood when you combine multiple Kernels. -This is not very relevant when you only use the built-in Advection kernels, but can be important when you are writing and combining your own Kernels! +This is not very relevant when you only use the built-in Advection Kernels, but can be important when you are writing and combining your own Kernels! ## Background When you run a Parcels simulation (i.e. a call to `pset.execute()`), the Kernel loop is the main part of the code that is executed. This part of the code loops through time and executes the Kernels for all particle. -In order to make sure that the displacements of a particle in the different Kernels can be summed, all Kernels add to a _change_ in position (`particles.dlon`, `particles.dlat`, and `particles.dz`). This is important, because there are situations where movement kernels would otherwise not commute. Take the example of advecting particles by currents _and_ winds. If the particle would first be moved by the currents and then by the winds, the result could be different from first moving by the winds and then by the currents. Instead, by summing the _changes_ in position, the ordering of the Kernels has no consequence on the particle displacement. +In order to make sure that the displacements of a particle in the different Kernels can be summed, all Kernels add to a _change_ in position (`particles.dlon`, `particles.dlat`, and `particles.dz`). This is important, because there are situations where movement Kernels would otherwise not commute. Take the example of advecting particles by currents _and_ winds. If the particle would first be moved by the currents and then by the winds, the result could be different from first moving by the winds and then by the currents. Instead, by summing the _changes_ in position, the ordering of the Kernels has no consequence on the particle displacement. ## Basic implementation @@ -91,7 +91,7 @@ windvector = parcels.VectorField( fieldset.add_field(windvector) ``` -Now we define a wind kernel that uses a forward Euler method to apply the wind forcing. Note that we update the `particles.dlon` and `particles.dlat` variables, rather than `particles.lon` and `particles.lat` directly. +Now we define a wind Kernel that uses a forward Euler method to apply the wind forcing. Note that we update the `particles.dlon` and `particles.dlat` variables, rather than `particles.lon` and `particles.lat` directly. ```{code-cell} def wind_kernel(particles, fieldset): @@ -100,7 +100,7 @@ def wind_kernel(particles, fieldset): particles.dlat += vwind * particles.dt ``` -First run a simulation where we apply kernels as `[AdvectionRK4, wind_kernel]` +First run a simulation where we apply Kernels as `[AdvectionRK2, wind_kernel]` ```{code-cell} :tags: [hide-output] @@ -114,14 +114,14 @@ output_file = parcels.ParticleFile( store="advection_then_wind.zarr", outputdt=np.timedelta64(6,'h') ) pset.execute( - [parcels.kernels.AdvectionRK4, wind_kernel], + [parcels.kernels.AdvectionRK2, wind_kernel], 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, AdvectionRK4]` +Then also run a simulation where we apply the Kernels in the reverse order as `[wind_kernel, AdvectionRK2]` ```{code-cell} :tags: [hide-output] @@ -132,7 +132,7 @@ output_file_reverse = parcels.ParticleFile( store="wind_then_advection.zarr", outputdt=np.timedelta64(6,"h") ) pset_reverse.execute( - [wind_kernel, parcels.kernels.AdvectionRK4], + [wind_kernel, 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_Argofloats.ipynb b/docs/user_guide/examples/tutorial_Argofloats.ipynb index 338541676d..0a37193ce7 100644 --- a/docs/user_guide/examples/tutorial_Argofloats.ipynb +++ b/docs/user_guide/examples/tutorial_Argofloats.ipynb @@ -15,7 +15,7 @@ "source": [ "This tutorial shows how simple it is to construct a Kernel in Parcels that mimics the [vertical movement of Argo floats](https://www.aoml.noaa.gov/phod/argo/images/argo_float_mission.jpg).\n", "\n", - "We first define the kernels for each phase of the Argo cycle." + "We first define the Kernel for the vertical movement cycle of the Argo float." ] }, { @@ -90,7 +90,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "And then we can run Parcels with this 'custom kernel'.\n", + "And then we can run Parcels with this `ArgoVerticalMovement` Kernel.\n", "\n", "Below we use the horizontal velocity fields of CopernicusMarine, which are provided as example_data with Parcels.\n" ] @@ -152,10 +152,10 @@ " z=[fieldset.mindepth],\n", ")\n", "\n", - "# combine Argo vertical movement kernel with built-in Advection kernel\n", + "# combine Argo vertical movement Kernel with built-in Advection Kernel\n", "kernels = [\n", " ArgoVerticalMovement,\n", - " parcels.kernels.AdvectionRK4,\n", + " parcels.kernels.AdvectionRK2,\n", "]\n", "\n", "# Create a ParticleFile object to store the output\n", @@ -165,7 +165,7 @@ " chunks=(1, 500), # setting to write in chunks of 500 observations\n", ")\n", "\n", - "# Now execute the kernels for 30 days, saving data every 30 minutes\n", + "# Now execute the Kernels for 30 days, saving data every 30 minutes\n", "pset.execute(\n", " kernels,\n", " runtime=timedelta(days=30),\n", @@ -248,7 +248,7 @@ ], "metadata": { "kernelspec": { - "display_name": "test-notebooks", + "display_name": "docs", "language": "python", "name": "python3" }, @@ -262,7 +262,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.0" + "version": "3.14.2" } }, "nbformat": 4, diff --git a/docs/user_guide/examples/tutorial_delaystart.ipynb b/docs/user_guide/examples/tutorial_delaystart.ipynb index 81ac797f4b..03bdb6a589 100644 --- a/docs/user_guide/examples/tutorial_delaystart.ipynb +++ b/docs/user_guide/examples/tutorial_delaystart.ipynb @@ -120,7 +120,7 @@ ")\n", "\n", "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", + " parcels.kernels.AdvectionRK2,\n", " runtime=np.timedelta64(24, \"h\"),\n", " dt=np.timedelta64(5, \"m\"),\n", " output_file=output_file,\n", @@ -259,7 +259,7 @@ ")\n", "\n", "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", + " parcels.kernels.AdvectionRK2,\n", " runtime=np.timedelta64(24, \"h\"),\n", " dt=np.timedelta64(5, \"h\"),\n", " output_file=output_file,\n", @@ -366,7 +366,7 @@ "\n", "output_file = parcels.ParticleFile(outfilepath, outputdt=np.timedelta64(2, \"h\"))\n", "pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", + " parcels.kernels.AdvectionRK2,\n", " endtime=ds_fields.time.values[0] + np.timedelta64(4, \"h\"),\n", " dt=np.timedelta64(1, \"h\"),\n", " output_file=output_file,\n", @@ -423,7 +423,7 @@ " )\n", " output_file = parcels.ParticleFile(outfilepath, outputdt=np.timedelta64(2, \"h\"))\n", " pset.execute(\n", - " parcels.kernels.AdvectionRK4,\n", + " parcels.kernels.AdvectionRK2,\n", " runtime=np.timedelta64(4, \"h\"),\n", " dt=np.timedelta64(1, \"h\"),\n", " output_file=output_file,\n", diff --git a/docs/user_guide/examples/tutorial_diffusion.ipynb b/docs/user_guide/examples/tutorial_diffusion.ipynb index eca4c35092..6e42b4ce91 100644 --- a/docs/user_guide/examples/tutorial_diffusion.ipynb +++ b/docs/user_guide/examples/tutorial_diffusion.ipynb @@ -545,7 +545,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "In the example, particles are released at one location and simulated for 2 days using advection (`AdvectionRK4`) and diffusion (`smagdiff`) kernels." + "In the example, particles are released at one location and simulated for 2 days using advection (`AdvectionRK2`) and diffusion (`smagdiff`) kernels." ] }, { @@ -578,7 +578,7 @@ "np.random.seed(1636) # Random seed for reproducibility\n", "\n", "pset.execute(\n", - " [parcels.kernels.AdvectionRK4, smagdiff],\n", + " [parcels.kernels.AdvectionRK2, smagdiff],\n", " runtime=np.timedelta64(2, \"D\"),\n", " dt=np.timedelta64(5, \"m\"),\n", " output_file=output_file,\n", diff --git a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb index f7b61b6f72..f5dc55571b 100644 --- a/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb +++ b/docs/user_guide/examples/tutorial_manipulating_field_data.ipynb @@ -135,7 +135,7 @@ " store=\"summed_advection_wind.zarr\", outputdt=np.timedelta64(6, \"h\")\n", ")\n", "pset.execute(\n", - " [parcels.kernels.AdvectionRK4],\n", + " [parcels.kernels.AdvectionRK2],\n", " runtime=np.timedelta64(5, \"D\"),\n", " dt=np.timedelta64(1, \"h\"),\n", " output_file=output_file,\n", diff --git a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb b/docs/user_guide/examples/tutorial_nemo.ipynb similarity index 59% rename from docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb rename to docs/user_guide/examples/tutorial_nemo.ipynb index ba26f3558e..d56aa5194f 100644 --- a/docs/user_guide/examples/tutorial_nemo_curvilinear.ipynb +++ b/docs/user_guide/examples/tutorial_nemo.ipynb @@ -5,7 +5,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# 🖥️ Curvilinear grids" + "# 🖥️ Nemo tutorial\n" ] }, { @@ -13,13 +13,31 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "One of the features of Parcels is that it can directly and natively work with `Field` data discretised on C-grids. These C grids are very popular in OGCMs, so velocity fields outputted by OGCMs are often provided on such grids, except if they have been firstly re-interpolated on a A grid.\n", + "\n", + "More information about C-grid interpolation can be found in [Delandmeter et al., 2019](https://www.geosci-model-dev-discuss.net/gmd-2018-339/).\n", + "An example of such a discretisation is the NEMO model, which is one of the models for which we provide support in Parcels. Other models are CROCO, ROMS, and MITgcm.\n", + "\n", + "```{note}\n", + "_How to know if your data is discretised on a C grid?_ The best way is to read the documentation that comes with the data. Alternatively, an easy check is to assess the coordinates of the U, V and W fields: \n", + "- for an **A grid**, U, V and W are distributed on the same nodes, such that the coordinates are the same. \n", + "- for a **C grid**, there is a shift of half a cell between the different variables.\n", + "```" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Curvilinear C-Grids\n", + "\n", "Parcels supports [curvilinear grids](https://www.nemo-ocean.eu/doc/node108.html) such as those used in the [NEMO models](https://www.nemo-ocean.eu/).\n", "\n", "```{note}\n", "TODO: make explicit how Parcels determines rotation\n", "```\n", "\n", - "We will be using the example dataset `NemoCurvilinear_data`. These fields are a purely zonal flow on an aqua-planet (so zonal-velocity is 1 m/s and meridional-velocity is 0 m/s everywhere, and no land). However, because of the curvilinear grid, the `U` and `V` fields vary for the rotated gridcells north of 20N.\n" + "We will be using the example dataset `NemoCurvilinear_data`. These fields are a purely zonal flow on an aqua-planet (so zonal-velocity is 1 m s-1 and meridional-velocity is 0 m s-1 everywhere, and no land). However, because of the curvilinear grid, the `U` and `V` fields vary for the rotated grid cells north of 20N." ] }, { @@ -36,25 +54,10 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "We can create a `FieldSet` just like we do for normal grids.\n", - "Note that NEMO is discretised on a C-grid. U and V velocities are not located on the same nodes (see https://www.nemo-ocean.eu/doc/node19.html ).\n", - "\n", - "```\n", - " __V1__\n", - "| |\n", - "U0 U1\n", - "|__V0__|\n", - "```\n", - "\n", - "To interpolate U, V velocities on the C-grid, Parcels needs to read the f-nodes, which are located on the corners of the cells (for indexing details: https://www.nemo-ocean.eu/doc/img360.png).\n", - "\n", - "```{note}\n", - "TODO: add link to grid indexing explanation once implemented in v4\n", - "```\n" + "We use the `parcels.convert.nemo_to_sgrid()` function to convert the NEMO curvilinear grid to a dataset that is SGrid-compliant. We then pass this new dataset into `parcels.FieldSet.from_sgrid_conventions()`. This function automatically detects the C-grid structure and sets up the grid accordingly." ] }, { @@ -80,11 +83,10 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "And we can plot the `U` field.\n" + "And we can plot the `U` field." ] }, { @@ -110,11 +112,10 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "As you see above, the `U` field indeed is 1 m/s south of 20N, but varies with longitude and latitude north of that. We can confirm that Parcels will take care to rotate the `U` and `V` fields by doing a field evaluation at (60N, 50E):" + "As you see above, the `U` field indeed is 1 m s-1 south of 20N, but varies with longitude and latitude north of that. We can confirm that Parcels will take care to rotate the `U` and `V` fields by doing a field evaluation at (60N, 50E):" ] }, { @@ -168,11 +169,10 @@ ] }, { - "attachments": {}, "cell_type": "markdown", "metadata": {}, "source": [ - "And then we can plot these trajectories. As expected, all trajectories go exactly zonal and due to the curvature of the earth, ones at higher latitude move more degrees eastward (even though the distance in km is equal for all particles). These particles at high latitudes cross the antimeridian (180 deg E) and keep going east.\n" + "And then we can plot these trajectories. As expected, all trajectories go exactly zonal and due to the curvature of the earth, ones at higher latitude move more degrees eastward (even though the distance in km is equal for all particles). These particles at high latitudes cross the antimeridian (180 deg E) and keep going east." ] }, { @@ -196,6 +196,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ + "### Handling longitude wrapping\n", "If we want the `particles.lon` to stay within `[-180,180]` (or `[0,360]`), we can either do this in post-processing, or add a periodic boundary Kernel:" ] }, @@ -267,11 +268,108 @@ "\n", "plt.show()" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Running 3D Nemo simulations\n", + "\n", + "Here, we focus on 3D fields. Basically, it is a straightforward extension of the 2D example above." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We again use the combination of the `parcels.convert.nemo_to_sgrid` function and the `parcels.FieldSet.from_sgrid_conventions()` to create a FieldSet object with the correct `Grid` information for the NEMO data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "data_folder = parcels.download_example_dataset(\"NemoNorthSeaORCA025-N006_data\")\n", + "ds_fields = xr.open_mfdataset(\n", + " data_folder.glob(\"ORCA*.nc\"),\n", + " data_vars=\"minimal\",\n", + " coords=\"minimal\",\n", + " compat=\"override\",\n", + ")\n", + "ds_coords = xr.open_dataset(data_folder / \"coordinates.nc\", decode_times=False)\n", + "ds_fset = parcels.convert.nemo_to_sgrid(\n", + " fields={\"U\": ds_fields[\"uo\"], \"V\": ds_fields[\"vo\"], \"W\": ds_fields[\"wo\"]},\n", + " coords=ds_coords,\n", + ")\n", + "fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "The code below is an example of how to then create a 3D simulation with particles, starting on a meridional line through the North Sea from the mouth of the river Rhine at 100m depth, and advecting them through the North Sea using the `AdvectionRK2_3D`" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "tags": [ + "hide-output" + ] + }, + "outputs": [], + "source": [ + "npart = 10\n", + "pset = parcels.ParticleSet(\n", + " fieldset=fieldset,\n", + " lon=np.linspace(1.9, 3.4, npart),\n", + " lat=np.linspace(65, 51.6, npart),\n", + " z=100 * np.ones(npart),\n", + ")\n", + "\n", + "pfile = parcels.ParticleFile(\n", + " store=\"output_nemo3D.zarr\", outputdt=np.timedelta64(1, \"D\")\n", + ")\n", + "\n", + "pset.execute(\n", + " parcels.kernels.AdvectionRK2_3D,\n", + " endtime=fieldset.time_interval.right,\n", + " dt=np.timedelta64(6, \"h\"),\n", + " output_file=pfile,\n", + ")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "We can then plot the trajectories on top of the surface U field" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "field = fieldset.U.data[0, 0, :, :]\n", + "field = field.where(field != 0, np.nan) # Mask land values for better plotting\n", + "plt.pcolormesh(fieldset.U.grid.lon, fieldset.U.grid.lat, field, cmap=\"RdBu\")\n", + "\n", + "ds_out = xr.open_zarr(\"output_nemo3D.zarr\")\n", + "plt.scatter(ds_out.lon.T, ds_out.lat.T, c=-ds_out.z.T, marker=\".\")\n", + "plt.colorbar(label=\"Depth (m)\")\n", + "plt.show()" + ] } ], "metadata": { "kernelspec": { - "display_name": "test-notebooks-latest", + "display_name": "parcels", "language": "python", "name": "python3" }, @@ -285,7 +383,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.13.9" + "version": "3.12.3" } }, "nbformat": 4, diff --git a/docs/user_guide/examples/tutorial_nemo_3D.ipynb b/docs/user_guide/examples/tutorial_nemo_3D.ipynb deleted file mode 100644 index 38d79a5bab..0000000000 --- a/docs/user_guide/examples/tutorial_nemo_3D.ipynb +++ /dev/null @@ -1,153 +0,0 @@ -{ - "cells": [ - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# 🖥️ Nemo 3D tutorial\n" - ] - }, - { - "attachments": {}, - "cell_type": "markdown", - "metadata": {}, - "source": [ - "One of the features of Parcels is that it can directly and natively work with `Field` data discretised on C-grids. These C grids are very popular in OGCMs, so velocity fields outputted by OGCMs are often provided on such grids, except if they have been firstly re-interpolated on a A grid.\n", - "\n", - "More information about C-grid interpolation can be found in [Delandmeter et al., 2019](https://www.geosci-model-dev-discuss.net/gmd-2018-339/).\n", - "An example of such a discretisation is the NEMO model, which is one of the models supported in Parcels. A tutorial teaching how to [interpolate 2D data on a NEMO grid](https://docs.oceanparcels.org/en/latest/examples/tutorial_nemo_curvilinear.html) is available within Parcels.\n", - "\n", - "```{note}\n", - "_How to know if your data is discretised on a C grid?_ The best way is to read the documentation that comes with the data. Alternatively, an easy check is to assess the coordinates of the U, V and W fields: for an A grid, U, V and W are distributed on the same nodes, such that the coordinates are the same. For a C grid, there is a shift of half a cell between the different variables.\n", - "```\n", - "\n", - "Here, we focus on 3D fields. Basically, it is a straightforward extension of the 2D example, but it is very easy to make a mistake in the setup of the vertical discretisation that would affect the interpolation scheme.\n", - "\n", - "For the C-grid interpolation in Parcels to work properly, it is important that `U`, `V` and `W` are on the same grid. All other tracers (e.g. `temperature`, `salinity`) should also be on this same grid. So even though these tracers are computed by NEMO on the T-points, Parcels expects them on the f-points (`glamf`, `gphif` and `depthw`). Parcels then under the hood makes sure the interpolation of these tracers is done correctly." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import numpy as np\n", - "import xarray as xr\n", - "\n", - "import parcels" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Parcels v4 comes with a `convert.nemo_to_sgrid` method that automatically sets up the correct `Grid` information for both 2D and 3D NEMO data." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "data_folder = parcels.download_example_dataset(\"NemoNorthSeaORCA025-N006_data\")\n", - "ds_fields = xr.open_mfdataset(\n", - " data_folder.glob(\"ORCA*.nc\"),\n", - " data_vars=\"minimal\",\n", - " coords=\"minimal\",\n", - " compat=\"override\",\n", - ")\n", - "ds_coords = xr.open_dataset(data_folder / \"coordinates.nc\", decode_times=False)\n", - "ds_fset = parcels.convert.nemo_to_sgrid(\n", - " fields={\"U\": ds_fields[\"uo\"], \"V\": ds_fields[\"vo\"], \"W\": ds_fields[\"wo\"]},\n", - " coords=ds_coords,\n", - ")\n", - "fieldset = parcels.FieldSet.from_sgrid_conventions(ds_fset)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The code below is an example of how to then create a 3D simulation with particles, starting on a meridional line through the North Sea from the mouth of the river Rhine at 100m depth, and advecting them through the North Sea using the `AdvectionRK2_3D`" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "tags": [ - "hide-output" - ] - }, - "outputs": [], - "source": [ - "npart = 10\n", - "pset = parcels.ParticleSet(\n", - " fieldset=fieldset,\n", - " lon=np.linspace(1.9, 3.4, npart),\n", - " lat=np.linspace(65, 51.6, npart),\n", - " z=100 * np.ones(npart),\n", - ")\n", - "\n", - "pfile = parcels.ParticleFile(\n", - " store=\"output_nemo3D.zarr\", outputdt=np.timedelta64(1, \"D\")\n", - ")\n", - "\n", - "pset.execute(\n", - " parcels.kernels.AdvectionRK2_3D,\n", - " endtime=fieldset.time_interval.right,\n", - " dt=np.timedelta64(6, \"h\"),\n", - " output_file=pfile,\n", - ")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "We can then plot the trajectories on top of the surface U field" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "field = fieldset.U.data[0, 0, :, :]\n", - "field = field.where(field != 0, np.nan) # Mask land values for better plotting\n", - "plt.pcolormesh(fieldset.U.grid.lon, fieldset.U.grid.lat, field, cmap=\"RdBu\")\n", - "\n", - "ds_out = xr.open_zarr(\"output_nemo3D.zarr\")\n", - "plt.scatter(ds_out.lon.T, ds_out.lat.T, c=-ds_out.z.T, marker=\".\")\n", - "plt.colorbar(label=\"Depth (m)\")\n", - "plt.show()" - ] - } - ], - "metadata": { - "kernelspec": { - "display_name": "parcels", - "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.12.3" - } - }, - "nbformat": 4, - "nbformat_minor": 4 -} diff --git a/docs/user_guide/examples/tutorial_sampling.ipynb b/docs/user_guide/examples/tutorial_sampling.ipynb index 28b9c3954f..27ff5af9de 100644 --- a/docs/user_guide/examples/tutorial_sampling.ipynb +++ b/docs/user_guide/examples/tutorial_sampling.ipynb @@ -140,7 +140,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We can then sample and Advect by combining the `SampleT` and `AdvectionRK4` kernels in a list. Note that the order does not matter.\n" + "We can then sample and advect by combining the `SampleT` and `AdvectionRK2` kernels in a list. Note that the order does not matter.\n" ] }, { @@ -160,7 +160,7 @@ "output_file = parcels.ParticleFile(\"sampletemp.zarr\", outputdt=timedelta(hours=1))\n", "\n", "pset.execute(\n", - " [parcels.kernels.AdvectionRK4, SampleT],\n", + " [parcels.kernels.AdvectionRK2, SampleT],\n", " runtime=timedelta(hours=30),\n", " dt=timedelta(minutes=5),\n", " output_file=output_file,\n", @@ -353,7 +353,7 @@ "output_file = parcels.ParticleFile(\"writeonce.zarr\", outputdt=timedelta(hours=1))\n", "\n", "pset.execute(\n", - " [parcels.kernels.AdvectionRK4, SampleT],\n", + " [parcels.kernels.AdvectionRK2, SampleT],\n", " runtime=timedelta(hours=24),\n", " dt=timedelta(minutes=5),\n", " output_file=output_file,\n", diff --git a/docs/user_guide/index.md b/docs/user_guide/index.md index 3e924ab6b1..35ecd9a7a8 100644 --- a/docs/user_guide/index.md +++ b/docs/user_guide/index.md @@ -25,11 +25,10 @@ The tutorials written for Parcels v3 are currently being updated for Parcels v4. :caption: Set up FieldSets :titlesonly: examples/explanation_grids.md -examples/tutorial_nemo_curvilinear.ipynb -examples/tutorial_nemo_3D.ipynb +examples/tutorial_nemo.ipynb examples/tutorial_velocityconversion.ipynb examples/tutorial_nestedgrids.ipynb -examples/tutorial_summingfields.ipynb +examples/tutorial_manipulating_field_data.ipynb ``` @@ -49,10 +48,6 @@ examples/tutorial_delaystart.ipynb examples/explanation_kernelloop.md examples/tutorial_sampling.ipynb examples/tutorial_statuscodes.ipynb -examples/tutorial_gsw_density.ipynb -examples/tutorial_Argofloats.ipynb -examples/tutorial_diffusion.ipynb -examples/tutorial_interaction.ipynb ``` ```{toctree} @@ -90,7 +85,12 @@ examples/tutorial_dt_integrators.ipynb ```{toctree} -:caption: Worked examples +:caption: Example Kernels +:titlesonly: +examples/tutorial_gsw_density.ipynb +examples/tutorial_Argofloats.ipynb +examples/tutorial_diffusion.ipynb +examples/tutorial_interaction.ipynb ``` diff --git a/docs/user_guide/v4-migration.md b/docs/user_guide/v4-migration.md index 6e2f4743df..263e006e8c 100644 --- a/docs/user_guide/v4-migration.md +++ b/docs/user_guide/v4-migration.md @@ -15,6 +15,7 @@ Version 4 of Parcels is unreleased at the moment. The information in this migrat - `math` functions should be replaced with array compatible equivalents (e.g., `math.sin` -> `np.sin`). Instead of `ParcelsRandom` you should use numpy's random functions. - `particle.depth` has been changed to `particles.z` to be consistent with the [CF conventions for trajectory data](https://cfconventions.org/cf-conventions/cf-conventions.html#trajectory-data), and to make Parcels also generalizable to atmospheric contexts. - The `InteractionKernel` class has been removed. Since normal Kernels now have access to _all_ particles, particle-particle interaction can be performed within normal Kernels. +- We added a new AdvectionRK2 Kernel. The AdvectionRK4 kernel is still available, but RK2 is now the recommended default advection scheme as it is faster while the accuracy is comparable for most applications. See also the Choosing an integration method tutorial. ## FieldSet