NestedModel abstraction + ERA5 regional hindcast example#233
Draft
ewquon wants to merge 82 commits into
Draft
Conversation
glwagner
reviewed
May 13, 2026
glwagner
reviewed
May 13, 2026
glwagner
reviewed
May 13, 2026
glwagner
reviewed
May 13, 2026
glwagner
reviewed
May 13, 2026
glwagner
reviewed
May 13, 2026
ewquon
force-pushed
the
eq/era5_breeze_land
branch
from
9d673f4 to
5463396
Compare
Codecov Report❌ Patch coverage is 📢 Thoughts on this report? Let us know! |
glwagner
reviewed
May 13, 2026
11 tasks
This reverts commit 5463396.
Collaborator
Author
|
Okay, the model state should be ready for dynamic initialization. A couple of important notes:
|
Member
|
@ewquon should we add terrain? Terrain following should work for the fully compressible formulation (just not substepping) |
Collaborator
Author
Hallelujah! I'm happy to add that in. Anything that gets us closer to reality sooner rather than later is a win in my book. |
glwagner
reviewed
May 14, 2026
glwagner
reviewed
May 14, 2026
The atmosphere→exchange state interpolation read dynamics.reference_state.density, which is `nothing` for CompressibleDynamics (prognostic density, no anelastic reference state), breaking AtmosphereLandModel coupling of a compressible Breeze atmosphere. Use Breeze's dynamics-agnostic dynamics_density / surface_pressure accessors instead. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Regression for the interpolate_state! compressible fix: an AtmosphereLandModel wrapping a Breeze CompressibleDynamics atmosphere constructs and wraps as a NestedSimulation child (exercises the coupled update_state! → interpolate_state! path). Construction-level — stepping the coupled child awaits a Breeze energy-flux/qᵛ fix for the compressible path. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Resolve prescribed_atmosphere.jl: combine the two_dimensional kwarg with main's tripolar velocity_boundary_conditions (#337) — the 2D surface velocities carry the north-fold sign-flip BC; the 3D NestedSimulation parent branch is unchanged. Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
Switch the study period from the 2016-09-10 HI-SCALE clear-sky day to the 20 May 2011 MCS squall line at ARM SGP — the Midlatitude Continental Convective Clouds Experiment (MC3E) case studied by Fan et al. (2017). Initialize at 0000 UTC and force for 18 h. Adds the Fan2017 citation. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Configure the limited-area grid as the inner (Domain 3) member of the 27 → 9 → 3 km telescoping WRF nest used by Fan et al. (2017): 300 × 270 ~3 km cells (their 301 × 271 staggered points − 1) and 50 vertical levels (~60 m near the surface stretching to ~490 m, matching their 51 staggered levels). The 1 km Domain 4 is the natural next NestedSimulation child. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
… ingest
Replace the sigma-z flat-bottom workaround with two coupled changes:
- Ingest ERA5 onto true geopotential heights z = Φ/g via NumericalEarth's
PressureLevelGrid (#241): `FieldTimeSeries(metadata, grid)` / `set!(field,
metadatum)` drive the regrid from the target node heights, deleting the
per-column `interp_z_masked` / `column_interp_z!` loop and the Φ₀
height-above-surface mapping. Pressure is regridded in log-p (the one
quantity for which linear-in-z and linear-in-log-p differ); the others are
linear-in-z, which equals linear-in-log-p under the affine layer map.
- Make the child grid terrain-following: `MutableVerticalDiscretization` +
`follow_terrain!(grid, regrid_topography(grid; ETOPO2022()))` (ETOPO 60",
~1.85 km). The parent stays a regular true-height grid; both regrid from
ERA5 by true Φ/g, so the Interpolated lateral BCs and Davies fringe sample
a consistent state. A local `set_topography!` method routes the ETOPO Field
through `set!` (bridge until Breeze ships it). `CompressibleDynamics` takes
the `terrain_metrics`; profile plots read true per-column heights.
Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add a map showing the ERA5 forcing region (parent) and the 3 km LAM — Fan (2017)'s Domain 3, the NestedSimulation child — as nested rectangles over ETOPO 2022 terrain, centered on ARM SGP. Drawn right after the grid setup so the domain geometry is written even if the run is cut short. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
`rnodes(::PressureLevelGrid)` (hit by the `grid == native_grid(...)` check in `FieldTimeSeries(metadata, grid)`) computed the column-mean z profile as `mean(selectdim(Φ, 3, k))` per level. On a GPU geopotential field that path scalar-indexes (`Statistics._mean` → `first`), which CUDA disallows. Reduce the non-vertical dims in one `mean(; dims)` call (a GPU kernel reduction) and copy the small result to host. CPU behavior is unchanged; surfaced by the era5_breeze terrain smoke test on a T4. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
`FieldTimeSeries(metadata, grid)` set `on_native_grid = grid == native_grid(...)`, and Oceananigans' grid `==` compares nodes regardless of type. For a PressureLevelGrid that means computing the column-mean geopotential profile (`mean_height_profile`) — discarded immediately whenever the target grid (e.g. an interpolation target like the LAM parent grid) differs in type from the native grid, which it always does there. Guard the comparison with a `typeof` check so the node reduction only runs when the grids could actually be equal. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
The prior `mean(Φi; dims)` still scalar-indexed on the GPU: `interior(Φ)` is a device `SubArray`, and `mean` of it routes through `Statistics._mean` → `first`, which CUDA disallows for `dims = :` and `dims = (1, 2)` alike. Copy the (small) geopotential to host before the reduction. The companion change to `FieldTimeSeries`' on_native_grid check means the era5_breeze ingest no longer calls this at all, but it must stay GPU-correct for plot recipes / Lz / the genuine native-grid comparison. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Add Natural Earth coastline, US state lines, and country borders to the nested-domains map (era5_breeze_domains.png) for geographic orientation — the SGP site and the 27→9→3 km nest context read at a glance. Borders are flattened from NaturalEarth (Multi)LineString feature collections to NaN-separated lon/lat via GeoInterface and drawn with plain CairoMakie; the axis is clipped to the map region. Declares NaturalEarth + GeoInterface in docs/Project.toml. Verified rendering standalone (no sim rerun). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…hysical!) `interpolate!` builds its target node via `_node`, whose vertical component is the reference `rnode`; on a MutableVerticalDiscretization (terrain-following) target that ignores the σ/η deformation, so `set!(field, metadatum)` / `Field(metadatum, grid)` placed ERA5 at the LAM's reference (~sea-level) heights. The lowest reference cells then fall below a PressureLevelGrid source's clipped surface, which clamps to finite-but-wrong on CPU and reads out of range on GPU (the 1e40/NaN seen in the era5_breeze terrain run). Add `interpolate_physical!`: a copy of the interpolate! kernel that samples the source at the target's deformed `znode`, dispatched only for mutable targets (regular targets fall back to `interpolate!`, unchanged). Route the metadata `set!`/`Field` regrids through it. CPU-verified: regridding a height field onto a terrain target now reproduces the terrain heights exactly. TODO in the docstring to drop this once Oceananigans' `interpolate!` resolves mutable target nodes from the physical `znode`. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ain map Overlay Natural Earth coastline, US state lines, and country borders on the nested-domains map (era5_breeze_domains.png) for geographic orientation — flattened from NaturalEarth (Multi)LineString feature collections to NaN-separated lon/lat via GeoInterface and drawn with plain CairoMakie; the axis is clipped to the map region. Widen the basemap to a 2.5° buffer around the ERA5 box (so the nest sits well inside the edge) and raise the elevation colorrange to 3000 m for the now-visible Rockies. Declares NaturalEarth + GeoInterface in docs/Project.toml. Verified rendering standalone (no sim rerun). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…xtent The ERA5 forcing region (the D3 parent) was the LAM box + 1° padding. Replace it with the spatial extent of Fan (2017)'s Domain 2 — 181×166 points → 180×165 cells at 9 km (3× the 3 km D3 step), SGP-centered like D3, snapped outward to ERA5's 0.25° grid: λ[-106.75, -88.25] × φ[29.75, 43.5] (~18.5°×13.75°). So the realized nest is ERA5 → D2 → 3 km D3, matching Fan's D2/D3 footprints (extent only; the parent stays at ERA5's ~0.25° resolution). Domain-map legend updated. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
column_fractional_z_index extrapolated the fractional level index outside [1, Nz] for target heights below a column's clipped surface (or above its top). The downstream @inbounds interpolator read then went out of bounds — finite on the CPU (halo), but uninitialized garbage on the GPU, which produced the NaN/Inf initial state in the terrain-following era5_breeze run (the target's lowest cells sit below ERA5's surface where ETOPO terrain is lower than ERA5's). Clamp to the nearest valid level, the intended out-of-range behavior. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…he domain plot Use the full ETOPO relief (regrid_bathymetry, negative over ocean) for the domain-map basemap and plot it with the :topo colormap over a symmetric range (-2500, 2500), so its sea-level break sits at z=0 — ocean depths in blues, land green→yellow→brown→white, with the deep Gulf / high Rockies clamped to keep the Plains gradient legible. Derive the land-sea mask from the relief's sign (is_ocean = bathy < 0) for the eventual SlabLand/ocean surface-BC split. Drop the redundant coastline overlay (the coloring renders it); keep state + country borders. The model terrain still uses regrid_topography (ocean clamped to sea level — correct for the terrain-following lower boundary). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Replace the symmetric ±2500 m colorrange with a two-sided (TwoSlopeNorm-style) mapping: the full bathymetry range fills :topo's lower [0,0.5] (blue) half and the full land range its upper [0.5,1] half, with z=0 pinned to the colormap's sea-level break. Implemented as a custom colormap resampled from :topo so a linear colorrange=(zmin,zmax) reproduces it — keeping the colorbar in physical metres with sea level at its (off-center) 0 tick. Ocean now uses its full blue range; land uses its full green→white range. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…ping (Breeze #712) Adopt the merged Breeze #712 API: - grid `z` → `TerrainFollowingVerticalDiscretization` + `materialize_terrain!`, replacing `MutableVerticalDiscretization` + `follow_terrain!` and dropping the in-example `set_topography!` bridge. - `CompressibleDynamics(SplitExplicitTimeDiscretization(); …)`: terrain physics auto-activates from the grid (no `terrain_metrics`), and the acoustic modes are substepped, so the outer Δt runs at the advection CFL (vertical-binding on the 60 m surface cells), ~25x the old acoustic-CFL Δt. - ERA5 ingest uses plain `interpolate!`; the merged `_node`->`znode` extension makes it terrain-aware, so the local `interpolate_physical!` workaround is unneeded here. Also lands the 80 m AGL u/v/w cut-plane 3x3 (ERA5 / 9 km / 3 km) diagnostic and its `cut_plane` / `era5_winds_on_grid` helpers. Smoke (100 iters, T4) runs end-to-end on the new API with an unchanged initial state. The cold-started (w = 0, hydrostatic) ERA5 IC is not in discrete nonhydrostatic balance; at the substepping Δt this surfaces as a large spurious adjustment, so a physical run needs dynamical initialization (`balance_adiabatically!`, #764) next. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…gid lid Add an `UpperSponge(damping_rate=1/5, depth=3e3)` to the production `SplitExplicitTimeDiscretization` — a 3 km-deep Rayleigh layer damping the vertical momentum (ρw)′ toward the 16.5 km rigid top at a 5 s timescale, so vertically-propagating modes don't reflect off the lid. Damps w only; the upper-level jet runaway and the near-surface terrain transient are separate, still-open issues. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…h sponge Raise the rigid lid from Lz ≈ 16.5 km (~100 hPa, inside the jet) to ~26.5 km (~20 hPa, well above it) at the same Nz = 50 — the stretch cap goes 490 → 1150 m so the upper cells coarsen while the 60 m surface layer is unchanged — and deepen the UpperSponge from 3 km to 5 km. This puts the lid and the Rayleigh layer in the quiescent lower stratosphere (standard NWP practice), so upward-propagating modes are absorbed above the troposphere instead of reflecting back into the jet. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Brings #341 (code alignment), #332 (SoilGrids dataset), #335 (OBC for ocean_simulation). Conflicts resolved: Project.toml Oceananigans compat -> "0.110.1, 0.111" (main's broader bound, for Oceananigans#main); NumericalEarth.bib -> keep both Fan2017 and poggio2021soilgrids. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
glwagner
changed the title
…ed vertical Reconfigure the LAM from the 3 km Domain 3 to the 9 km Domain 2 at a uniform 1/12° step (ERA5 0.25° / 3), so the nest telescopes cleanly 3:1 at each level: ERA5 (D1 role) → D2 (1/12°) → D3 (1/36°, the next nest-down). The ERA5 forcing region is now the child footprint padded 1° and snapped to 0.25°, so the parent strictly encloses the Davies fringe (previously it matched the D2 extent, which the child now occupies). Match Fan (2017)'s vertical grid: 60 m surface cell, 490 m maximum spacing, Nz=50. Fan publishes no stretching ratio or model top (WRF η-coordinate), so a 1.15× stretch lands the top at ~20 km (≈50 hPa), above the ~16 km jet; the Rayleigh sponge is 3 km deep (~17–20 km), in the lower stratosphere. Update the domain map and log lines to the D2 framing, and reduce the cut-plane comparison to ERA5 vs the 9 km child (the 3 km row returns with the Domain 3 nest-down). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Tie τ_relax to the outer step (10·Δt) instead of the advective-crossing estimate 5·Δx/U, which worked out to O(50–700)·Δt — far too weak to hold the lateral boundaries toward the parent. Move the Δt computation up to just after the grid is built so the Davies block can reference it. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
Collaborator
Author
|
Netcdf output (if desired) pending CliMA/Oceananigans.jl#5689 |
The ERA5 analysis cold-starts Cartesian w = 0, but on a terrain-following grid the contravariant velocity w̃ = w − ∇h·u_h should be ≈ 0 at the surface (flow follows the ground). With w = 0 the IC instead carries w̃ = −∇h·u_h, which the surface kinematic BC discharges on the first step. Graft ρw ← ρw − ρw̃ after the first update_state! (which, with ρw = 0, sets ρw̃ = −∇h·ρu_h), giving ρw = ∇h·ρu_h — the contravariant-zeroing momentum, reusing Breeze's own discrete slope so w̃ → 0 to machine precision. This is a correct terrain-consistency improvement to the IC; it does not by itself cure the cold-start blow-up (which is the vertical-advection CFL on the 60 m cells). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…he fringe The lateral BCs feed the child the smooth ERA5 parent state, which assumes the surface sits at the parent orography. At the west inflow edge the child ETOPO is up to +713 m above it, so the boundary-supplied hydrostatic pressure is grossly inconsistent with the child terrain-following surface and discharges as a spurious near-surface horizontal pressure-gradient force — the cold-start blow-up. Blend ETOPO toward the ERA5 parent orography over the first N_taper cells of every lateral edge so the boundary terrain matches the parent and ramps to full ETOPO by the inner fringe edge. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
…_pressure Add `hydrostatic_pressure_from_surface` (src/Atmospheres): a per-column hydrostatic integration of d(ln p)/dz = −g/(RᵐT) upward from the ERA5 surface pressure, anchored at each column's grid bottom face (terrain surface for the terrain-following child, domain bottom for the regular parent) with the moist mixture gas constant Rᵐ = qᵈRᵈ + qᵛRᵛ. This replaces the log-p interpolation `regrid_pressure`. Over high terrain the sub-surface ERA5 pressure levels are clamped to the surface, so the interpolated near-surface pressure (and the EOS density derived from it) was spuriously too dense — a near-sea-level density at ~2.6 km. Injected through the parent-driven lateral BCs, that drove the LAM inflow blow-up. Anchoring on the surface pressure and integrating up keeps the IC in discrete hydrostatic balance and the near-surface density physical. Use the helper for the parent forcing state, the child initial condition, and the ERA5-winds cut-plane diagnostic; delete `regrid_pressure`. Add a unit test (dry isothermal analytic, moist Rᵐ, orography-offset anchor). Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
5 participants
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Summary
@glwagner @kaiyuan-cheng
Adds
examples/era5_breeze.jl— a building block for a regional modeling example that will eventually couple Breeze (compressible solver, in development) to forthcoming SlabLand and SlabOcean components.We download ERA5 reanalysis over an SGP-centered LAM bounding box (HI-SCALE 2016-09-10 case day) and interpolate onto a
LatitudeLongitudeGridsized for ~3 km horizontal cells at the domain center latitude.Tᵥis computed as a derived field usingBreeze.ThermodynamicConstantsfor theRᵥ/Rd − 1coefficient.Punch list
Notes
Nz=1rather than a 2-D(Bounded, Bounded, Flat)grid, sidesteppingCliMA/Oceananigans.jl#5473(Flat↔non-Flatinterpolate!errors with a crypticBoundsError; #5474 will convert that to a clearArgumentErrorbut does not lift the restriction). Mirrors the pattern inexamples/ERA5_hourly_data.jl.LatitudeLongitudeGrid. A future variant on a projected Cartesian grid would benefit from theset!projection support proposed in Add map projection support toset!forRectilinearGridtargets #232.Test plan
~/.cdsapircconfigured.