From df59982e1e965551b0489af924e450f660ef7c7d Mon Sep 17 00:00:00 2001 From: viraj ravani Date: Tue, 10 Nov 2026 15:06:22 +0530 Subject: [PATCH 1/3] Default C-grid tracer fields to CGrid_Tracer interpolator When constructing fields from SGRID conventions in StructuredModelData, default tracer fields defined at cell centers (face dimensions) on staggered C-grids to CGrid_Tracer rather than hardcoding all scalar fields to XLinear. --- src/parcels/_core/model.py | 48 ++++++++++++- tests/test_fieldset.py | 142 ++++++++++++++++++++++++++++++++++++- 2 files changed, 188 insertions(+), 2 deletions(-) diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index e1c44de30..b70028be7 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -25,6 +25,7 @@ from parcels._logger import logger from parcels._python import NOTSET, NotSetType from parcels.interpolators import ( + CGrid_Tracer, CGrid_Velocity, Ux_Velocity, UxConstantEdgeConstantZC, @@ -275,7 +276,13 @@ def from_sgrid_conventions( model._fields = model.construct_fields() for f in model._fields: if isinstance(f, Field): - f.interp_method = XLinear() + interp_cls = _select_xinterpolator( + model.data[f.name], + model.grid.sgrid_metadata, + model.vector_field_components, + ) + if interp_cls is not None: + f.interp_method = interp_cls() return model @@ -512,6 +519,45 @@ def _discover_ux_U_and_V(ds: ux.UxDataset) -> ux.UxDataset: return ds +def _select_xinterpolator( + da: xr.DataArray, + sgrid_metadata: sgrid.SGrid2DMetadata, + vector_fields: ptyping.VectorFields | None = None, +) -> type[ScalarInterpolator]: + """Selects the appropriate interpolator for a given xarray DataArray on an SGRID grid. + + Following Delandmeter and Van Sebille (2019), tracer fields on an Arakawa C-grid + are located at cell centers (face dimensions) and should be interpolated constant + over the grid cell using CGrid_Tracer. Non-tracer fields, velocity components, + node-located variables, or fields on unstaggered grids default to XLinear. + """ + fnp_x, fnp_y = sgrid_metadata.face_dimensions + if fnp_x.face == fnp_x.node or fnp_y.face == fnp_y.node: + return XLinear + + # Tracers on a C-grid must be defined at the face (cell center) dimensions + face_dims = {fnp_x.face, fnp_y.face} + if not face_dims.issubset(da.dims): + return XLinear + + # Velocity components (U, V, W or components of vector fields) are not tracers + velocity_names = {"U", "V", "W"} + if vector_fields: + for comps in vector_fields.values(): + velocity_names.update(comps) + if da.name is not None and da.name in velocity_names: + return XLinear + + # If vertical dimensions are present and staggered, interface-defined variables + # (e.g. vertical velocity W defined on node) are not cell-centered tracers + if sgrid_metadata.vertical_dimensions: + fnp_z = sgrid_metadata.vertical_dimensions[0] + if fnp_z.face != fnp_z.node and fnp_z.node in da.dims: + return XLinear + + return CGrid_Tracer + + def _select_uxinterpolator(da: ux.UxDataArray): """Selects the appropriate uxarray interpolator for a given uxarray UxDataArray""" supported_uxinterp_mapping = { diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 7fd63f5c9..98f75abc5 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -15,7 +15,11 @@ from parcels import Particle, ParticleFile, ParticleSet, Variable, convert from parcels._core.fieldset import CalendarError, FieldSet, IncompatibleMeshesException, _datetime_to_msg from parcels._core.mesh import SphericalMesh -from parcels._core.model import _default_vector_field_components, _infer_mesh_from_sgrid_data +from parcels._core.model import ( + _default_vector_field_components, + _infer_mesh_from_sgrid_data, + _select_xinterpolator, +) from parcels._core.utils.time import TimeInterval from parcels._datasets.structured.generated import simple_UV_dataset, simple_UV_sigma_dataset from parcels._datasets.structured.generic import datasets as datasets_structured @@ -24,6 +28,7 @@ from parcels._datasets.unstructured.generated import simple_small_delaunay from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels._sgrid.accessor import SGridDatasetInconsistency +from parcels.interpolators import CGrid_Tracer, CGrid_Velocity, XLinear ds = datasets_structured["ds_2d_left"] @@ -722,3 +727,138 @@ def test_fieldset_from_simple_UV_sigma_dataset(n_snapshots): assert depth.dims == ("time", "ZF", "YC", "XC") assert depth.sizes["time"] == n_snapshots assert fieldset.U.grid.get_axis_dim("Z") == nz + + +@pytest.mark.parametrize( + "dataset_name", + [ + "ds_2d_left", + "ds_2d_right", + "ds_2d_inner", + "ds_2d_outer", + "2d_left_rotated", + "2d_left_unrolled_cone", + ], +) +def test_cgrid_tracer_default_interpolator(dataset_name): + """Verify that tracer fields on C-grid default to CGrid_Tracer while non-tracers default to XLinear.""" + ds = datasets_structured[dataset_name] + fieldset = FieldSet.from_sgrid_conventions( + ds, + vector_fields={"UV": ("U_C_grid", "V_C_grid")}, + ) + + # data_c is defined on cell centers (face dimensions) -> CGrid_Tracer + assert isinstance(fieldset.data_c.interp_method, CGrid_Tracer) + + # data_g is defined on cell vertices (node dimensions) -> XLinear + assert isinstance(fieldset.data_g.interp_method, XLinear) + + # Scalar velocity components default to XLinear + assert isinstance(fieldset.U_C_grid.interp_method, XLinear) + assert isinstance(fieldset.V_C_grid.interp_method, XLinear) + assert isinstance(fieldset.U_A_grid.interp_method, XLinear) + assert isinstance(fieldset.V_A_grid.interp_method, XLinear) + + # Vector velocity field defaults to CGrid_Velocity + assert isinstance(fieldset.UV.interp_method, CGrid_Velocity) + + +def test_cgrid_tracer_piecewise_constant_evaluation(): + """Verify that CGrid_Tracer evaluates to a constant value everywhere inside a C-grid cell.""" + ds = datasets_structured["ds_2d_left"] + fieldset = FieldSet.from_sgrid_conventions(ds) + + # Both points (0.5, 0.5, 0.5) and (0.5, 0.51, 0.51) lie within the same C-grid cell + val1 = fieldset.data_c.eval(0.0, 0.5, 0.5, 0.5) + val2 = fieldset.data_c.eval(0.0, 0.5, 0.51, 0.51) + np.testing.assert_allclose(val1, val2) + + # Confirm that XLinear would produce differing values at the two locations + fieldset.data_c.interp_method = XLinear() + val_lin1 = fieldset.data_c.eval(0.0, 0.5, 0.5, 0.5) + val_lin2 = fieldset.data_c.eval(0.0, 0.5, 0.51, 0.51) + assert not np.allclose(val_lin1, val_lin2) + + +def test_select_xinterpolator_edge_cases(): + """Directly test _select_xinterpolator behavior across various dimensional configurations.""" + metadata = datasets_structured["ds_2d_left"].sgrid.metadata + da_3d_tracer = datasets_structured["ds_2d_left"]["data_c"] + da_node = datasets_structured["ds_2d_left"]["data_g"] + da_edge_u = datasets_structured["ds_2d_left"]["U_C_grid"] + + # 3D cell-center tracer -> CGrid_Tracer + assert _select_xinterpolator(da_3d_tracer, metadata) is CGrid_Tracer + + # 2D cell-center tracer (no vertical dim) -> CGrid_Tracer + da_2d_tracer = da_3d_tracer.isel(ZC=0, drop=True) + assert _select_xinterpolator(da_2d_tracer, metadata) is CGrid_Tracer + + # 3D vertical interface variable (node along Z) -> XLinear + da_z_interface = da_3d_tracer.rename({"ZC": "ZG"}) + assert _select_xinterpolator(da_z_interface, metadata) is XLinear + + # Node-defined variable -> XLinear + assert _select_xinterpolator(da_node, metadata) is XLinear + + # Edge-defined velocity variable -> XLinear + assert _select_xinterpolator(da_edge_u, metadata) is XLinear + + # Velocity variable names U, V, W -> XLinear even if defined on face dims + da_u_named = da_3d_tracer.rename("U") + assert _select_xinterpolator(da_u_named, metadata) is XLinear + + # Vector field component names -> XLinear + assert ( + _select_xinterpolator( + da_3d_tracer, + metadata, + vector_fields={"UV": ("data_c", "other")}, + ) + is XLinear + ) + + # Unstaggered grid metadata (face == node) -> XLinear + from parcels._sgrid.core import FaceNodePadding, Padding, SGrid2DMetadata + + unstaggered_metadata = SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("X", "Y"), + face_dimensions=( + FaceNodePadding("X", "X", Padding.NONE), + FaceNodePadding("Y", "Y", Padding.NONE), + ), + ) + da_unstaggered = xr.DataArray(np.zeros((3, 3)), dims=["Y", "X"], name="tracer") + assert _select_xinterpolator(da_unstaggered, unstaggered_metadata) is XLinear + + +def test_cgrid_tracer_particle_sampling(): + """Verify that particles can sample a CGrid_Tracer field during simulation.""" + ds = datasets_structured["ds_2d_left"].copy() + fieldset = FieldSet.from_sgrid_conventions( + ds, + vector_fields={"UV": ("U_C_grid", "V_C_grid")}, + ) + + SampleParticle = Particle.add_variable(Variable("c", dtype=np.float32)) + + def SampleKernel(particles, fieldset): + particles.c = fieldset.data_c[particles.t, particles.z, particles.y, particles.x] + + # Initialize particles within the same C-grid cell + pset = ParticleSet( + fieldset, + pclass=SampleParticle, + x=[0.5, 0.51], + y=[0.5, 0.51], + z=[0.5, 0.5], + t=np.datetime64("2000-01-01"), + ) + pset.execute(SampleKernel, dt=timedelta(hours=1), runtime=timedelta(hours=1)) + + # Both particles started in the same cell, so their sampled tracer value should match + assert pset.c[0] == pset.c[1] + assert not np.isnan(pset.c[0]) From 15fe1d175f4e32282a7d3020eaa046406250e63c Mon Sep 17 00:00:00 2001 From: viraj ravani Date: Tue, 10 Nov 2026 22:47:28 +0530 Subject: [PATCH 2/3] Fix CI failures: handle A-grid scalar defaults and align peninsula dataset P coordinates --- src/parcels/_core/model.py | 12 +++++++- src/parcels/_datasets/structured/generated.py | 2 +- tests/test_fieldset.py | 29 ++++++++++++++++++- 3 files changed, 40 insertions(+), 3 deletions(-) diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index b70028be7..e62c4daae 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -280,6 +280,7 @@ def from_sgrid_conventions( model.data[f.name], model.grid.sgrid_metadata, model.vector_field_components, + dataset=model.data, ) if interp_cls is not None: f.interp_method = interp_cls() @@ -523,14 +524,23 @@ def _select_xinterpolator( da: xr.DataArray, sgrid_metadata: sgrid.SGrid2DMetadata, vector_fields: ptyping.VectorFields | None = None, + dataset: xr.Dataset | None = None, ) -> type[ScalarInterpolator]: """Selects the appropriate interpolator for a given xarray DataArray on an SGRID grid. Following Delandmeter and Van Sebille (2019), tracer fields on an Arakawa C-grid are located at cell centers (face dimensions) and should be interpolated constant over the grid cell using CGrid_Tracer. Non-tracer fields, velocity components, - node-located variables, or fields on unstaggered grids default to XLinear. + node-located variables, or fields on unstaggered/A-grids default to XLinear. """ + # If the model grid is an A-grid (velocity components defined on identical dimensions), + # all fields are collocated; default to XLinear. + if dataset is not None and vector_fields: + for comps in vector_fields.values(): + if len(comps) >= 2 and comps[0] in dataset and comps[1] in dataset: + if _is_agrid(dataset, u=comps[0], v=comps[1]): + return XLinear + fnp_x, fnp_y = sgrid_metadata.face_dimensions if fnp_x.face == fnp_x.node or fnp_y.face == fnp_y.node: return XLinear diff --git a/src/parcels/_datasets/structured/generated.py b/src/parcels/_datasets/structured/generated.py index 04c6fe473..8c696bb43 100644 --- a/src/parcels/_datasets/structured/generated.py +++ b/src/parcels/_datasets/structured/generated.py @@ -310,7 +310,7 @@ def peninsula_dataset(xdim=100, ydim=50, mesh_str="flat", grid_type="A"): { "U": (Udims, U), "V": (Vdims, V), - "P": (["YC", "XC"], P), + "P": (["YG", "XG"], P), }, coords={ "YC": (["YC"], np.arange(ydim) - 0.5, {"axis": "Y"}), diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 98f75abc5..3407dff58 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -28,7 +28,7 @@ from parcels._datasets.unstructured.generated import simple_small_delaunay from parcels._datasets.unstructured.generic import datasets as datasets_unstructured from parcels._sgrid.accessor import SGridDatasetInconsistency -from parcels.interpolators import CGrid_Tracer, CGrid_Velocity, XLinear +from parcels.interpolators import CGrid_Tracer, CGrid_Velocity, XLinear, XLinear_Velocity ds = datasets_structured["ds_2d_left"] @@ -834,6 +834,33 @@ def test_select_xinterpolator_edge_cases(): da_unstaggered = xr.DataArray(np.zeros((3, 3)), dims=["Y", "X"], name="tracer") assert _select_xinterpolator(da_unstaggered, unstaggered_metadata) is XLinear + # A-grid dataset with collocated velocity components -> XLinear for all scalar fields + ds_agrid = datasets_structured["ds_2d_left"] + assert ( + _select_xinterpolator( + da_3d_tracer, + metadata, + vector_fields={"UV": ("U_A_grid", "V_A_grid")}, + dataset=ds_agrid, + ) + is XLinear + ) + + +def test_agrid_all_fields_default_to_xlinear(): + """Verify that when a dataset is configured as an A-grid, all scalar fields default to XLinear.""" + ds = datasets_structured["ds_2d_left"] + fieldset = FieldSet.from_sgrid_conventions( + ds, + vector_fields={"UV": ("U_A_grid", "V_A_grid")}, + ) + # Even though data_c has face dimensions, on an A-grid it must default to XLinear + assert isinstance(fieldset.data_c.interp_method, XLinear) + assert isinstance(fieldset.data_g.interp_method, XLinear) + assert isinstance(fieldset.U_A_grid.interp_method, XLinear) + assert isinstance(fieldset.V_A_grid.interp_method, XLinear) + assert isinstance(fieldset.UV.interp_method, XLinear_Velocity) + def test_cgrid_tracer_particle_sampling(): """Verify that particles can sample a CGrid_Tracer field during simulation.""" From 996ee7b26c774dff3705f14c191a5f783bb75b26 Mon Sep 17 00:00:00 2001 From: viraj ravani Date: Tue, 10 Nov 2026 23:17:00 +0530 Subject: [PATCH 3/3] Address review comments: simplify interpolator selection and scope tests --- src/parcels/_core/model.py | 23 +++------- tests/test_fieldset.py | 90 +------------------------------------- 2 files changed, 7 insertions(+), 106 deletions(-) diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index e62c4daae..a2d91512a 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -276,14 +276,12 @@ def from_sgrid_conventions( model._fields = model.construct_fields() for f in model._fields: if isinstance(f, Field): - interp_cls = _select_xinterpolator( + f.interp_method = _select_xinterpolator( model.data[f.name], model.grid.sgrid_metadata, model.vector_field_components, dataset=model.data, - ) - if interp_cls is not None: - f.interp_method = interp_cls() + )() return model @@ -530,8 +528,8 @@ def _select_xinterpolator( Following Delandmeter and Van Sebille (2019), tracer fields on an Arakawa C-grid are located at cell centers (face dimensions) and should be interpolated constant - over the grid cell using CGrid_Tracer. Non-tracer fields, velocity components, - node-located variables, or fields on unstaggered/A-grids default to XLinear. + over the grid cell using CGrid_Tracer. Non-tracer fields, node-located variables, + or fields on unstaggered/A-grids default to XLinear. """ # If the model grid is an A-grid (velocity components defined on identical dimensions), # all fields are collocated; default to XLinear. @@ -541,23 +539,12 @@ def _select_xinterpolator( if _is_agrid(dataset, u=comps[0], v=comps[1]): return XLinear - fnp_x, fnp_y = sgrid_metadata.face_dimensions - if fnp_x.face == fnp_x.node or fnp_y.face == fnp_y.node: - return XLinear - # Tracers on a C-grid must be defined at the face (cell center) dimensions + fnp_x, fnp_y = sgrid_metadata.face_dimensions face_dims = {fnp_x.face, fnp_y.face} if not face_dims.issubset(da.dims): return XLinear - # Velocity components (U, V, W or components of vector fields) are not tracers - velocity_names = {"U", "V", "W"} - if vector_fields: - for comps in vector_fields.values(): - velocity_names.update(comps) - if da.name is not None and da.name in velocity_names: - return XLinear - # If vertical dimensions are present and staggered, interface-defined variables # (e.g. vertical velocity W defined on node) are not cell-centered tracers if sgrid_metadata.vertical_dimensions: diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 3407dff58..221559f7e 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -729,20 +729,9 @@ def test_fieldset_from_simple_UV_sigma_dataset(n_snapshots): assert fieldset.U.grid.get_axis_dim("Z") == nz -@pytest.mark.parametrize( - "dataset_name", - [ - "ds_2d_left", - "ds_2d_right", - "ds_2d_inner", - "ds_2d_outer", - "2d_left_rotated", - "2d_left_unrolled_cone", - ], -) -def test_cgrid_tracer_default_interpolator(dataset_name): +def test_cgrid_tracer_default_interpolator(): """Verify that tracer fields on C-grid default to CGrid_Tracer while non-tracers default to XLinear.""" - ds = datasets_structured[dataset_name] + ds = datasets_structured["ds_2d_left"] fieldset = FieldSet.from_sgrid_conventions( ds, vector_fields={"UV": ("U_C_grid", "V_C_grid")}, @@ -764,23 +753,6 @@ def test_cgrid_tracer_default_interpolator(dataset_name): assert isinstance(fieldset.UV.interp_method, CGrid_Velocity) -def test_cgrid_tracer_piecewise_constant_evaluation(): - """Verify that CGrid_Tracer evaluates to a constant value everywhere inside a C-grid cell.""" - ds = datasets_structured["ds_2d_left"] - fieldset = FieldSet.from_sgrid_conventions(ds) - - # Both points (0.5, 0.5, 0.5) and (0.5, 0.51, 0.51) lie within the same C-grid cell - val1 = fieldset.data_c.eval(0.0, 0.5, 0.5, 0.5) - val2 = fieldset.data_c.eval(0.0, 0.5, 0.51, 0.51) - np.testing.assert_allclose(val1, val2) - - # Confirm that XLinear would produce differing values at the two locations - fieldset.data_c.interp_method = XLinear() - val_lin1 = fieldset.data_c.eval(0.0, 0.5, 0.5, 0.5) - val_lin2 = fieldset.data_c.eval(0.0, 0.5, 0.51, 0.51) - assert not np.allclose(val_lin1, val_lin2) - - def test_select_xinterpolator_edge_cases(): """Directly test _select_xinterpolator behavior across various dimensional configurations.""" metadata = datasets_structured["ds_2d_left"].sgrid.metadata @@ -805,35 +777,6 @@ def test_select_xinterpolator_edge_cases(): # Edge-defined velocity variable -> XLinear assert _select_xinterpolator(da_edge_u, metadata) is XLinear - # Velocity variable names U, V, W -> XLinear even if defined on face dims - da_u_named = da_3d_tracer.rename("U") - assert _select_xinterpolator(da_u_named, metadata) is XLinear - - # Vector field component names -> XLinear - assert ( - _select_xinterpolator( - da_3d_tracer, - metadata, - vector_fields={"UV": ("data_c", "other")}, - ) - is XLinear - ) - - # Unstaggered grid metadata (face == node) -> XLinear - from parcels._sgrid.core import FaceNodePadding, Padding, SGrid2DMetadata - - unstaggered_metadata = SGrid2DMetadata( - cf_role="grid_topology", - topology_dimension=2, - node_dimensions=("X", "Y"), - face_dimensions=( - FaceNodePadding("X", "X", Padding.NONE), - FaceNodePadding("Y", "Y", Padding.NONE), - ), - ) - da_unstaggered = xr.DataArray(np.zeros((3, 3)), dims=["Y", "X"], name="tracer") - assert _select_xinterpolator(da_unstaggered, unstaggered_metadata) is XLinear - # A-grid dataset with collocated velocity components -> XLinear for all scalar fields ds_agrid = datasets_structured["ds_2d_left"] assert ( @@ -860,32 +803,3 @@ def test_agrid_all_fields_default_to_xlinear(): assert isinstance(fieldset.U_A_grid.interp_method, XLinear) assert isinstance(fieldset.V_A_grid.interp_method, XLinear) assert isinstance(fieldset.UV.interp_method, XLinear_Velocity) - - -def test_cgrid_tracer_particle_sampling(): - """Verify that particles can sample a CGrid_Tracer field during simulation.""" - ds = datasets_structured["ds_2d_left"].copy() - fieldset = FieldSet.from_sgrid_conventions( - ds, - vector_fields={"UV": ("U_C_grid", "V_C_grid")}, - ) - - SampleParticle = Particle.add_variable(Variable("c", dtype=np.float32)) - - def SampleKernel(particles, fieldset): - particles.c = fieldset.data_c[particles.t, particles.z, particles.y, particles.x] - - # Initialize particles within the same C-grid cell - pset = ParticleSet( - fieldset, - pclass=SampleParticle, - x=[0.5, 0.51], - y=[0.5, 0.51], - z=[0.5, 0.5], - t=np.datetime64("2000-01-01"), - ) - pset.execute(SampleKernel, dt=timedelta(hours=1), runtime=timedelta(hours=1)) - - # Both particles started in the same cell, so their sampled tracer value should match - assert pset.c[0] == pset.c[1] - assert not np.isnan(pset.c[0])