Default C-grid tracer fields to CGrid_Tracer interpolator by v3ravani · Pull Request #3053 · Parcels-code/Parcels

v3ravani · 2026-11-10T09:37:19Z

Description

Problem

In StructuredModelData.from_sgrid_conventions, all scalar fields defaulted to XLinear, even for tracer fields located at cell centers on staggered C-grids. Following Delandmeter and Van Sebille (2019), C-grid tracers should default to CGrid_Tracer to be interpolated constant over the grid cell.

Solution

Added _select_xinterpolator in src/parcels/_core/model.py to classify scalar fields on SGRID grids. When the model is a staggered C-grid, fields defined at cell centers (lateral face dimensions) default to CGrid_Tracer; all other fields default to XLinear.
Integrated _is_cgrid (aligning with PR Rename _is_agrid to _is_cgrid #3059) to identify C-grid topology using vector field components, default velocity discovery, or SGRID edge dimensions. This reliably identifies C-grids even when a FieldSet has no VectorFields, while ensuring A-grids and B-grids default to XLinear.
Node-located fields, edge-located fields, and unstaggered grids default to XLinear.
In peninsula_dataset, aligned streamfunction P dimensions to ["YG", "XG"] to match its underlying node coordinates lat (YG) and lon (XG) (consistent with stommel_gyre_dataset).

Testing

Added regression tests in tests/test_fieldset.py:
- test_cgrid_tracer_default_interpolator: verifies default interpolators on C-grids.
- test_select_xinterpolator_edge_cases: verifies dimensional edge cases across 2D/3D tracers, nodes, edges, and A-grid configurations.
- test_cgrid_tracer_without_vector_fields: verifies that cell-centred tracers on C-grids default to CGrid_Tracer even when no VectorFields are present.
- test_agrid_and_bgrid_without_vector_fields: verifies that A-grid and B-grid fields default to XLinear without VectorFields.
Verified test suite passes locally (tests/test_fieldset.py, tests/test_advection.py, tests/test_interpolation.py, tests/test_field.py, tests/test_xgrid.py).
Verified with pixi run lint and pixi run typing.

Checklist

Closes CGrid Tracers do not default to CGrid_Tracer Interpolator #3049
Tests added
This PR targets the correct branch (main for normal development, v3-support for v3 support)

AI Disclosure

This PR contains AI-generated content.
- I have tested any AI-generated content in my PR.
- I take responsibility for any AI-generated content in my PR.
- Describe how you used it (e.g., by pasting your prompt): I used an AI coding assistant to help draft _select_xinterpolator, align dataset coordinates, and author regression tests. The code was developed at my direction, and I have read, tested, and validated all changes locally.

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.

erikvansebille · 2026-11-10T11:26:06Z

Thanks for this PR< but you have not filled out the AI statement in our PR template. As per our AI-usage policy, that is very important for us

v3ravani · 2026-11-10T16:10:17Z

Thanks for this PR< but you have not filled out the AI statement in our PR template. As per our AI-usage policy, that is very important for us

Thanks, I have updated the PR description with the AI disclosure.

…terpolator

…taset P coordinates

v3ravani · 2026-11-10T17:18:46Z

Updated with root-cause fix and synced with latest main:

Handled A-grid datasets in _select_xinterpolator so collocated scalar fields default to XLinear.
In peninsula_dataset, aligned streamfunction P dimensions to ['YG', 'XG'] to match its underlying node coordinates lat and lon (consistent with stommel_gyre_dataset), ensuring continuous linear interpolation along particle streamlines.
Verified all unit tests (test_peninsula_fieldset, test_fieldset, test_interpolation, test_field, test_xgrid), typing (mypy), and linting (pre-commit) pass cleanly.

wyatt-fluidnumerics

Hello, thanks for this PR! See comments.

Broadly I think some of the testing introduced in this PR is out of scope since it focuses more on the correctness of the CGrid_Tracer interpolator rather than the interpolator selection.

It also broke CI but that seemed to be resolved as I finalized this review.

wyatt-fluidnumerics · 2026-11-10T17:06:05Z

+                    model.grid.sgrid_metadata,
+                    model.vector_field_components,
+                )
+                if interp_cls is not None:


_select_xinterpolator never returns None

Fixed, directly instantiated without the redundant None check.

wyatt-fluidnumerics · 2026-11-10T17:09:03Z

+        return XLinear
+
+    # Velocity components (U, V, W or components of vector fields) are not tracers
+    velocity_names = {"U", "V", "W"}


I don't think this is needed. CGrid horizontal velocities will be caught in the second check above since they are not on the face, and the vertical velocity gets caught below.

Agreed, removed the velocity names check.

wyatt-fluidnumerics · 2026-11-10T17:10:05Z

+@pytest.mark.parametrize(
+    "dataset_name",
+    [
+        "ds_2d_left",


These datasets all store their velocities/tracers on CGrids, so they all check the exact same logic, only one should be used.

Updated to test only \ds_2d_left.

wyatt-fluidnumerics · 2026-11-10T17:11:50Z

+    assert isinstance(fieldset.UV.interp_method, CGrid_Velocity)
+
+
+def test_cgrid_tracer_piecewise_constant_evaluation():


I think this is out of scope for this PR.

Removed as out of scope.

wyatt-fluidnumerics · 2026-11-10T17:12:52Z

+    )
+
+    # Unstaggered grid metadata (face == node) -> XLinear
+    from parcels._sgrid.core import FaceNodePadding, Padding, SGrid2DMetadata


This import should be moved to the top of the file.

Moved to the top of the file.

wyatt-fluidnumerics · 2026-11-10T17:13:16Z

+    assert _select_xinterpolator(da_unstaggered, unstaggered_metadata) is XLinear
+
+
+def test_cgrid_tracer_particle_sampling():


Out of scope.

Removed as out of scope.

wyatt-fluidnumerics · 2026-11-10T17:18:41Z

+    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


I am confused on what this is checking for, my understanding is that if a single dimension is assigned to both the faces and nodes it would break assertions during Fieldset creation, so I believe this is unneeded.

Good catch, removed.

v3ravani · 2026-11-10T17:51:09Z

Thanks for the review and guidance @wyatt-fluidnumerics!

I have addressed all the feedback in commit 996ee7b2:

Cleaned up _select_xinterpolator:
- Removed the redundant if interp_cls is not None: check in model.py.
- Removed the redundant face == node check and hardcoded velocity name sets; selection now relies cleanly on dimension placement and SGRID face/node coordinates.
Scoped tests strictly to interpolator selection:
- De-parametrized test_cgrid_tracer_default_interpolator to use only ds_2d_left.
- Removed the out-of-scope piecewise constant evaluation and particle sampling tests.
- Moved all SGRID imports to the module top in tests/test_fieldset.py.
- Retained targeted tests for interpolator selection edge cases and A-grid defaults.
Verification:
- All tests in tests/test_fieldset.py, tests/test_advection.py, tests/test_interpolation.py, tests/test_field.py, and tests/test_xgrid.py pass.
- Code formatting, linting (pre-commit), and typing (mypy) pass with zero issues.

wyatt-fluidnumerics

Looks good! Just one question regarding testing. I'll approve once its resolved.

wyatt-fluidnumerics · 2026-11-12T20:34:07Z

+    )
+
+
+def test_agrid_all_fields_default_to_xlinear():


Is this not already tested by the assert above at line 782?

Good catch, you're right. I've removed est_agrid_all_fields_default_to_xlinear as est_select_xinterpolator_edge_cases already verifies that A-grid scalar fields evaluate to XLinear.

erikvansebille · 2026-11-12T21:01:16Z

+    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.
+    """
+    # 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


Can this not be made much simpler by only catching the case where we expect a CGrid_Tracer (a non-vector-field that resides on a c-grid) and in all other cases return XLinear? This seems overly-complicated? Or am I missing something?

Done! Simplified _select_xinterpolator to positively check for the C-grid tracer condition (cell-centered on a C-grid) and default to XLinear in all other cases, removing the extra dataset/vector_fields inspection from the function.

v3ravani · 2026-11-13T10:47:00Z

Thanks @erikvansebille and @wyatt-fluidnumerics for the feedback! I've simplified the interpolator selection logic and addressed the redundant test coverage. Pushed the follow-up changes in 155d6b2 and updated the PR description. Appreciate the guidance!

erikvansebille · 2026-11-13T14:43:45Z

+        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


Thanks for simplifying, @v3ravani - but I still think this logic is too complex. Why the is_vertical_interface check? What does that do? Which cases does it remove? Can you perhaps link to some sketches or background material that show how this check exactly works?

The is_vertical_interface check was originally intended for 3D Lorenz-staggered C-grids (e.g., Fig. 1 in Delandmeter & Van Sebille 2019, https://doi.org/10.5194/gmd-12-3571-2019), where vertical velocity w and vertical mixing reside at cell interfaces in the vertical while lying at cell centers laterally (XC, YC). It guarded against such interface variables defaulting to CGrid_Tracer.

However, you are completely right that this is unnecessary complexity here: velocities are handled via VectorField (CGrid_Velocity), and SGRID metadata specifies 2D lateral topology where {face_x, face_y} cleanly identifies cell-centered tracers.

I have removed the vertical interface logic entirely in commit 7d893552. _select_xinterpolator now simply checks {fnp_x.face, fnp_y.face}.issubset(da.dims) on C-grids and defaults all other fields to XLinear.

erikvansebille · 2026-11-13T15:27:11Z

            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))


I now see that this logic wouldn't work for C-grid FieldSets that have no VectorFields at all.

Can we use the _is_agrid or the new _is_cgrid function from #3059 here - instead of making yet another logic?

Updated! We now adopt the canonical _is_cgrid logic from PR #3059 and expose StructuredModelData.is_cgrid. When vector fields are defined, it delegates to _is_cgrid(ds, u, v). If no vector fields are specified (e.g. vector_fields={}), it inspects the dataset for default velocity components or SGRID edge-staggered variables to distinguish C-grids from A-grids and B-grids.

Added regression tests verifying that cell-centred tracers on C-grid FieldSets without vector fields default to CGrid_Tracer, while A-grid and B-grid fields default to XLinear.

v3ravani · 2026-11-13T16:07:21Z

Updated C-grid detection logic in commit d0ef2545 to address review feedback:

Reused and extended _is_cgrid:
- Replaced redundant inline vector field inspection with StructuredModelData.is_cgrid, delegating to _is_cgrid (from PR Rename _is_agrid to _is_cgrid #3059).
- _is_cgrid supports explicit component pairs (u, v), auto-discovered velocity components, and fallback to SGRID edge dimensions for datasets without vector fields.
- B-grids and A-grids continue to correctly identify as non-C-grids.
Added regression tests:
- test_cgrid_tracer_without_vector_fields: verifies that cell-centred tracers on a C-grid default to CGrid_Tracer even when vector_fields={}.
- test_agrid_and_bgrid_without_vector_fields: verifies that A-grid and B-grid fields default to XLinear without vector fields.
Verification:
- All tests in tests/test_fieldset.py, tests/test_interpolation.py, and tests/test_advection.py pass.
- pixi run lint (pre-commit) and pixi run typing (mypy) pass with no issues.

Default C-grid tracer fields to CGrid_Tracer interpolator

df59982

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.

github-project-automation Bot added this to Parcels development Nov 10, 2026

github-project-automation Bot moved this to Backlog in Parcels development Nov 10, 2026

erikvansebille requested a review from wyatt-fluidnumerics November 10, 2026 16:13

v3ravani added 2 commits November 10, 2026 22:41

Merge remote-tracking branch 'upstream/main' into fix-cgrid-tracer-in…

9109734

…terpolator

Fix CI failures: handle A-grid scalar defaults and align peninsula da…

15fe1d1

…taset P coordinates

wyatt-fluidnumerics requested changes Nov 10, 2026

View reviewed changes

github-project-automation Bot moved this from Backlog to Ready in Parcels development Nov 10, 2026

Address review comments: simplify interpolator selection and scope tests

996ee7b

v3ravani requested a review from wyatt-fluidnumerics November 11, 2026 13:33

wyatt-fluidnumerics reviewed Nov 12, 2026

View reviewed changes

erikvansebille requested changes Nov 12, 2026

View reviewed changes

v3ravani added 2 commits November 13, 2026 15:36

Merge branch 'upstream/main' into fix-cgrid-tracer-interpolator

8f0d33a

Simplify C-grid tracer interpolator selection

155d6b2

erikvansebille reviewed Nov 13, 2026

View reviewed changes

Simplify C-grid tracer selection to face dimensions check

7d89355

erikvansebille reviewed Nov 13, 2026

View reviewed changes

Support C-grid detection without vector fields using _is_cgrid

d0ef254

		assert isinstance(fieldset.UV.interp_method, CGrid_Velocity)


		def test_cgrid_tracer_piecewise_constant_evaluation():

		assert _select_xinterpolator(da_unstaggered, unstaggered_metadata) is XLinear


		def test_cgrid_tracer_particle_sampling():

Conversation

v3ravani commented Nov 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Description

Problem

Solution

Testing

Checklist

AI Disclosure

Uh oh!

erikvansebille commented Nov 10, 2026

Uh oh!

v3ravani commented Nov 10, 2026

Uh oh!

v3ravani commented Nov 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

wyatt-fluidnumerics left a comment • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

wyatt-fluidnumerics Nov 10, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

v3ravani commented Nov 10, 2026

Uh oh!

wyatt-fluidnumerics left a comment

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

v3ravani commented Nov 13, 2026

Uh oh!

Choose a reason for hiding this comment

Uh oh!

v3ravani Nov 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

v3ravani Nov 13, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Choose a reason for hiding this comment

Uh oh!

v3ravani commented Nov 13, 2026

Uh oh!

v3ravani commented Nov 10, 2026 •

edited

Loading

v3ravani commented Nov 10, 2026 •

edited

Loading

wyatt-fluidnumerics left a comment •

edited

Loading

wyatt-fluidnumerics Nov 10, 2026 •

edited

Loading

v3ravani Nov 13, 2026 •

edited

Loading

v3ravani Nov 13, 2026 •

edited

Loading