From df59982e1e965551b0489af924e450f660ef7c7d Mon Sep 17 00:00:00 2001 From: viraj ravani Date: Tue, 10 Nov 2026 15:06:22 +0530 Subject: [PATCH 1/6] 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/6] 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/6] 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]) From 155d6b27de5ddd967bfba80bce0c2245cb6e9a33 Mon Sep 17 00:00:00 2001 From: viraj ravani Date: Fri, 13 Nov 2026 16:15:13 +0530 Subject: [PATCH 4/6] Simplify C-grid tracer interpolator selection --- src/parcels/_core/model.py | 47 +++++++++++++++----------------------- tests/test_fieldset.py | 30 +++--------------------- 2 files changed, 21 insertions(+), 56 deletions(-) diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index a2d91512a..6e4d8fd8e 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -274,13 +274,13 @@ def from_sgrid_conventions( skip_field_data_validation=skip_field_data_validation, ) model._fields = model.construct_fields() + is_cgrid = any(isinstance(f.interp_method, CGrid_Velocity) for f in model._fields if isinstance(f, VectorField)) for f in model._fields: if isinstance(f, Field): f.interp_method = _select_xinterpolator( model.data[f.name], model.grid.sgrid_metadata, - model.vector_field_components, - dataset=model.data, + is_cgrid=is_cgrid, )() return model @@ -521,38 +521,27 @@ def _discover_ux_U_and_V(ds: ux.UxDataset) -> ux.UxDataset: def _select_xinterpolator( da: xr.DataArray, sgrid_metadata: sgrid.SGrid2DMetadata, - vector_fields: ptyping.VectorFields | None = None, - dataset: xr.Dataset | None = None, + *, + is_cgrid: bool = True, ) -> 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, node-located variables, - or fields on unstaggered/A-grids default to XLinear. + are located at cell centers (face dimensions) and default to CGrid_Tracer. + All other fields 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 - - # 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 - - # 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 + if is_cgrid: + fnp_x, fnp_y = sgrid_metadata.face_dimensions + is_face_centered = {fnp_x.face, fnp_y.face}.issubset(da.dims) + is_vertical_interface = False + if sgrid_metadata.vertical_dimensions: + fnp_z = sgrid_metadata.vertical_dimensions[0] + is_vertical_interface = fnp_z.node in da.dims and fnp_z.face != fnp_z.node + + if is_face_centered and not is_vertical_interface: + return CGrid_Tracer + + return XLinear def _select_uxinterpolator(da: ux.UxDataArray): diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 221559f7e..664250c21 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, XLinear_Velocity +from parcels.interpolators import CGrid_Tracer, CGrid_Velocity, XLinear ds = datasets_structured["ds_2d_left"] @@ -777,29 +777,5 @@ def test_select_xinterpolator_edge_cases(): # Edge-defined velocity variable -> XLinear assert _select_xinterpolator(da_edge_u, 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) + # A-grid configuration -> XLinear + assert _select_xinterpolator(da_3d_tracer, metadata, is_cgrid=False) is XLinear From 7d893552ab8a1da1293e85f1b58bd03d3c8ba962 Mon Sep 17 00:00:00 2001 From: viraj ravani Date: Fri, 13 Nov 2026 20:51:09 +0530 Subject: [PATCH 5/6] Simplify C-grid tracer selection to face dimensions check --- src/parcels/_core/model.py | 8 +------- tests/test_fieldset.py | 4 ---- 2 files changed, 1 insertion(+), 11 deletions(-) diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index 6e4d8fd8e..a102526cc 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -532,13 +532,7 @@ def _select_xinterpolator( """ if is_cgrid: fnp_x, fnp_y = sgrid_metadata.face_dimensions - is_face_centered = {fnp_x.face, fnp_y.face}.issubset(da.dims) - is_vertical_interface = False - if sgrid_metadata.vertical_dimensions: - fnp_z = sgrid_metadata.vertical_dimensions[0] - is_vertical_interface = fnp_z.node in da.dims and fnp_z.face != fnp_z.node - - if is_face_centered and not is_vertical_interface: + if {fnp_x.face, fnp_y.face}.issubset(da.dims): return CGrid_Tracer return XLinear diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index 664250c21..b6c307046 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -767,10 +767,6 @@ def test_select_xinterpolator_edge_cases(): 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 From d0ef254548dcbe2541dda70c46ef857637fb359e Mon Sep 17 00:00:00 2001 From: viraj ravani Date: Fri, 13 Nov 2026 21:35:35 +0530 Subject: [PATCH 6/6] Support C-grid detection without vector fields using _is_cgrid --- src/parcels/_core/model.py | 63 +++++++++++++++++++++++++++++++-- tests/test_fieldset.py | 72 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 130 insertions(+), 5 deletions(-) diff --git a/src/parcels/_core/model.py b/src/parcels/_core/model.py index a102526cc..f72b48146 100644 --- a/src/parcels/_core/model.py +++ b/src/parcels/_core/model.py @@ -217,7 +217,7 @@ def construct_fields(self) -> list[Field | VectorField]: for vfield_name, components in self.vector_field_components.items(): interp_method = ( - XLinear_Velocity() if _is_agrid(self.data, u=components[0], v=components[1]) else CGrid_Velocity() + CGrid_Velocity() if _is_cgrid(self.data, u=components[0], v=components[1]) else XLinear_Velocity() ) component_fields = [single_fields[name] for name in components] @@ -226,6 +226,17 @@ def construct_fields(self) -> list[Field | VectorField]: fields: dict[str, Field | VectorField] = {**single_fields, **vector_fields} return list(fields.values()) + @property + def is_cgrid(self) -> bool: + """Whether the model grid is a staggered C-grid.""" + if self.vector_field_components: + return any( + _is_cgrid(self.data, u=comps[0], v=comps[1]) + for comps in self.vector_field_components.values() + if len(comps) >= 2 and comps[0] in self.data and comps[1] in self.data + ) + return _is_cgrid(self.data, sgrid_metadata=self.grid.sgrid_metadata) + @classmethod def from_sgrid_conventions( cls, @@ -274,7 +285,7 @@ def from_sgrid_conventions( skip_field_data_validation=skip_field_data_validation, ) model._fields = model.construct_fields() - is_cgrid = any(isinstance(f.interp_method, CGrid_Velocity) for f in model._fields if isinstance(f, VectorField)) + is_cgrid = model.is_cgrid for f in model._fields: if isinstance(f, Field): f.interp_method = _select_xinterpolator( @@ -579,10 +590,56 @@ def _select_uxinterpolator(da: ux.UxDataArray): return None +def _is_cgrid( + ds: xr.Dataset, + u: str | None = None, + v: str | None = None, + sgrid_metadata: sgrid.SGrid2DMetadata | None = None, +) -> bool: + """Check if a dataset or specified (u, v) components represent a staggered C-grid. + + Following Arakawa C-grid conventions, horizontal velocity components U and V + are defined on different dimensions (edges). If u and v are defined on identical + dimensions, the grid is an A-grid or B-grid. + + When u and v are not specified (e.g. for scalar-only FieldSets), the dataset + is checked for default velocity components or edge-staggered variables. + """ + if u is not None and v is not None: + # check if U and V are defined on different dimensions + # if yes, interpret as C grid (if no, it could be an A or B grid) + return set(ds[u].dims) != set(ds[v].dims) + + # If dataset has default U and V variables, check their dimensions + default_vfields = _default_vector_field_components(list(ds.data_vars)) + for comps in default_vfields.values(): + if len(comps) >= 2 and comps[0] in ds and comps[1] in ds: + return _is_cgrid(ds, u=comps[0], v=comps[1]) + + # If SGRID metadata is available, check for edge-staggered variables + if sgrid_metadata is None and hasattr(ds, "sgrid"): + try: + sgrid_metadata = ds.sgrid.metadata + except Exception: + pass + + if sgrid_metadata is not None: + fnp_x, fnp_y = sgrid_metadata.face_dimensions + if fnp_x.face != fnp_x.node and fnp_y.face != fnp_y.node: + edge1 = {fnp_x.face, fnp_y.node} + edge2 = {fnp_x.node, fnp_y.face} + has_edge1 = any(edge1.issubset(da.dims) for da in ds.data_vars.values()) + has_edge2 = any(edge2.issubset(da.dims) for da in ds.data_vars.values()) + if has_edge1 or has_edge2: + return True + + return False + + def _is_agrid(ds: xr.Dataset, u: str, v: str) -> bool: # check if U and V are defined on the same dimensions # if yes, interpret as A grid - return set(ds[u].dims) == set(ds[v].dims) + return not _is_cgrid(ds, u=u, v=v) def _get_time_interval(data: xr.DataArray | ux.UxDataArray) -> TimeInterval | None: diff --git a/tests/test_fieldset.py b/tests/test_fieldset.py index b6c307046..598e47484 100644 --- a/tests/test_fieldset.py +++ b/tests/test_fieldset.py @@ -12,9 +12,9 @@ import parcels.tutorial import tests -from parcels import Particle, ParticleFile, ParticleSet, Variable, convert +from parcels import Particle, ParticleFile, ParticleSet, Variable, VectorField, convert from parcels._core.fieldset import CalendarError, FieldSet, IncompatibleMeshesException, _datetime_to_msg -from parcels._core.mesh import SphericalMesh +from parcels._core.mesh import FlatMesh, SphericalMesh from parcels._core.model import ( _default_vector_field_components, _infer_mesh_from_sgrid_data, @@ -775,3 +775,71 @@ def test_select_xinterpolator_edge_cases(): # A-grid configuration -> XLinear assert _select_xinterpolator(da_3d_tracer, metadata, is_cgrid=False) is XLinear + + +def test_cgrid_tracer_without_vector_fields(): + """Verify that cell-centred tracers on C-grid default to CGrid_Tracer even when no VectorFields are present.""" + ds = datasets_structured["ds_2d_left"] + fieldset = FieldSet.from_sgrid_conventions( + ds, + vector_fields={}, + ) + + # 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 retain XLinear + assert isinstance(fieldset.U_C_grid.interp_method, XLinear) + assert isinstance(fieldset.V_C_grid.interp_method, XLinear) + + # FieldSet contains no VectorFields + assert not any(isinstance(f, VectorField) for f in fieldset.fields.values()) + + +def test_agrid_and_bgrid_without_vector_fields(): + """Verify that A-grid and B-grid fields default to XLinear when no VectorFields are present.""" + from parcels._datasets.structured.generated import peninsula_dataset + from parcels._sgrid.core import FaceNodePadding, Padding, SGrid2DMetadata, _attach_sgrid_metadata + + # A-grid dataset with no VectorFields + ds_a = peninsula_dataset(grid_type="A") + fieldset_a = FieldSet.from_sgrid_conventions(ds_a, vector_fields={}) + assert isinstance(fieldset_a.U.interp_method, XLinear) + assert isinstance(fieldset_a.V.interp_method, XLinear) + assert isinstance(fieldset_a.P.interp_method, XLinear) + + # B-grid dataset where U and V are both collocated on nodes (YG, XG) + ds_b = xr.Dataset( + { + "tracer": (["YC", "XC"], np.ones((5, 5))), + "U": (["YG", "XG"], np.ones((6, 6))), + "V": (["YG", "XG"], np.ones((6, 6))), + }, + coords={ + "YC": (["YC"], np.arange(5) + 0.5, {"axis": "Y"}), + "YG": (["YG"], np.arange(6), {"axis": "Y"}), + "XC": (["XC"], np.arange(5) + 0.5, {"axis": "X"}), + "XG": (["XG"], np.arange(6), {"axis": "X"}), + "lon": (["XG"], np.linspace(0, 10, 6), {"axis": "X", "units": "m"}), + "lat": (["YG"], np.linspace(0, 10, 6), {"axis": "Y", "units": "m"}), + }, + ).pipe( + _attach_sgrid_metadata, + SGrid2DMetadata( + cf_role="grid_topology", + topology_dimension=2, + node_dimensions=("XG", "YG"), + node_coordinates=("lon", "lat"), + face_dimensions=( + FaceNodePadding("XC", "XG", Padding.NONE), + FaceNodePadding("YC", "YG", Padding.NONE), + ), + ), + ) + fieldset_b = FieldSet.from_sgrid_conventions(ds_b, mesh=FlatMesh(), vector_fields={}) + assert isinstance(fieldset_b.tracer.interp_method, XLinear) + assert isinstance(fieldset_b.U.interp_method, XLinear) + assert isinstance(fieldset_b.V.interp_method, XLinear)