Skip to content

Commit 182688a

Browse files
fluidnumericsJoefluidnumericsJoe
committed
Draft unstructured interpolation tutorial (wip)
1 parent 54fbc40 commit 182688a

2 files changed

Lines changed: 261 additions & 2 deletions

File tree

‎docs/user_guide/examples/tutorial_interpolation.ipynb‎

Lines changed: 120 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -251,6 +251,124 @@
251251
"- UXPiecewiseLinearNode"
252252
]
253253
},
254+
{
255+
"cell_type": "code",
256+
"execution_count": null,
257+
"metadata": {},
258+
"outputs": [],
259+
"source": [
260+
"from parcels._datasets.unstructured.generated import simple_small_delaunay\n",
261+
"\n",
262+
"ds = simple_small_delaunay(nx=4, ny=4)\n",
263+
"ds"
264+
]
265+
},
266+
{
267+
"cell_type": "code",
268+
"execution_count": null,
269+
"metadata": {},
270+
"outputs": [],
271+
"source": [
272+
"import parcels\n",
273+
"\n",
274+
"grid = parcels.UxGrid(ds.uxgrid, ds.nz, mesh=\"flat\")\n",
275+
"U = parcels.Field(\n",
276+
" \"U\", ds[\"U\"], grid, interp_method=parcels.interpolators.UXPiecewiseConstantFace\n",
277+
")\n",
278+
"V = parcels.Field(\n",
279+
" \"V\", ds[\"V\"], grid, interp_method=parcels.interpolators.UXPiecewiseConstantFace\n",
280+
")\n",
281+
"UV = parcels.VectorField(\"UV\", U, V)\n",
282+
"Tnode = parcels.Field(\n",
283+
" \"T_node\",\n",
284+
" ds[\"T_node\"],\n",
285+
" grid,\n",
286+
" interp_method=parcels.interpolators.UXPiecewiseLinearNode,\n",
287+
")\n",
288+
"Tface = parcels.Field(\n",
289+
" \"T_face\",\n",
290+
" ds[\"T_face\"],\n",
291+
" grid,\n",
292+
" interp_method=parcels.interpolators.UXPiecewiseConstantFace,\n",
293+
")\n",
294+
"fieldset = parcels.FieldSet([U, V, UV, Tnode, Tface])"
295+
]
296+
},
297+
{
298+
"cell_type": "code",
299+
"execution_count": null,
300+
"metadata": {},
301+
"outputs": [],
302+
"source": [
303+
"import matplotlib.pyplot as plt\n",
304+
"import matplotlib.tri as mtri\n",
305+
"import numpy as np\n",
306+
"\n",
307+
"fig, ax = plt.subplots(1, 2, figsize=(10, 5), constrained_layout=True)\n",
308+
"\n",
309+
"x = ds.uxgrid.node_lon.values\n",
310+
"y = ds.uxgrid.node_lat.values\n",
311+
"tri = ds.uxgrid.face_node_connectivity.values\n",
312+
"triang = mtri.Triangulation(x, y, triangles=tri)\n",
313+
"\n",
314+
"\n",
315+
"# Shade the triangle faces using smooth shading between the node registered data\n",
316+
"T_node = np.squeeze(ds[\"T_node\"][0, 0, :].values)\n",
317+
"tpc = ax[0].tripcolor(\n",
318+
" triang,\n",
319+
" T_node,\n",
320+
" shading=\"gouraud\", # smooth shading\n",
321+
" cmap=\"viridis\",\n",
322+
" vmin=-1.0,\n",
323+
" vmax=1.0,\n",
324+
")\n",
325+
"# Plot the element edges\n",
326+
"ax[0].triplot(triang, color=\"black\", linewidth=0.5)\n",
327+
"\n",
328+
"# Plot the node registered data\n",
329+
"sc1 = ax[0].scatter(\n",
330+
" x, y, c=T_node, cmap=\"viridis\", s=150, edgecolors=\"black\", vmin=-1.0, vmax=1.0\n",
331+
")\n",
332+
"\n",
333+
"cbar = fig.colorbar(sc1, ax=ax[0])\n",
334+
"cbar.set_label(\"T\")\n",
335+
"\n",
336+
"ax[0].set_aspect(\"equal\", \"box\")\n",
337+
"ax[0].set_xlabel(\"x\")\n",
338+
"ax[0].set_ylabel(\"y\")\n",
339+
"ax[0].set_title(\"Node Registered Data\")\n",
340+
"\n",
341+
"# # Plot the face registered data\n",
342+
"T_face = np.squeeze(ds[\"T_face\"][0, 0, :].values)\n",
343+
"tpc = ax[1].tripcolor(\n",
344+
" triang,\n",
345+
" T_face,\n",
346+
" shading=\"flat\",\n",
347+
" cmap=\"viridis\",\n",
348+
" edgecolors=\"k\",\n",
349+
" linewidth=0.5,\n",
350+
" vmin=-1.0,\n",
351+
" vmax=1.0,\n",
352+
")\n",
353+
"# Plot the element edges\n",
354+
"ax[1].triplot(triang, color=\"black\", linewidth=0.5)\n",
355+
"\n",
356+
"# Plot the node registered data\n",
357+
"xf = ds.uxgrid.face_lon.values\n",
358+
"yf = ds.uxgrid.face_lat.values\n",
359+
"sc2 = ax[1].scatter(\n",
360+
" xf, yf, c=T_face, cmap=\"viridis\", s=150, edgecolors=\"black\", vmin=-1.0, vmax=1.0\n",
361+
")\n",
362+
"\n",
363+
"cbar = fig.colorbar(sc2, ax=ax[1])\n",
364+
"cbar.set_label(\"T\")\n",
365+
"\n",
366+
"ax[1].set_aspect(\"equal\", \"box\")\n",
367+
"ax[1].set_xlabel(\"x\")\n",
368+
"ax[1].set_ylabel(\"y\")\n",
369+
"ax[1].set_title(\"Face Registered Data\")"
370+
]
371+
},
254372
{
255373
"cell_type": "markdown",
256374
"metadata": {},
@@ -272,7 +390,7 @@
272390
],
273391
"metadata": {
274392
"kernelspec": {
275-
"display_name": "test-notebooks",
393+
"display_name": "Python 3 (ipykernel)",
276394
"language": "python",
277395
"name": "python3"
278396
},
@@ -286,7 +404,7 @@
286404
"name": "python",
287405
"nbconvert_exporter": "python",
288406
"pygments_lexer": "ipython3",
289-
"version": "3.11.0"
407+
"version": "3.12.11"
290408
}
291409
},
292410
"nbformat": 4,

‎src/parcels/_datasets/unstructured/generated.py‎

Lines changed: 141 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,141 @@
1+
import math
2+
3+
import numpy as np
4+
import uxarray as ux
5+
import xarray as xr
6+
7+
T = 2
8+
vmax = 1.0
9+
delta = 0.1
10+
TIME = xr.date_range("2000", "2001", T)
11+
12+
13+
def simple_small_delaunay(nx=10, ny=10):
14+
"""
15+
Data on a small Delaunay grid. The naming convention of the dataset and grid is consistent with what is
16+
provided by UXArray when reading in FESOM2 datasets.
17+
"""
18+
lon, lat = np.meshgrid(np.linspace(0, 1.0, nx, dtype=np.float32), np.linspace(0, 1.0, ny, dtype=np.float32))
19+
lon_flat = lon.ravel()
20+
lat_flat = lat.ravel()
21+
zf = np.linspace(0.0, 1000.0, 2, endpoint=True, dtype=np.float32) # Vertical element faces
22+
zc = 0.5 * (zf[:-1] + zf[1:]) # Vertical element centers
23+
nz = zf.size
24+
nz1 = zc.size
25+
26+
# mask any point on one of the boundaries
27+
mask = np.isclose(lon_flat, 0.0) | np.isclose(lon_flat, 1.0) | np.isclose(lat_flat, 0.0) | np.isclose(lat_flat, 1.0)
28+
29+
boundary_points = np.flatnonzero(mask)
30+
31+
uxgrid = ux.Grid.from_points(
32+
(lon_flat, lat_flat),
33+
method="regional_delaunay",
34+
boundary_points=boundary_points,
35+
)
36+
uxgrid.attrs["Conventions"] = "UGRID-1.0"
37+
38+
# Define arrays U (zonal), V (meridional) and P (sea surface height)
39+
U = np.zeros((1, nz1, uxgrid.n_face), dtype=np.float64)
40+
V = np.zeros((1, nz1, uxgrid.n_face), dtype=np.float64)
41+
W = np.zeros((1, nz, uxgrid.n_node), dtype=np.float64)
42+
P = np.zeros((1, nz1, uxgrid.n_face), dtype=np.float64)
43+
Tface = np.zeros((1, nz1, uxgrid.n_face), dtype=np.float64)
44+
45+
for i, (x, y) in enumerate(zip(uxgrid.face_lon, uxgrid.face_lat, strict=False)):
46+
P[0, 0, i] = -vmax * delta * (1 - x) * (math.exp(-x / delta) - 1) * np.sin(math.pi * y)
47+
U[0, 0, i] = -vmax * (1 - math.exp(-x / delta) - x) * np.cos(math.pi * y)
48+
V[0, 0, i] = vmax * ((2.0 - x) * math.exp(-x / delta) - 1) * np.sin(math.pi * y)
49+
Tface[0, 0, i] = np.sin(math.pi * y) * np.cos(math.pi * x)
50+
51+
Tnode = np.zeros((1, nz1, uxgrid.n_node), dtype=np.float64)
52+
for i, (x, y) in enumerate(zip(uxgrid.node_lon, uxgrid.node_lat, strict=False)):
53+
Tnode[0, 0, i] = np.sin(math.pi * y) * np.cos(math.pi * x)
54+
55+
u = ux.UxDataArray(
56+
data=U,
57+
name="U",
58+
uxgrid=uxgrid,
59+
dims=["time", "nz1", "n_face"],
60+
coords=dict(
61+
time=(["time"], [TIME[0]]),
62+
nz1=(["nz1"], zc),
63+
),
64+
attrs=dict(
65+
description="zonal velocity", units="m/s", location="face", mesh="delaunay", Conventions="UGRID-1.0"
66+
),
67+
)
68+
v = ux.UxDataArray(
69+
data=V,
70+
name="V",
71+
uxgrid=uxgrid,
72+
dims=["time", "nz1", "n_face"],
73+
coords=dict(
74+
time=(["time"], [TIME[0]]),
75+
nz1=(["nz1"], zc),
76+
),
77+
attrs=dict(
78+
description="meridional velocity", units="m/s", location="face", mesh="delaunay", Conventions="UGRID-1.0"
79+
),
80+
)
81+
w = ux.UxDataArray(
82+
data=W,
83+
name="W",
84+
uxgrid=uxgrid,
85+
dims=["time", "nz", "n_node"],
86+
coords=dict(
87+
time=(["time"], [TIME[0]]),
88+
nz=(["nz"], zf),
89+
),
90+
attrs=dict(
91+
description="meridional velocity", units="m/s", location="node", mesh="delaunay", Conventions="UGRID-1.0"
92+
),
93+
)
94+
p = ux.UxDataArray(
95+
data=P,
96+
name="p",
97+
uxgrid=uxgrid,
98+
dims=["time", "nz1", "n_face"],
99+
coords=dict(
100+
time=(["time"], [TIME[0]]),
101+
nz1=(["nz1"], zc),
102+
),
103+
attrs=dict(description="pressure", units="N/m^2", location="face", mesh="delaunay", Conventions="UGRID-1.0"),
104+
)
105+
106+
tface = ux.UxDataArray(
107+
data=Tface,
108+
name="T_face",
109+
uxgrid=uxgrid,
110+
dims=["time", "nz1", "n_face"],
111+
coords=dict(
112+
time=(["time"], [TIME[0]]),
113+
nz1=(["nz1"], zc),
114+
),
115+
attrs=dict(
116+
description="Tracer field sampled on face centers",
117+
units="None",
118+
location="face",
119+
mesh="delaunay",
120+
Conventions="UGRID-1.0",
121+
),
122+
)
123+
tnode = ux.UxDataArray(
124+
data=Tnode,
125+
name="T_node",
126+
uxgrid=uxgrid,
127+
dims=["time", "nz1", "n_node"],
128+
coords=dict(
129+
time=(["time"], [TIME[0]]),
130+
nz1=(["nz1"], zc),
131+
),
132+
attrs=dict(
133+
description="Tracer field sampled on face vertices",
134+
units="None",
135+
location="node",
136+
mesh="delaunay",
137+
Conventions="UGRID-1.0",
138+
),
139+
)
140+
141+
return ux.UxDataset({"U": u, "V": v, "W": w, "p": p, "T_face": tface, "T_node": tnode}, uxgrid=uxgrid)

0 commit comments

Comments
 (0)