Computed kinetic energy dissipation rate seems too small

[SamDeAbreu]

Hi everyone!

I'm currently working on a GPU simulation of the ice-ocean boundary layer, and in my analysis I am computing all the terms in the turbulent kinetic energy budget:

where nu_e is the eddy viscosity (denote the dissipation term as epsilon hereafter). My Reynolds averaging is just a horizontal spatial average. I've also neglected the TKE advection as it is small. I'm also computing the energy budget:

where the LHS is the sum of the change in kinetic and potential energy as well as the dissipation (plus the change in TKE, but this is small) and the RHS is the wind work (tau is the imparted surface stress, in this case just u_star_0^2 in the x direction).

Now I have two issues. First, the TKE budget doesn't close (see below for an example simulation with nonzero residual averaged over half an inertial period). I expect near the surface that the shear production would be balanced by the dissipation, but it seems that the shear production is more than two times larger. Note that terms seem to sort of blow up near the surface, but I assume that is due to insufficient resolution there.

The second is that the second energy budget doesn't close (see below, the same simulation as above).

There is a strange spike roughly before 0.1 in the second energy budget, but that is associated with spin-up. Interestingly, if I multiply my dissipation by ~2.8, it seems to roughly close both budgets. I was wondering if anyone had any ideas as to why this is happening, or if someone could spot a stupid mistake I am making. Below is the simplified code that reproduces this issue (very similar to the ocean mixing example):

#========================= Imports =======================#
using Random
using Printf
using Oceananigans
using Oceananigans.Units: minute, minutes, hour
using NPZ
using Adapt, CUDA
using Statistics
using Oceanostics: IsotropicKineticEnergyDissipationRate, TurbulentKineticEnergy, ZShearProductionRate, KineticEnergyDissipationRate

#========================= Constants =======================#
# Simulation config
Nx = 64          # number of points in x direction
Ny = 64          # number of points in y direction
Nz = 64          # number of points in the vertical direction

Lx = Ly = 64     # (m) domain horizontal extents
Lz = 50          # (m) domain depth
duration = 12hour

# Sponge Layer constants
const tau = 60 # Damping timescale (s)
const sigma = 0.05*Lz # Sponge layer extent

# Forcing constants
const u_star_0 = 0.028

const f = 4*π/(60*60*24)*sin(75*π/180) # Coriolis parameter (1/s)
const alpha = 1.00834e-5 # Thermal expansion coefficient (1/K)
const beta = 7.93e-4 # Saline Contraction coefficient (1/psu)

# Read files
#direc = ARGS[1]
#filename = ARGS[2]
direc = "."
filename = "energy_reynolds17_git"
b = npzread(direc*"//conditions.npz")
time_steps = CuArray(b["time_steps"])
#========================= Stretching Functions and Grid =======================#
grid = RectilinearGrid(GPU(), size = (Nx, Ny, Nz), x = (0, Lx), y = (0, Ly), z = (-Lz, 0))

#========================= Boundary Conditions =======================#

# Formulate boundary conditions
u_bcs = FieldBoundaryConditions(top = FluxBoundaryCondition(-u_star_0^2 ))

#========================= Forcing =======================#
bottom_mask = GaussianMask{:z}(center=-grid.Lz, width=sigma) # Mask that decays as a Gaussian from the bottom
uvw_sponge = Relaxation(rate=1/tau, mask=bottom_mask) # Sponge to dampen internal wave reflections

#========================= Model =======================#
buoyancy = SeawaterBuoyancy(equation_of_state=LinearEquationOfState(thermal_expansion=alpha, haline_contraction=beta))
model = NonhydrostaticModel(; grid, buoyancy,
                            advection = WENO(grid; order=9),
                            timestepper = :RungeKutta3,
                            tracers = (:T, :S),
                            coriolis = FPlane(f=f),
                            closure = SmagorinskyLilly(),
                            boundary_conditions = (u=u_bcs,),
                            forcing=(u=uvw_sponge, v=uvw_sponge, w=uvw_sponge))

#========================= Initial Conditions =======================#
@inline Ξ(z) = randn() * z / model.grid.Lz * (1 + z / model.grid.Lz) # Random noise damped at top and bottom
@inline w_int(x, y, z) = 1e-3 * u_star_0 * Ξ(z) # Initial vertical velocity is random noise, velocity scale comes from the 2% rule of thumb 

# Import initial conditions (vertical profiles)
z_values_S = b["S"]
z_values_T = b["T"]
z_values_u = b["u"]
z_values_v = b["v"] 
# Construct arrays with initial conditions
arr_S = zeros(Float64, Nx, Ny, Nz)
arr_T = zeros(Float64, Nx, Ny, Nz)
arr_u = zeros(Float64, Nx, Ny, Nz)
arr_v = zeros(Float64, Nx, Ny, Nz)
for k in 1:Nz
    arr_S[:, :, k] .= z_values_S[k]
    arr_T[:, :, k] .= z_values_T[k] + 273
    arr_u[:, :, k] .= z_values_u[k]
    arr_v[:, :, k] .= z_values_v[k]
end

#set!(model, u=CuArray(arr_u), v=CuArray(arr_v), w=w_int, S=CuArray(arr_S), T=CuArray(arr_T))
set!(model, u=w_int, v=w_int, w=w_int, S=CuArray(arr_S), T=CuArray(arr_T))

#========================= Simulation Init =======================#
simulation = Simulation(model, Δt=3.0, stop_time=duration)

wizard = TimeStepWizard(cfl=1.0, max_change=1.1, max_Δt=1minute)
simulation.callbacks[:wizard] = Callback(wizard, IterationInterval(10))

# Print a progress message
progress_message(sim) = @printf("Iteration: %04d, time: %s, Δt: %s, max(|u|) = %.1e ms⁻¹, max(|w|) = %.1e ms⁻¹, wall time: %s\n",

                                iteration(sim), prettytime(sim), prettytime(sim.Δt),
                                maximum(abs, sim.model.velocities.u), maximum(abs, sim.model.velocities.w), prettytime(sim.run_wall_time))

# Add callbacks
add_callback!(simulation, progress_message, IterationInterval(20))

#========================= Diagnostics =======================#
# Compute the TKE budget terms
# Actual velocities
u = @at (Center, Center, Center) model.velocities.u
v = @at (Center, Center, Center) model.velocities.v
w = @at (Center, Center, Center) model.velocities.w
# Horizontally-Mean velocities
U = Field(Average(u, dims=(1,2)))
V = Field(Average(v, dims=(1,2)))
W = Field(Average(w, dims=(1,2)))
# Primed velocities (fluctuations)
u_prime = @at (Center, Center, Center) Field(u - U)
v_prime = @at (Center, Center, Center) Field(v - V)
w_prime = @at (Center, Center, Center) Field(w - W)
# Pressure calculations (only use nonhydrostatic pressure as primed hydrostatic is 0 for this advection scheme)
p = model.pressures.pNHS
P = Field(Average(p, dims=(1, 2)))
p_prime = @at (Center, Center, Center) (p-P)

# TKE 
TKE_pointwise = Field(0.5*(u_prime^2 + v_prime^2 + w_prime^2))
TKE = Average(TKE_pointwise, dims=(1,2))

# TKE Advection
TKE_adv_1 = Field(U*Field(Average(∂x(TKE_pointwise), dims=(1,2))))
TKE_adv_2 = Field(V*Field(Average(∂y(TKE_pointwise), dims=(1,2))))
TKE_adv_3 = Field(W*Field(Average(∂z(TKE_pointwise), dims=(1,2))))
TKE_adv = TKE_adv_1 + TKE_adv_2 + TKE_adv_3

# Shear Production (using Oceanostics)
uw = Field(Average(u_prime*w_prime, dims=(1,2)))
vw = Field(Average(v_prime*w_prime, dims=(1,2)))
ww = Field(Average(w_prime*w_prime, dims=(1,2)))
U_shear = Field(∂z(U))
V_shear = Field(∂z(V))
W_shear = Field(∂z(W))
P_s = Average(ZShearProductionRate(model, u, v, w, U, V, W), dims=(1,2))

# Buoyancy Flux
buoyancy_field = 9.8*(beta*model.tracers.S - alpha*model.tracers.T)
b_prime = buoyancy_field - Field(Average(buoyancy_field, dims=(1,2)))
B_prod = Field(Average(-b_prime*w_prime, dims=(1,2))) 

# Viscous Dissipation (includes both SGS and molecular contributions)
ν = Field(model.diffusivity_fields.νₑ)
S_11 = Field(∂x(u_prime)^2)
S_22 = Field(∂y(v_prime)^2)
S_33 = Field(∂z(w_prime)^2)
S_12_1 = Field(∂y(u_prime) + ∂x(v_prime))
S_12 = Field(0.25*S_12_1^2)
S_13_1 = Field(∂z(u_prime) + ∂x(w_prime))
S_13 = Field(0.25*S_13_1^2)
S_23_1 = Field((∂y(w_prime) + ∂z(v_prime)))
S_23 = Field(0.25*S_23_1^2)
term1 = Field(S_11 + S_22 + S_33)
term2 = Field(S_12 + S_13 + S_23)
epsilon = Field(Average(2*ν*(term1 + 2*term2), dims=(1,2)))
epsilon2 = Field(Average(IsotropicKineticEnergyDissipationRate(model; U, V, W), dims=(1,2)))

# Transport terms
# Turbulent transport
trans_1_1 = Field(TKE_pointwise*w_prime)
trans_1_d = @at (Center, Center, Center) Field(-∂z(trans_1_1))
turbulent_transport = Field(Average(trans_1_d, dims=(1,2)))

# Pressure diffusion
trans_2 = Field(p_prime*w_prime)
trans_2_d = @at (Center, Center, Center) Field(-∂z(trans_2))
pressure_transport = Field(Average(trans_2_d, dims=(1,2)))

# Viscous transport (includes both molecular and SGS contributions)
xu =  @at (Center, Center, Center) Field(∂z(u_prime))
yv =  @at (Center, Center, Center) Field(∂y(v_prime))
zw =  @at (Center, Center, Center) Field(∂z(w_prime))
yu =  @at (Center, Center, Center) Field(∂y(u_prime))
xv =  @at (Center, Center, Center) Field(∂x(v_prime))
xw =  @at (Center, Center, Center) Field(∂x(w_prime))
zu =  @at (Center, Center, Center) Field(∂z(u_prime))
zv =  @at (Center, Center, Center) Field(∂z(v_prime))
yw =  @at (Center, Center, Center) Field(∂y(w_prime))
tau13_SGS = Field(ν*(zu+xw)) 
tau23_SGS = Field(ν*(yw+zv))
tau33_SGS = Field(ν*(zw+zw))
trans_3 = Field(-tau13_SGS*u_prime - tau23_SGS*v_prime - tau33_SGS*
w_prime)
trans_3_d = @at (Center, Center, Center) Field(-∂z(trans_3))
viscous_transport = Field(Average(trans_3_d,dims=(1,2)))

# Other diagnostics
S_mean = Field(Average(model.tracers.S, dims=(1,2)))
T_mean = Field(Average(model.tracers.T, dims=(1,2)))
T_mean_C = Field(T_mean - 273)
Tz = Field(∂z(model.tracers.T))
turbulent_heat_flux = Field(Average(Tz*model.diffusivity_fields.νₑ, dims=(1,2)))
T_prime = @at (Center, Center, Center) Field(model.tracers.T - T_mean)
S_prime = @at (Center, Center, Center) Field(model.tracers.S - S_mean)
heat_full_flux = Field(Average(w*model.tracers.T, dims=(1,2)))
salt_full_flux = Field(Average(w*model.tracers.S, dims=(1,2)))
heat_flux = Field(Average(T_prime*w_prime, dims=(1,2)))
salt_flux = Field(Average(S_prime*w_prime, dims=(1,2)))
eddy_diffusivity = Field(Average(model.diffusivity_fields.νₑ, dims=(1,2)))
IW_flux = Field(Average(p_prime*w_prime/1024, dims=(1,2)))

#========================= Simulation Output =======================#
xC, yC, zC = xnodes(grid, Center()), ynodes(grid, Center()), znodes(grid, Center())

S_mean_out = (; S_mean=S_mean)
T_mean_out = (; T_mean=T_mean_C)
U_out = (; U=(u^2 + v^2 + w^2)^0.5)

output = Dict("U"=>U, "V"=>V, "W"=>W, "heat_flux"=>heat_flux, "turbulent_heat_flux"=>turbulent_heat_flux, "heat_full_flux"=>heat_full_flux, "salt_flux"=>salt_flux, "salt_full_flux"=>salt_full_flux, "eddy_diffusivity"=>eddy_diffusivity, "S_mean"=>S_mean, "T_mean"=>T_mean_C, "TKE"=>TKE, "TKE_adv"=> TKE_adv, "P_s"=>P_s, "B_prod"=>B_prod, "epsilon"=>epsilon, "pressure_transport"=>pressure_transport, "turbulent_transport"=>turbulent_transport, "viscous_transport"=>viscous_transport, "epsilon2"=>epsilon2, "IW_flux"=>IW_flux)
dims = Dict("U"=>("zC",), "V"=>("zC",), "W"=>("zC",), "heat_flux"=>("zC",), "turbulent_heat_flux"=>("zC",), "heat_full_flux"=>("zC",), "salt_flux"=>("zC",), "salt_full_flux"=>("zC",), "eddy_diffusivity"=>("zC",), "S_mean"=>("zC",), "T_mean"=>("zC",
), "TKE"=>("zC",), "TKE_adv"=>("zC,"), "P_s"=>("zC",), "B_prod"=>("zC",), "epsilon"=>("zC",), "pressure_transport"=>("zC",), "turbulent_transport"=>("zC",), "viscous_transport"=>("zC",), "epsilon2"=>("zC,"), "IW_flux"=>("zC,"))

#JLD2 output writer (for animating the simulation run)
simulation.output_writers[:slices] =
    JLD2OutputWriter(model, merge(model.velocities, model.tracers, S_mean_out, T_mean_out, U_out),
                     filename = filename * ".jld2",
                     schedule = TimeInterval(10minute),
                     indices = (:, grid.Ny/2, :),
                     overwrite_existing = true)

# NetCDF output writer (for data analysis)
simulation.output_writers[:data_analysis] =
NetCDFOutputWriter(model, filename=filename * ".nc", output,
                   dimensions=dims,
                   schedule=TimeInterval(1minute), 
                   overwrite_existing=true)

#========================= Run =======================#
run!(simulation)

Thanks!

[glwagner]

It looks like you are using the WENO advection scheme which means that there will be numerical kinetic energy dissipation. Might this be the culprit for the unclosed budget? You can try switching to Centered(order=2) for the advection scheme. Generally it's probably a good idea to either use upwinded advection to provided numerical dissipation or use an LES closure like SmagorinskyLilly, but not both at the same time (an explicit closure may be better for you if you would like to close the TKE budget).

There's been similar discussions recently over on Oceanostics, eg tomchor/Oceanostics.jl#170. Maybe @tomchor and @zhihua-zheng have comments to make.

Your code looks great and I appreciate how neat it is -- that makes it easy to read and understand. Here are just a few notes I got from a quick read:

• Are you able to use Oceanostics to compute viscous transport? It looks like you are using the package but not for all terms... if terms are missing, it may make sense to open an issue over there.
• Very minor considering the code as a whole, but I don't think you need to wrap trans_3_d in a Field before Average -- it will be more efficient to average the operator directly (if that will compile). That's true more broadly, I notice that essentially all operations are wrapped in Field.
• Your IW_flux is divided by 1024 which looks like a reference density, but note that p is kinematic pressure and there is no reference density in Oceananigans. So you may want to remove the division by 1024.
• You may be able to write FPlane(latitude=75) which could be more understandable than the manual Coriolis calculation.

[SamDeAbreu]

I have tried 10th order center differencing, but the simulation became quickly filled with error/oscillations from what I believe the creation and maintenance of sharp gradients.

After playing around with it, it seems AMD with a potentially modified model constant and 2nd order center differencing seems to perform the best. Thank you everyone for your help!

[francispoulin]

Very good to know. I was thinking 4th or 6th but glad to know that 10th is easily accessible. But I can believe that the dispersion is much more intense at that high order.

[glwagner]

High order central differencing is more appropriate for DNS-type situations where gradients tend to be smooth at the grid scale. Note though that the "high-order" is not necessarily indicative of accuracy, since the reconstructions are 1D rather than 3D as they would need to be (plus, even with high-order advection, every other difference / reconstruction would still be 2nd order).

[parfenyev]

High order central differencing is more appropriate for DNS-type situations where gradients tend to be smooth at the grid scale. Note though that the "high-order" is not necessarily indicative of accuracy, since the reconstructions are 1D rather than 3D as they would need to be (plus, even with high-order advection, every other difference / reconstruction would still be 2nd order).

Maybe you can recommend a source where we can read about different advection and timestepper schemes, their pros and cons, and when it is best to use them? I am particularly interested in what scheme is best to use for DNS of 2D turbulence in a domain with no-slip walls in NonhydrostaticModel. Now I use timestepper = :RungeKutta3 and advection = UpwindBiasedFifthOrder(), but I am not sure if this is optimal.

[glwagner]

For DNS you should use advection=Centered(order=2). This way, ringing (and instability) will signal the need to increase resolution; when the problem is resolved diffusion dominates at the grid scale and the oscillatory errors of second order advection are no issue. I also think timestepper=:RungeKutta3 is going to give you fastest time to solution because of the diffusive CFL it allows (in DNS your time step will be limited by diffusion) but I should find a reference for the exact stability critieria. @jm-c might know...

[zhihua-zheng]

1. To calculate TKE budget terms using existing functions in Oceanostic, you can remove the mean (https://github.com/tomchor/Oceanostics.jl/blob/51f477ed1e3d9d2ab0683c47609f2f8e2afd6243/src/Oceanostics.jl#L82-L98) and pass the perturbation fields only. However, I tried this but it doesn't seem to give the right subgrid stress perturbation ($\tau^{\prime}_{ij}$, which is the full subgrid stress with horizontal mean removed). My guess is that the viscosities are not being computed correctly if the kernel (e.g., viscous_flux_uz) only has perturbation fields as input.
   What I end up doing is just passing perturbation velocities and taking advantage of horizontal average. For example (here), the TKE dissipation term is calculated as:

$$\Big\langle \tau_{ij} \frac{\partial u^{\prime}_i}{\partial x_j} \Big\rangle = \Big\langle \langle \tau_{ij} \rangle \frac{\partial u^{\prime}_i}{\partial x_j} + \tau^{\prime}_{ij} \frac{\partial u^{\prime}_i}{\partial x_j} \Big\rangle = \Big\langle \tau^{\prime}_{ij} \frac{\partial u^{\prime}_i}{\partial x_j} \Big\rangle$$

2. Viscous transport can be calculated as the difference between KineticEnergyStressTerm and KineticEnergyDissipationRate in Oceanostics. Likewise, you can use the same trick to get this term in TKE budget. The turbulent transport term can be calculated from AdvectionTerm using perturbation velocities.

3. As @glwagner pointed out, the surface momentum flux is not included in the sub-grid stress for FluxBoundaryCondition, so you have to add that separately (see below). I recall that was the cause for my issue with the topmost cell. I think this is also maybe missing in the KineticEnergyTendency in Oceanostics, as the flux is added after the tendency calculation (not sure)?

@inline kernel_getbc(i, j, k, grid, boundary_condition, clock, fields) =
    getbc(boundary_condition, i, j, grid, clock, fields)

@inline function SurfaceMomentumFlux(model::NonhydrostaticModel, momentum_name)
    model_fields = fields(model)
    momentum = model.velocities[momentum_name]
    mom_bc = momentum.boundary_conditions.top
    LX = location(momentum, Int32(1))
    LY = location(momentum, Int32(2))
    return KernelFunctionOperation{LX, LY, Nothing}(kernel_getbc, model.grid, mom_bc, model.clock, model_fields)
end

# Surface fluxes
Qu = Field(Average(SurfaceMomentumFlux(model, :u), dims=(1,2)))
Qv = Field(Average(SurfaceMomentumFlux(model, :v), dims=(1,2)))

[zhihua-zheng]

In this case, I want to calculate $\tau^{\prime}_{ij} = -2 \nu_e \partial_j u^{\prime}_i$. Is there already decomposition of this kind?

[glwagner]

Yes. Every computations (every single one!) is an operator that can be called with KernelFunctionOperation. For example the ones you're referring to are here:

Oceananigans.jl/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl

Lines 153 to 157 in 41209a2

	# Horizontal viscous fluxes for isotropic diffusivities
	@inline ν_σᶜᶜᶜ(i, j, k, grid, closure, K, clock, fields, σᶜᶜᶜ, args...) = νᶜᶜᶜ(i, j, k, grid, closure, K, clock, fields) * σᶜᶜᶜ(i, j, k, grid, args...)
	@inline ν_σᶠᶠᶜ(i, j, k, grid, closure, K, clock, fields, σᶠᶠᶜ, args...) = νᶠᶠᶜ(i, j, k, grid, closure, K, clock, fields) * σᶠᶠᶜ(i, j, k, grid, args...)
	@inline ν_σᶠᶜᶠ(i, j, k, grid, closure, K, clock, fields, σᶠᶜᶠ, args...) = νᶠᶜᶠ(i, j, k, grid, closure, K, clock, fields) * σᶠᶜᶠ(i, j, k, grid, args...)
	@inline ν_σᶜᶠᶠ(i, j, k, grid, closure, K, clock, fields, σᶜᶠᶠ, args...) = νᶜᶠᶠ(i, j, k, grid, closure, K, clock, fields) * σᶜᶠᶠ(i, j, k, grid, args...)

The object K is model.diffusivity_fields. This also makes me realize that your guess about the viscosities isn't right --- the eddy viscosities are precomputed in a separate step. So if you put in eddy viscosities computed with the full velocity field, then calling for example viscous_flux_ux(i, j, k, grid, closure, diffusivity_fields, clock, perturbation_model_fields, buoyancy) will use the precomputed eddy viscosity in diffusivity_fields but the perturbation_model_fields for ux.

[glwagner]

Here's the implementation of SmagorinskyLilly for example: https://github.com/CliMA/Oceananigans.jl/blob/main/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

The eddy viscosity is computed here:

Oceananigans.jl/src/TurbulenceClosures/turbulence_closure_implementations/smagorinsky_lilly.jl

Lines 123 to 134 in 41209a2

	function compute_diffusivities!(diffusivity_fields, closure::SmagorinskyLilly, model; parameters = :xyz)
	arch = model.architecture
	grid = model.grid
	buoyancy = model.buoyancy
	velocities = model.velocities
	tracers = model.tracers
	launch!(arch, grid, parameters, _compute_smagorinsky_viscosity!,
	diffusivity_fields.νₑ, grid, closure, buoyancy, velocities, tracers)
	return nothing
	end

Note that compute_diffusivities is called in update_state! at the end of every time step, before diagnostics / output is computed.

[zhihua-zheng]

calling for example viscous_flux_ux(i, j, k, grid, closure, diffusivity_fields, clock, perturbation_model_fields, buoyancy) will use the precomputed eddy viscosity in diffusivity_fields but the perturbation_model_fields for ux.

Then I don't understand why using perturbation_model_fields for the dissipation gave me much smaller values back then. I need to go back and confirm it. Thank you for pointing out that!

[glwagner]

I also suspect that the eddy viscosity would not be too different when calculated from perturbation fields, because it seems likely that gradients are dominated by the perturbation... but not positive about that of course...