Skip to content

Materialize a time-varying pressure-level geopotential once per state update - #624

Open
ewquon wants to merge 7 commits into
mainfrom
eq/geopotential-materialization
Open

ewquon wants to merge 7 commits into
mainfrom
eq/geopotential-materialization

Conversation

@ewquon

@ewquon ewquon commented Aug 30, 2026 •

Copy link
Copy Markdown
Collaborator

Closes #623.

The discretization now always stores a plain Field: a TimeSeriesInterpolation is materialized at construction and recomputed by materialize_geopotential! from update_state!(::PrescribedAtmosphere), which NestedModel.time_step! reaches on every nested step. Each probe of the column bisection is a single array read again.

This is a Field type alias rather than a new struct, so rnode, ColumnView, first_above_surface_level, Adapt and on_architecture keep working unchanged. geopotential_data_for_extrema forwards through .operand so Lz and the column-mean profile still span every time slice, not just the resident snapshot, and clipping runs before materialization so every snapshot inherits the clipped columns.

End-to-end on the ERA5 → 12 km nested hindcast of examples/breeze_downscaling_era5.jl (150×136×50, Float32), A100-SXM4-40GB, exclusive node, 400 iterations, arms alternated A/B/A/B:

median s/step
main 0.0478, 0.0495
this branch 0.0368, 0.0359 1.34×

Within-arm spread is 0.0017 and 0.0009 s against a between-arm gap of 0.0123 s, and the arms' ranges do not overlap. Physics is unchanged: ρ ∈ [0.0947, 1.1857] and max|u| 98.6–101.9 m/s across all four runs.

The accuracy cost is zero rather than small. The parent clock advances only inside time_step!(atmos, Δt), which calls update_state! immediately, and the child reads znode between parent ticks, so the snapshot always holds Φ at the parent clock time — over 500 steps of ERA5 data, max |materialized − interpolated| is 0 exactly.

One behavior change: reading rnode after hand-setting clock.time returns the previous snapshot until a refresh runs, and the existing clock-following test is updated to assert that. Separately, FieldTimeSeries defaults to Clamp() time indexing, so past the end of the window the heights freeze on the last snapshot rather than extrapolating; that is now pinned by a test.

🤖 Generated with Claude Code

… update

`column_fractional_z_index` bisects a column of the geopotential for every cell it
interpolates. When the geopotential is a `TimeSeriesInterpolation` — the default
for `ERA5HourlyPressureLevels()`, which passes `z = nothing` — each of those
probes was a `FieldTimeSeries[i, j, k, Time(t)]`, carrying its own binary search
over `times` and a two-slice blend. Within a step every probe asks for the same
Φ(t), so that work was redundant.

Materialize it instead. The discretization now always stores a plain `Field`; a
`TimeSeriesInterpolation` is materialized at construction and recomputed by
`materialize_geopotential!` from `update_state!(::PrescribedAtmosphere)`, which
`NestedModel.time_step!` reaches on every nested step. Each probe is a single
array read again.

Implemented as a `Field` type alias rather than a new struct, so `rnode`,
`ColumnView`, `first_above_surface_level`, `Adapt` and `on_architecture` keep
working unchanged. `geopotential_data_for_extrema` forwards through `.operand` so
`Lz` and the column-mean profile still span every time slice, not just the
resident snapshot. Clipping runs before materialization, so every snapshot
inherits the clipped columns.

Microbenchmark, paired and interleaved in one process, parent 40x40x37 -> child
150x136x50, output bit-identical:

    interpolate! end-to-end      121.5 ms -> 57.4 ms   2.12x
    column_fractional_z_index    109.3 ms -> 20.4 ms   5.36x
    materialize_geopotential!                0.53 ms/tick

End-to-end on the ERA5 -> 12 km nested hindcast
(`examples/breeze_downscaling_era5.jl`, 150x136x50, Float32), A100-SXM4-40GB,
exclusive, 400 iterations, arms alternated A/B/A/B against 2eec75b:

    main       0.0478, 0.0495 s/step
    this       0.0368, 0.0359 s/step        1.34x

Within-arm spread is 0.0017 and 0.0009 s against a between-arm gap of 0.0123 s,
and the two arms' ranges do not overlap. Physics envelopes agree: rho in
[0.0947, 1.1857], max|u| 98.6-101.9 m/s across all four runs.

The accuracy cost is zero rather than small. The parent clock advances only
inside `time_step!(atmos, Δt)`, which calls `update_state!` immediately, and the
child reads `znode` between parent ticks, so the snapshot always holds Φ at the
parent clock time. Over 500 steps of ERA5 data, max |materialized - interpolated|
is 0 exactly.

One behavior change: reading `rnode` after hand-setting `clock.time` returns the
previous snapshot until a refresh runs. The existing clock-following test is
updated to assert that. Note also that `FieldTimeSeries` defaults to `Clamp()`
time indexing, so past the end of the window the heights freeze on the last
snapshot rather than extrapolating; that is now pinned by a test.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01YHo11y7uLQQfPZjo3e12tx
@codecov

codecov Bot commented Aug 30, 2026 •

Copy link
Copy Markdown

Codecov Report

✅ All modified and coverable lines are covered by tests.

📢 Thoughts on this report? Let us know!

ewquon and others added 2 commits September 1, 2026 16:02
The one conflict was in test/test_pressure_level_grid.jl, and it was
structural rather than textual: main parameterized the whole file over
`test_architectures`, while this branch had added testsets written against
the old CPU-only structure.

Resolution: main's harness is kept whole, and this branch's testsets are
ported into it. Every added testset now runs under every architecture in
`test_architectures`, not just CPU:

  - materialized Phi reproduces the interpolated-in-time Phi
  - materialize_geopotential! is a no-op on a static-Field Phi
  - update_state! refreshes the materialized Phi
  - materialized Phi follows the series' time extrapolation
  - single-snapshot TimeSeriesInterpolation
  - materialize_geopotential! refills the snapshot's halos

`LatitudeLongitudeGrid(CPU(); ...)` becomes `LatitudeLongitudeGrid(arch; ...)`,
`make_plg()` becomes `make_plg(arch)`, scalar writes (`interior(f)[i, j, k] = `,
`fts[n][i, j, k] = `) become `set!(field, comprehension)`, host comparisons of
device data go through `Array(...)`, and host-side `rnode` /
`column_fractional_z_index` reads are wrapped in `@allowscalar`, matching main's
idiom: both index the geopotential field directly, so reading them from the host
on a device grid is a scalar read by construction.

Both sides had edited "TimeSeriesInterpolation-backed Phi ignores halo zeros"
and "... Phi heights follow the clock". These keep main's `@allowscalar`
structure and this branch's stale-then-refresh assertions, which pin that
`rnode` holds still until `materialize_geopotential!` runs.

No assertion had to be restricted to CPU. `(@allocated
materialize_geopotential!(grid)) == 0` on the static-Field path dispatches to
the `nothing` fallback with no device work, so it is architecture-independent;
it is verified on CPU here and will first run on GPU in CI.

src/Grids/pressure_level_vertical_discretization.jl auto-merged correctly:
main's only change there shortened a comment above the `column_fractional_z_index`
clamp, and it does not touch the materialization this branch added.

Verified on CPU (CUDA_VISIBLE_DEVICES=-1, Julia 1.12.7, Oceananigans 0.111.0):
test_pressure_level_grid 1471/1471, test_prescribed_atmosphere 18/18,
test_clock_consistency 22/22, test_nested_simulation 263/263. The GPU pass of
the loop is unverified locally and will first run in CI.

Co-Authored-By: Claude Opus 5 (1M context) <noreply@anthropic.com>
@bischtob bischtob changed the title Materialize a time-varying pressure-level geopotential once per state update Materialize a time-varying pressure-level geopotential once per state update [NUM-232] Sep 24, 2026
@bischtob bischtob changed the title Materialize a time-varying pressure-level geopotential once per state update [NUM-232] Materialize a time-varying pressure-level geopotential once per state update Sep 24, 2026

Copy link
Copy Markdown

NUM-232

This branch has not been deployed

No deployments
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.

Time-varying pressure-level geopotential is re-interpolated at every probe of the column bisection

2 participants