From 6a3f4418d2a1869c009da7f01fcc814f84265c8c Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 9 Jun 2025 16:41:15 -0400 Subject: [PATCH 01/12] Fix call to grid.search to remove reference to field --- parcels/field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/parcels/field.py b/parcels/field.py index 81116dee5a..34d1a842aa 100644 --- a/parcels/field.py +++ b/parcels/field.py @@ -327,7 +327,7 @@ def eval(self, time: datetime, z, y, x, particle=None, applyConversion=True): try: tau, ti = _search_time_index(self, time) - bcoords, _ei = self.grid.search(self, z, y, x, ei=_ei) + bcoords, _ei = self.grid.search(z, y, x, ei=_ei) value = self._interp_method(self, ti, _ei, bcoords, tau, time, z, y, x) if np.isnan(value): From 851e131ca5a229b981a58f3c0a0f0d36d7f7f6df Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 9 Jun 2025 17:05:37 -0400 Subject: [PATCH 02/12] [#2031] Add vertical coordinate search --- parcels/uxgrid.py | 38 ++++++++++++++++++++++++++++++-------- 1 file changed, 30 insertions(+), 8 deletions(-) diff --git a/parcels/uxgrid.py b/parcels/uxgrid.py index abdca45931..023b2eac99 100644 --- a/parcels/uxgrid.py +++ b/parcels/uxgrid.py @@ -15,8 +15,19 @@ class UxGrid(BaseGrid): for interpolation on unstructured grids. """ - def __init__(self, grid: ux.grid.Grid) -> UxGrid: + def __init__(self, grid: ux.grid.Grid, z: np.ndarray) -> UxGrid: + """ + Initializes the UxGrid with a uxarray grid and vertical coordinate array. + + Parameters + ---------- + grid : ux.grid.Grid + The uxarray grid object containing the unstructured grid data. + z : np.ndarray + A 1D array of vertical coordinates (depths) corresponding to the vertical position of layer faces of the grid + """ self.uxgrid = grid + self.z = z def search( self, z: float, y: float, x: float, ei: int | None = None, search2D: bool = False @@ -24,25 +35,36 @@ def search( tol = 1e-10 def try_face(fid): - # TODO : Vertical search is not implemented yet, so we assume z is not used. bcoords, err = self.uxgrid._get_barycentric_coordinates(y, x, fid) if (bcoords >= 0).all() and (bcoords <= 1).all() and err < tol: - return bcoords, self.ravel_index(0, fid) # Z and time indices are 0 for now + return bcoords, fid return None, None + def find_vertical_index() -> int: + if len(self.z) == 1: + return 0 + zi = np.searchsorted(self.z, z, side="right") - 1 # Search assumes that z is positive and increasing with i + if zi < 0 or zi >= len(self.z) - 1: + raise FieldOutOfBoundError(z, y, x) + return zi + + if not search2D: + zi = find_vertical_index() # Find the vertical cell center nearest to z + else: + zi = 0 + if ei is not None: zi, fi = self.unravel_index(ei) - bcoords, ei_new = try_face(fi) + bcoords, fi_new = try_face(fi) if bcoords is not None: - return bcoords, ei_new - + return bcoords, self.ravel_index(zi, fi_new) # Try neighbors of current face for neighbor in self.uxgrid.face_face_connectivity[fi, :]: if neighbor == -1: continue - bcoords, ei_new = try_face(neighbor) + bcoords, fi_new = try_face(neighbor) if bcoords is not None: - return bcoords, ei_new + return bcoords, self.ravel_index(zi, fi_new) # Global fallback using spatial hash fi, bcoords = self.uxgrid.get_spatial_hash().query([[x, y]]) From dc3c071a8dc1a5ceb92dce73a48f9cd985e90256 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 9 Jun 2025 17:06:19 -0400 Subject: [PATCH 03/12] [#2031] Change call to unravel_index via field.grid --- parcels/application_kernels/interpolation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/parcels/application_kernels/interpolation.py b/parcels/application_kernels/interpolation.py index a1abf14064..2801f46faa 100644 --- a/parcels/application_kernels/interpolation.py +++ b/parcels/application_kernels/interpolation.py @@ -27,7 +27,7 @@ def UXPiecewiseConstantFace( face registered, such as u,v in FESOM. """ # TODO joe : handle vertical interpolation - zi, fi = field.unravel_index(ei) + zi, fi = field.grid.unravel_index(ei) return field.data[ti, zi, fi] @@ -48,6 +48,6 @@ def UXPiecewiseLinearNode( such as the vertical velocity w in FESOM. """ # TODO joe : handle vertical interpolation - zi, fi = field.unravel_index(ei) - node_ids = field.data.uxgrid.face_node_connectivity[fi, :] + zi, fi = field.grid.unravel_index(ei) + node_ids = field.grid.uxgrid.face_node_connectivity[fi, :] return np.dot(field.data[ti, zi, node_ids], bcoords) From 3f5afe7d239434e1ea63cd592c4a88023479510c Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Mon, 9 Jun 2025 17:07:20 -0400 Subject: [PATCH 04/12] [#2031] Update tests to pass vertical face locations to UxGrid constructor --- tests/v4/test_uxarray_fieldset.py | 23 +++++++++++++---------- 1 file changed, 13 insertions(+), 10 deletions(-) diff --git a/tests/v4/test_uxarray_fieldset.py b/tests/v4/test_uxarray_fieldset.py index 7b0e64b2a3..975c21b413 100644 --- a/tests/v4/test_uxarray_fieldset.py +++ b/tests/v4/test_uxarray_fieldset.py @@ -37,13 +37,13 @@ def uv_fesom_channel(ds_fesom_channel) -> VectorField: U=Field( name="U", data=ds_fesom_channel.U, - grid=UxGrid(ds_fesom_channel.uxgrid), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), interp_method=UXPiecewiseConstantFace, ), V=Field( name="V", data=ds_fesom_channel.V, - grid=UxGrid(ds_fesom_channel.uxgrid), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), interp_method=UXPiecewiseConstantFace, ), ) @@ -57,17 +57,20 @@ def uvw_fesom_channel(ds_fesom_channel) -> VectorField: U=Field( name="U", data=ds_fesom_channel.U, - grid=UxGrid(ds_fesom_channel.uxgrid), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), interp_method=UXPiecewiseConstantFace, ), V=Field( name="V", data=ds_fesom_channel.V, - grid=UxGrid(ds_fesom_channel.uxgrid), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), interp_method=UXPiecewiseConstantFace, ), W=Field( - name="W", data=ds_fesom_channel.W, grid=UxGrid(ds_fesom_channel.uxgrid), interp_method=UXPiecewiseLinearNode + name="W", + data=ds_fesom_channel.W, + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), + interp_method=UXPiecewiseLinearNode, ), ) return UVW @@ -123,14 +126,14 @@ def test_fesom2_square_delaunay_uniform_z_coordinate_eval(): ds = datasets_unstructured["fesom2_square_delaunay_uniform_z_coordinate"] UVW = VectorField( name="UVW", - U=Field(name="U", data=ds.U, grid=UxGrid(ds.uxgrid), interp_method=UXPiecewiseConstantFace), - V=Field(name="V", data=ds.V, grid=UxGrid(ds.uxgrid), interp_method=UXPiecewiseConstantFace), - W=Field(name="W", data=ds.W, grid=UxGrid(ds.uxgrid), interp_method=UXPiecewiseLinearNode), + U=Field(name="U", data=ds.U, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseConstantFace), + V=Field(name="V", data=ds.V, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseConstantFace), + W=Field(name="W", data=ds.W, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseLinearNode), ) - P = Field(name="p", data=ds.p, grid=UxGrid(ds.uxgrid), interp_method=UXPiecewiseConstantFace) + P = Field(name="p", data=ds.p, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseConstantFace) fieldset = FieldSet([UVW, P, UVW.U, UVW.V, UVW.W]) assert fieldset.U.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 assert fieldset.V.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 - # assert fieldset.W.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 0.0 + assert fieldset.W.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 0.0 assert fieldset.P.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 From bb37185b26db6b6b491327c321db8fb9f96b0361 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Tue, 10 Jun 2025 08:58:30 -0400 Subject: [PATCH 05/12] [#2031] Change z to be ux.UxDataArray; require 1d vertical coordinate --- parcels/uxgrid.py | 18 +++++++++++++----- 1 file changed, 13 insertions(+), 5 deletions(-) diff --git a/parcels/uxgrid.py b/parcels/uxgrid.py index 023b2eac99..780078b235 100644 --- a/parcels/uxgrid.py +++ b/parcels/uxgrid.py @@ -15,7 +15,7 @@ class UxGrid(BaseGrid): for interpolation on unstructured grids. """ - def __init__(self, grid: ux.grid.Grid, z: np.ndarray) -> UxGrid: + def __init__(self, grid: ux.grid.Grid, z: ux.UxDataArray) -> UxGrid: """ Initializes the UxGrid with a uxarray grid and vertical coordinate array. @@ -23,10 +23,16 @@ def __init__(self, grid: ux.grid.Grid, z: np.ndarray) -> UxGrid: ---------- grid : ux.grid.Grid The uxarray grid object containing the unstructured grid data. - z : np.ndarray + z : ux.UxDataArray A 1D array of vertical coordinates (depths) corresponding to the vertical position of layer faces of the grid + While uxarray allows nz to be spatially and temporally varying, the parcels.UxGrid class considers the case where + the vertical coordinate is constant in time and space. This implies flat bottom topography and no moving ALE vertical grid. """ self.uxgrid = grid + if not isinstance(z, ux.UxDataArray): + raise TypeError("z must be an instance of ux.UxDataArray") + if z.ndim != 1: + raise ValueError("z must be a 1D array of vertical coordinates") self.z = z def search( @@ -41,10 +47,12 @@ def try_face(fid): return None, None def find_vertical_index() -> int: - if len(self.z) == 1: + nz = self.z.shape[0] + if nz == 1: return 0 - zi = np.searchsorted(self.z, z, side="right") - 1 # Search assumes that z is positive and increasing with i - if zi < 0 or zi >= len(self.z) - 1: + zf = self.z.values + zi = np.searchsorted(zf, z, side="right") - 1 # Search assumes that z is positive and increasing with i + if zi < 0 or zi >= nz - 1: raise FieldOutOfBoundError(z, y, x) return zi From ddfa6bad78c5628f392adc22556c87a4ff98ac23 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Tue, 10 Jun 2025 09:00:46 -0400 Subject: [PATCH 06/12] [#2031] Update tests to use nz coord for z in tests --- tests/v4/test_uxarray_fieldset.py | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/tests/v4/test_uxarray_fieldset.py b/tests/v4/test_uxarray_fieldset.py index 975c21b413..c62df25a77 100644 --- a/tests/v4/test_uxarray_fieldset.py +++ b/tests/v4/test_uxarray_fieldset.py @@ -37,13 +37,13 @@ def uv_fesom_channel(ds_fesom_channel) -> VectorField: U=Field( name="U", data=ds_fesom_channel.U, - grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"]), interp_method=UXPiecewiseConstantFace, ), V=Field( name="V", data=ds_fesom_channel.V, - grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"]), interp_method=UXPiecewiseConstantFace, ), ) @@ -57,19 +57,19 @@ def uvw_fesom_channel(ds_fesom_channel) -> VectorField: U=Field( name="U", data=ds_fesom_channel.U, - grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"]), interp_method=UXPiecewiseConstantFace, ), V=Field( name="V", data=ds_fesom_channel.V, - grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"]), interp_method=UXPiecewiseConstantFace, ), W=Field( name="W", data=ds_fesom_channel.W, - grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.nz), + grid=UxGrid(ds_fesom_channel.uxgrid, z=ds_fesom_channel.coords["nz"]), interp_method=UXPiecewiseLinearNode, ), ) @@ -126,11 +126,11 @@ def test_fesom2_square_delaunay_uniform_z_coordinate_eval(): ds = datasets_unstructured["fesom2_square_delaunay_uniform_z_coordinate"] UVW = VectorField( name="UVW", - U=Field(name="U", data=ds.U, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseConstantFace), - V=Field(name="V", data=ds.V, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseConstantFace), - W=Field(name="W", data=ds.W, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseLinearNode), + U=Field(name="U", data=ds.U, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseConstantFace), + V=Field(name="V", data=ds.V, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseConstantFace), + W=Field(name="W", data=ds.W, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseLinearNode), ) - P = Field(name="p", data=ds.p, grid=UxGrid(ds.uxgrid, z=ds.nz), interp_method=UXPiecewiseConstantFace) + P = Field(name="p", data=ds.p, grid=UxGrid(ds.uxgrid, z=ds.coords["nz"]), interp_method=UXPiecewiseConstantFace) fieldset = FieldSet([UVW, P, UVW.U, UVW.V, UVW.W]) assert fieldset.U.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 From 90e73df2952fc6afe942c33d8e6c7461db68bcd1 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Tue, 10 Jun 2025 09:10:11 -0400 Subject: [PATCH 07/12] [#2031] Add `nz` coordinate and vertical velocity to stommel gyre example `nz` coordinate, which is the vertical face locations of the vertical grid are required for the `UxGrid` vertical search --- parcels/_datasets/unstructured/generic.py | 30 ++++++++++++++++++----- 1 file changed, 24 insertions(+), 6 deletions(-) diff --git a/parcels/_datasets/unstructured/generic.py b/parcels/_datasets/unstructured/generic.py index 91704adeb2..ac9542638d 100644 --- a/parcels/_datasets/unstructured/generic.py +++ b/parcels/_datasets/unstructured/generic.py @@ -24,6 +24,10 @@ def _stommel_gyre_delaunay(): lon, lat = np.meshgrid(np.linspace(0, 60.0, Nx, dtype=np.float32), np.linspace(0, 60.0, Nx, dtype=np.float32)) lon_flat = lon.ravel() lat_flat = lat.ravel() + zf = np.linspace(0.0, 1000.0, 2, endpoint=True, dtype=np.float32) # Vertical element faces + zc = 0.5 * (zf[:-1] + zf[1:]) # Vertical element centers + nz = zf.size + nz1 = zc.size # mask any point on one of the boundaries mask = ( @@ -40,9 +44,10 @@ def _stommel_gyre_delaunay(): uxgrid.attrs["Conventions"] = "UGRID-1.0" # Define arrays U (zonal), V (meridional) and P (sea surface height) - U = np.zeros((1, 1, lat.size), dtype=np.float64) - V = np.zeros((1, 1, lat.size), dtype=np.float64) - P = np.zeros((1, 1, lat.size), dtype=np.float64) + U = np.zeros((1, nz1, lat.size), dtype=np.float64) + V = np.zeros((1, nz1, lat.size), dtype=np.float64) + W = np.zeros((1, nz, lat.size), dtype=np.float64) + P = np.zeros((1, nz1, lat.size), dtype=np.float64) for i, (x, y) in enumerate(zip(lon_flat, lat_flat, strict=False)): xi = x / 60.0 @@ -72,7 +77,20 @@ def _stommel_gyre_delaunay(): dims=["time", "nz1", "n_node"], coords=dict( time=(["time"], [TIME[0]]), - nz1=(["nz1"], [0]), + nz1=(["nz1"], zc), + ), + attrs=dict( + description="meridional velocity", units="m/s", location="node", mesh="delaunay", Conventions="UGRID-1.0" + ), + ) + w = ux.UxDataArray( + data=W, + name="W", + uxgrid=uxgrid, + dims=["time", "nz", "n_node"], + coords=dict( + time=(["time"], [TIME[0]]), + nz=(["nz"], zf), ), attrs=dict( description="meridional velocity", units="m/s", location="node", mesh="delaunay", Conventions="UGRID-1.0" @@ -85,12 +103,12 @@ def _stommel_gyre_delaunay(): dims=["time", "nz1", "n_node"], coords=dict( time=(["time"], [TIME[0]]), - nz1=(["nz1"], [0]), + nz1=(["nz1"], zc), ), attrs=dict(description="pressure", units="N/m^2", location="node", mesh="delaunay", Conventions="UGRID-1.0"), ) - return ux.UxDataset({"U": u, "V": v, "p": p}, uxgrid=uxgrid) + return ux.UxDataset({"U": u, "V": v, "W": w, "p": p}, uxgrid=uxgrid) def _fesom2_square_delaunay_uniform_z_coordinate(): From b1f5f1cabec858bef6ac273fdc7e11018797a704 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 11 Jun 2025 09:20:28 -0400 Subject: [PATCH 08/12] [#2031] Add vertical search for layer index --- parcels/uxgrid.py | 31 ++++++++++++++++--------------- 1 file changed, 16 insertions(+), 15 deletions(-) diff --git a/parcels/uxgrid.py b/parcels/uxgrid.py index 780078b235..eccd51aa58 100644 --- a/parcels/uxgrid.py +++ b/parcels/uxgrid.py @@ -24,7 +24,7 @@ def __init__(self, grid: ux.grid.Grid, z: ux.UxDataArray) -> UxGrid: grid : ux.grid.Grid The uxarray grid object containing the unstructured grid data. z : ux.UxDataArray - A 1D array of vertical coordinates (depths) corresponding to the vertical position of layer faces of the grid + A 1D array of vertical coordinates (depths) associated with the layer interface heights (not the mid-layer depths). While uxarray allows nz to be spatially and temporally varying, the parcels.UxGrid class considers the case where the vertical coordinate is constant in time and space. This implies flat bottom topography and no moving ALE vertical grid. """ @@ -47,22 +47,23 @@ def try_face(fid): return None, None def find_vertical_index() -> int: - nz = self.z.shape[0] - if nz == 1: + if search2D: return 0 - zf = self.z.values - zi = np.searchsorted(zf, z, side="right") - 1 # Search assumes that z is positive and increasing with i - if zi < 0 or zi >= nz - 1: - raise FieldOutOfBoundError(z, y, x) - return zi - - if not search2D: - zi = find_vertical_index() # Find the vertical cell center nearest to z - else: - zi = 0 + else: + nz = self.z.shape[0] + if nz == 1: + return 0 + zf = self.z.values + # Return zi such that zf[zi] <= z < zf[zi+1] + zi = np.searchsorted(zf, z, side="right") - 1 # Search assumes that z is positive and increasing with i + if zi < 0 or zi >= nz - 1: + raise FieldOutOfBoundError(z, y, x) + return zi + + zi = find_vertical_index() # Find the vertical cell center nearest to z if ei is not None: - zi, fi = self.unravel_index(ei) + _, fi = self.unravel_index(ei) bcoords, fi_new = try_face(fi) if bcoords is not None: return bcoords, self.ravel_index(zi, fi_new) @@ -79,7 +80,7 @@ def find_vertical_index() -> int: if fi == -1: raise FieldOutOfBoundError(z, y, x) - return bcoords, self.ravel_index(zi, fi) + return bcoords[0], self.ravel_index(zi, fi[0]) def _get_barycentric_coordinates(self, y, x, fi): """Checks if a point is inside a given face id on a UxGrid.""" From 86f39be41dc1dc8760f7a50fed6d59613e142ad9 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 11 Jun 2025 10:24:15 -0400 Subject: [PATCH 09/12] [#2031] Add linear vertical interpolation for uxgrid This commit adds a test case to verify we can exactly interpolate fields that depend linearly on z. Additionally, this commit fixes the previously failing `test_fesom2_square_delaunay_uniform_z_coordinate_eval`. The fix for the apparent memory bug is to reference `field.values` --- parcels/application_kernels/interpolation.py | 24 +++++++++----- tests/v4/test_field.py | 35 ++++++++++++++++++-- tests/v4/test_uxarray_fieldset.py | 3 +- 3 files changed, 50 insertions(+), 12 deletions(-) diff --git a/parcels/application_kernels/interpolation.py b/parcels/application_kernels/interpolation.py index 2801f46faa..55d2abd5a7 100644 --- a/parcels/application_kernels/interpolation.py +++ b/parcels/application_kernels/interpolation.py @@ -26,9 +26,8 @@ def UXPiecewiseConstantFace( This interpolation method is appropriate for fields that are face registered, such as u,v in FESOM. """ - # TODO joe : handle vertical interpolation zi, fi = field.grid.unravel_index(ei) - return field.data[ti, zi, fi] + return field.values[ti, zi, fi] def UXPiecewiseLinearNode( @@ -43,11 +42,20 @@ def UXPiecewiseLinearNode( x: np.float32 | np.float64, ): """ - Piecewise linear interpolation kernel for node registered data. This - interpolation method is appropriate for fields that are node registered - such as the vertical velocity w in FESOM. + Piecewise linear interpolation kernel for node registered data located at vertical interface levels. + This interpolation method is appropriate for fields that are node registered such as the vertical + velocity W in FESOM2. Effectively, it applies barycentric interpolation in the lateral direction + and piecewise linear interpolation in the vertical direction. """ - # TODO joe : handle vertical interpolation - zi, fi = field.grid.unravel_index(ei) + k, fi = field.grid.unravel_index(ei) node_ids = field.grid.uxgrid.face_node_connectivity[fi, :] - return np.dot(field.data[ti, zi, node_ids], bcoords) + # The zi refers to the vertical layer index. The field in this routine are assumed to be defined at the vertical interface levels. + # For interface zi, the interface indices are [zi, zi+1], so we need to use the values at zi and zi+1. + # First, do barycentric interpolation in the lateral direction for each interface level + fzk = np.dot(field.values[ti, k, node_ids], bcoords) + fzkp1 = np.dot(field.values[ti, k + 1, node_ids], bcoords) + + # Then, do piecewise linear interpolation in the vertical direction + zk = field.grid.z.values[k] + zkp1 = field.grid.z.values[k + 1] + return (fzk * (zkp1 - z) + fzkp1 * (z - zk)) / (zkp1 - zk) # Linear interpolation in the vertical direction diff --git a/tests/v4/test_field.py b/tests/v4/test_field.py index 5dcad3dccd..12c619e6fc 100644 --- a/tests/v4/test_field.py +++ b/tests/v4/test_field.py @@ -3,7 +3,7 @@ import uxarray as ux import xarray as xr -from parcels import Field +from parcels import Field, UXPiecewiseConstantFace, UXPiecewiseLinearNode from parcels._datasets.structured.generic import T as T_structured from parcels._datasets.structured.generic import datasets as datasets_structured from parcels._datasets.unstructured.generic import datasets as datasets_unstructured @@ -38,7 +38,10 @@ def test_field_init_param_types(): pytest.param(ux.UxDataArray(), Grid(xr.Dataset()), id="uxdata-grid"), pytest.param( xr.DataArray(), - UxGrid(datasets_unstructured["stommel_gyre_delaunay"].uxgrid), + UxGrid( + datasets_unstructured["stommel_gyre_delaunay"].uxgrid, + z=datasets_unstructured["stommel_gyre_delaunay"].coords["nz"], + ), id="xarray-uxgrid", ), ], @@ -121,3 +124,31 @@ def test_field_interpolation_out_of_spatial_bounds(): ... def test_field_interpolation_out_of_time_bounds(): ... + + +def test_field_unstructured_z_linear(): + ds = datasets_unstructured["fesom2_square_delaunay_uniform_z_coordinate"] + + # Change the pressure values to be linearly dependent on the vertical coordinate + for k, z in enumerate(ds.coords["nz1"]): + ds["p"].values[:, k, :] = z + + # Change the vertical velocity values to be linearly dependent on the vertical coordinate + for k, z in enumerate(ds.coords["nz"]): + ds["W"].values[:, k, :] = z + + grid = UxGrid(ds.uxgrid, z=ds.coords["nz"]) + # Note that the vertical coordinate is required to be the position of the layer interfaces ("nz"), not the mid-layers ("nz1") + P = Field(name="p", data=ds.p, grid=grid, interp_method=UXPiecewiseConstantFace) + + # Test above first cell center - for piecewise constant, should return the depth of the first cell center + assert np.isclose(P.eval(time=ds.time[0].values, z=10.0, y=30.0, x=30.0, applyConversion=False), 55.555557) + # Test below first cell center, but in the first layer - for piecewise constant, should return the depth of the first cell center + assert np.isclose(P.eval(time=ds.time[0].values, z=65.0, y=30.0, x=30.0, applyConversion=False), 55.555557) + # Test bottom layer - for piecewise constant, should return the depth of the of the bottom layer cell center + assert np.isclose(P.eval(time=ds.time[0].values, z=900.0, y=30.0, x=30.0, applyConversion=False), 944.44445801) + + W = Field(name="W", data=ds.W, grid=grid, interp_method=UXPiecewiseLinearNode) + assert np.isclose(W.eval(time=ds.time[0].values, z=10.0, y=30.0, x=30.0, applyConversion=False), 10.0) + assert np.isclose(W.eval(time=ds.time[0].values, z=65.0, y=30.0, x=30.0, applyConversion=False), 65.0) + assert np.isclose(W.eval(time=ds.time[0].values, z=900.0, y=30.0, x=30.0, applyConversion=False), 900.0) diff --git a/tests/v4/test_uxarray_fieldset.py b/tests/v4/test_uxarray_fieldset.py index b804a7750b..fbb2cc3840 100644 --- a/tests/v4/test_uxarray_fieldset.py +++ b/tests/v4/test_uxarray_fieldset.py @@ -117,7 +117,6 @@ def test_fesom_channel(ds_fesom_channel, uvw_fesom_channel): pset.execute(endtime=timedelta(days=1), dt=timedelta(hours=1)) -@pytest.mark.xfail(reason="https://github.com/OceanParcels/Parcels/pull/2026#issuecomment-2945609874") # TODO: Fix def test_fesom2_square_delaunay_uniform_z_coordinate_eval(): """ Test the evaluation of a fieldset with a FESOM2 square Delaunay grid and uniform z-coordinate. @@ -137,4 +136,4 @@ def test_fesom2_square_delaunay_uniform_z_coordinate_eval(): assert fieldset.U.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 assert fieldset.V.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 assert fieldset.W.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 0.0 - assert fieldset.P.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 + assert fieldset.p.eval(time=ds.time[0].values, z=1.0, y=30.0, x=30.0, applyConversion=False) == 1.0 From 7e5d33d8b8e27bd7db0cff86da434e83d2491087 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Wed, 11 Jun 2025 10:31:27 -0400 Subject: [PATCH 10/12] [#2031] Add docstring for new test --- tests/v4/test_field.py | 28 ++++++++++++++++------------ 1 file changed, 16 insertions(+), 12 deletions(-) diff --git a/tests/v4/test_field.py b/tests/v4/test_field.py index 0d914c8d17..f92af2ff6e 100644 --- a/tests/v4/test_field.py +++ b/tests/v4/test_field.py @@ -113,19 +113,11 @@ def test_vectorfield_init_different_time_intervals(): ... -def test_field_unstructured_grid_creation(): ... - - -def test_field_interpolation(): ... - - -def test_field_interpolation_out_of_spatial_bounds(): ... - - -def test_field_interpolation_out_of_time_bounds(): ... - - def test_field_unstructured_z_linear(): + """Tests correctness of piecewise constant and piecewise linear interpolation methods on an unstructured grid with a vertical coordinate. + The example dataset is a FESOM2 square Delaunay grid with uniform z-coordinate. Cell centered and layer registered data are defined to be + linear functions of the vertical coordinate. This allows for testing of exactness of the interpolation methods. + """ ds = datasets_unstructured["fesom2_square_delaunay_uniform_z_coordinate"] # Change the pressure values to be linearly dependent on the vertical coordinate @@ -151,3 +143,15 @@ def test_field_unstructured_z_linear(): assert np.isclose(W.eval(time=ds.time[0].values, z=10.0, y=30.0, x=30.0, applyConversion=False), 10.0) assert np.isclose(W.eval(time=ds.time[0].values, z=65.0, y=30.0, x=30.0, applyConversion=False), 65.0) assert np.isclose(W.eval(time=ds.time[0].values, z=900.0, y=30.0, x=30.0, applyConversion=False), 900.0) + + +def test_field_unstructured_grid_creation(): ... + + +def test_field_interpolation(): ... + + +def test_field_interpolation_out_of_spatial_bounds(): ... + + +def test_field_interpolation_out_of_time_bounds(): ... From 5aea20dad9c329ea4600c809a809891583430462 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Thu, 12 Jun 2025 16:07:48 -0400 Subject: [PATCH 11/12] Deep copy the dataset for modifying values to avoid test pollution --- tests/v4/test_field.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/tests/v4/test_field.py b/tests/v4/test_field.py index f92af2ff6e..6dbe77f85b 100644 --- a/tests/v4/test_field.py +++ b/tests/v4/test_field.py @@ -118,7 +118,7 @@ def test_field_unstructured_z_linear(): The example dataset is a FESOM2 square Delaunay grid with uniform z-coordinate. Cell centered and layer registered data are defined to be linear functions of the vertical coordinate. This allows for testing of exactness of the interpolation methods. """ - ds = datasets_unstructured["fesom2_square_delaunay_uniform_z_coordinate"] + ds = datasets_unstructured["fesom2_square_delaunay_uniform_z_coordinate"].copy(deep=True) # Change the pressure values to be linearly dependent on the vertical coordinate for k, z in enumerate(ds.coords["nz1"]): From c8923532171a7e03fe37c24a67795255593b8f57 Mon Sep 17 00:00:00 2001 From: Joe Schoonover Date: Thu, 12 Jun 2025 16:10:06 -0400 Subject: [PATCH 12/12] Use field.data.values to be explicit about what we're referencing --- parcels/application_kernels/interpolation.py | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/parcels/application_kernels/interpolation.py b/parcels/application_kernels/interpolation.py index 55d2abd5a7..18efd02763 100644 --- a/parcels/application_kernels/interpolation.py +++ b/parcels/application_kernels/interpolation.py @@ -27,7 +27,7 @@ def UXPiecewiseConstantFace( face registered, such as u,v in FESOM. """ zi, fi = field.grid.unravel_index(ei) - return field.values[ti, zi, fi] + return field.data.values[ti, zi, fi] def UXPiecewiseLinearNode( @@ -52,8 +52,8 @@ def UXPiecewiseLinearNode( # The zi refers to the vertical layer index. The field in this routine are assumed to be defined at the vertical interface levels. # For interface zi, the interface indices are [zi, zi+1], so we need to use the values at zi and zi+1. # First, do barycentric interpolation in the lateral direction for each interface level - fzk = np.dot(field.values[ti, k, node_ids], bcoords) - fzkp1 = np.dot(field.values[ti, k + 1, node_ids], bcoords) + fzk = np.dot(field.data.values[ti, k, node_ids], bcoords) + fzkp1 = np.dot(field.data.values[ti, k + 1, node_ids], bcoords) # Then, do piecewise linear interpolation in the vertical direction zk = field.grid.z.values[k]