Skip to content

Commit 481414e

Browse files
reint-fischerreint-fischer
authored
Merge branch 'v4-dev' into structure_overview
2 parents 8b64522 + db0f8fe commit 481414e

File tree

10 files changed

+486
-87
lines changed

10 files changed

+486
-87
lines changed

‎docs/user_guide/examples/tutorial_interpolation.ipynb‎

Lines changed: 234 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -131,7 +131,10 @@
131131
" fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten()\n",
132132
" )\n",
133133
" pset[p_interp.__name__].execute(\n",
134-
" SampleP, runtime=np.timedelta64(1, \"D\"), dt=np.timedelta64(1, \"D\")\n",
134+
" SampleP,\n",
135+
" runtime=np.timedelta64(1, \"D\"),\n",
136+
" dt=np.timedelta64(1, \"D\"),\n",
137+
" verbose_progress=False,\n",
135138
" )"
136139
]
137140
},
@@ -239,26 +242,245 @@
239242
"```"
240243
]
241244
},
245+
{
246+
"cell_type": "markdown",
247+
"metadata": {},
248+
"source": [
249+
"### Interpolation at boundaries\n",
250+
"In some cases, we need to implement specific boundary conditions, for example to prevent particles from \n",
251+
"getting \"stuck\" near land. [This guide](../examples_v3/documentation_unstuck_Agrid.ipynb) describes \n",
252+
"how to implement this in parcels using `parcels.interpolators.XFreeslip` and `parcels.interpolators.XPartialslip`."
253+
]
254+
},
242255
{
243256
"cell_type": "markdown",
244257
"metadata": {},
245258
"source": [
246259
"## Interpolation on unstructured grids\n",
247-
"```{note}\n",
248-
"TODO: add example on simple unstructured fieldset\n",
249-
"```\n",
250-
"- UXPiecewiseConstantFace\n",
251-
"- UXPiecewiseLinearNode"
260+
"Parcels v4 supports the use of general circulation model output that is defined on unstructured grids. We include basic interpolators to help you get started, including\n",
261+
"- `UxPiecewiseConstantFace` - this interpolator implements piecewise constant interpolation and is appropriate for data that is registered to the face centers of the unstructured grid\n",
262+
"- `UxPiecewiseLinearNode` - this interpolator implements barycentric interpolation and is appropriate for data that is registered to the corner vertices of the unstructured grid faces\n",
263+
"\n",
264+
"To get started, we use a very simple generated `UxArray.UxDataset` that is included with Parcels."
265+
]
266+
},
267+
{
268+
"cell_type": "code",
269+
"execution_count": null,
270+
"metadata": {},
271+
"outputs": [],
272+
"source": [
273+
"from parcels._datasets.unstructured.generated import simple_small_delaunay\n",
274+
"\n",
275+
"ds = simple_small_delaunay(nx=4, ny=4)\n",
276+
"ds"
252277
]
253278
},
254279
{
255280
"cell_type": "markdown",
256281
"metadata": {},
257282
"source": [
258-
"## Interpolation at boundaries\n",
259-
"In some cases, we need to implement specific boundary conditions, for example to prevent particles from \n",
260-
"getting \"stuck\" near land. [This guide](../examples_v3/documentation_unstuck_Agrid.ipynb) describes \n",
261-
"how to implement this in parcels using `parcels.interpolators.XFreeslip` and `parcels.interpolators.XPartialslip`."
283+
"Next, we create the `Field` and `Fieldset` objects that will be used later in advancing particles. When creating a `Field` or `VectorField` object for unstructured grid data, we attach a `parcels.UxGrid` object and attach an `interp_method` to each object. For data that is defined on face centers, we use the `UxPiecewiseConstantFace` interpolator and for data that is defined on the face vertices, we use the `UxPiecewiseLinearNode` interpolator. In this example, we will look specifically at interpolating a tracer field that is defined by the same underlying analytical function, but is defined on both faces and vertices as separate fields."
284+
]
285+
},
286+
{
287+
"cell_type": "code",
288+
"execution_count": null,
289+
"metadata": {},
290+
"outputs": [],
291+
"source": [
292+
"import numpy as np\n",
293+
"\n",
294+
"import parcels\n",
295+
"\n",
296+
"grid = parcels.UxGrid(ds.uxgrid, ds.nz, mesh=\"flat\")\n",
297+
"\n",
298+
"Tnode = parcels.Field(\n",
299+
" \"T_node\",\n",
300+
" ds[\"T_node\"],\n",
301+
" grid,\n",
302+
" interp_method=parcels.interpolators.UxPiecewiseLinearNode,\n",
303+
")\n",
304+
"Tface = parcels.Field(\n",
305+
" \"T_face\",\n",
306+
" ds[\"T_face\"],\n",
307+
" grid,\n",
308+
" interp_method=parcels.interpolators.UxPiecewiseConstantFace,\n",
309+
")\n",
310+
"fieldset = parcels.FieldSet([Tnode, Tface])"
311+
]
312+
},
313+
{
314+
"cell_type": "markdown",
315+
"metadata": {},
316+
"source": [
317+
"We can now create particles to sample data using either the node registered or face registered data."
318+
]
319+
},
320+
{
321+
"cell_type": "code",
322+
"execution_count": null,
323+
"metadata": {},
324+
"outputs": [],
325+
"source": [
326+
"SampleParticle = parcels.Particle.add_variable(\n",
327+
" parcels.Variable(\"tracer\", dtype=np.float32, initial=np.nan)\n",
328+
")\n",
329+
"\n",
330+
"\n",
331+
"def SampleTracer_Node(particles, fieldset):\n",
332+
" particles.tracer = fieldset.T_node[particles]\n",
333+
"\n",
334+
"\n",
335+
"def SampleTracer_Face(particles, fieldset):\n",
336+
" particles.tracer = fieldset.T_face[particles]"
337+
]
338+
},
339+
{
340+
"cell_type": "markdown",
341+
"metadata": {},
342+
"source": [
343+
"Now, we can create a `ParticleSet` for each interpolation experiment. In one of the `ParticleSet` objects we use the `SampleTracer_Node` kernel to showcase interpolating with the `UxPiecewiseLinearNode` interpolator for node-registered data. In the other, we use the `SampleTracer_Face` to showcase interpolating with the `UxPiecewiseConstantFace` interpolator for face-registered data."
344+
]
345+
},
346+
{
347+
"cell_type": "code",
348+
"execution_count": null,
349+
"metadata": {},
350+
"outputs": [],
351+
"source": [
352+
"pset = {}\n",
353+
"\n",
354+
"xv, yv = np.meshgrid(np.linspace(0.2, 0.8, 8), np.linspace(0.2, 0.9, 8))\n",
355+
"\n",
356+
"pset[\"node\"] = parcels.ParticleSet(\n",
357+
" fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten()\n",
358+
")\n",
359+
"pset[\"node\"].execute(\n",
360+
" SampleTracer_Node,\n",
361+
" runtime=np.timedelta64(1, \"D\"),\n",
362+
" dt=np.timedelta64(1, \"D\"),\n",
363+
" verbose_progress=False,\n",
364+
")\n",
365+
"\n",
366+
"pset[\"face\"] = parcels.ParticleSet(\n",
367+
" fieldset, pclass=SampleParticle, lon=xv.flatten(), lat=yv.flatten()\n",
368+
")\n",
369+
"pset[\"face\"].execute(\n",
370+
" SampleTracer_Face,\n",
371+
" runtime=np.timedelta64(1, \"D\"),\n",
372+
" dt=np.timedelta64(1, \"D\"),\n",
373+
" verbose_progress=False,\n",
374+
")"
375+
]
376+
},
377+
{
378+
"cell_type": "markdown",
379+
"metadata": {},
380+
"source": [
381+
"Now, we can plot the node-registered (`T_node`) and face-registered (`T_face`) fields with the values interpolated onto the particles for each field. "
382+
]
383+
},
384+
{
385+
"cell_type": "code",
386+
"execution_count": null,
387+
"metadata": {},
388+
"outputs": [],
389+
"source": [
390+
"import matplotlib.pyplot as plt\n",
391+
"import matplotlib.tri as mtri\n",
392+
"import numpy as np\n",
393+
"\n",
394+
"fig, ax = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)\n",
395+
"\n",
396+
"x = ds.uxgrid.node_lon.values\n",
397+
"y = ds.uxgrid.node_lat.values\n",
398+
"tri = ds.uxgrid.face_node_connectivity.values\n",
399+
"triang = mtri.Triangulation(x, y, triangles=tri)\n",
400+
"\n",
401+
"\n",
402+
"# Shade the triangle faces using smooth shading between the node registered data\n",
403+
"T_node = np.squeeze(ds[\"T_node\"][0, 0, :].values)\n",
404+
"tpc = ax[0].tripcolor(\n",
405+
" triang,\n",
406+
" T_node,\n",
407+
" shading=\"gouraud\", # smooth shading\n",
408+
" cmap=\"viridis\",\n",
409+
" vmin=-1.0,\n",
410+
" vmax=1.0,\n",
411+
")\n",
412+
"# Plot the element edges\n",
413+
"ax[0].triplot(triang, color=\"black\", linewidth=0.5)\n",
414+
"\n",
415+
"# Plot the node registered data\n",
416+
"sc1 = ax[0].scatter(\n",
417+
" x, y, c=T_node, cmap=\"viridis\", s=150, edgecolors=\"black\", vmin=-1.0, vmax=1.0\n",
418+
")\n",
419+
"\n",
420+
"ax[0].scatter(\n",
421+
" pset[\"node\"].lon,\n",
422+
" pset[\"node\"].lat,\n",
423+
" c=pset[\"node\"].tracer,\n",
424+
" cmap=\"viridis\",\n",
425+
" edgecolors=\"k\",\n",
426+
" s=50,\n",
427+
" vmin=-1.0,\n",
428+
" vmax=1.0,\n",
429+
")\n",
430+
"\n",
431+
"cbar = fig.colorbar(sc1, ax=ax[0])\n",
432+
"cbar.set_label(\"T\")\n",
433+
"\n",
434+
"ax[0].set_aspect(\"equal\", \"box\")\n",
435+
"ax[0].set_xlabel(\"x\")\n",
436+
"ax[0].set_ylabel(\"y\")\n",
437+
"ax[0].set_title(\"Node Registered Data\")\n",
438+
"\n",
439+
"# # Plot the face registered data\n",
440+
"T_face = np.squeeze(ds[\"T_face\"][0, 0, :].values)\n",
441+
"assert triang.triangles.shape[0] == T_face.shape[0]\n",
442+
"tpc2 = ax[1].tripcolor(\n",
443+
" triang,\n",
444+
" facecolors=T_face,\n",
445+
" shading=\"flat\",\n",
446+
" cmap=\"viridis\",\n",
447+
" edgecolors=\"k\",\n",
448+
" linewidth=0.5,\n",
449+
" vmin=-1.0,\n",
450+
" vmax=1.0,\n",
451+
")\n",
452+
"# Plot the element edges\n",
453+
"ax[1].triplot(triang, color=\"black\", linewidth=0.5)\n",
454+
"\n",
455+
"# Plot the face registered data\n",
456+
"xf = ds.uxgrid.face_lon.values\n",
457+
"yf = ds.uxgrid.face_lat.values\n",
458+
"\n",
459+
"ax[1].scatter(\n",
460+
" pset[\"face\"].lon,\n",
461+
" pset[\"face\"].lat,\n",
462+
" c=pset[\"face\"].tracer,\n",
463+
" cmap=\"viridis\",\n",
464+
" edgecolors=\"k\",\n",
465+
" s=50,\n",
466+
" vmin=-1.0,\n",
467+
" vmax=1.0,\n",
468+
")\n",
469+
"\n",
470+
"cbar = fig.colorbar(tpc, ax=ax[1])\n",
471+
"cbar.set_label(\"T\")\n",
472+
"\n",
473+
"ax[1].set_aspect(\"equal\", \"box\")\n",
474+
"ax[1].set_xlabel(\"x\")\n",
475+
"ax[1].set_ylabel(\"y\")\n",
476+
"ax[1].set_title(\"Face Registered Data\")"
477+
]
478+
},
479+
{
480+
"cell_type": "markdown",
481+
"metadata": {},
482+
"source": [
483+
"In both plots the black lines show the edges that define the boundaries between the faces in the unstructured grid. For the node registered field, the background coloring is done using smooth shading based on the value of `T_node` at the corner nodes. For the face registered fields, each face is colored according to the value of `T_face` at the face center. The color of the particles is the value of the function that the particles take on via the corresponding interpolation method."
262484
]
263485
},
264486
{
@@ -272,7 +494,7 @@
272494
],
273495
"metadata": {
274496
"kernelspec": {
275-
"display_name": "test-notebooks",
497+
"display_name": "Python 3 (ipykernel)",
276498
"language": "python",
277499
"name": "python3"
278500
},
@@ -286,7 +508,7 @@
286508
"name": "python",
287509
"nbconvert_exporter": "python",
288510
"pygments_lexer": "ipython3",
289-
"version": "3.11.0"
511+
"version": "3.12.11"
290512
}
291513
},
292514
"nbformat": 4,

‎docs/user_guide/examples_v3/tutorial_stommel_uxarray.ipynb‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -112,7 +112,7 @@
112112
"\n",
113113
"In Parcels, grid searching is conducted with respect to the faces. In other words, when a grid index `ei` is provided to an interpolation method, this refers the face index `fi` at vertical layer `zi` (when unraveled). Within the interpolation method, the `field.grid.uxgrid.face_node_connectivity` attribute can be used to obtain the node indices that surround the face. Using these connectivity tables is necessary for properly indexing node registered data.\n",
114114
"\n",
115-
"For the example Stommel gyre dataset in this tutorial, the `u` and `v` velocity components are face registered (similar to FESOM). Parcels includes a nearest neighbor interpolator for face registered unstructured grid data through `Parcels.application_kernels.interpolation.UXPiecewiseConstantFace`. Below, we create the `Field`s `U` and `V` and associate them with the `UxGrid` we created previously and this interpolation method."
115+
"For the example Stommel gyre dataset in this tutorial, the `u` and `v` velocity components are face registered (similar to FESOM). Parcels includes a nearest neighbor interpolator for face registered unstructured grid data through `Parcels.application_kernels.interpolation.UxPiecewiseConstantFace`. Below, we create the `Field`s `U` and `V` and associate them with the `UxGrid` we created previously and this interpolation method."
116116
]
117117
},
118118
{
@@ -121,26 +121,26 @@
121121
"metadata": {},
122122
"outputs": [],
123123
"source": [
124-
"from parcels.application_kernels.interpolation import UXPiecewiseConstantFace\n",
124+
"from parcels.application_kernels.interpolation import UxPiecewiseConstantFace\n",
125125
"from parcels.field import Field\n",
126126
"\n",
127127
"U = Field(\n",
128128
" name=\"U\",\n",
129129
" data=ds.U,\n",
130130
" grid=grid,\n",
131-
" interp_method=UXPiecewiseConstantFace,\n",
131+
" interp_method=UxPiecewiseConstantFace,\n",
132132
")\n",
133133
"V = Field(\n",
134134
" name=\"V\",\n",
135135
" data=ds.V,\n",
136136
" grid=grid,\n",
137-
" interp_method=UXPiecewiseConstantFace,\n",
137+
" interp_method=UxPiecewiseConstantFace,\n",
138138
")\n",
139139
"P = Field(\n",
140140
" name=\"P\",\n",
141141
" data=ds.p,\n",
142142
" grid=grid,\n",
143-
" interp_method=UXPiecewiseConstantFace,\n",
143+
" interp_method=UxPiecewiseConstantFace,\n",
144144
")"
145145
]
146146
},

‎src/parcels/_core/fieldset.py‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -18,7 +18,7 @@
1818
from parcels._core.xgrid import _DEFAULT_XGCM_KWARGS, XGrid
1919
from parcels._logger import logger
2020
from parcels._typing import Mesh
21-
from parcels.interpolators import UXPiecewiseConstantFace, UXPiecewiseLinearNode, XConstantField, XLinear
21+
from parcels.interpolators import UxPiecewiseConstantFace, UxPiecewiseLinearNode, XConstantField, XLinear
2222

2323
if TYPE_CHECKING:
2424
from parcels._core.basegrid import BaseGrid
@@ -459,9 +459,9 @@ def _select_uxinterpolator(da: ux.UxDataArray):
459459
"""Selects the appropriate uxarray interpolator for a given uxarray UxDataArray"""
460460
supported_uxinterp_mapping = {
461461
# (nz1,n_face): face-center laterally, layer centers vertically — piecewise constant
462-
"nz1,n_face": UXPiecewiseConstantFace,
462+
"nz1,n_face": UxPiecewiseConstantFace,
463463
# (nz,n_node): node/corner laterally, layer interfaces vertically — barycentric lateral & linear vertical
464-
"nz,n_node": UXPiecewiseLinearNode,
464+
"nz,n_node": UxPiecewiseLinearNode,
465465
}
466466
# Extract only spatial dimensions, neglecting time
467467
da_spatial_dims = tuple(d for d in da.dims if d not in ("time",))

0 commit comments

Comments
 (0)