Skip to content

Compute reference density with derivative #185

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
glwagner wants to merge 16 commits into
base: main
Choose a base branch
from
+52 −30
Open

Compute reference density with derivative #185

Changes from 1 commit
Commits
Show all changes
16 commits
Select commit Hold shift + click to select a range
6d3678c
compute reference pressure with integral
glwagner Nov 21, 2025
d33613f
compute density from pressure
glwagner Nov 21, 2025
312e3b3
clean up?
glwagner Nov 21, 2025
c2a44ba
Merge branch 'main' into glw/ref-integral
glwagner Nov 25, 2025
df38360
Merge branch 'main' into glw/ref-integral
glwagner Jan 8, 2026
0a49c6a
Add partial derivative operator to imports
glwagner Jan 8, 2026
bf9dac5
Merge branch 'main' into glw/ref-integral
glwagner Jan 8, 2026
fae0a01
Apply suggestions from code review
giordano Jan 12, 2026
a09cbcf
Merge branch 'main' into glw/ref-integral
giordano Jan 12, 2026
b0ecdd5
Merge branch 'main' into glw/ref-integral
glwagner Jan 19, 2026
7972fca
rm whitespace
glwagner Jan 19, 2026
1e3feca
Merge branch 'main' into glw/ref-integral
glwagner Jan 19, 2026
b2c781e
Fix discrete hydrostatic balance at top boundary in reference state
glwagner Jan 19, 2026
31542cc
Fix discrete hydrostatic balance in reference state and update tests
glwagner Jan 19, 2026
4c6ac15
Fix discrete hydrostatic balance in reference state
glwagner Jan 19, 2026
cec296b
Resolve merge conflicts: use GradientBoundaryCondition at both bounda…
glwagner Jan 19, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
33 changes: 30 additions & 3 deletions src/Thermodynamics/reference_states.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,5 +1,9 @@
using Oceananigans: Oceananigans, Center, Field, set!, fill_halo_regions!
using Oceananigans.Utils: launch!
using Oceananigans.Operators: Δzᶜᶜᶠ, ℑzᵃᵃᶠ

using Adapt: Adapt, adapt
using KernelAbstractions: @kernel, @index

#####
##### Reference state computations for Boussinesq and Anelastic models
Expand Down Expand Up @@ -91,12 +95,35 @@ function ReferenceState(grid, thermo=ThermodynamicConstants(eltype(grid));
p₀ = convert(FT, base_pressure)
θ₀ = convert(FT, potential_temperature)

pᵣ = Field{Nothing, Nothing, Center}(grid)
ρᵣ = Field{Nothing, Nothing, Center}(grid)
set!(pᵣ, z -> adiabatic_hydrostatic_pressure(z, p₀, θ₀, thermo))
set!(ρᵣ, z -> adiabatic_hydrostatic_density(z, p₀, θ₀, thermo))
fill_halo_regions!(pᵣ)
fill_halo_regions!(ρᵣ)

pᵣ = Field{Nothing, Nothing, Center}(grid)
compute_reference_pressure!(pᵣ, grid, p₀, ρᵣ, thermo)
# set!(pᵣ, z -> adiabatic_hydrostatic_pressure(z, p₀, θ₀, thermo))
fill_halo_regions!(pᵣ)

return ReferenceState(p₀, θ₀, pᵣ, ρᵣ)
end

@kernel function _compute_reference_pressure!(pᵣ, grid, p₀, ρᵣ, thermo)
i, j = @index(Global, NTuple)

# ∂z p = - ρᵣ g
# (p⁺ - pᵏ) / Δz = - ℑzᶠ(ρᵣ) * g
# p⁺ = pᵏ - ℑzᶠ(ρᵣ) * g Δz

@inbounds pᵣ[i, j, 1] = p₀
g = thermo.gravitational_acceleration

for k = 2:grid.Nz
@inbounds pᵣ[i, j, k] = pᵣ[i, j, k-1] - ℑzᵃᵃᶠ(i, j, k, grid, ρᵣ) * g * Δzᶜᶜᶠ(i, j, k, grid)
Copy link
Member

@navidcy navidcy Nov 21, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

are we sure it's not

... - ℑzᵃᵃᶠ(i, j, k-1, grid, ρᵣ) * g * Δzᶜᶜᶠ(i, j, k-1, grid)

?

Copy link
Member

@navidcy navidcy Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We want to interpolate the density on the face between cell center k-1 and k and add it to the integral, right? But that's face k-1, no?

Copy link
Member

@navidcy navidcy Nov 21, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more. 

using Oceananigans
using CairoMakie
using Oceananigans.Operators: Δzᶜᶜᶠ, ℑzᵃᵃᶠ

grid = RectilinearGrid(size=20, z=(0, 5), topology=(Flat, Flat, Bounded))

ρ = CenterField(grid)

p₁ = CenterField(grid)
p₂ = CenterField(grid)
p_analytic = CenterField(grid)

p₀ = 123.2

set!(ρ, z -> exp(-2z))

# dp/dz = - ρ, p(0) = p₀
set!(p_analytic, z -> p₀ + 1/2 * (exp(-2z) - 1))

for f in (p₁, p₂, p_analytic)
    f[1] = p₀
end

i, j = 1, 1
for k = 2:grid.Nz
    @inbounds p₁[i, j, k] = p₁[i, j, k-1] - ℑzᵃᵃᶠ(i, j, k, grid, ρ) * Δzᶜᶜᶠ(i, j, k, grid)
    @inbounds p₂[i, j, k] = p₂[i, j, k-1] - ℑzᵃᵃᶠ(i, j, k-1, grid, ρ) * Δzᶜᶜᶠ(i, j, k-1, grid)
end

fig = Figure()
axρ = Axis(fig[1, 1], ylabel="ρ(z)")
axp = Axis(fig[1, 2], ylabel="p(z)")
axe = Axis(fig[1, 3], ylabel="error")

lines!(axρ, ρ)

scatterlines!(axp, p₁)
scatterlines!(axp, p₂)
scatterlines!(axp, p_analytic)

scatterlines!(axe, abs(p_analytic - p₁), label="|p_analytic - p₁|")
scatterlines!(axe, abs(p_analytic - p₂), label="|p_analytic - p₂|")
axislegend(axe, position=:rt)

fig
display

Copy link
Member

@navidcy navidcy Nov 21, 2025 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Not sure which is better; given that it's only derivatives of pressure that matter the blue might be better since it's more "parallel" to the analytic? If abs value matters, seems like yellow is better as it doesn't have that offset error...

Copy link
Member Author

@glwagner glwagner Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

oh you are right!

Copy link
Member Author

@glwagner glwagner Nov 21, 2025

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

pretty nice demo btw

end
end

function compute_reference_pressure!(pᵣ, grid, p₀, ρᵣ, thermo)
arch = grid.architecture
launch!(arch, grid, (1, 1), _compute_reference_pressure!, pᵣ, grid, p₀, ρᵣ, thermo)
return nothing
end
Loading