Compute reference density with derivative #185
Conversation
|
For the test, it may make sense to have a still atmosphere as an initial condition to see if any perturbation is generated during the initialization. |
| 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) |
There was a problem hiding this comment.
are we sure it's not
... - ℑzᵃᵃᶠ(i, j, k-1, grid, ρᵣ) * g * Δzᶜᶜᶠ(i, j, k-1, grid)?
There was a problem hiding this comment.
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?
There was a problem hiding this comment.
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
There was a problem hiding this comment.
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...
|
We can have a test to compare with analytic? |
well, there is a test that fails actually, which involves the potential temperature - energy equivalence. I also realized that it will be simpler to compute the density from pressure (not the other way around) by computing the derivative. What do you think about that? |
|
I'm bit puzzled for the use of computing numerically since there are analytic expressions... |
that's the idea. but we don't know if it matters. |
|
Gotcha... |
we'd have to read Pauluis 2008 more closely (we also should put in a total energy test) |
|
what do you think about my method using a derivative instead of integral? |
Yes, I'm not sure. Either methods use the analytical for p or ρ and compute the other numerically... But how do we pick which is better to be computed numerically vs analytically? |
|
I think the key property is that the vertical momentum equation is discretely satisfied (whereas it was not before). I believe this could enter into considerations about conservation of total discrete energy (although I am not sure, just speculating). If that is what matters, then you can achieve it either way. If there is another consideration, then I'm not sure. Using the analytical profiles might be fine too. |
glwagner
changed the title
|
The reason to do this has become clearer. If we want to use the reference state as a state with CompressibleDynamics, it is important that the reference state is in fact a discrete solution of the system. Thus, while it is not clear whether being a discrete solution matters for AnelasticDynamics, it may matter for CompressibleDynamics. So I would like to move forward with this. |
Use GradientBoundaryCondition for pressure at the top boundary to ensure the discrete pressure gradient correctly extrapolates into halo regions. This fixes a 50% error in density at the top cell. At the bottom boundary, keep ValueBoundaryCondition with the known surface pressure p₀ to ensure exact interpolation to the surface (required by atmosphere_model_construction tests). Also fixes a bug where surface_density was called with θ₀ (potential temperature) instead of the correct 4-argument form with pˢᵗ.
Reference state changes: - Use GradientBoundaryCondition at BOTH top and bottom boundaries for pressure to ensure correct discrete hydrostatic balance (ρ = -∂z(p)/g) regardless of whether the grid extends above or below z=0 - Remove stale ValueBoundaryCondition import Test updates: - diagnostics.jl: Use explicit z=(0, 1000) instead of extent=(100, 100, 1000) which was creating a grid with z ∈ [-1000, 0] and the surface at the top - atmosphere_model_construction.jl: Use rtol=1e-4 for surface pressure/density interpolation tests since GradientBoundaryCondition introduces small errors - dynamics.jl: Exclude Float32 from vertical momentum conservation test because integrated roundoff errors over the large domain (4e12 m³) make the test meaningless for Float32 (O(100) vs O(1e-7) for Float64)
Codecov Report✅ All modified and coverable lines are covered by tests. 📢 Thoughts on this report? Let us know! |
| z_bottom = first(z) | ||
| z_top = last(z) |
There was a problem hiding this comment.
These lines are accessing individual elements of z, which causes the "Scalar indexing is disallowed" error on GPU.
| # Uses rtol=1e-4 because GradientBoundaryCondition is used for discrete hydrostatic balance, | ||
| # which introduces small interpolation errors at the surface. |
There was a problem hiding this comment.
Thanks for the comment, it preempted my question! 😃
This PR computes the reference state density using a derivative, from the reference pressure, rather than evaluating the analytical formula. This ensures that
discretely exactly rather than just approximately.
Might make sense to polish this PR off with a test.