How do I calculate the total pressure contribution in a turbulent kinetic energy budget using the Oceananigans kernels?

[tomchor]

I'm trying to calculate a kinetic energy budget for my nonhydrostatic simulation (that's part of what #3423 was about), and I figured the easiest way to do this is to use Oceananigans' own kernels. For most terms that's relatively simple. To calculate the contribution from the shear stress term, for example, I can do something like

u * ∂ⱼ_τ₁ⱼ(args...) + v * ∂ⱼ_τ₂ⱼ(args...) + w * ∂ⱼ_τ₃ⱼ(args...)

...(along with the proper interpolations, of course, as it is done in Oceanostics.)

However, I'm having a bit of trouble with the pressure term. I understand pressure is used to correct the velocity in order to enforce incompressibility, but that makes the calculation of its role (as it appears in the Navier-Stokes equations: $\nabla p$) kinda confusing to me.

If I understand correctly, when run!ning a simulation, the run!() function will call time_step!, which, for a Runge-Kutta substep as an example, looks like:

Oceananigans.jl/src/TimeSteppers/runge_kutta_3.jl

Lines 116 to 124 in 1c2a6f8

	rk3_substep!(model, Δt, γ², ζ²)
	calculate_pressure_correction!(model, second_stage_Δt)
	pressure_correct_velocities!(model, second_stage_Δt)
	tick!(model.clock, second_stage_Δt; stage=true)
	store_tendencies!(model)
	update_state!(model, callbacks)
	step_lagrangian_particles!(model, second_stage_Δt)

So after the substep (where we get the predictor velocities) the velocities are corrected using the nonhydrostatic pressure as such

Oceananigans.jl/src/Models/NonhydrostaticModels/pressure_correction.jl

Lines 31 to 37 in 1c2a6f8

	@kernel function _pressure_correct_velocities!(U, grid, Δt, pNHS)
	i, j, k = @index(Global, NTuple)
	@inbounds U.u[i, j, k] -= ∂xᶠᶜᶜ(i, j, k, grid, pNHS) * Δt
	@inbounds U.v[i, j, k] -= ∂yᶜᶠᶜ(i, j, k, grid, pNHS) * Δt
	@inbounds U.w[i, j, k] -= ∂zᶜᶜᶠ(i, j, k, grid, pNHS) * Δt
	end

Meaning the (nonhydrostatic) pressure gradient contribution is kinda hidden in the new velocities, and only then we get a call to update_state!(), where things that don't clearly relate to the nonhydrostatic pressure happen; for example the hydrostatic pressure will be updated.

Then there will be a call to compute_tendencies!(), which will eventually launch u_velocity_tendency! (and v, etc.). Taking the u case as an example you can clearly see the contributions from terms like advection, Coriolis, hydrostatic pressure greadient, etc., but there's no clear nonhydrostatic pressure gradient term:

Oceananigans.jl/src/Models/NonhydrostaticModels/nonhydrostatic_tendency_kernel_functions.jl

Lines 71 to 80 in 1c2a6f8

	return ( - div_𝐯u(i, j, k, grid, advection, total_velocities, velocities.u)
	- div_𝐯u(i, j, k, grid, advection, velocities, background_fields.velocities.u)
	- x_f_cross_U(i, j, k, grid, coriolis, velocities)
	- hydrostatic_pressure_gradient_x(i, j, k, grid, hydrostatic_pressure)
	- ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, diffusivities, clock, model_fields, buoyancy)
	- immersed_∂ⱼ_τ₁ⱼ(i, j, k, grid, velocities, u_immersed_bc, closure, diffusivities, clock, model_fields)
	+ x_curl_Uˢ_cross_U(i, j, k, grid, stokes_drift, velocities, clock.time)
	+ ∂t_uˢ(i, j, k, grid, stokes_drift, clock.time)
	+ x_dot_g_bᶠᶜᶜ(i, j, k, grid, buoyancy, tracers)
	+ forcings.u(i, j, k, grid, clock, model_fields))

So my question is, if I wanted to calculate (and "close") a kinetic energy budget accurately (really I want to calculate a turbulent kinetic energy budget, but let's start with the simple stuff) how would I calculate the contribution from the pressure gradient term?

CC @zhihua-zheng

[glwagner]

I believe that computing the exact, pointwise kinetic energy budget is challenging. This is equivalent to being able to compute the numerical dissipation (the residual of such an exact budget), which I am not sure is possible in general or has been done. But here are some references:

• https://agupubs.onlinelibrary.wiley.com/doi/epdf/10.1029/2022MS003380
• https://www.sciencedirect.com/science/article/abs/pii/S1463500314000754

Note that I think your shear production term is not exact, but approximate. To be exact, I think it would have to combine information at two different time-steps (because that is necessary to derive an equation for the discrete variance). If you're seeking an approximate budget that respects the C-grid but possibly contains temporal approximations, then I would use the continuous analogue of the pressure contribution to kinetic energy, u . grad(p).

@simone-silvestri might have more to say because he has been working on quantifying numerical dissipation.

[tomchor]

Note that I think your shear production term is not exact, but approximate. To be exact, I think it would have to combine information at two different time-steps (because that is necessary to derive an equation for the discrete variance).

Yes, agreed. It's okay for my purposes that the budget is only approximate in time, but I'd like it to be as exact as possible in space (like you said, respecting the C-grid relations, preferably in conservative form). I'm aware it will only be as accurate as the numerical/spurious mixing will allow (as your references illustrate), but I'm okay with that.

then I would use the continuous analogue of the pressure contribution to kinetic energy, u . grad(p).

So you're suggesting using an AbstractOperation like u * ∂x(sum(model.pressures)) + v * ∂y(...)...?

My worry about that approach is the treatment of the boundaries, since I suspect it can non-negligibly impact the form drag estimation, although I don't fully understand the pressure correction well enough to be sure. Is there somewhere in the docs where we go into that in a bit more detail? The only thing I could find is this section. I also skimmed Schumann and Sweet (1988), but it appears to cover only the numerics, not how the NS equations boundary conditions map to the Poisson equation.

[simone-silvestri]

You can treat pressure like any other term, there should be no ambiguity contingent on the fact you use the pressure at n+1. In fact, the logical steps that we are taking are:

uⁿ⁺¹  = uⁿ +  Δt Gⁿ
pⁿ⁺¹  = poisson_solution(∇⋅uⁿ⁺¹...)
uⁿ⁺¹  = uⁿ⁺¹ + Δt ∂ₓpⁿ⁺¹

so the non-hydrostatic pressure at n+1 is a part of the tendency as much as all the other terms. Also, as Greg said, you need the information at n+1 to calculate the budgets because the correct way to assess the discrete kinetic energy is to do

(uⁿ⁺¹ + uⁿ) Gⁿ / 2

For the pressure this means that the discrete contribution of the pressure term to the discrete kinetic energy is

(uⁿ⁺¹ + uⁿ) ∂ₓpⁿ⁺¹ / 2

If you want to be extremely precise you have to account for the time-stepping scheme used.

[tomchor]

Thanks for the clarification. Being precise in time isn't very important to me to be honest. What I was having a hard time understanding is the fact that if I add:

- div_𝐯u(i, j, k, grid, advection, total_velocities, velocities.u) 
- div_𝐯u(i, j, k, grid, advection, velocities, background_fields.velocities.u) 
- x_f_cross_U(i, j, k, grid, coriolis, velocities) 
- hydrostatic_pressure_gradient_x(i, j, k, grid, hydrostatic_pressure) 
- ∂ⱼ_τ₁ⱼ(i, j, k, grid, closure, diffusivities, clock, model_fields, buoyancy) 
- immersed_∂ⱼ_τ₁ⱼ(i, j, k, grid, velocities, u_immersed_bc, closure, diffusivities, clock, model_fields) 
+ x_curl_Uˢ_cross_U(i, j, k, grid, stokes_drift, velocities, clock.time) 
+ ∂t_uˢ(i, j, k, grid, stokes_drift, clock.time) 
+ x_dot_g_bᶠᶜᶜ(i, j, k, grid, buoyancy, tracers) 
+ forcings.u(i, j, k, grid, clock, model_fields)

PLUS add to that ∂x(model.pressures.pNHS), then the result (obviously) won't equal u_velocity_tendency(), which I was associating with ∂t(u).

But what you're saying is that u_velocity_tendency() + ∂x(model.pressures.pNHS) is what should approximate ∂t(u), not u_velocity_tendency(), correct? And that comes from the fact that in time-stepping we correct the velocities using ∂x(model.pressures.pNHS) which it's "hidden" in the calculated velocities.

Like you said, there's some nuance regarding which time is used for what (n+1 or n), but in an approximated sense you can say that ∂t(u) ≈ u_velocity_tendency() + ∂x(model.pressures.pNHS), correct?

[simone-silvestri]

depends, if you want to have to close the budget you probably have to be precise with the temporal discretization. At least this is what we noticed when trying to close tracer variance budgets in the hydrostatic model. In terms of correctness, we noticed that there was quite a difference using n or the correct 0.5 * ((n+1) + (n)) variable. If you just want a qualitative answer then maybe you can make do with quantities at n

[tomchor]

Thanks @simone-silvestri this clarifies a lot and I think it answers my question :)