From 70c848c50cdcbbe3352e4a59fe3917503c8a965a Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 6 May 2026 11:43:31 -0400 Subject: [PATCH] Add unstructured uniform translation and solid body rotation exactness tests --- .../_datasets/unstructured/generated.py | 203 ++++++++++++++++++ tests/uxvalidation.py | 133 ++++++++++++ 2 files changed, 336 insertions(+) create mode 100644 tests/uxvalidation.py diff --git a/src/parcels/_datasets/unstructured/generated.py b/src/parcels/_datasets/unstructured/generated.py index 69b425ea23..6e789095a3 100644 --- a/src/parcels/_datasets/unstructured/generated.py +++ b/src/parcels/_datasets/unstructured/generated.py @@ -139,3 +139,206 @@ def simple_small_delaunay(nx=10, ny=10): ) return ux.UxDataset({"U": u, "V": v, "W": w, "p": p, "T_face": tface, "T_node": tnode}, uxgrid=uxgrid) + + +def _build_delaunay_grid(nx, lon_range, lat_range): + """Build a Delaunay-triangulated UxGrid from a regular nx-by-nx point lattice.""" + lon, lat = np.meshgrid( + np.linspace(lon_range[0], lon_range[1], nx, dtype=np.float64), + np.linspace(lat_range[0], lat_range[1], nx, dtype=np.float64), + ) + lon_flat = lon.ravel() + lat_flat = lat.ravel() + mask = ( + np.isclose(lon_flat, lon_range[0]) + | np.isclose(lon_flat, lon_range[1]) + | np.isclose(lat_flat, lat_range[0]) + | np.isclose(lat_flat, lat_range[1]) + ) + boundary_points = np.flatnonzero(mask) + uxgrid = ux.Grid.from_points( + (lon_flat, lat_flat), + method="regional_delaunay", + boundary_points=boundary_points, + ) + uxgrid.attrs["Conventions"] = "UGRID-1.0" + return uxgrid + + +def _wrap_uvw_dataset(uxgrid, U, V, W, zc, zf, time, uv_dim, uv_location, description): + """Wrap (U, V, W) numpy arrays into a UxDataset following Parcels' UGRID conventions. + + U, V are placed on ``uv_dim`` (``"n_face"`` or ``"n_node"``) at vertical centres ``zc``. + W is always on ``n_node`` at vertical interfaces ``zf``. + """ + u = ux.UxDataArray( + data=U, + name="U", + uxgrid=uxgrid, + dims=["time", "zc", uv_dim], + coords=dict(time=(["time"], time), zc=(["zc"], zc)), + attrs=dict( + description=f"zonal velocity ({description})", + units="degrees/s", + location=uv_location, + mesh="delaunay", + Conventions="UGRID-1.0", + ), + ) + v = ux.UxDataArray( + data=V, + name="V", + uxgrid=uxgrid, + dims=["time", "zc", uv_dim], + coords=dict(time=(["time"], time), zc=(["zc"], zc)), + attrs=dict( + description=f"meridional velocity ({description})", + units="degrees/s", + location=uv_location, + mesh="delaunay", + Conventions="UGRID-1.0", + ), + ) + w = ux.UxDataArray( + data=W, + name="W", + uxgrid=uxgrid, + dims=["time", "zf", "n_node"], + coords=dict(time=(["time"], time), zf=(["zf"], zf)), + attrs=dict( + description=f"vertical velocity ({description})", + units="degrees/s", + location="node", + mesh="delaunay", + Conventions="UGRID-1.0", + ), + ) + return ux.UxDataset({"U": u, "V": v, "W": w}, uxgrid=uxgrid) + + +def uniform_translation_face_centered(nx=20, U0=0.001, V0=0.0005): + """T1-1 uniform translation, face-centered (U, V on ``n_face``). + + Verification field `u = U_0`,`v = V_0` on a Delaunay triangulation of a + regular ``nx``-by-``nx`` point lattice over ``[0, 10] x [0, 10]``. Single vertical layer + spans ``z in [0, 1]``. ``W`` is identically zero on the corner nodes. + """ + uxgrid = _build_delaunay_grid(nx, (0.0, 10.0), (0.0, 10.0)) + + zf = np.array([0.0, 1.0], dtype=np.float64) + zc = np.array([0.5], dtype=np.float64) + time = xr.date_range("2000-01-01", periods=1, freq="2h") + + U = np.full((1, zc.size, uxgrid.n_face), U0, dtype=np.float64) + V = np.full((1, zc.size, uxgrid.n_face), V0, dtype=np.float64) + W = np.zeros((1, zf.size, uxgrid.n_node), dtype=np.float64) + + return _wrap_uvw_dataset(uxgrid, U, V, W, zc, zf, time, "n_face", "face", "uniform translation") + + +def uniform_translation_node_centered(nx=20, U0=0.001, V0=0.0005): + """T1-1 uniform translation, node-centered (U, V on ``n_node``). + + Same field as :func:`uniform_translation_face_centered` but with horizontal velocity + sampled at corner nodes for use with barycentric (linear) interpolators. + """ + uxgrid = _build_delaunay_grid(nx, (0.0, 10.0), (0.0, 10.0)) + + zf = np.array([0.0, 1.0], dtype=np.float64) + zc = np.array([0.5], dtype=np.float64) + time = xr.date_range("2000-01-01", periods=1, freq="2h") + + U = np.full((1, zc.size, uxgrid.n_node), U0, dtype=np.float64) + V = np.full((1, zc.size, uxgrid.n_node), V0, dtype=np.float64) + W = np.zeros((1, zf.size, uxgrid.n_node), dtype=np.float64) + + return _wrap_uvw_dataset(uxgrid, U, V, W, zc, zf, time, "n_node", "node", "uniform translation") + + +def solid_body_rotation_face_centered(nx=40, omega=2.0 * math.pi / 3600.0): + """T1-2 2D solid-body rotation, face-centered (U, V on ``n_face``). + + Verification field :math:`u = -\\Omega y`, :math:`v = \\Omega x` on a Delaunay + triangulation of a regular ``nx``-by-``nx`` point lattice over ``[-5, 5] x [-5, 5]`` + centred on the rotation axis. Single vertical layer; ``W = 0``. + """ + uxgrid = _build_delaunay_grid(nx, (-5.0, 5.0), (-5.0, 5.0)) + + zf = np.array([0.0, 1.0], dtype=np.float64) + zc = np.array([0.5], dtype=np.float64) + time = xr.date_range("2000-01-01", periods=1, freq="2h") + + u_vals = -omega * uxgrid.face_lat.values + v_vals = omega * uxgrid.face_lon.values + U = np.broadcast_to(u_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_face)).copy() + V = np.broadcast_to(v_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_face)).copy() + W = np.zeros((1, zf.size, uxgrid.n_node), dtype=np.float64) + + return _wrap_uvw_dataset(uxgrid, U, V, W, zc, zf, time, "n_face", "face", "solid-body rotation") + + +def solid_body_rotation_node_centered(nx=40, omega=2.0 * math.pi / 3600.0): + """T1-2 2D solid-body rotation, node-centered (U, V on ``n_node``). + + Same field as :func:`solid_body_rotation_face_centered` but sampled at corner + nodes. Because the field is linear in space, barycentric interpolation reproduces it + exactly and any error in particle trajectories is attributable to the time integrator. + """ + uxgrid = _build_delaunay_grid(nx, (-5.0, 5.0), (-5.0, 5.0)) + + zf = np.array([0.0, 1.0], dtype=np.float64) + zc = np.array([0.5], dtype=np.float64) + time = xr.date_range("2000-01-01", periods=1, freq="2h") + + u_vals = -omega * uxgrid.node_lat.values + v_vals = omega * uxgrid.node_lon.values + U = np.broadcast_to(u_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_node)).copy() + V = np.broadcast_to(v_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_node)).copy() + W = np.zeros((1, zf.size, uxgrid.n_node), dtype=np.float64) + + return _wrap_uvw_dataset(uxgrid, U, V, W, zc, zf, time, "n_node", "node", "solid-body rotation") + + +def solid_body_rotation_3d_face_centered(nx=40, nz=10, omega=2.0 * math.pi / 3600.0, W0=0.005): + """T1-3 3D helical motion, face-centered horizontal velocity. + + Verification field :math:`u = -\\Omega y`, :math:`v = \\Omega x`, :math:`w = W_0` on a + Delaunay triangulation of a regular ``nx``-by-``nx`` point lattice over + ``[-5, 5] x [-5, 5]`` with ``nz`` vertical layers spanning ``[0, 100]``. ``U`` and ``V`` + are sampled at face centres; ``W`` is sampled at corner nodes on the layer interfaces. + """ + uxgrid = _build_delaunay_grid(nx, (-5.0, 5.0), (-5.0, 5.0)) + + zf = np.linspace(0.0, 100.0, nz + 1, dtype=np.float64) + zc = 0.5 * (zf[:-1] + zf[1:]) + time = xr.date_range("2000-01-01", periods=1, freq="2h") + + u_vals = -omega * uxgrid.face_lat.values + v_vals = omega * uxgrid.face_lon.values + U = np.broadcast_to(u_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_face)).copy() + V = np.broadcast_to(v_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_face)).copy() + W = np.full((1, zf.size, uxgrid.n_node), W0, dtype=np.float64) + + return _wrap_uvw_dataset(uxgrid, U, V, W, zc, zf, time, "n_face", "face", "3D solid-body rotation") + + +def solid_body_rotation_3d_node_centered(nx=40, nz=10, omega=2.0 * math.pi / 3600.0, W0=0.005): + """T1-3 3D helical motion, node-centered (U, V, W on nodes). + + Same field as :func:`solid_body_rotation_3d_face_centered` with horizontal velocity + sampled at corner nodes. The horizontal field is linear and the vertical field is + constant, so barycentric horizontal + linear vertical interpolation are both exact. + """ + uxgrid = _build_delaunay_grid(nx, (-5.0, 5.0), (-5.0, 5.0)) + + zf = np.linspace(0.0, 100.0, nz + 1, dtype=np.float64) + zc = 0.5 * (zf[:-1] + zf[1:]) + time = xr.date_range("2000-01-01", periods=1, freq="2h") + + u_vals = -omega * uxgrid.node_lat.values + v_vals = omega * uxgrid.node_lon.values + U = np.broadcast_to(u_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_node)).copy() + V = np.broadcast_to(v_vals[np.newaxis, np.newaxis, :], (1, zc.size, uxgrid.n_node)).copy() + W = np.full((1, zf.size, uxgrid.n_node), W0, dtype=np.float64) + + return _wrap_uvw_dataset(uxgrid, U, V, W, zc, zf, time, "n_node", "node", "3D solid-body rotation") diff --git a/tests/uxvalidation.py b/tests/uxvalidation.py new file mode 100644 index 0000000000..b16ee019fb --- /dev/null +++ b/tests/uxvalidation.py @@ -0,0 +1,133 @@ +"""Validation tests for the analytical unstructured datasets in +``parcels._datasets.unstructured.generated``. + +Each test runs a particle through a generated field and asserts the result against the +closed-form trajectory at the strongest accuracy the configuration permits: + +* **T1-1** (uniform translation): the field is constant in space and time, so every + interpolator and every Runge-Kutta stage evaluates the same value. All four + combinations of face/node placement and EE/RK4 must hit machine precision. +* **T1-2** (2D solid-body rotation): the field is linear in space, so only the + node-centered grid (barycentric interpolation) reproduces it exactly. Trajectory + error is then attributable purely to the time integrator. RK4 on the node-centered + field must return a particle to its starting point after one full orbit. +* **T1-3** (3D helix): horizontal field as T1-2 plus a constant vertical velocity. + The vertical ODE has constant RHS so depth advection is exact for both face- and + node-centered datasets under any 3D integrator. +""" + +import math + +import numpy as np +import pytest + +import parcels +from parcels._datasets.unstructured.generated import ( + solid_body_rotation_3d_face_centered, + solid_body_rotation_3d_node_centered, + solid_body_rotation_face_centered, + solid_body_rotation_node_centered, + uniform_translation_face_centered, + uniform_translation_node_centered, +) +from parcels.kernels import ( + AdvectionEE, + AdvectionRK4, + AdvectionRK4_3D, +) + +# Uniform translation parameters +T1_1_U0 = 0.001 +T1_1_V0 = 0.0005 +T1_1_RUNTIME = np.timedelta64(3600, "s") +T1_1_DT = np.timedelta64(300, "s") +T1_1_LON0 = 3.0 +T1_1_LAT0 = 3.0 +T1_1_TOL = 1e-5 # particle positions are float32; ~6 decimal digits of precision + +# Solid body rotation - 2D parameters +T1_2_OMEGA = 2.0 * math.pi / 3600.0 +T1_2_PERIOD = 2.0 * math.pi / T1_2_OMEGA +T1_2_RUNTIME = np.timedelta64(int(T1_2_PERIOD), "s") +T1_2_DT = np.timedelta64(15, "s") +T1_2_LON0 = 2.0 +T1_2_LAT0 = 0.0 + +# Solid body rotation - 3D parameters +T1_3_W0 = 0.005 +T1_3_Z0 = 50.0 +T1_3_DT = np.timedelta64(15, "s") +T1_3_RUNTIME = T1_2_RUNTIME + + +@pytest.mark.parametrize( + "dataset_fn", + [uniform_translation_face_centered, uniform_translation_node_centered], + ids=["face-centered", "node-centered"], +) +@pytest.mark.parametrize("integrator", [AdvectionEE, AdvectionRK4], ids=["EE", "RK4"]) +def test_uniform_translation_exact(dataset_fn, integrator): + ds = dataset_fn(nx=20, U0=T1_1_U0, V0=T1_1_V0) + fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") + + pset = parcels.ParticleSet(fieldset, lon=[T1_1_LON0], lat=[T1_1_LAT0]) + pset.execute(integrator, runtime=T1_1_RUNTIME, dt=T1_1_DT, verbose_progress=False) + + t = T1_1_RUNTIME / np.timedelta64(1, "s") + np.testing.assert_allclose(pset.lon, T1_1_LON0 + T1_1_U0 * t, atol=T1_1_TOL) + np.testing.assert_allclose(pset.lat, T1_1_LAT0 + T1_1_V0 * t, atol=T1_1_TOL) + + +def test_solid_body_rotation_node_centered_rk4_returns_to_start(): + ds = solid_body_rotation_node_centered(nx=40, omega=T1_2_OMEGA) + fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") + + pset = parcels.ParticleSet(fieldset, lon=[T1_2_LON0], lat=[T1_2_LAT0]) + pset.execute(AdvectionRK4, runtime=T1_2_RUNTIME, dt=T1_2_DT, verbose_progress=False) + + np.testing.assert_allclose(pset.lon, T1_2_LON0, atol=1e-6) + np.testing.assert_allclose(pset.lat, T1_2_LAT0, atol=1e-6) + + +@pytest.mark.parametrize("integrator", [AdvectionEE, AdvectionRK4], ids=["EE", "RK4"]) +def test_solid_body_rotation_face_centered_runs_bounded(integrator): + # Piecewise-constant interpolation has spatial truncation error on a linear field, + # so we only assert the trajectory stays inside the [-5, 5] mesh. + ds = solid_body_rotation_face_centered(nx=40, omega=T1_2_OMEGA) + fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") + + pset = parcels.ParticleSet(fieldset, lon=[T1_2_LON0], lat=[T1_2_LAT0]) + pset.execute(integrator, runtime=T1_2_RUNTIME, dt=T1_2_DT, verbose_progress=False) + assert abs(float(pset.lon[0])) < 5.0 + assert abs(float(pset.lat[0])) < 5.0 + + +def test_helix_node_centered_rk4_3d_returns_to_start_with_exact_depth(): + ds = solid_body_rotation_3d_node_centered(nx=40, nz=10, omega=T1_2_OMEGA, W0=T1_3_W0) + fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") + + pset = parcels.ParticleSet(fieldset, lon=[T1_2_LON0], lat=[T1_2_LAT0], z=[T1_3_Z0]) + pset.execute(AdvectionRK4_3D, runtime=T1_3_RUNTIME, dt=T1_3_DT, verbose_progress=False) + + t = T1_3_RUNTIME / np.timedelta64(1, "s") + np.testing.assert_allclose(pset.lon, T1_2_LON0, atol=1e-6) + np.testing.assert_allclose(pset.lat, T1_2_LAT0, atol=1e-6) + np.testing.assert_allclose(pset.z, T1_3_Z0 + T1_3_W0 * t, atol=1e-4) + + +@pytest.mark.parametrize( + "dataset_fn", + [solid_body_rotation_3d_face_centered, solid_body_rotation_3d_node_centered], + ids=["face-centered", "node-centered"], +) +def test_helix_constant_vertical_velocity_exact_depth(dataset_fn): + # Constant w implies linear-in-z interpolation and constant-RHS depth ODE are both exact, + # independent of horizontal grid placement. + ds = dataset_fn(nx=40, nz=10, omega=T1_2_OMEGA, W0=T1_3_W0) + fieldset = parcels.FieldSet.from_ugrid_conventions(ds, mesh="flat") + + pset = parcels.ParticleSet(fieldset, lon=[T1_2_LON0], lat=[T1_2_LAT0], z=[T1_3_Z0]) + pset.execute(AdvectionRK4_3D, runtime=T1_3_RUNTIME, dt=T1_3_DT, verbose_progress=False) + + t = T1_3_RUNTIME / np.timedelta64(1, "s") + np.testing.assert_allclose(pset.z, T1_3_Z0 + T1_3_W0 * t, atol=1e-4)