zhihua-zheng
commented
2 months ago On the second point, I think it is because KineticEnergyTendency does not include boundary tendency from prescribed surface fluxes?
Open zhihua-zheng opened 2 months ago
zhihua-zheng
commented
2 months ago On the second point, I think it is because KineticEnergyTendency does not include boundary tendency from prescribed surface fluxes?
tomchor
commented
2 months ago @zhihua-zheng thanks for opening this! Here are a couple of comments:
These 5 terms correspond to the
KineticEnergyTendency,AdvectionTerm,PressureRedistributionTerm,BuoyancyProductionTerm,KineticEnergyStressTermin Oceanostics'TKEBudgetTerms, respectively.
This is not important here, but just to clarify, the KETendency isn't exactly $\partial k / \partial t$. It is that, minus the contribution from the nonhydrostatic pressure. This is said in the docstring, but I just realized it's incomplete!:
In this case the integrated contributions from the nonhydrostatic pressure are zero, so once you integrate the tendency there's no difference.
- The nonzero integral of
AdvectionTermconflicts with what I imagined from the equation.
This indeed is weird. My first guess is that this might come from the momentum boundary condition, although I thought that was implemented through the stress term. (I'm actually having a hard time finding where exactly that BC is applied, but I'll look more into it later.)
- The
KineticEnergyTendencybeing mostly negative conflicts with the fact that the integral of KE grows with time.
That is also weird. The first thing I'd compare is the KineticEnergyTendency with the actual calculated $\partial k / \partial t$ from KineticEnergy. I'm guessing there might be a sign wrong somewhere...?
The way I'd move forward is to remove the buoyancy forcing, keep only the momentum forcing for simplicity, and add Coriolis. That way the BL depth will equilibrate (at $\sim u_\star / f$ if I recall correctly) and the whole system reaches a steady state. Then you can plot the relevant terms as a function of $z$. This should give us insight into where the advection contribution is coming from. Then we can deal with the tendency later.
Also a note: once we understand this better, I think it would be great to convert this code into an example of KE budgeting for the docs!
glwagner
commented
2 months ago Can you link to the code where these things are implemented? AdvectionTerm must somehow be a numerical approximation, so perhaps there is a discrete discrepany. Ultimately to get discrete correctness we have to derive every term from the discretely exact kinetic energy equation. This can be a little onerous though...
As a sanity check and to shed some light on whether discreteness is playing a role --- you should try repeating the problem with different time-steps (perhaps 10x smaller) and also on different resolution grids.
zhihua-zheng
commented
2 months ago Can you link to the code where these things are implemented?
It is using the velocity multiplied by the advective flux divergence kernel: https://github.com/tomchor/Oceanostics.jl/blob/ffa54aab5df7978aaf92368dd5db106599e591f5/src/TKEBudgetTerms.jl#L161C1-L166C4
zhihua-zheng
commented
2 months ago That is also weird. The first thing I'd compare is the
KineticEnergyTendencywith the actual calculated $\partial k/ \partial t$ fromKineticEnergy. I'm guessing there might be a sign wrong somewhere...?
They don't match, the blue line shows the KE grows with time.
zhihua-zheng
commented
2 months ago Following your suggestions, I have done two more experiments, both with finer stretched vertical grid (~ 0.5 m at surface), Coriolis and wind forcing only. Different advection schemes are used, and results are sensitive to this choice (sign of numerical dissipation?).
Below I show two figures for each experiment: evolution of vertically integrated horizontal mean budget term; time averaged horizontal mean of budget term as a function of depth. Note that different than the figure posted above, I have included the sign of each RHS term here to indicate the contribution to local KE tendency. The top most grid point is neglected due to outstandingly large value, and this appears to bring the KineticEnergyTendency inline with actual change of KineticEnergy with time.
The first one uses CenteredSecondOrder().
With WENO(order=5):
tomchor
commented
2 months ago Okay, just to make I'm understanding the main takeaways here:
CenteredSecondOrder() the main balance is between the advection term creating KE and the stress term removing it. With WENO() the balance seems similar, but the sign of the contributions is flipped, no?AdvectionTerm is unreasonably large and ignoring it improved results (specifically dk/dt matches the tendency term). @glwagner can this be related to the halos not being filled? I don't really understand very well how halo filling works in the code, but I can't think of anything else.Also can you pelase share the new code along with the code used to plot the figures? I wanna do some investigating of my own.
tomchor
commented
2 months ago It boggles me the most that the contribution from the AdvectionTerm isn't zero even far from the surface. Inspecting the code, I don't see anything wrong with it. Like @zhihua-zheng mentioned, this is the relevant kernel:
And that's being called from here: https://github.com/tomchor/Oceanostics.jl/blob/ffa54aab5df7978aaf92368dd5db106599e591f5/src/TKEBudgetTerms.jl#L196-L199
What the code is doing, using index notation, is
\begin{align}
u_i \partial_j(u_j u_i) &= \partial_j(u_j u_i^2) - u_j u_i \partial_j u_i\\
&= \partial_j(u_j u_i^2) - u_j \partial_j u_i^2/2\\
&= \partial_j(u_j u_i^2)/2 - u_i^2 \partial_j u_j/2
\end{align}I think this is highly unlikely, but is it possible that $\partial_j u_j$ is not zero? i.e. can the flow be erroneously divergent here?
But also, once you remove the first grid point, it seems that the budget is approximately closing, indicating that the diagnostics are (apart from the first point) quantifying things correctly... or at least with compensating errors...
zhihua-zheng
commented
2 months ago Also can you pelase share the new code along with the code used to plot the figures? I wanna do some investigating of my own.
Sure. Code for simulation:
using Oceananigans
using Oceananigans.Units: minutes, days
using Oceananigans.Advection: div_Uc
using Oceananigans.AbstractOperations: KernelFunctionOperation
using Oceanostics
using Oceanostics.TKEBudgetTerms: KineticEnergyTendency, AdvectionTerm, KineticEnergyStressTerm,
PressureRedistributionTerm, BuoyancyProductionTerm
function KineticEnergyAdvectionTerm(model::NonhydrostaticModel, ke; velocities = model.velocities, location = (Center, Center, Center))
return KernelFunctionOperation{Center, Center, Center}(div_Uc, model.grid, model.advection, velocities, ke)
end
const Nz = 64
const Lz = 50
const refinement = 1.5 # controls spacing near surface (higher means finer spaced)
const stretching = 10 # controls rate of stretching at bottom
h(k) = (k - 1) / Nz
ζ₀(k) = 1 + (h(k) - 1) / refinement
Σ(k) = (1 - exp(-stretching * h(k))) / (1 - exp(-stretching))
z_faces(k) = Lz * (ζ₀(k) * Σ(k) - 1)
grid = RectilinearGrid(GPU(), size=(32, 32, Nz), x=(0, 100), y=(0, 100), z=z_faces,
topology=(Periodic, Periodic, Bounded))
const hᵢ = 20
const Qᵘ = -2e-4
const Qᵇ = 0#1.3e-7
const N² = 6.4e-5 # s⁻², initial and bottom buoyancy gradient
b_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵇ),
bottom = GradientBoundaryCondition(N²))
u_boundary_conditions = FieldBoundaryConditions(top = FluxBoundaryCondition(Qᵘ))
model = NonhydrostaticModel(; grid, coriolis = FPlane(f=1e-4),
advection = WENO(),#CenteredSecondOrder(),#WENO(),
timestepper = :RungeKutta3,
tracers = :b,
buoyancy = BuoyancyTracer(),
closure = AnisotropicMinimumDissipation(),
boundary_conditions = (b=b_boundary_conditions, u=u_boundary_conditions))
@inline Ξ(z) = randn() * exp(z / 4)
@inline stratification(z) = z < - hᵢ ? N² * z : N² * (-hᵢ)
@inline bᵢ(x, y, z) = stratification(z) + 1e-1 * Ξ(z) * N² * model.grid.Lz
@inline uᵢ(x, y, z) = 1e-1 * Ξ(z)
@inline wᵢ(x, y, z) = 1e-1 * Ξ(z)
set!(model, u=uᵢ, w=wᵢ, b=bᵢ)
simulation = Simulation(model, Δt=5, stop_time=4days)
wizard = TimeStepWizard(cfl=0.65, diffusive_cfl=0.65)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(5))
# KE budget using resolved velocities
ke = Field(KineticEnergy(model))
KE = Field(Average(ke, dims=(1,2)))
KE_tendency = Field(Average(KineticEnergyTendency(model), dims=(1,2)))
KE_advection = Field(Average(AdvectionTerm(model), dims=(1,2)))
#KE_adv = Field(Average(KineticEnergyAdvectionTerm(model, ke), dims=(1,2)))
KE_stress = Field(Average(KineticEnergyStressTerm(model), dims=(1,2)))
KE_pressure = Field(Average(PressureRedistributionTerm(model), dims=(1,2)))
KE_buoyancy = Field(Average(BuoyancyProductionTerm(model), dims=(1,2)))
iKE = Field(Integral(KE))
iKE_tendency = Field(Integral(KE_tendency))
iKE_advection = Field(Integral(KE_advection))
#iKE_adv = Field(Integral(KE_adv))
iKE_stress = Field(Integral(KE_stress))
iKE_pressure = Field(Integral(KE_pressure))
iKE_buoyancy = Field(Integral(KE_buoyancy))
outputs = (; KE, KE_tendency, KE_advection, KE_stress, KE_pressure, KE_buoyancy,
iKE, iKE_tendency, iKE_advection, KE_stress, iKE_pressure, iKE_buoyancy)
output_interval = 5minutes
simulation.output_writers[:profile] =
NetCDFOutputWriter(model, outputs,
schedule = TimeInterval(output_interval),
filename = "mwe_KE_budget.nc",
overwrite_existing = true)
run!(simulation)Code for plotting in python:
import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
ds = xr.open_dataset('mwe_KE_budget.nc')
ds.close()
ds['timeTf'] = ds.time/np.timedelta64(int(np.around(2*np.pi/1e-4)), 's')
### Time series
plt.close()
plt.figure(figsize=(8,4), constrained_layout=True)
deltaz = np.diff(ds.zF)
ds['IKE_tendency'] = (ds.KE_tendency*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_advection'] = (ds.KE_advection*deltaz).isel(zC=slice(None,-1)).sum('zC')
#ds['IKE_adv'] = (ds.KE_adv*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_pressure'] = (ds.KE_pressure*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_buoyancy'] = (ds.KE_buoyancy*deltaz).isel(zC=slice(None,-1)).sum('zC')
ds['IKE_stress'] = (ds.KE_stress*deltaz).isel(zC=slice(None,-1)).sum('zC')
plt.plot(ds.timeTf, ds.IKE_tendency, 'k');
plt.plot(ds.timeTf, -ds.IKE_advection, 'C1');
plt.plot(ds.timeTf, -ds.IKE_pressure, 'C2');
plt.plot(ds.timeTf, ds.IKE_buoyancy, 'C3');
plt.plot(ds.timeTf, -ds.IKE_stress, 'C4');
plt.plot(ds.timeTf, (ds.IKE_tendency + ds.IKE_advection + ds.IKE_pressure - ds.IKE_buoyancy + ds.IKE_stress), ':C5', lw=1)
plt.legend(['tendency term', 'advection term', 'pressure term', 'buoyancy term', 'stress term', 'residual'],
frameon=False, loc='lower right', ncol=3)
#plt.plot(ds.timeTf, -ds.IKE_adv, ':b', lw=1);
plt.xlabel(r'$T_{inertial}$')
plt.ylabel('Integrated KE budgets')
plt.ylim(-3.6e-5, 4.6e-5)
plt.ticklabel_format(axis='y', scilimits=(0,0))
plt.xlim(0,5.5)
ax1 = plt.gca().twinx()
color = 'C0'
ax1.set_ylabel(r'$\int k dV$', color=color)
ax1.plot(ds.timeTf, ds.iKE, color=color)
ax1.set_yticks([0, 1, 2, 3, 4])
ax1.tick_params(axis='y', labelcolor=color)
ax1.set_ylim(-3.6, 4.6);
# plt.savefig('mwe_iKE_budgets', dpi=250)
### Profiles
plt.close()
plt.figure(figsize=(7.3,5), constrained_layout=True)
plt.subplot(121)
time_interval = (ds.timeTf >= 4) & (ds.timeTf < 4.5)
ke = ds.KE.where(time_interval, drop=True).isel(time=[0,-1]).T
dkedt = ke.diff('time')/ke.time.diff('time').dt.seconds
plt.plot(ds.KE_tendency.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'k')
plt.plot(dkedt, ds.zC, '--k')
plt.plot(-ds.KE_advection.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C1')
plt.plot(-ds.KE_pressure.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C2')
plt.plot(ds.KE_buoyancy.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C3')
plt.plot(-ds.KE_stress.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C4')
plt.plot((ds.KE_tendency + ds.KE_advection + ds.KE_pressure - ds.KE_buoyancy +
ds.KE_stress).where(time_interval).mean('time')[:-1], ds.zC[:-1], ':C5', lw=1)
#plt.plot(-ds.KE_adv.where(time_interval).mean('time')[:-1], ds.zC[:-1], ':b', lw=1)
plt.text(0.02,0.42, r'$T_{inertial} \in [4,4.5)$', fontsize=12, transform=plt.gca().transAxes)
plt.grid('on', ls='--', lw=0.4)
plt.xlim(-3e-6, 3e-6)
plt.ylim(-40,0)
plt.xlabel('Time averaged KE budget term')
ax1 = plt.gca().twiny()
color = 'C0'
ax1.set_xlabel(r'$k$', color=color)
ax1.plot(ds.KE.where(time_interval).mean('time'), ds.zC, color=color)
ax1.set_xticks(np.array([0, 0.5, 1, 1.5])*1e-2)
ax1.tick_params(axis='x', labelcolor=color)
plt.ticklabel_format(axis='x', scilimits=(0,0))
ax1.set_xlim(-0.019, 0.019)
plt.subplot(122)
time_interval = (ds.timeTf >= 4.5) & (ds.timeTf < 5)
ke = ds.KE.where(time_interval, drop=True).isel(time=[0,-1]).T
dkedt = ke.diff('time')/ke.time.diff('time').dt.seconds
plt.plot(ds.KE_tendency.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'k')
plt.plot(dkedt, ds.zC, '--k')
plt.plot(-ds.KE_advection.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C1')
plt.plot(-ds.KE_pressure.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C2')
plt.plot(ds.KE_buoyancy.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C3')
plt.plot(-ds.KE_stress.where(time_interval).mean('time')[:-1], ds.zC[:-1], 'C4')
plt.plot((ds.KE_tendency + ds.KE_advection + ds.KE_pressure - ds.KE_buoyancy +
ds.KE_stress).where(time_interval).mean('time')[:-1], ds.zC[:-1], ':C5', lw=1)
plt.legend(['tendency term', r'$\Delta k/\Delta t$', 'advection term', 'pressure term', 'buoyancy term',
'stress term', 'residual'], frameon=False, loc='lower left', borderpad=0.1)
#plt.plot(-ds.KE_adv.where(time_interval).mean('time')[:-1], ds.zC[:-1], ':b', lw=1)
plt.text(0.02,0.42, r'$T_{inertial} \in [4.5,5)$', fontsize=12, transform=plt.gca().transAxes)
plt.gca().set_yticklabels([])
plt.grid('on', ls='--', lw=0.4)
plt.xlim(-3e-6, 3e-6)
plt.ylim(-40,0)
plt.xlabel('Time averaged KE budget term')
ax2 = plt.gca().twiny()
color = 'C0'
ax2.set_xlabel(r'$k$', color=color)
ax2.plot(ds.KE.where(time_interval).mean('time'), ds.zC, color=color)
ax2.set_xticks(np.array([0, 0.5, 1, 1.5])*1e-2)
ax2.tick_params(axis='x', labelcolor=color)
plt.ticklabel_format(axis='x', scilimits=(0,0))
ax2.set_xlim(-0.019, 0.019);
# plt.savefig('mwe_KE_budgets', dpi=250)This is not discretely correct. The discrete kinetic energy equation is obtained by multiplying the momentum equation by
$$ \frac{u^{n} + u^{n+1}}{2} $$
To understand this consider the time-discretized "momentum" equation (considering only $u$ for simplicity)
$$ u^{n+1} - u^{n} = \Delta t \int_{tn}^{t{n+1}} G_u \mathrm{d} t $$
If you multiply the above equation by $u^n$ you're not going to get a kinetic energy equation. You have to multiply it by $(u^{n} + u^{n+1}) / 2$, which yields
$$ \frac{1}{2} \left ( u^{n+1} \right )^2 - \frac{1}{2} \left ( u^n \right )^2 = \frac{1}{2} (u^{n} + u^{n+1}) \Delta t \int_{tn}^{t{n+1}} G_u \mathrm{d} t $$
The formulas in the code don't incorporate the velocity field from two time-steps so they just can't be correct. It also doesn't look like it's consistent with the time-stepping scheme, which involves tendencies from both $n$ and time-step $n-1$... ?
What you can try to show to get around this is that the budget doesn't change when you change the time-step. This would indicate that the error you are making wrt temporal discretization is negligible. But it's definitely not going to be negligible in general. On the contrary, we usually use time-steps close to the stability boundary.
glwagner
commented
2 months ago @glwagner can this be related to the halos not being filled? I don't really understand very well how halo filling works in the code, but I can't think of anything else.
The halos are filled within update_state!:
which is called just before completing a time-step or substep.
Therefore whenever diagnostics are performed the halos for all prognostic variables should be filled.
zhihua-zheng
commented
2 months ago
- You observe qualitative differences when you change advection scheme. For example with the
CenteredSecondOrder()the main balance is between the advection term creating KE and the stress term removing it. WithWENO()the balance seems similar, but the sign of the contributions is flipped, no?
With WENO() it is more like the time-averaged contribution from AdvectionTerm is ~ 0, while KineticEnergyStressTerm is positive, reflecting the wind work wins over dissipation in the mean? In this case, dissipation is indeed much smaller than the one with CenteredSecondOrder().
Ideally just based on the KE equation, I would only expect KineticEnergyStressTerm as the sole driver for KE tendency. In other words, the partition between advection and stress term does not seem right to me in either case.
jhdong2016
commented
2 months ago Hi Zhihua,
I am also doing some energy budget analysis and also have trouble in closing the budget. Have you ever closed the surface momentum flux to check the balance? I found that the result just becomes much harder to understand if there is only sea surface buoyancy flux.
zhihua-zheng
commented
2 months ago Have you ever closed the surface momentum flux to check the balance? I found that the result just becomes much harder to understand if there is only sea surface buoyancy flux.
Haven't tried a purely buoyancy forced case for KE budget yet, but it's on my list. I have tried that for a TKE budget though, the dissipation is similarly smaller with WENO().
jhdong2016
commented
2 months ago Have you ever closed the surface momentum flux to check the balance? I found that the result just becomes much harder to understand if there is only sea surface buoyancy flux.
Haven't tried a purely buoyancy forced case for KE budget yet, but it's on my list. I have tried that for a TKE budget though, the dissipation is similarly smaller with
WENO().
I am Jihai Dong. Jacob and Baylor told me about you. Do you have wechat? I think it will be convenient for us to discuss. You can send it to my email, jihai_dong@nuist.edu.cn.
tomchor
commented
1 month ago This is not discretely correct. The discrete kinetic energy equation is obtained by multiplying the momentum equation by
un+un+12
Agreed, but I think for steady-state problems this error should be small, no? This is why I suggested working with a steady-state problem.
Also, even if the tendency error is not small, the fact that contributions from advective (more generally, divergent) terms should be zero should hold independently for each time step as long as the flow is divergence-free, no? The pressure contributions, for example, seem to be zero here, as expected.
glwagner
commented
1 month ago Oh of course if the solution isn't changing in time. But if its only statistically steady (and thus everything is pointwise changing in time) I'm not so sure that we get a magical time-mean cancellation of errors...
I'm thinking about the second point...
glwagner
commented
1 month ago A few more things...
On the mystery about the near-surface value and where momentum fluxes are added: if using a FluxBoundaryCondition, the momentum fluxes are added directly to the tendency:
As far as I can tell, this is currently missing from the way that the "total" tendency is computed. That is probably a significant source of discrepancy. Adding fluxes would be a little annoying. There is a workaround, however, which is to use a GradientBoundaryCondition instead --- in that case, any boundary fluxes do indeed enter through the stress term (and thus also involve diffusivities / viscosities). This is a little more work if you're using an SGS closure though (ie not constant diffusivity). But it can be done. I think there are also probably some source code changes we could make to streamline this way of forcing a simulation.
Note that one can also access the tendencies directly from model.timestepper, and those would have fluxes in them. They don't necessarily have to recomputed (maybe there is another reason for this, I am not sure and cannot remember...) Right now the tendencies are computed in update_state!, which means they should always be coincident in time with the velocity field.
As for this point:
Also, even if the tendency error is not small, the fact that contributions from advective (more generally, divergent) terms should be zero should hold independently for each time step as long as the flow is divergence-free, no?
I am not sure, given the discretization and high-order advection. Can we work this out at least for second-order advection to understand whether the discretization does not introduce anything new into the continuous identities?
There is a more serious issue as well, which is that the true, temporally correct kinetic energy equation involves the term
$$ u^{n+1}_i \partial_j \left ( u^n_i u^n_j \right ) $$
I don't think this term can be expressed as a pure divergence. So perhaps one really is forced to rely on the temporally approximate form of this budget. In that case, we should always be confirming results by halving the time-step...
zhihua-zheng
commented
1 month ago I can confirm that the volume integrated advection term goes to ~ 0 if the first grid point is used, so the flow is divergence free. The nonzero integral of advection in the initial post is likely due to the coarse resolution near surface.
zhihua-zheng
commented
1 month ago in that case, any boundary fluxes do indeed enter through the stress term (and thus also involve diffusivities / viscosities)
So if using a FluxBoundaryCondition for momentum, the wind work would not show up in the KineticEnergyStressTerm? KineticEnergyStressTerm is just velocity multiplied by stress divergence.
glwagner
commented
1 month ago I can confirm that the volume integrated advection term goes to ~ 0 if the first grid point is used, so the flow is divergence free. The nonzero integral of advection in the initial post is likely due to the coarse resolution near surface.
And just to clarify, I think this shows quite a bit more than just that the flow is divergence free. It also proves there is somehow a way to form the discrete equivalent of the continuous identities that make this term vanish, which is non-trivial because of the staggered grid.
glwagner
commented
1 month ago in that case, any boundary fluxes do indeed enter through the stress term (and thus also involve diffusivities / viscosities)
So if using a
FluxBoundaryConditionfor momentum, the wind work would not show up in theKineticEnergyStressTerm?KineticEnergyStressTermis just velocity multiplied by stress divergence.
I doubt it, but if show me the definition of this term, I can try to confirm. Presumably this term is associated with the turbulence closure / internal stresses. However, one can have a flux of momentum even with zero interior viscosity (eg closure=nothing for Oceananigans) in a finite volume system. So if the definitions are all consistent, then directly prescribed fluxes should not enter via closure-associated stresses.
If you instead prescribe the velocity gradient at the boundary, and also have viscosity, then there will indeed be a closure-associated stress across the surface and that would be included, probably.
zhihua-zheng
commented
1 month ago This is kernel for KineticEnergyStressTerm:
https://github.com/tomchor/Oceanostics.jl/blob/ffa54aab5df7978aaf92368dd5db106599e591f5/src/TKEBudgetTerms.jl#L297-L308
And it is being called from here:
https://github.com/tomchor/Oceanostics.jl/blob/ffa54aab5df7978aaf92368dd5db106599e591f5/src/TKEBudgetTerms.jl#L322-L335
zhihua-zheng
commented
1 month ago Let's say we don't care about calculating the full KE tendency for now, and the goal is to form a budget that explains the evolution of observed KE. I think it would help if I explicitly add contribution of top boundary flux to the stress term. I tried something like this:
@inline ψ∂₃Qtopᶠᶜᶜ(i, j, grid, ψ, args...) = @inbounds ψ[i, j, grid.Nz] * getbc(ψ.boundary_conditions.top, i, j, grid, args...) *
Azᶠᶜᶠ(i, j, grid.Nz+1, grid) / Vᶠᶜᶜ(i, j, grid.Nz, grid)
@inline ψ∂₃Qtopᶜᶠᶜ(i, j, grid, ψ, args...) = @inbounds ψ[i, j, grid.Nz] * getbc(ψ.boundary_conditions.top, i, j, grid, args...) *
Azᶜᶠᶠ(i, j, grid.Nz+1, grid) / Vᶜᶠᶜ(i, j, grid.Nz, grid)
@inline function uₕtop_stress_divergenceᶜᶜᶜ(i, j, grid, clock, model_fields)
u∂₃_τ₀ = ℑxᶜᵃᵃ(i, j, grid.Nz, grid, ψ∂₃Qtopᶠᶜᶜ, model_fields.u, clock, model_fields)
v∂₃_τ₀ = ℑyᵃᶜᵃ(i, j, grid.Nz, grid, ψ∂₃Qtopᶜᶠᶜ, model_fields.v, clock, model_fields)
return u∂₃_τ₀ + v∂₃_τ₀
end
function KineticEnergyBoundaryStressTerm(model::NonhydrostaticModel)
return KernelFunctionOperation{Center, Center, Nothing}(uₕtop_stress_divergenceᶜᶜᶜ, model.grid, model.clock, fields(model))
endObviously it doesn't work as there is probably many errors in this code snippet (e.g., I'm using KernelFunctionOperation to calculate a 2D field). I don't have much experience with this type of kernel calculation, can you help me to make this work?
zhihua-zheng
commented
1 month ago By adding the contribution from surface momentum flux offline, the integrated KE budget is almost closed.
glwagner
commented
1 month ago Hmm, yes, I think it is a good idea to compute the boundary term separately. You need to make sure you return 0 except in the upper grid point though, if it is going to be combined with other contributions and summed in a 3D kernel.
Another possibility is to start with an example that either doesn't have a boundary condition, or uses GradientBoundaryCondition with a constant viscosity (for example). That might be easier for debugging. There's a lot going on here.
I personally think the beginning point should be to use the tendency directly and confirm that we can recover dk/dt. That already has to confront (1) the decomposition of the pressure contribution from the tendency and (2) the inexact time-discretization of the diagnostics. Only after that test works (I'd also be interesting in characterizing the error made by the inexact time-discretization) do I think we can feel confident delving into the decomposition of the tendency. And once we do that, build things up one by one, considering incrementally more complex examples. For example, you can consider a problem that only has advection, with no closure or buoyancy.
tomchor
commented
1 month ago A few more things...
On the mystery about the near-surface value and where momentum fluxes are added: if using a
FluxBoundaryCondition, the momentum fluxes are added directly to the tendency:As far as I can tell, this is currently missing from the way that the "total" tendency is computed. That is probably a significant source of discrepancy. Adding fluxes would be a little annoying. There is a workaround, however, which is to use a
GradientBoundaryConditioninstead --- in that case, any boundary fluxes do indeed enter through the stress term (and thus also involve diffusivities / viscosities). This is a little more work if you're using an SGS closure though (ie not constant diffusivity). But it can be done. I think there are also probably some source code changes we could make to streamline this way of forcing a simulation.
This is interesting. Here's a question: is there a way to implement FluxBoundaryConditions so that they would always show up when calculating the kernel for the proper term? So for example, in this case they'd have to show up in the ∂ⱼ_τ₃ⱼ kernel. I think that'd be extremely useful for budgets, although I do recognize that it would be hard and might require some non-trivial code refactoring in Oceananigans.
tomchor
commented
1 month ago By adding the contribution from surface momentum flux offline, the integrated KE budget is almost closed.
This is looking pretty good. Do you think this is as precise as we can get given the approximations we've discussed here?
zhihua-zheng
commented
1 month ago Do you think this is as precise as we can get given the approximations we've discussed here?
No, the residual does shrink if I reduce the time step to 1 minute. I don't know when the error would converge.
tomchor
commented
1 month ago Okay, it's good that we see convergence with time then. Should we consider this issue solved?
If so, I'd like to propose that we turn this whole exercise into a docs example. This whole discussion was super helpful for me and I bet it'll be for other people as well since I think lots of people try to to close budgets in similar configurations. Also I think an example in docs would be best, instead of a discussion for example, since a docs example:
glwagner
commented
1 month ago Do you think this is as precise as we can get given the approximations we've discussed here?
No, the residual does shrink if I reduce the time step to 1 minute. I don't know when the error would converge.
Given the time discretization error, you will always observe improvements when you reduce the time-step.
glwagner
commented
1 month ago I know I'm a broken record but I really think it would be a service to the world to do the time discretization correctly. I fear the misunderstandings that will continue to spread...
zhihua-zheng
commented
1 month ago Do we know the discrete form of the KE equation? I don't think I know enough about the numerical details to do it correctly.
tomchor
commented
1 month ago I know I'm a broken record but I really think it would be a service to the world to do the time discretization correctly. I fear the misunderstandings that will continue to spread...
I'm okay with doing an approximation, which we can see here works well, and making it absolutely clear that it is an approximation, and explaining what are the approximations.
But I obviously also agree that it would be more valuable to do it properly with $u^n$ and $u^{n+1}$. I just think it'll be an order of magnitude more work and I'm guessing none of us has the time to dedicate to that. But maybe I'm wrong...? @glwagner do you feel the extra work would be minor compared to what @zhihua-zheng already did?
glwagner
commented
1 month ago Order of magnitude can't possibly be correct given the huge amount of effort that's already gone into this... can it? 😱
At the very least, as a starting point, the functions we're implementing here could have some kwarg like
cool_diagnostics(model; previous_velocities=model.velocities, ...)with a note that if not provided, then previous_velocities defaults to the current velocities --- and we're somehow extrapolating in time to make this computation, accurate within O(dt), hopefully...
This presumably would not require much work, at least nothing behind developing a precise understanding of the approximation that is being made (which is currently acknowledged as important). Then each user can decide whether they want to bother with a callback that stores the previous velocities (this seems quite easy but to each their own!)
Consider also the intellectual cost of trying to evaluate whether this approximation is causing a bug / problem with closing the budget every time you approach a new problem with this software package... in complex systems, every inch of rationality we can get is a huge victory I think.
It's important to note that the trade-offs are substantially different when one sits down to write a software package, versus writing something for a paper. In a paper, we get limited gains from investment because its only used one time and we can verify approximations on a one-off basis. But with software like this, our investments pay dividends and hopefully make us rich eventually.
This is meant to be a discussion. I need some help to interpret my KE ($
k = u_i u_i /2$) budget calculations, which follow this equation:These 5 terms correspond to the
KineticEnergyTendency,AdvectionTerm,PressureRedistributionTerm,BuoyancyProductionTerm,KineticEnergyStressTermin Oceanostics'TKEBudgetTerms, respectively.When integrated in volume for the entire domain, this becomes
For horizontally periodic, vertically bounded domain, the integrated
AdvectionTermandPressureRedistributionTermshould be zero.This is the code I used:
The result from the simulation is shown above. The residual is nearly 0. I mainly have two points of confusion:
AdvectionTermconflicts with what I imagined from the equation.KineticEnergyTendencybeing mostly negative conflicts with the fact that the integral of KE grows with time.I could have overlooked something here, please let me know what you think.