Acoustic substepping cleanup: typed AcousticDampingStrategy + new default + BW example by glwagner · Pull Request #622 · NumericalEarth/Breeze.jl

glwagner · 2026-04-08T18:45:53Z

Summary

Substepper damping infrastructure: typed AcousticDampingStrategy, two new
damping mechanisms (UpperSponge, VorticityDamping), an alternative
vector-invariant momentum-advection scheme (CompressibleVectorInvariant),
and a constructor-time check that rejects HydrostaticSphericalCoriolis on
non-hydrostatic dynamics.

What's new in the substepper

field on `SplitExplicitTimeDiscretization`	new types
`damping`	`NoDivergenceDamping`, `ThermalDivergenceDamping`, `HyperdiffusiveDivergenceDamping`, `VorticityDamping` (new), or a `Tuple` of these
`sponge`	new — `UpperSponge(damping_rate, depth, ramp=CubicRamp())` for implicit Rayleigh damping at the lid, folded into the column tridiag (Klemp–Dudhia–Hassiotis 2008 form, as in WRF / MPAS-A)

sponge = UpperSponge(damping_rate = 0.2, depth = 5kilometers)
damping = (ThermalDivergenceDamping(coefficient = 0.1),
           VorticityDamping(coefficient = 0.05))
SplitExplicitTimeDiscretization(; damping, sponge)

AbstractRamp family for the sponge profile: CubicRamp (Hermite
smoothstep, default — GPU-cheap), Sin2Ramp (WRF/MPAS literature form),
LinearRamp. Custom shapes via <:AbstractRamp.

Vector-invariant momentum advection (FV3/MPAS-style)

CompressibleVectorInvariant <: AbstractAdvectionScheme decomposes the
flux-form momentum advection as

∇·(ρ𝐯⊗u) = ρ (𝐯·∇)u + u (∇·ρ𝐯)

with (𝐯·∇)u computed in vector-invariant form (vorticity flux +
Bernoulli + vertical advection) by an inner WENOVectorInvariant. Pair
with scalar_advection = WENO():

AtmosphereModel(grid;
                momentum_advection = CompressibleVectorInvariant(),
                scalar_advection   = WENO(),
                ...)

Coriolis safety check

AtmosphereModel now errors at construction if HydrostaticSphericalCoriolis
is passed: it drops the 2Ω cos(φ) cross-terms that couple horizontal
momentum to w, and Breeze always carries prognostic ρw. The
SphericalCoriolis(rotation_rate=Ω) non-hydrostatic form is the only
self-consistent choice on the sphere.

Example

examples/baroclinic_wave.jl drops the explicit timestepper kwarg —
AtmosphereModel already auto-selects :AcousticRungeKutta3 for
CompressibleDynamics(SplitExplicitTimeDiscretization()).

Tests

test/substepper_structural.jl — 17/17 pass (sponge added to
get_coefficient signature; S1 updated to pass nothing).
test/latitude_longitude_grid.jl — switched from
HydrostaticSphericalCoriolis to SphericalCoriolis.
test/acoustic_substepping.jl — passes with new damping API.

Note on branch scope

This branch (glw/hevi-imex-docs) carries a long line of acoustic-substepper
cleanup work plus an unrelated older HEVI/IMEX track that's not the focus of
this PR. Reviewers should focus on the recent commits — the typed damping
API, UpperSponge, VorticityDamping, CompressibleVectorInvariant, and
the HydrostaticSphericalCoriolis guard.

🤖 Generated with Claude Code

numterra-bot

⚠️ Performance Alert ⚠️

Possible performance regression was detected for benchmark 'Breeze.jl Benchmarks'.
Benchmark result of this commit is worse than the previous benchmark result exceeding threshold 1.10.

Benchmark suite	Current: `6ef1a8b`	Previous: `b0e3e87`	Ratio
`CBL; Dynamics: anelastic; Microphysics: nothing [Float32]/Compare advections/NVIDIA L4/WENO5 [256, 256, 128]`	`106590562.60741785` points/s	`121640055.53819701` points/s	`1.14`
`CBL; Dynamics: anelastic; Microphysics: nothing [Float32]/Advection: WENO5/NVIDIA L4/256x256x128`	`106590562.60741785` points/s	`121640055.53819701` points/s	`1.14`

This comment was automatically generated by workflow using github-action-benchmark.

codecov · 2026-04-23T23:45:38Z

Codecov Report

❌ Patch coverage is 50.93633% with 131 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/Thermodynamics/reference_states.jl	34.52%	110 Missing ⚠️
src/CompressibleEquations/compressible_dynamics.jl	50.00%	12 Missing ⚠️
src/TimeSteppers/ssp_runge_kutta_3.jl	44.44%	5 Missing ⚠️
src/CompressibleEquations/time_discretizations.jl	83.33%	2 Missing ⚠️
...ressibleEquations/compressible_density_tendency.jl	80.00%	1 Missing ⚠️
src/TimeSteppers/acoustic_runge_kutta_3.jl	90.90%	1 Missing ⚠️

📢 Thoughts on this report? Let us know!

giordano · 2026-04-25T00:32:39Z

This file isn't added to docs/make.jl, is that intentional?

giordano · 2026-04-28T23:37:47Z

+                                      nothing, ZDirection(),
+                                      Π⁰, θ⁰, γRᵈ, g, δτ_new)
+
+    @info @sprintf("[S1] b[1] = %.6f, c[1] = %.6e (must be 1.0 and 0.0)", b₁, c₁)


Are these prints really necessary? I'm not a fan of printing during tests, if there's any useful information that'd should be tested, not printed. Same applies to all other new tests which print stuff

probably not. I agree, printing during testing is annoying. Learned that lesson with Oceananigans.

giordano · 2026-04-29T13:50:03Z

-    g = constants.gravitational_acceleration
-    Nz = size(grid, 3)
-    loc = (nothing, nothing, Center())
+    g Tᵣstants.gravitational_acceleration
+    NzTᵣze(grid, 3)
+Tᵣ


Something went wrong here

argh. tried to search and replace

I reverted 511a875 entirely, I think many things went wrong with that change

@sprintf

Per Giordano's review on PR #622: tests shouldn't print. The diagnostic @info @sprintf("[Sx] ...", value) lines added for debugging substepper behavior had no role beyond noise — every quantity printed has an accompanying @test that does the actual check. Cleared from: - test/substepper_structural.jl (15 prints, plus a dead growth-rate regression block whose only output was an @info) - test/substepper_rest_state.jl (6 prints) - test/atmosphere_model_construction.jl (1 "Testing NaN Checker..." banner) `using Printf` removed from both substepper test files (no remaining sprintf consumers). Co-Authored-By: Claude Opus 4.7 (1M context) <noreply@anthropic.com>

giordano · 2026-05-07T00:56:59Z

+function SplitExplicitTimeDiscretization(; substeps = nothing,
+                                           forward_weight = 0.65,
+                                           damping = ThermalDivergenceDamping(coefficient = 0.1),
+                                           sponge = nothing,
+                                           substep_distribution = ProportionalSubsteps())
+
+    damping isa Union{AcousticDampingStrategy, Tuple} ||
+        throw(ArgumentError("`damping` must be an `AcousticDampingStrategy` or a tuple of them"))
+
+    sponge isa Union{Nothing, UpperSponge} ||
+        throw(ArgumentError("`sponge` must be `nothing` or an `UpperSponge`"))
+
    return SplitExplicitTimeDiscretization(substeps,
                                           forward_weight,
-                                           divergence_damping_coefficient,
-                                           acoustic_damping_coefficient)
+                                           damping,
+                                           sponge,
+                                           substep_distribution)
 end


I just realised #678 conflicts with this PR, but this code has the same problem as what that PR was trying to solve: it's hard to enforce the same float type everywhere. For example running on Metal GPU it's hard to consistently use Float32 everywhere.

mm yes! We need to add an FT option to the constructor. Is that what you mean?

Yes, and make sure that all fields have the same floating point type.

giordano · 2026-05-07T10:51:17Z

+@inline function stage_substep_count_and_size(::ProportionalSubsteps, β_stage, Δt, N)
+    Δτ = Δt / N
+    Nτ = max(1, round(Int, β_stage * N))
+    return Nτ, Δτ
+end


I think it'd be better to take also the grid as argument and do

FT = eltype(grid) Δτ = FT(Δt) / N

because it's generally more accurate to do the conversion before the division. Same for the other methods below. And then you can replace all the FT(Δτ) -> Δτ

ewquon · 2026-05-12T15:59:48Z

-[`CompressibleDynamics`](@ref) solves the fully compressible Euler equations with prognostic density ``ρ``.
-This formulation retains acoustic waves and is suitable for problems where full compressibility is important.
+[`CompressibleDynamics`](@ref) solves the fully compressible Euler equations with prognostic
+density ``ρ``. The formulation retains acoustic waves and is suitable for problems where full


Suggested change

density ``ρ``. The formulation retains acoustic waves and is suitable for problems where full

total density ``ρ`` (including dry air, vapor, and condensate). The formulation retains acoustic waves and is suitable for problems where full

ewquon · 2026-05-12T16:16:10Z

-    double-counts the ``\theta`` change from earlier stages (since ``\theta(U^*) = \theta^n +
-    \beta_1 \Delta t \, \dot{\theta}^s`` already includes the stage 1 contribution), producing
-    an ``O(\Delta t)`` error per time step that causes instability at larger ``\Delta t``.
+Breeze can also fold the vertical component into the column tridiag by setting


Suggested change

Breeze can also fold the vertical component into the column tridiag by setting

This horizontal divergence damping is applied by default. Breeze can also fold the vertical component into the column tridiag by setting

ewquon · 2026-05-12T16:18:42Z

+Δ(ρu)' &= - γ_h \, ∂_x D_τ , \\
+Δ(ρv)' &= - γ_h \, ∂_y D_τ .


Suggested change

Δ(ρu)' &= - γ_h \, ∂_x D_τ , \\

Δ(ρv)' &= - γ_h \, ∂_y D_τ .

Δ(ρu)' &= - γ_x \, ∂_x D_τ , \\

Δ(ρv)' &= - γ_y \, ∂_y D_τ .

ewquon · 2026-05-12T16:27:24Z

+The off-centering parameter ``ω = 1/2`` is classical centered Crank–Nicolson — neutrally
+stable for the linearized inviscid system but susceptible to amplification of distributed
+floating-point noise through the non-normal substep operator (see
+[Stability analysis](@ref stability-analysis)). The default ``ω = 0.65`` adds modest
+dissipation; the dimensionless parameter ``ε = 2ω - 1 = 0.3`` quantifies the deviation from
+centered.


Suggested change

The off-centering parameter ``ω = 1/2`` is classical centered Crank–Nicolson — neutrally

stable for the linearized inviscid system but susceptible to amplification of distributed

floating-point noise through the non-normal substep operator (see

[Stability analysis](@ref stability-analysis)). The default ``ω = 0.65`` adds modest

dissipation; the dimensionless parameter ``ε = 2ω - 1 = 0.3`` quantifies the deviation from

centered.

The off-centering parameter ``ω = 1/2`` is classical centered Crank–Nicolson — neutrally

stable for the linearized inviscid system but susceptible to amplification of distributed

floating-point noise through the non-normal substep operator (see

[Stability analysis](@ref stability-analysis)).

A fully implicit backward Euler scheme is obtained with ``ω = 1`` and offers the most dissipation.

The default ``ω = 0.65`` adds modest dissipation; the dimensionless parameter ``ε = 2ω - 1 = 0.3`` quantifies the deviation from centered.

glwagner marked this pull request as draft April 8, 2026 19:48

glwagner added the build all examples 🏗️ PRs which should build all examples label Apr 9, 2026

numterra-bot Bot reviewed Apr 9, 2026

View reviewed changes

giordano reviewed Apr 25, 2026

View reviewed changes

giordano force-pushed the glw/hevi-imex-docs branch 2 times, most recently from 3744152 to 2bc4670 Compare April 25, 2026 11:01

giordano reviewed Apr 28, 2026

View reviewed changes

giordano reviewed Apr 29, 2026

View reviewed changes

glwagner commented Apr 29, 2026

View reviewed changes

Comment thread examples/baroclinic_wave.jl Outdated

This was referenced Apr 30, 2026

Microphysics scheduler #664

Draft

Add LatitudeLongitudeGrid support with CompressibleDynamics #520

Closed

glwagner marked this pull request as ready for review May 5, 2026 18:49

glwagner requested review from ewquon and kaiyuan-cheng May 7, 2026 00:17

glwagner commented May 7, 2026

View reviewed changes

Comment thread ext/BreezeReactantExt/Timesteppers.jl Outdated

giordano reviewed May 7, 2026

View reviewed changes

glwagner mentioned this pull request May 7, 2026

_update_pressure_and_average! kernel always called with Float64 argument #679

Closed

giordano reviewed May 7, 2026

View reviewed changes

Rebuild acoustic substepper branch without validation outputs

d92a092

glwagner force-pushed the glw/hevi-imex-docs branch from 6de3601 to d92a092 Compare May 9, 2026 12:03

glwagner and others added 8 commits May 10, 2026 02:52

Fix vertical damping notation in substepper docs

d8c71fb

Merge remote-tracking branch 'origin/main' into glw/hevi-imex-docs

8a3add7

Avoid duplication

e92669c

Match MPAS first acoustic substep sequencing

a8d0729

Allow direct GPU sponge coefficient test

9e43bf7

Document split-explicit outer step cap

93971dc

Document LES timestep limits for split explicit runs

face5b3

Mask normal momentum tendencies on bounded faces

6a24128

glwagner added 3 commits May 11, 2026 18:37

Revert bounded-face tendency masking

3e409c5

Update substepper stability documentation

c603174

Test acoustic parameter float conversion

bb58eeb

dkytezab mentioned this pull request May 12, 2026

Add / Re-factor AcousticSubstepper #567

Closed

glwagner added 7 commits May 12, 2026 02:06

Test acoustic vertical tridiagonal coefficients

c4b7c68

Tighten latitude-longitude metric drift test

2966e52

Apply frozen pressure gradient on first acoustic substep

fc7ebe1

Clarify first substep pressure-gradient comment

62a8e77

Clarify upper sponge ramp documentation

c9bd8cb

Merge branch 'main' into glw/hevi-imex-docs

eb16591

Allow traditional spherical Coriolis

48fdb6c

ewquon reviewed May 12, 2026

View reviewed changes

	density ``ρ``. The formulation retains acoustic waves and is suitable for problems where full
	total density ``ρ`` (including dry air, vapor, and condensate). The formulation retains acoustic waves and is suitable for problems where full

	Breeze can also fold the vertical component into the column tridiag by setting
	This horizontal divergence damping is applied by default. Breeze can also fold the vertical component into the column tridiag by setting

		Δ(ρu)' &= - γ_h \, ∂_x D_τ , \\
		Δ(ρv)' &= - γ_h \, ∂_y D_τ .

Conversation

glwagner commented Apr 8, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Summary

What's new in the substepper

Vector-invariant momentum advection (FV3/MPAS-style)

Coriolis safety check

Example

Tests

Note on branch scope

Uh oh!

numterra-bot Bot 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

⚠️ Performance Alert ⚠️

Uh oh!

codecov Bot commented Apr 23, 2026 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Codecov Report

Uh oh!

Uh oh!

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!

Uh oh!

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!

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!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

3 participants

glwagner commented Apr 8, 2026 •

edited

Loading

numterra-bot Bot left a comment •

edited

Loading

codecov Bot commented Apr 23, 2026 •

edited

Loading