Acoustic refraction in wind shear — a differentiable example

This example does two things on top of the same compressible forward model. First, we simulate an acoustic pulse propagating through a wind shear layer using the fully compressible Euler equations, and observe how the shear refracts the wave: waves traveling with the wind bend downward (trapped near the surface), while waves traveling against the wind bend upward. The effective propagation speed for a wave in direction $\hat{\boldsymbol{n}}$ is

\[\mathbb{C}^{ac} + \boldsymbol{u} \cdot \hat{\boldsymbol{n}}\]

where $ℂᵃᶜ$ is the acoustic sound speed and $\boldsymbol{u}$ is the wind velocity. Wavefronts tilt toward regions of lower effective propagation speed, "ducting" sound energy along the surface — which is why distant sounds are often heard more clearly downwind. For more on this topic, see

Ostashev, V. E. and Wilson, D. K. (2015). Acoustics in moving inhomogeneous media (CRC Press).
Pierce, A. D. (2019). Acoustics: An introduction to its physical principles and applications (Springer Cham).

Second, we use this setup as a minimal introduction to differentiable atmospheric simulation in Breeze. After running the forward problem, we take a gradient through the entire compressible time-stepping — asking which parts of the wind profile control how much acoustic energy ends up trapped near the surface? — using Enzyme.jl for reverse-mode automatic differentiation and Reactant.jl to compile the model down to XLA. The result is a 2D sensitivity field obtained in a single backward pass; the same answer via finite differences would cost one model rerun per grid cell.

We use stable stratification to suppress Kelvin-Helmholtz instability and a logarithmic wind profile consistent with the atmospheric surface layer.

using Breeze
using Oceananigans: Oceananigans
using Oceananigans.Units
using Printf
using CairoMakie

Grid and model setup

Nx, Nz = 128, 64
Lx, Lz = 1000, 200  # (m)

grid = RectilinearGrid(size = (Nx, Nz), x = (-Lx/2, Lx/2), z = (0, Lz),
                       topology = (Periodic, Flat, Bounded))

128×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── Periodic x ∈ [-500.0, 500.0) regularly spaced with Δx=7.8125
├── Flat y                       
└── Bounded  z ∈ [0.0, 200.0]    regularly spaced with Δz=3.125

This example is dry, so the θˡⁱ→T inversion has an exact closed form and needs no Newton steps. temperature_tolerance = 0 selects the fixed-trip (unrolled) inversion and temperature_maxiter = 2 keeps the trip count tiny: the differentiable pass at the end of this example takes a reverse-mode gradient through the compressible time step with Reactant/Enzyme, and the default tolerance-based while inversion compiles to an XLA while op that does not differentiate cheaply, while a high iteration count needlessly inflates the traced graph (NumericalEarth/Breeze.jl#767). The forward and adjoint models use identical dynamics.

model = AtmosphereModel(grid; dynamics = CompressibleDynamics(ExplicitTimeStepping(); temperature_tolerance = 0, temperature_maxiter = 2))

AtmosphereModel{CPU, RectilinearGrid}(time = 0 seconds, iteration = 0)
├── grid: 128×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on CPU with 3×0×3 halo
├── dynamics: CompressibleDynamics{ExplicitTimeStepping}
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵛ: Centered(order=2)
├── forcing: @NamedTuple{ρ::Returns{Float64}, ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵛ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Background state

We build a hydrostatically balanced reference state using ReferenceState. This provides the background density and pressure profiles.

constants = model.thermodynamic_constants

θ₀ = 300      # Reference potential temperature (K)
p₀ = 101325   # Surface pressure (Pa)
pˢᵗ = 1e5     # Standard pressure (Pa)

reference = ReferenceState(grid, constants; surface_pressure=p₀, potential_temperature=θ₀, standard_pressure=pˢᵗ)

ReferenceState{Float64}(p₀=101325.0, θ₀=300.0, pˢᵗ=100000.0)

The sound speed at the surface determines the acoustic wave propagation speed.

Rᵈ = constants.molar_gas_constant / constants.dry_air.molar_mass
cᵖᵈ = constants.dry_air.heat_capacity
γ = cᵖᵈ / (cᵖᵈ - Rᵈ)
ℂᵃᶜ = sqrt(γ * Rᵈ * θ₀)

347.15629052210863

The wind profile follows the classic log-law of the atmospheric surface layer.

U₀ = 20 # Surface velocity (m/s), u★/κ
ℓ = 1   # Roughness length (m), like, shrubs and stuff

Uᵢ(z) = U₀ * log((z + ℓ) / ℓ)

Uᵢ (generic function with 1 method)

Initial conditions

We initialize a localized Gaussian density pulse representing an acoustic disturbance. For a rightward-propagating acoustic wave, the velocity perturbation is in phase with the density perturbation: $u' = (ℂᵃᶜ / ρ₀) ρ'$.

δρ = 0.01    # Density perturbation amplitude (kg/m³)
σ = 20       # Pulse width (m)

gaussian(x, z) = exp(-(x^2 + z^2) / 2σ^2)
ρ₀ = interior(reference.density, 1, 1, 1)[]

ρᵢ_func(x, z) = adiabatic_hydrostatic_density(z, p₀, θ₀, pˢᵗ, constants) + δρ * gaussian(x, z)
uᵢ_func(x, z) = Uᵢ(z) # + (ℂᵃᶜ / ρ₀) * δρ * gaussian(x, z)

set!(model, ρ=ρᵢ_func, θ=θ₀, u=uᵢ_func)

Simulation setup

Acoustic waves travel fast ($ℂᵃᶜ ≈ 347$ m/s), so we need a small time step. The Courant–Friedrichs–Lewy (CFL) condition is based on the effective propagation speed $ℂᵃᶜ + \mathrm{max}(U)$.

Δx, Δz = Lx / Nx, Lz / Nz
Δt = 0.5 * min(Δx, Δz) / (ℂᵃᶜ + Uᵢ(Lz))
stop_time = 0.5 # (s) — long enough for the wave to traverse the domain and for refraction to bend rays visibly

simulation = Simulation(model; Δt, stop_time)
Oceananigans.Diagnostics.erroring_NaNChecker!(simulation)

function progress(sim)
    u, v, w = sim.model.velocities
    msg = @sprintf("Iter: %d, t: %s, max|u|: %.2f m/s, max|w|: %.2f m/s",
                   iteration(sim), prettytime(sim),
                   maximum(abs, u), maximum(abs, w))
    @info msg
end

add_callback!(simulation, progress, IterationInterval(100))

Output

We perturbation fields for density and x-velocity for visualization.

ρ = model.dynamics.density
u, v, w = model.velocities

ρᵇᵍ = CenterField(grid)
uᵇᵍ = XFaceField(grid)

set!(ρᵇᵍ, (x, z) -> adiabatic_hydrostatic_density(z, p₀, θ₀, pˢᵗ, constants))
set!(uᵇᵍ, (x, z) -> Uᵢ(z))

ρ′ = Field(ρ - ρᵇᵍ)
u′ = Field(u - uᵇᵍ)

U = Average(u, dims = 1)
R = Average(ρ, dims = 1)
W² = Average(w^2, dims = 1)

filename = "acoustic_wave.jld2"
outputs = (; ρ′, u′, w, U, R, W²)

simulation.output_writers[:jld2] = JLD2Writer(model, outputs; filename,
                                              schedule = TimeInterval(0.01),
                                              overwrite_existing = true)

run!(simulation)

[ Info: Initializing simulation...
[ Info: Iter: 0, t: 0 seconds, max|u|: 105.91 m/s, max|w|: 0.00 m/s
[ Info:     ... simulation initialization complete (25.860 seconds)
[ Info: Executing initial time step...
[ Info:     ... initial time step complete (1.903 seconds).
[ Info: Iter: 100, t: 333.448 ms, max|u|: 105.91 m/s, max|w|: 0.48 m/s
[ Info: Simulation is stopping after running for 30.731 seconds.
[ Info: Simulation time 500 ms equals or exceeds stop time 500 ms.

Visualization

Load the saved perturbation fields and create a snapshot.

ρ′ts = FieldTimeSeries(filename, "ρ′")
u′ts = FieldTimeSeries(filename, "u′")
wts = FieldTimeSeries(filename, "w")
Uts = FieldTimeSeries(filename, "U")
Rts = FieldTimeSeries(filename, "R")
W²ts = FieldTimeSeries(filename, "W²")

times = ρ′ts.times
Nt = length(times)

fig = Figure(size = (900, 600), fontsize = 12)

axρ = Axis(fig[1, 2]; aspect = 5, ylabel = "z (m)")
axw = Axis(fig[2, 2]; aspect = 5, ylabel = "z (m)")
axu = Axis(fig[3, 2]; aspect = 5, xlabel = "x (m)", ylabel = "z (m)")
axR = Axis(fig[1, 1]; xlabel = "⟨ρ⟩ (kg/m³)")
axW = Axis(fig[2, 1]; xlabel = "⟨w²⟩ (m²/s²)", limits = (extrema(W²ts), nothing))
axU = Axis(fig[3, 1]; xlabel = "⟨u⟩ (m/s)")

hidexdecorations!(axρ)
hidexdecorations!(axw)
colsize!(fig.layout, 1, Relative(0.2))

n = Observable(Nt)
ρ′n = @lift ρ′ts[$n]
u′n = @lift u′ts[$n]
wn = @lift wts[$n]
Un = @lift Uts[$n]
Rn = @lift Rts[$n]
W²n = @lift W²ts[$n]

ρlim = δρ / 4
ulim = 1

hmρ = heatmap!(axρ, ρ′n; colormap = :balance, colorrange = (-ρlim, ρlim))
hmw = heatmap!(axw, wn; colormap = :balance, colorrange = (-ulim, ulim))
hmu = heatmap!(axu, u′n; colormap = :balance, colorrange = (-ulim, ulim))

lines!(axR, Rn)
lines!(axW, W²n)
lines!(axU, Un)

Colorbar(fig[1, 3], hmρ; label = "ρ′ (kg/m³)")
Colorbar(fig[2, 3], hmw; label = "w (m/s)")
Colorbar(fig[3, 3], hmu; label = "u′ (m/s)")

title = @lift "Acoustic wave in log-layer shear — t = $(prettytime(times[$n]))"
fig[0, :] = Label(fig, title, fontsize = 16, tellwidth = false)

CairoMakie.record(fig, "acoustic_wave.mp4", 1:Nt, framerate = 18) do nn
    n[] = nn
end

A differentiable workflow

The forward simulation above gives us the physics; the rest of the example treats that same forward model as a function and takes its gradient. This is the minimal pattern you'd reach for whenever you want to do data-assimilation, parameter calibration, or sensitivity analysis with a Breeze atmosphere — wrap the time-stepping in a scalar-valued loss, compile it with Reactant.@compile, and differentiate it with Enzyme.autodiff.

The pattern is:

Rebuild the model on a ReactantState grid so all arrays are XLA buffers.
Choose a differentiated input (here, the initial wind field).
Define a scalar loss that re-initializes the model from that input, runs nsteps of time_step! inside a @trace loop, and reduces to a scalar diagnostic.
Wrap loss in a grad_loss that calls Enzyme.autodiff(...) with the input as Duplicated and everything else as Const.
Compile once with Reactant.@compile raise=true raise_first=true; run many times.

Why Reactant?

Reactant traces Julia code into an intermediate representation (StableHLO) that XLA can optimize and Enzyme can differentiate. The key requirement is that the model lives on ReactantState — Reactant's architecture in Oceananigans — so that all arrays are XLA buffers. We therefore rebuild the same physical setup on a new grid whose architecture is ReactantState().

using Reactant, CUDA    # CUDA is required for loading the Reactant extension
using Enzyme
using Statistics: mean
using Oceananigans.Architectures: ReactantState
using Reactant: @trace

Reactant.set_default_backend("cpu")

Rebuild the grid and model on ReactantState.

grid_ad = RectilinearGrid(ReactantState(); size = (Nx, Nz),
                          x = (-Lx/2, Lx/2), z = (0, Lz),
                          topology = (Periodic, Flat, Bounded))

model_ad = AtmosphereModel(grid_ad; dynamics = CompressibleDynamics(ExplicitTimeStepping(); temperature_tolerance = 0, temperature_maxiter = 2)) # fixed-trip, low-iteration EOS inversion so Enzyme can differentiate it cheaply (see forward model)

AtmosphereModel{ReactantState, RectilinearGrid}(time = 0 seconds, iteration = Reactant.ConcretePJRTNumber{Int64, 1}(0))
├── grid: 128×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on ReactantState with 3×0×3 halo
├── dynamics: CompressibleDynamics{ExplicitTimeStepping}
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵛ: Centered(order=2)
├── forcing: @NamedTuple{ρ::Returns{Float64}, ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵛ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Fixed and varying fields

In this experiment the initial density pulse and hydrostatic background are held fixed; only the wind profile varies. We therefore precompute the total initial density once (background + Gaussian pulse) and use it as a Const. We also keep the standalone background $\bar\rho(z)$ around so the loss can subtract it from the model density to isolate the acoustic perturbation.

ρᵇᵍ = CenterField(grid_ad)
set!(ρᵇᵍ, (x, z) -> adiabatic_hydrostatic_density(z, p₀, θ₀, pˢᵗ, constants))

ρ_total = CenterField(grid_ad)
set!(ρ_total, ρᵢ_func)

The initial wind field is the quantity we differentiate with respect to. Enzyme accumulates $∂J / ∂u_i$ into the shadow buffer $du_i$.

uᵢ  = XFaceField(grid_ad)
duᵢ = XFaceField(grid_ad)
set!(uᵢ, (x, z) -> Uᵢ(z))
set!(duᵢ, 0)

128×1×64 Field{Oceananigans.Grids.Face, Oceananigans.Grids.Center, Oceananigans.Grids.Center} on Oceananigans.Grids.RectilinearGrid on ReactantState
├── grid: 128×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on ReactantState with 3×0×3 halo
├── boundary conditions: FieldBoundaryConditions
│   └── west: Periodic, east: Periodic, south: Nothing, north: Nothing, bottom: ZeroFlux, top: ZeroFlux, immersed: Nothing
└── data: 134×1×70 OffsetArray(::Reactant.ConcretePJRTArray{Float64,3}, -2:131, 1:1, -2:67) with eltype Float64 with indices -2:131×1:1×-2:67
    └── max=0.0, min=0.0, mean=0.0

The shadow model stores accumulated adjoints for every prognostic field.

dmodel_ad = Enzyme.make_zero(model_ad)

AtmosphereModel{ReactantState, RectilinearGrid}(time = 0 seconds, iteration = Reactant.ConcretePJRTNumber{Int64, 1}(0))
├── grid: 128×1×64 RectilinearGrid{Float64, Periodic, Flat, Bounded} on ReactantState with 3×0×3 halo
├── dynamics: CompressibleDynamics{ExplicitTimeStepping}
├── formulation: LiquidIcePotentialTemperatureFormulation
├── thermodynamic_constants: ThermodynamicConstants{Float64}
├── timestepper: SSPRungeKutta3
├── advection scheme: 
│   ├── momentum: Centered(order=2)
│   ├── ρθ: Centered(order=2)
│   └── ρqᵛ: Centered(order=2)
├── forcing: @NamedTuple{ρ::Returns{Float64}, ρu::Returns{Float64}, ρv::Returns{Float64}, ρw::Returns{Float64}, ρθ::Returns{Float64}, ρqᵛ::Returns{Float64}, ρe::Returns{Float64}}
├── tracers: ()
├── coriolis: Nothing
└── microphysics: Nothing

Time step and integration length

We reuse the CFL-based time step and the exact number of iterations from the forward simulation above. The Simulation API is not used here because Reactant compiles a fixed-length traced loop instead. Gradient checkpointing requires a perfect-square step count, so we round up to the next perfect square.

Nt = simulation.model.clock.iteration
Nsteps = (isqrt(Nt - 1) + 1)^2

Defining the objective

We measure the mean squared acoustic density anomaly along the bottom of the domain:

\[J \;=\; \frac{1}{N_x}\sum_{i}\bigl[\rho(x_i, z_1) - \bar\rho(x_i, z_1)\bigr]^2\]

This is a global measure of how much acoustic energy ends up trapped near the surface. Averaging along the whole bottom row gives a sensitivity field that lights up wherever the wind affects any part of the surface response, making the ducting pattern visible across the entire domain.

The set! inside loss is what re-initializes the model from the current wind field on every backward evaluation. Without it, AD would differentiate a stale trajectory.

function loss(model, uᵢ, ρ_total, ρᵇᵍ, θ₀, Δt, nsteps)
    set!(model; ρ = ρ_total, θ = θ₀, u = uᵢ)
    @trace mincut=true checkpointing=true track_numbers=false for _ in 1:nsteps
        time_step!(model, Δt)
    end
    ρ₀  = interior(model.dynamics.density, :, :, 1)
    ρᵇ₀ = interior(ρᵇᵍ, :, :, 1)
    return mean((ρ₀ .- ρᵇ₀).^2)
end

loss (generic function with 1 method)

The gradient wrapper

grad_loss zeroes the adjoint buffer and calls Enzyme.autodiff in reverse mode. The model and the initial wind are Duplicated (primal + shadow); everything else is Const.

function grad_loss(model, dmodel, uᵢ, duᵢ, ρ_total, ρᵇᵍ, θ₀, Δt, nsteps)
    parent(duᵢ) .= 0
    _, J = Enzyme.autodiff(
        Enzyme.set_strong_zero(Enzyme.ReverseWithPrimal),
        loss, Enzyme.Active,
        Enzyme.Duplicated(model, dmodel),
        Enzyme.Duplicated(uᵢ, duᵢ),
        Enzyme.Const(ρ_total),
        Enzyme.Const(ρᵇᵍ),
        Enzyme.Const(θ₀),
        Enzyme.Const(Δt),
        Enzyme.Const(nsteps))
    return duᵢ, J
end

grad_loss (generic function with 1 method)

Compilation and execution

Reactant.@compile traces the function once to build an XLA executable. The flags raise=true and raise_first=true ensure that every KernelAbstractions kernel is "raised" to StableHLO before Enzyme differentiates through it — a requirement for the backward pass.

@info "Compiling differentiated model — this may take a minute..."
compiled_grad = Reactant.@compile raise=true raise_first=true sync=true grad_loss(
    model_ad, dmodel_ad, uᵢ, duᵢ,
    ρ_total, ρᵇᵍ, θ₀, Δt, Nsteps)

@info "Running gradient..."
du, J = compiled_grad(model_ad, dmodel_ad, uᵢ, duᵢ, ρ_total, ρᵇᵍ, θ₀, Δt, Nsteps)

xs_u = xnodes(grid_ad, Face())
zs   = znodes(grid_ad, Center())

@info @sprintf("Surface-mean (ρ - ρ̄)² = %.6e after %d steps", Float64(only(J)), Nsteps)

[ Info: Compiling differentiated model — this may take a minute...
[ Info: Running gradient...
┌ Error: Exception while generating log record in module Main.var"##277" at /__w/Breeze.jl/Breeze.jl/docs/src/literated/acoustic_wave.md:12
│   exception =
│    MethodError: no method matching writeexp(::Vector{UInt8}, ::Int64, ::Reactant.ConcretePJRTNumber{Float64, 1}, ::Int64, ::Bool, ::Bool, ::Bool, ::Char, ::UInt8)
│    The function `writeexp` exists, but no method is defined for this combination of argument types.
│    
│    Closest candidates are:
│      writeexp(::Any, ::Any, !Matched::T, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any, !Matched::Any) where T<:Union{Float16, Float32, Float64}
│       @ Base ryu/exp.jl:1
│      writeexp(::Any, ::Any, !Matched::T, ::Any, ::Any, ::Any, ::Any, ::Any, ::Any) where T<:Union{Float16, Float32, Float64}
│       @ Base ryu/exp.jl:1
│      writeexp(::Any, ::Any, !Matched::T, ::Any, ::Any, ::Any, ::Any, ::Any) where T<:Union{Float16, Float32, Float64}
│       @ Base ryu/exp.jl:1
│      ...
│    
│    Stacktrace:
│      [1] fmt(buf::Vector{UInt8}, pos::Int64, arg::Reactant.ConcretePJRTNumber{Float64, 1}, spec::Printf.Spec{Val{'e'}})
│        @ Printf /usr/local/julia/share/julia/stdlib/v1.12/Printf/src/Printf.jl:523
│      [2] fmt
│        @ /usr/local/julia/share/julia/stdlib/v1.12/Printf/src/Printf.jl:315 [inlined]
│      [3] format
│        @ /usr/local/julia/share/julia/stdlib/v1.12/Printf/src/Printf.jl:836 [inlined]
│      [4] format(::Printf.Format{Base.CodeUnits{UInt8, String}, Tuple{Printf.Spec{Val{'e'}}, Printf.Spec{Val{'d'}}}}, ::Reactant.ConcretePJRTNumber{Float64, 1}, ::Int64)
│        @ Printf /usr/local/julia/share/julia/stdlib/v1.12/Printf/src/Printf.jl:946
│      [5] top-level scope
│        @ /__w/Breeze.jl/Breeze.jl/docs/src/literated/acoustic_wave.md:397
│      [6] macro expansion
│        @ logging/logging.jl:387 [inlined]
│      [7] eval(m::Module, e::Any)
│        @ Core ./boot.jl:489
│      [8] include_string(mapexpr::typeof(identity), mod::Module, code::String, filename::String)
│        @ Base ./loading.jl:2874
│      [9] include_string
│        @ ./loading.jl:2884 [inlined]
│     [10] #61
│        @ /usr/local/share/julia/packages/Literate/On3EJ/src/Literate.jl:959 [inlined]
│     [11] task_local_storage(body::Literate.var"#61#62"{String, Bool, Module, String}, key::Symbol, val::String)
│        @ Base ./task.jl:298
│     [12] #59
│        @ /usr/local/share/julia/packages/Literate/On3EJ/src/Literate.jl:955 [inlined]
│     [13] (::IOCapture.var"#12#13"{Type{Union{}}, Literate.var"#59#60"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}})()
│        @ IOCapture /usr/local/share/julia/packages/IOCapture/MR051/src/IOCapture.jl:170
│     [14] with_logstate(f::IOCapture.var"#12#13"{Type{Union{}}, Literate.var"#59#60"{String, Bool, Module, String}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}, IOContext{Base.PipeEndpoint}}, logstate::Base.CoreLogging.LogState)
│        @ Base.CoreLogging ./logging/logging.jl:542
│     [15] with_logger(f::Function, logger::Base.CoreLogging.ConsoleLogger)
│        @ Base.CoreLogging ./logging/logging.jl:653
│     [16] capture(f::Literate.var"#59#60"{String, Bool, Module, String}; rethrow::Type, color::Bool, passthrough::Bool, capture_buffer::IOBuffer, io_context::Vector{Any})
│        @ IOCapture /usr/local/share/julia/packages/IOCapture/MR051/src/IOCapture.jl:167
│     [17] capture
│        @ /usr/local/share/julia/packages/IOCapture/MR051/src/IOCapture.jl:100 [inlined]
│     [18] execute_block(sb::Module, block::String; inputfile::String, fake_source::String, softscope::Bool, continue_on_error::Bool)
│        @ Literate /usr/local/share/julia/packages/Literate/On3EJ/src/Literate.jl:953
│     [19] execute_block
│        @ /usr/local/share/julia/packages/Literate/On3EJ/src/Literate.jl:935 [inlined]
│     [20] execute_markdown!(io::IOBuffer, sb::Module, block::String, outputdir::String; inputfile::String, fake_source::String, flavor::Literate.DocumenterFlavor, image_formats::Vector{Tuple{MIME, String}}, file_prefix::String, softscope::Bool, continue_on_error::Bool)
│        @ Literate /usr/local/share/julia/packages/Literate/On3EJ/src/Literate.jl:683
│     [21] (::Literate.var"#42#43"{Dict{String, Any}, IOBuffer, Module, Literate.CodeChunk, Int64})()
│        @ Literate /usr/local/share/julia/packages/Literate/On3EJ/src/Literate.jl:651
│     [22] cd(f::Literate.var"#42#43"{Dict{String, Any}, IOBuffer, Module, Literate.CodeChunk, Int64}, dir::String)
│        @ Base.Filesystem ./file.jl:112
│     [23] markdown(inputfile::String, outputdir::String; config::Dict{Any, Any}, kwargs::@Kwargs{flavor::Literate.DocumenterFlavor, preprocess::var"#2#3", execute::Bool})
│        @ Literate /usr/local/share/julia/packages/Literate/On3EJ/src/Literate.jl:650
│     [24] macro expansion
│        @ ./timing.jl:697 [inlined]
│     [25] top-level scope
│        @ /__w/Breeze.jl/Breeze.jl/docs/literate.jl:358
│     [26] include(mod::Module, _path::String)
│        @ Base ./Base.jl:306
│     [27] exec_options(opts::Base.JLOptions)
│        @ Base ./client.jl:317
│     [28] _start()
│        @ Base ./client.jl:550
└ @ Main.var"##277" /__w/Breeze.jl/Breeze.jl/docs/src/literated/acoustic_wave.md:12

Sensitivity visualization

The heatmap shows $\partial J / \partial u_i(x, z)$: positive values are wind perturbations that would increase surface acoustic energy, negative values would decrease it. Because $J$ integrates along the entire bottom, the pattern reveals which parts of the wind profile feed energy into the surface duct from anywhere along it.

sensitivity = Array(interior(du, :, 1, :))
abs_max     = maximum(abs, sensitivity)

fig_sens = Figure(size = (800, 350), fontsize = 12)
Label(fig_sens[0, :],
      "∂J / ∂uᵢ  (J = ⟨(ρ - ρ̄)²⟩ at surface,  t=$(prettytime(Nsteps * Δt))",
      fontsize = 14, tellwidth = false)
ax_sens = Axis(fig_sens[1, 1]; xlabel = "x (m)", ylabel = "z (m)")
hm = heatmap!(ax_sens, xs_u, zs, sensitivity; colormap = :balance,
              colorrange = (-abs_max, abs_max))
Colorbar(fig_sens[1, 2], hm; label = "∂J / ∂uᵢ")

fig_sens

Julia version and environment information

This example was executed with the following version of Julia:

using InteractiveUtils: versioninfo
versioninfo()

Julia Version 1.12.6
Commit 15346901f00 (2026-04-09 19:20 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × AMD EPYC 7R13 Processor
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, znver3)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_GPG = 3673DF529D9049477F76B37566E3C7DC03D6E495
  JULIA_LOAD_PATH = :@breeze
  JULIA_PKG_SERVER_REGISTRY_PREFERENCE = eager
  JULIA_VERSION = 1.12.6
  JULIA_DEPOT_PATH = /usr/local/share/julia:
  JULIA_PATH = /usr/local/julia
  JULIA_PROJECT = @breeze

These were the top-level packages installed in the environment:

import Pkg
Pkg.status()

Status `/__w/Breeze.jl/Breeze.jl/docs/Project.toml`
  [86bc3604] AtmosphericProfilesLibrary v0.1.7
  [660aa2fb] Breeze v0.5.3 `.`
⌅ [052768ef] CUDA v6.0.0
  [13f3f980] CairoMakie v0.15.11
⌅ [6a9e3e04] CloudMicrophysics v0.35.0
  [e30172f5] Documenter v1.17.0
  [daee34ce] DocumenterCitations v1.4.1
  [7da242da] Enzyme v0.13.154
  [46192b85] GPUArraysCore v0.2.0
  [63c18a36] KernelAbstractions v0.9.41
  [98b081ad] Literate v2.21.0
  [85f8d34a] NCDatasets v0.14.15
  [9e8cae18] Oceananigans v0.110.1
  [a01a1ee8] RRTMGP v0.21.9
  [3c362404] Reactant v0.2.264
  [b77e0a4c] InteractiveUtils v1.11.0
  [44cfe95a] Pkg v1.12.1
  [9a3f8284] Random v1.11.0
Info Packages marked with ⌅ have new versions available but compatibility constraints restrict them from upgrading. To see why use `status --outdated`

This page was generated using Literate.jl.