Implement a KinematicModel based on PrescribedVelocityFields #394

glwagner · 2026-01-08T06:00:04Z

This "kinematic driver" can be used to test microphysics schemes.

kaiyuan-cheng · 2026-01-08T13:28:25Z

The code appears to work as expected, except for the boundary conditions, see

The desired behavior is an upward translation of the profile. I think we may need an inlet boundary condition.

kaiyuan-cheng · 2026-01-08T13:51:23Z

src/AtmosphereModels/set_atmosphere_model.jl

+Set the velocity component `name` (`:u`, `:v`, or `:w`) to `value`.
+Also updates the corresponding momentum field.
+"""
+function set_velocity!(model::AtmosphereModel, name::Symbol, value)


Can we set velocity as a function of time? I believe so, but I’m not certain.

Yes! For time-dependent velocity fields, you need to build them explicitly in advance of the model with PrescribedVelocityFields and then pass them into the model:

model = AtmosphereModel(grid; velocities=prescribed_velocities)

Initially I was thinking taht we wouldn't support set! at all, but it is actually simply enough if the velocities are ordinary Field.

For boundary conditinons, do you think that we should accept boundary conditions on the velocity field? eg something like

w_bcs = FieldBoundaryConditions(bottom=OpenBoundaryCondition(w_inlet)) boundary_conditions = (; w = w_bcs) model = AtmosphereModel(grid; boundary_conditions)

I believe it should also be easy to have inlet boundary conditions with PrescribedVelocityFields (in that case no special thing is needed except defining functions w(z, t) that are non-zero at z=0)

Initially I was thinking taht we wouldn't support set! at all, but it is actually simply enough if the velocities are ordinary Field.

For boundary conditinons, do you think that we should accept boundary conditions on the velocity field? eg something like

w_bcs = FieldBoundaryConditions(bottom=OpenBoundaryCondition(w_inlet)) boundary_conditions = (; w = w_bcs) model = AtmosphereModel(grid; boundary_conditions)

Looks good to me!

codecov · 2026-01-08T14:54:28Z

Codecov Report

❌ Patch coverage is 77.60417% with 43 lines in your changes missing coverage. Please review.

Files with missing lines	Patch %	Lines
src/KinematicDriver/prescribed_dynamics.jl	78.78%	21 Missing ⚠️
.../KinematicDriver/kinematic_driver_time_stepping.jl	37.93%	18 Missing ⚠️
src/Thermodynamics/reference_states.jl	0.00%	4 Missing ⚠️

📢 Thoughts on this report? Let us know!

glwagner · 2026-01-08T15:20:11Z

The code appears to work as expected, except for the boundary conditions, see

The desired behavior is an upward translation of the profile. I think we may need an inlet boundary condition.

nice! can you paste your code that produced this in this PR?

kaiyuan-cheng · 2026-01-08T15:56:27Z

nice! can you paste your code that produced this in this PR?

using Breeze
using CairoMakie
using Printf

# Grid
Nz = 100
Lz = 20000
grid = RectilinearGrid(CPU(); size=Nz, x=0, y=0, z=(0, Lz),
                       topology=(Flat, Flat, Bounded))

# Reference state
constants = ThermodynamicConstants()
θ₀ = 300
p₀ = 100000
reference_state = ReferenceState(grid, constants;
                                 surface_pressure = p₀,
                                 potential_temperature = θ₀)

# Model setup
ρ_surface = reference_state.density[1, 1, 1]
ρθ_bcs = FieldBoundaryConditions(bottom=ValueBoundaryCondition(ρ_surface * θ₀))

dynamics = PrescribedDynamics(reference_state)
model = AtmosphereModel(grid; dynamics,
                        advection = UpwindBiased(order=1),
                        boundary_conditions = (ρθ=ρθ_bcs,),
                        thermodynamic_constants = constants)

# Initial conditions
θₜᵣ = 343
zₜᵣ = 12000
Tₜᵣ = 213
g = constants.gravitational_acceleration
cᵖᵈ = constants.dry_air.heat_capacity

function θᵇᵍ(z)
    θ_troposphere = θ₀ + (θₜᵣ - θ₀) * (z / zₜᵣ)^(5/4)
    θ_stratosphere = θₜᵣ * exp(g / (cᵖᵈ * Tₜᵣ) * (z - zₜᵣ))
    return ifelse(z <= zₜᵣ, θ_troposphere, θ_stratosphere)
end

set!(model; θ=θᵇᵍ, qᵗ=0, w=10.0)

# Simulation
simulation = Simulation(model; Δt=5.0, stop_time=1800, verbose=false)

times = Float64[]
θ_profiles = []

function record_profiles(sim)
    push!(times, time(sim))
    push!(θ_profiles, Array(interior(sim.model.formulation.potential_temperature, 1, 1, :)))
end
add_callback!(simulation, record_profiles, TimeInterval(300))

function progress(sim)
    θ = sim.model.formulation.potential_temperature
    @info "t = $(prettytime(sim)), θ_surface = $(@sprintf("%.1f", θ[1, 1, 1])) K"
end
add_callback!(simulation, progress, TimeInterval(300))

@info "Running kinematic advection test..."
run!(simulation)

# Visualization
z_km = znodes(grid, Center()) ./ 1000
fig = Figure(size=(500, 600))
ax_θ = Axis(fig[1, 1]; xlabel="θ (K)", ylabel="z (km)", title="Potential temperature")

colors = cgrad(:viridis, length(times), categorical=true)
for (n, t) in enumerate(times)
    lines!(ax_θ, θ_profiles[n], z_km; color=colors[n], label="t = $(Int(t/60)) min")
end

hlines!(ax_θ, [zₜᵣ/1000]; color=:gray, linestyle=:dash, label="Tropopause")
axislegend(ax_θ; position=:rb)
Label(fig[0, 1], "Kinematic Driver: dry atmosphere, w = 10 m/s"; fontsize=16, tellwidth=false)

fig

glwagner · 2026-01-09T23:17:40Z

So, after adding support for velocity boundary conditions we can in principle use the kinematic driver. Note however that flux form advection requires a velocity field that satisfied the anelastic continuity equation, eg we have to use something like

using Oceananigans.Operators: ℑzᵃᵃᶠ
ρ = reference_state.density
ρ₀ = ℑzᵃᵃᶠ(1, 1, 1, grid, ρ)
ρᵀ = ℑzᵃᵃᶠ(1, 1, grid.Nz+1, grid, ρ)
w = ZFaceField(grid)
W₀ = 10
ρᶠ = @at (Center, Center, Face) 1 * ρ
set!(w, ρ₀ * W₀ / ρᶠ)

we could use a different advection operator for the kinematic driver, although that would mean we can't use WENO.

Here's what I get

with

using Breeze
using CairoMakie
using Printf

# Grid
Nz = 100
Lz = 20000
grid = RectilinearGrid(CPU(); size=Nz, x=0, y=0, z=(0, Lz),
                       topology=(Flat, Flat, Bounded))

# Reference state
constants = ThermodynamicConstants()
θ₀ = 300
p₀ = 100000
reference_state = ReferenceState(grid, constants;
                                 surface_pressure = p₀,
                                 potential_temperature = θ₀)

using Oceananigans.Operators: ℑzᵃᵃᶠ
ρ = reference_state.density
ρ₀ = ℑzᵃᵃᶠ(1, 1, 1, grid, ρ)
ρᵀ = ℑzᵃᵃᶠ(1, 1, grid.Nz+1, grid, ρ)
w = ZFaceField(grid)
W₀ = 10
ρᶠ = @at (Center, Center, Face) 1 * ρ
set!(w, ρ₀ * W₀ / ρᶠ)

# Model setup
ρθ_bcs = FieldBoundaryConditions(bottom=ValueBoundaryCondition(ρ₀ * θ₀))

top_w = OpenBoundaryCondition(W₀ * ρ₀ / ρᵀ)
bottom_w = OpenBoundaryCondition(W₀)
w_bcs = FieldBoundaryConditions(bottom=bottom_w, top=top_w)

dynamics = PrescribedDynamics(reference_state)
model = AtmosphereModel(grid; dynamics,
                        advection = UpwindBiased(order=1),
                        boundary_conditions = (ρθ=ρθ_bcs, w=w_bcs),
                        thermodynamic_constants = constants)

# Initial conditions
θₜᵣ = 343
zₜᵣ = 12000
Tₜᵣ = 213
g = constants.gravitational_acceleration
cᵖᵈ = constants.dry_air.heat_capacity

function θᵇᵍ(z)
    θ_troposphere = θ₀ + (θₜᵣ - θ₀) * (z / zₜᵣ)^(5/4)
    θ_stratosphere = θₜᵣ * exp(g / (cᵖᵈ * Tₜᵣ) * (z - zₜᵣ))
    return ifelse(z <= zₜᵣ, θ_troposphere, θ_stratosphere)
end

set!(model; θ=θᵇᵍ, qᵗ=0, w=w)

# Simulation
simulation = Simulation(model; Δt=1, stop_time=180, verbose=false)

times = Float64[]
θ_profiles = []

function record_profiles(sim)
    push!(times, time(sim))
    push!(θ_profiles, Array(interior(sim.model.formulation.potential_temperature, 1, 1, :)))
end
add_callback!(simulation, record_profiles, TimeInterval(60))

function progress(sim)
    θ = sim.model.formulation.potential_temperature
    @info "t = $(prettytime(sim)), θ_surface = $(@sprintf("%.1f", θ[1, 1, 1])) K"
end
add_callback!(simulation, progress, TimeInterval(60))

@info "Running kinematic advection test..."
run!(simulation)

# Visualization
z_km = znodes(grid, Center()) ./ 1000
fig = Figure(size=(500, 600))
ax_θ = Axis(fig[1, 1]; xlabel="θ (K)", ylabel="z (km)", title="Potential temperature")

colors = cgrad(:viridis, length(times), categorical=true)
for (n, t) in enumerate(times)
    lines!(ax_θ, θ_profiles[n], z_km; color=colors[n], label="t = $(Int(t/60)) min")
end

hlines!(ax_θ, [zₜᵣ/1000]; color=:gray, linestyle=:dash, label="Tropopause")
axislegend(ax_θ; position=:rb)
Label(fig[0, 1], "Kinematic Driver: dry atmosphere, w = 10 m/s"; fontsize=16, tellwidth=false)

fig

kaiyuan-cheng · 2026-01-09T23:53:58Z

Great! We now have a working kinematic model. While WENO is desirable, I think this should be good for now.

glwagner · 2026-01-10T00:01:27Z

Great! We now have a working kinematic model. While WENO is desirable, I think this should be good for now.

@kaiyuan-cheng these results were generated by using a velocity field that satisfies the anelastic continuity equation.

I looked into it and I think there are two cases that we want:

PrescribedDensity and add the divergence term, so that advection looks like this:

$$ \partial_t ( \rho c ) = - \nabla \cdot \left ( \rho u c \right ) + c \nabla \cdot \left ( \rho u \right ) $$

Also evolve the density (compressible kinematic formulation)

The first case is needed for P3 validation I think.

Here's another reference for kinematic drivers:

https://github.com/Adehill/KiD-A/blob/master/docs/KiD_2.3.2625.pdf

note, our kinematic driver will support 1D, 2D, and 3D.

glwagner · 2026-01-10T00:02:14Z

@kaiyuan-cheng let me know if you agree! The divergence term is easy to add, and with that form we can still use WENO for $\nabla \cdot \left ( \rho u c \right )$

kaiyuan-cheng · 2026-01-10T01:53:12Z

@kaiyuan-cheng let me know if you agree! The divergence term is easy to add, and with that form we can still use WENO for ∇ ⋅ ( ρ u c )

The divergence term at the boundaries is necessary if we want to enforce mass continuity. Admittedly, I don’t understand why WENO works in this case but not in the previous one.

…gnece correction

glwagner · 2026-01-11T03:46:26Z

@kaiyuan-cheng let me know if you agree! The divergence term is easy to add, and with that form we can still use WENO for ∇ ⋅ ( ρ u c )

The divergence term at the boundaries is necessary if we want to enforce mass continuity. Admittedly, I don’t understand why WENO works in this case but not in the previous one.

Here's an explanation. The advection term is

$$ \rho u \cdot \nabla c $$

We don't have a WENO method to evaluate this term as is.

Also note that on a staggered grid, this term must be evaluated at cell centers. However, neither $u$ nor $\nabla c$ are at cell centers on a staggered grid.

However, we can rewrite this term because

$$ \nabla \cdot \left ( \rho u c \right ) = \rho u \cdot \nabla c + c \nabla \cdot \left ( \rho u \right ) $$

therefore

$$ \rho u \cdot \nabla c = \nabla \cdot \left ( \rho u c \right ) - c \nabla \cdot \left ( \rho u \right ) $$

The expression on the right has different numerical properties than the expression on the left. The term $\nabla \cdot \left ( \rho u c \right )$ is conservative (eg summing over the volume returns 0). Moreover, the divergence can be evaluated on cells, if the flux $\rho u c$ can be evaluated on the cell interfaces. This is natural on the C-grid, using WENO reconstruction for $c$.

The second term is $c \nabla \cdot \left ( \rho u \right )$. This too is naturally evaluated on the C-grid, because $c$ is already at the right location ans the mass divergence $\nabla \cdot \left ( \rho u \right )$ naturally. Note, we need to evaluate $\nabla \cdot \left ( \rho u \right )$ to step forward density with CompressibleDynamics or to compute the pressure correction with AnelasticDynamics. For that reason, we already have operators for it.

When we have CompressibleDynamics or AnelasticDynamics, this second term combines with the tendency term or vanishes (respectively).

However with PrescribedDynamics, we want to keep this second term.

kaiyuan-cheng · 2026-01-11T11:41:14Z

This is the desired result! Is WENO being used in this case?

src/KinematicDriver/prescribed_dynamics.jl

giordano · 2026-01-12T02:12:31Z

src/KinematicDriver/prescribed_dynamics.jl

+# Convenient method for letting the user specify only the divergence parameter
+# and infer the others.
+function PrescribedDynamics{Div}(density::D,
+                                 pressure::P,
+                                 surface_pressure::FT,
+                                 standard_pressure::FT) where {Div, D, P, FT}
+    return PrescribedDynamics{Div, D, P, FT}(density, pressure, surface_pressure, standard_pressure)
+end


I added this method for resolving the method error when trying to run

Breeze.jl/src/KinematicDriver/prescribed_dynamics.jl

Line 81 in 2ccb1fb

return PrescribedDynamics{divergence_correction}(density, pressure, p₀, pˢᵗ)

You expected to be able to call PrescribedDynamics with a single type parameter, but in reality it requires 3 more and there was no method to actually do it. I presume you wanted something like this (and tests pass for me now), but feel free to adapt it if you were thinking something else.

do we need a inner constructor then?

docs/make.jl

giordano · 2026-01-12T02:27:45Z

docs/make.jl

+example_pages = [] # Empty for fast builds
+


This is all green now, can this be removed? Do you want to add the new example to the list?

yes and yes

glwagner · 2026-01-13T05:42:34Z

This is the desired result! Is WENO being used in this case?

yes

Co-authored-by: Mosè Giordano <[email protected]>

examples/kinematic_driver.jl

test/kinematic_driver.jl

KinematicModel based on PrescribedVelocityFields

d90600f

glwagner requested a review from kaiyuan-cheng January 8, 2026 06:00

glwagner added 12 commits January 7, 2026 22:08

refactor

37ffad5

cleaner implementation

fef5e4a

cleaner design

f42faed

fix

e60b296

fix show

7412ffe

fix show

9453a0c

fix filenames

117b735

fix filenames

d1faab5

better errror

c0e152c

fixes

6eefe5b

fixes

c082155

fix

3d190c2

kaiyuan-cheng reviewed Jan 8, 2026

View reviewed changes

rm stale imports

1a68967

Merge branch 'main' into glw/kinematic-driver

bcfc332

support velocity boundary conditions with PrescribedDynamics

ffb3fff

kaiyuan-cheng approved these changes Jan 9, 2026

View reviewed changes

glwagner mentioned this pull request Jan 10, 2026

Implement PredictedParticlePropertyMicrophysics #395

Draft

revamp kinematic driver to support non-prescribed density, plus diver…

cbe2404

…gnece correction

glwagner added 2 commits January 10, 2026 20:58

add adapt_structure for PrescribedDensity

e069e96

delete some code

ed78fdd

giordano added 2 commits January 12, 2026 01:47

Add Adapt to test environment

af1126d

Add missing method

2ccb1fb

giordano reviewed Jan 12, 2026

View reviewed changes

src/KinematicDriver/prescribed_dynamics.jl Outdated Show resolved Hide resolved

giordano reviewed Jan 12, 2026

View reviewed changes

glwagner and others added 8 commits January 12, 2026 22:42

Update docs/make.jl

b304385

Co-authored-by: Mosè Giordano <[email protected]>

Update docs/make.jl

f34855d

Co-authored-by: Mosè Giordano <[email protected]>

Update docs/make.jl

c8d625e

Co-authored-by: Mosè Giordano <[email protected]>

Update src/KinematicDriver/prescribed_dynamics.jl

c63cf5e

Co-authored-by: Mosè Giordano <[email protected]>

Merge branch 'main' into glw/kinematic-driver

ede602b

Update KinematicDriver.jl

2c77331

add kinematic driver to make.jl

98cec51

export surface density

35d2ce5

giordano reviewed Jan 14, 2026

View reviewed changes

examples/kinematic_driver.jl Outdated Show resolved Hide resolved

glwagner added 2 commits January 13, 2026 18:21

rm reference to stationary parcel model

bc31766

Merge branch 'main' into glw/kinematic-driver

8080d3d

giordano reviewed Jan 14, 2026

View reviewed changes

examples/kinematic_driver.jl Outdated Show resolved Hide resolved

giordano and others added 2 commits January 14, 2026 01:26

Apply suggestion from @giordano

6311e0b

fix materialize dynamics

75d51b8

giordano reviewed Jan 14, 2026

View reviewed changes

test/kinematic_driver.jl Show resolved Hide resolved

fix constructor for AtmosphereModel

c2d1caf

glwagner merged commit 8f62cd7 into main Jan 14, 2026
10 of 11 checks passed

glwagner deleted the glw/kinematic-driver branch January 14, 2026 04:01

giordano added the breaking change 💔 Concerning a change which breaks the API label Jan 14, 2026

giordano mentioned this pull request Jan 14, 2026

Bump to 0.3.0 #403

Merged

Implement a KinematicModel based on PrescribedVelocityFields #394

Implement a KinematicModel based on PrescribedVelocityFields #394

Uh oh!

Conversation

glwagner commented Jan 8, 2026

Uh oh!

kaiyuan-cheng commented Jan 8, 2026

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!

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

Uh oh!

Codecov Report

Uh oh!

glwagner commented Jan 8, 2026

Uh oh!

kaiyuan-cheng commented Jan 8, 2026

Uh oh!

glwagner commented Jan 9, 2026

Uh oh!

kaiyuan-cheng commented Jan 9, 2026

Uh oh!

glwagner commented Jan 10, 2026

Uh oh!

glwagner commented Jan 10, 2026

Uh oh!

kaiyuan-cheng commented Jan 10, 2026

Uh oh!

glwagner commented Jan 11, 2026

Uh oh!

kaiyuan-cheng commented Jan 11, 2026

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Choose a reason for hiding this comment

Uh oh!

Choose a reason for hiding this comment

Uh oh!

glwagner commented Jan 13, 2026

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

codecov bot commented Jan 8, 2026 •

edited

Loading