Skip to content

Use gt4py staggered dimensions for KHalfDim - #1429

Draft
havogt wants to merge 30 commits into
mainfrom
khalfdim-staggering
Draft

Use gt4py staggered dimensions for KHalfDim#1429
havogt wants to merge 30 commits into
mainfrom
khalfdim-staggering

Conversation

@havogt

@havogt havogt commented Aug 11, 2026

Copy link
Copy Markdown
Contributor

Introduces real staggered vertical dimensions using gt4py 1.2.0's
gtx.flip_staggered (ADR next/0026, GridTools/gt4py#2667), replacing the
convention where half-level fields were KDim arrays allocated with
extend={KDim: 1}.

KHalfDim = gtx.flip_staggered(KDim)   # Dimension("_StaggeredK", VERTICAL)

Grid.size[KHalfDim] == num_levels + 1. Static vertical shifts become
half-integer cartesian offsets instead of a FieldOffset.

Index convention

Half level k bounds model level k from above, matching ICON:

read from expression resolves to
a K output index a(KDim - 0.5) half[k]
a K output index a(KDim + 0.5) half[k+1]
a KHalf output index a(KHalfDim + 0.5) full[k]
a KHalf output index a(KHalfDim - 0.5) full[k-1]

The practical rule: an integer ±1 offset that used to cross the two grids
collapses into the ±0.5.

Scope

All packages that allocate vertical fields: common, diffusion, dycore,
testing, bindings, muphys, microphysics, tracer_advection,
standalone_driver. tools and physics_driver have no vertical allocations.

Points worth review:

  • Mixed-grid fused operators use a tuple of per-output domain= dicts, one
    per output field, so a single fused out=(...) can write fields on different
    vertical grids.
  • pressure_ifc is computed by keeping the scan on KDim and staggering its
    result in a separate plain field operator with
    concat_where(KHalfDim == nlev, ...). A scan whose axis differs from its
    output's vertical dimension is rejected by the embedded backend, and embedded
    also cannot take tuple domains together with a scan.
  • Half↔full copy stencils are split per direction, since a single program
    can no longer serve both.
  • no dimension generics in gt4py means vertical-agnostic horizontal helpers
    need one copy per grid (_init_cell_khalf_field_with_zero_wp and friends).
  • as_offset on a staggered field: stagger first, then apply the dynamic
    shift — the offset field must carry the dimension being shifted.

Verification

Against serialized reference data on gtfn_cpu:

suite result
test_velocity_advection.py 22/22
test_solve_nonhydro.py, driver subset 8/8
test_advection.py 4/4
metrics datatests 111

Plus the non-datatest suites (common 606, diffusion 25, dycore stencil 76,
tracer_advection 30, muphys 58, standalone_driver 45), and pre-commit.

Worth flagging for reviewers: the dycore stencil suite's 73 xfails are
backend-conditional.
embedded_remap_error and uses_concat_where xfail only
on embedded, which is the backend that suite is normally verified on, so those
fused stencils were never type-checked. Running them on gtfn_cpu showed their
fixtures still allocated half-level fields on KDim and their numpy references
assigned nlev+1-wide results into nlev-wide slices — several of those
assignments are broken on main too, as the shared *_numpy helpers always
return nlev+1. Those are fixed here; 14 of the 17 affected files are now green
on gtfn_cpu.

Draft — open items

  • 21 solver-test failures fixed (136a4fce2). Root cause: the numpy helper
    solve_tridiagonal_matrix_for_w_back_substitution_numpy was rewritten to take a
    KHalfDim-wide w, but both composite solver tests still handed it the
    KDim-truncated slice next_w[:, :n_lev]. Its nlev = w.shape[1] - 1 then came
    out one short, so the backward sweep pinned half level nlev-2 instead of nlev-1
    and dropped that level's recurrence term; the error propagated upward through
    w_new[k+1], decaying by the z_q factor per level. That is exactly the
    measured signature (deviations at k=35..38 only, peak k=37/38, k=39 and k=40
    clean, magnitudes in the ratio 1 : 1 : z_q : z_q^2).
    test_vertically_implicit_dycore_solver_at_{predictor,corrector}_step.py now
    pass 9/9 and 12/12 on gtfn_cpu with the fixtures' usual unseeded random data.

  • Not yet run: bindings datatests, MPI, diffusion integration,
    test_solve_nonhydro.py in full, and any GPU backend.

  • compute_hydrostatic_correction_term keeps the pre-staggering formulation and
    is handed theta_v_ic relabelled onto KDim via a zero-copy view, because
    dace cannot lower as_offset applied to a staggered access (standalone
    reproducer filed separately). It is the one place in the dycore not using
    staggering, marked with a TODO.

🤖 Generated with Claude Code

havogt and others added 30 commits August 7, 2026 08:14
Make KHalfDim a real staggered dimension (gtx.flip_staggered(KDim)) instead of
a metadata-only marker, and carry it through common, diffusion and dycore.

Half-level fields are now allocated with nlev+1 levels on KHalfDim rather than
over-allocated on KDim via extend={KDim: 1}, and the KHalfDim -> KDim mapping
in states.factory is removed.

Reads that cross the two grids use half-integer cartesian shifts, which follow
ICON's convention that half level k bounds model level k from above:

    div_ic = wgtfac_c * div(KHalfDim + 0.5) + (1 - wgtfac_c) * div(KHalfDim - 0.5)

Programs whose outputs straddle both grids pass a tuple of per-output domains
rather than a single domain dict, keeping the fused field operators intact.

Two helpers needed a per-grid twin because gt4py has no dimension generics:
_cell_2_edge_interpolation and _init_cell_kdim_field_with_zero_wp.

Verified on the embedded backend: dycore stencil tests 75 passed, model/common
606 passed, diffusion 25 passed, and the metrics datatests reproduce the
serialized ICON reference exactly. The metrics -> dycore integration tests and
the gtfn backend are not yet exercised.

Left on KDim for now: pressure_ifc (produced by a KDim scan, consumed only by
muphys/IO), and the tracer_advection, muphys and standalone_driver packages.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Upstream 4c858a6 reworks compute_perturbed_quantities_and_interpolation and
the surrounding stencils: it merges interpolate_cell_field_to_half_levels_{vp,wp}
into one module, renames interpolate_to_surface to
extrapolate_quadratically_to_surface, makes exner_at_cells_on_half_levels and
perturbed_theta_v_at_cells_on_half_levels program local, and gives the outputs
individual domains.

Rather than resolve the conflicts line by line, the conflicted files are taken
from upstream unchanged and the KHalfDim conversion is re-applied on top. That
keeps upstream's restructuring intact and avoids reconciling two rewrites of the
same code by hand.

Upstream's per-output domains use the same tuple-of-dicts idiom this branch
already relies on, so the half-level outputs only needed their dimension
switched.

Verified on the embedded backend after the merge: dycore stencil tests 76
passed, model/common 606 passed, diffusion 25 passed, and the metrics datatests
still reproduce the serialized ICON reference exactly. pre-commit is clean.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Drop direct imports of Dimension objects from icon4py.model.common.dimension
and reference them through the `dims` namespace, matching upstream #1400.
FieldOffsets keep their direct imports.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The gtfn integration tests exposed three fields with inconsistent vertical
grids across the metrics/dycore seam: vn_on_half_levels was still declared
EdgeKField in compute_advection_in_horizontal_momentum_equation, while
horizontal_advection_of_w_at_edges_on_half_levels and vertical_cfl were still
allocated on KDim.

z_v_grad_w now spans nlev+1 levels; ICON stores only nlev, so the test compares
against the reference extent.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
py2fgen takes the array extents from Fortran but the dimension identity from
the wrapper annotation, so a KDim-annotated parameter feeding a KHalfDim state
field is rejected when the program is called.

MetricStateNonHydro.rayleigh_w was still annotated KField while the metrics
factory produces and compute_rayleigh_damping_factor consumes KHalfField; the
stale annotation made the wrapper look consistent.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
VerticalGrid is a dataclass and accepts the mislabelled field, but _vct_a
reaches _compute_rayleigh_w and _compute_scaling_factor_for_3d_divdamp, which
declare fa.KHalfField and reject it.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The gtfn integration tests rejected three fields still allocated on KDim:
nonhydro_buoy_at_cells_on_half_levels, rayleigh_damping_factor and
tangential_wind_on_half_levels, the last also carrying a stale EdgeKField
annotation on IntermediateFields.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Converts p_mflx_contra_v, p_mflx_tracer_v, p_upflux, p_face, vfl_tracer and the
z_cfl/z_face buffers, resolving the TODO(dastrm) markers. This also closes the
dycore seam: the dycore already produced
dynamical_vertical_mass_flux_at_cells_on_half_levels as CellKHalfField into
AdvectionPrepAdvState.mass_flx_ic, which was still declared CellKField.

Still on KDim: the ppm4gpu stencils and the _k_field index fields.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The rayleigh chain is declared fa.KHalfField in dycore_utils and the implicit
solver, so the savepoint readers and the stencil test allocations had to follow.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
dycore_states already declared it CellKHalfField; the test allocations lagged.
Also clamp the nonhydro_buoy comparison to the reference extent: ICON's
z_th_ddz_exner_c stores only nlev levels.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The corrector zeroes prep_adv.dynamical_vertical_mass_flux_at_cells_on_half_levels,
now a CellKHalfField, and it already wrote num_levels + 1 levels. Add the missing
program wrapper around the existing _init_cell_khalf_field_with_zero_wp and use it.

Also clamp the tangential_wind_on_half_levels comparison: ICON's z_vt_ie stores
only nlev levels.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Separates the recurrence from the staggering, which is what makes this work on
every backend: the scan stays on KDim, where its axis matches its output
dimension, and a plain field_operator moves the interface pressure onto
KHalfDim. The surface level comes from a concat_where on an equality predicate;
an open predicate leaves the range unsizable on the embedded backend when the
other branch does not span KHalfDim.

This also converts pressure_on_cells_half_levels in muphys and driver_io, the
pressure_ifc savepoint readers, and DiagnosticState.pressure_ifc.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Reusing _new_te for pressure_ifc_on_model_levels happened to be safe, since it is
rewritten in scatter_to_prognostic before it is read, but a field named after the
new temperature is a confusing place to keep a pressure temporary.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The three ppm4gpu programs write half-level outputs, so their full-level operands
(p_cc, p_cellmass_now, z_delta_q, z_a1) are read at KHalfDim +- n.5. The
_sum_neighbor_contributions helpers are file-local and only ever called from
half-level outputs, so they convert in place.

The PiecewiseParabolicMethod index field had to be split: the three ppm4gpu
programs index half levels while integrate_tracer_vertically indexes model levels.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
mass_flx_ic feeds AdvectionPrepAdvState, now a CellKHalfField, so the two
branches of _make_prep_adv agree again: the other one forwards the dycore's
dynamical_vertical_mass_flux_at_cells_on_half_levels, which has been KHalfDim
since the dycore conversion.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The PPM path zeroes and copies between the full and half vertical grids.
Splitting the copy stencils per direction turns each integer offset that
crossed grids into a half-integer one.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The half-level fixtures were still allocated on KDim, and the numpy
references assigned nlev+1 wide results into nlev wide slices.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…level grid

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
…l grid

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
wgtfac_c moved to KHalfDim, so domain[KDim] raises at granule init.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
dace cannot lower as_offset applied to a staggered access, so this one
stencil keeps the pre-staggering formulation and the caller hands it a
zero-copy view of the half-level buffer.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
The numpy helper was rewritten to take a KHalfDim-wide w, but both composite
solver tests still passed it the KDim-truncated slice, so its sweep stopped one
half level early and dropped the recurrence term at nlev-2.

Co-Authored-By: Claude Opus 5 <noreply@anthropic.com>
@github-actions

Copy link
Copy Markdown

When developing, you can test your changes on CSCS CI before merge with the default pipeline: cscs-ci run default. This will run a default subset of tests.

You can pass options to override pipeline variables, for example:

  • cscs-ci run default;BACKENDS=gtfn_cpu;LEVELS=unit
  • cscs-ci run default;MODEL_SUBPACKAGES=common:driver;SESSIONS=model
    Avoid running the pipeline for all tests when you are developing.

Available options are:

  • SESSIONS: model, model_mpi, or tools (correspond to nox sessions)
  • MODEL_SUBSETS: datatest, basic, or stencils (correspond to nox session selections)
  • MODEL_SUBPACKAGES: subpackages for non-MPI tests (last component, e.g. diffusion, driver)
  • MODEL_MPI_SUBPACKAGES: subpackages for MPI tests (as above)
  • BACKENDS: backends
  • GRIDS: grids for stencil tests (simple, icon_regional, or icon_global)
  • LEVELS: testing level for non-stencil tests (unit or integration)

For each option, all can be used as a shorthand for all possible values of that variable, e.g. LEVELS=all.

See scripts/python/generate_ci_pipeline.py and noxfile.py for available values for each option.

The all pipeline can be run with cscs-ci run all. This will run all icon4py tests in CSCS CI which can be expensive. This pipeline runs on a schedule on main, and can be run when extensive validation is needed (e.g. before releases).

Merging

Once your PR is approved and ready for merging, add it to the merge queue. The merge CSCS CI pipeline will run automatically on the merge-queue branch and must pass before the PR is merged. A dummy merge check will be triggered on the PR itself since it's required to add a PR to the merge queue.

Optional Tests

To run benchmarks you can use:

  • cscs-ci run benchmark-bencher

For more detailed information please look at CI in the EXCLAIM universe.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

1 participant