Expands `TracerVarianceDissipationRate` to work with both AMD and Smagorinsky closure

[tomchor]

At the moment TracerVarianceDissipationRate only works with AMD since the function _calc_κccc isn't defined for any other closure in the Oceananigans side (see https://github.com/CliMA/Oceananigans.jl/issues/2751).

[tomchor]

@glwagner if you have an easier way to support closures other than AMD please feel free to chime in. I couldn't come up with a better way than this one, although I'm sure it exists.

[glwagner]

I suggest using the diffusivity interface to extract the diffusion coefficient.

[glwagner]

Another possibility is to use the lower-level functions κᶠᶜᶜ, κᶜᶠᶜ, κᶜᶜᶠ. These are called with the function signature

κᶠᶜᶜ(i, j, k, grid, closure, diffusivity_fields, id, clock, model_fields)

but thinking even more about this, I think an even better way would be to write the variance dissipation in terms of the fluxes directly, and use diffusive_flux_x, etc.

[glwagner]

Ok, here's one way. It should be tested, but will work (in principle) with any turbulence closure. It may only return expected results for second-order operators though.

# χ is variance dissipation
#
# There are three discrete contributions to χ, each with a different location:
#
# χx at fcc
# χy at cfc
# χz at ccf

@inline χxᶠᶜᶜ(i, j, k, grid, c, args...) = - δxᶠᵃᵃ(i, j, k, grid, c) * diffusive_flux_x(i, j, k, grid, args...) 
@inline χyᶜᶠᶜ(i, j, k, grid, c, args...) = - δyᵃᶠᵃ(i, j, k, grid, c) * diffusive_flux_y(i, j, k, grid, args...)
@inline χzᶜᶜᶠ(i, j, k, grid, c, args...) = - δzᵃᵃᶠ(i, j, k, grid, c) * diffusive_flux_z(i, j, k, grid, args...)

# We can reconstruct this at cell centers:

@inline χᶜᶜᶜ(i, j, k, grid, c, args...) = (ℑxᶜᵃᵃ(i, j, k, grid, χxᶠᶜᶜ, c, args...) +
                                           ℑyᵃᶜᵃ(i, j, k, grid, χyᶜᶠᶜ, c, args...) +
                                           ℑzᵃᵃᶜ(i, j, k, grid, χzᶜᶜᶠ, c, args...)) / Vᶜᶜᶜ(i, j, k, grid)

We'll have to figure out what args... needs to be, I don't know off the top of my head.

EDIT: args can be found here:

https://github.com/CliMA/Oceananigans.jl/blob/4b4bef76d82fb573778c3321b9b4a49503ed0d9c/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl#L201

so its

args = (closure, diffusivity_fields, id, c, clock, model_fields, buoyancy)

[glwagner]

Hmm maybe there's some wrong thinking there. I'm going to go back and try to get the discrete dissipation more carefully.

[glwagner]

Ok, amazingly I think what I wrote above is the correct tracer dissipation on the staggered grid. It takes a bit of algebra to figure out, and it's best to consider the one-dimensional case first. We can write this up in the docs if you like.

[glwagner]

Hmm, looks like there are no docs, but we should probably change that

[tomchor]

Thanks for the help @glwagner !

Ok, amazingly I think what I wrote above is the correct tracer dissipation on the staggered grid. It takes a bit of algebra to figure out, and it's best to consider the one-dimensional case first. We can write this up in the docs if you like.

Assuming you went through the derivation to arrive at this, do you have it handy somewhere already? If so can you just paste here and then I'll transcribe it to the docs once those are set up?

It may only return expected results for second-order operators though.

Do you mind elaborating on this comment, please?

[glwagner]

It may only return expected results for second-order operators though.

Do you mind elaborating on this comment, please?

For higher-order operators, there are further identities (algebraic identities on the discrete staggered grid) that one can do to isolate the part of dissipation that contains no exact integrals.

For example, just multiplying the tracer equation by c produces a term

c * k * lap(c)

For this second-order operator, we can separate this into two terms,

c * k * lap(c) = div(k * c * grad(c)) - k * |grad(c)|^2

The second term is the dissipation, while the first term integrates to zero. For higher-order operators, we need more steps to extract things like the first term that integrate to zero.

Assuming you went through the derivation to arrive at this, do you have it handy somewhere already?

I scribbled it on a piece of paper. Probably best to write it up when we have docs. It involves writing the discrete versions of the above. The important identity can be proven in 1D; it's that

c * δ(q) = δ(ℑ(c) * q) - ℑ(q * δ(c))

where δ is a difference and ℑ is interpolation/reconstruction. We can prove this simply by writing out the terms and doing the algebra.

[tomchor]

@glwagner I just tested this on GPUs and CPUs for Smag and AMD and it seems to be working. The new formulation produces results that are about 2 times smaller than the previous formulation for my simulations. That discrepancy seems a bit high to me, but I checked the the calculations and couldn't find any errors, so I'm tempted to attribute this to the fact that we were approximating a highly fluctuating quantity before.

Any thoughts?

[glwagner]

A factor of two does seem large, but plausible...

It may be possible to test this. For example, consider the diffusion of a random tracer field with advection=nothing (so there's no numerical diffusion from an advection scheme). Then we can compare time-series of Average(c^2/2) to another time-series generated by integrating the variance equation with the RHS calculated from (1) the "exact" variance we're using here and (2) the "approximate" variance (perhaps just with forward Euler to start, which should be ok if we keep the time steps very small I think).

[tomchor]

Okay FYI I ran some not-so-controlled tests in one of my simulations and this new formulation seems to be farther from what I expect to be the correct value of variance dissipation, indicating it isn't correct. So we should not merge this before running some rigorous tests (like the one you mentioned).

The model option you mentioned for the tests (advection=nothing) only works for the hydrostatic model, correct?

[glwagner]

If there's no velocity field then there is no difference between the nonhydrostatic and hydrostatic model.

[glwagner]

But what is the test that verifies that some other formula for the tracer variance dissipation is correct? I'm not sure I understand the reasoning.

[tomchor]

Just recording some progress to come back to it later. The plotted quantities (which should be identifiable from the legend and are integrated in the domain) are from a NonhydrostaticModel simulation with zero velocity and random initial condition for a tracer c. Each panel shows two different ways of quantifying the tracer variance dissipation rate using the tools available in this PR (the top panel is offset so that both quantities start at zero). These are pretty different so something's not right:

Although I should note that the ratio of both estimates is almost exactly 16, which probably means there is a factor of 2 (or more) missing somewhere...

PS: These plots ignore the molecular diffusion of tracer variance, but that term is pretty small.

[glwagner]

You should be comparing c^2 / 2 right? Not c^2. That might be 1 factor of 2...

As for the other factor of 8, well thats 2^3... hm...

Can you confirm the identity

c * δ(q) = δ(ℑ(c) * q) - ℑ(q * δ(c))

?

Maybe just turn molecular diffusion off? You may also be able to use the closure tuple directly with this formulation if you have to include it.

[tomchor]

You should be comparing c^2 / 2 right? Not c^2. That might be 1 factor of 2...

I'll work on the identity later, but for reference, as of right now I think the comparisons are correct, no? From Stull's book:

(For folks who are more familiar with Wyngaard's book, that's equation (5.7) there...)

The factor of 2 is included in the TracerVarDissipRate definition: https://github.com/tomchor/Oceanostics.jl/blob/fdad22b7deba9f12427f4405e8f14ca3c1032999/src/FlowDiagnostics.jl#L364-L369

So unless I'm missing something the correct comparison is indeed with $c^2$.

[tomchor]

Can you confirm the identity

c * δ(q) = δ(ℑ(c) * q) - ℑ(q * δ(c))

?

Assuming δ is the difference operator and ℑ is the interpolation operator, I was able to get the same identity.

I couldn't find it online, but doesn't MITgcm have these diagnostics figured out already? Maybe we could double-check.

[glwagner]

I guess it depends whether you define variance as c^2 or 1/2 c^2. Since kinetic energy is 1/2 u^2, I've used 1/2 c^2 by analogy. Stull also defines "epsilon" as the dissipation rate of 1/2 * c^2.

[tomchor]

I spent some time today investigating this and I think we're forgetting some unit factor in the definition of the tracer variance.

To exemplify, here's a new example that uses a resolved initial noise (to minimize issues related to grid-size noise):

using Oceananigans
using Oceanostics

N = 4
κ = 1

grid = RectilinearGrid(topology=(Periodic, Periodic, Periodic),
                       size=(N,N,N), extent=(1,1,1))
model = NonhydrostaticModel(grid=grid, tracers=:c, closure=ScalarDiffusivity(κ=κ))

# A kind of convoluted way to create x-periodic, resolved initial noise
σx = 2grid.Δxᶜᵃᵃ # x length scale of the noise
σy = 2grid.Δyᵃᶜᵃ # x length scale of the noise
σz = 2grid.Δzᵃᵃᶜ # z length scale of the noise

N = 2^4 # How many Gaussians do we want sprinkled throughout the domain?
x₀ = grid.Lx * rand(N); y₀ = grid.Ly * rand(N); z₀ = -grid.Lz * rand(N) # Locations of the Gaussians

xₚ = x₀ .+ (grid.Lx .* [-2;;-1;;0;;1;;2]) # Make that noise periodic by "infinite" horizontal reflection
yₚ = y₀ .+ (grid.Ly .* [-2;;-1;;0;;1;;2]) # Make that noise periodic by "infinite" horizontal reflection
zₚ = z₀ .+ (grid.Lz .* [-2;;-1;;0;;1;;2]) # Make that noise periodic by "infinite" horizontal reflection

resolved_noise(x, y, z) = sum(@. exp(-(x-xₚ)^2/σx^2 -(y-yₚ)^2/σy^2 -(z-zₚ)^2/σz^2))
set!(model, c=resolved_noise)

simulation = Simulation(model; Δt=grid.Δxᶜᵃᵃ^2/κ/20, stop_time=0.05)

c = model.tracers.c
χ  = Oceanostics.FlowDiagnostics.IsotropicTracerVarianceDissipationRate(model, :c)
χ2 = @at (Center, Center, Center) 2 * κ * (∂x(c)^2 + ∂y(c)^2 + ∂z(c)^2) # Same as χ but non-conservative
ε_q = @at (Center, Center, Center) κ * (∂x(∂x(c^2)) + ∂y(∂y(c^2)) + ∂z(∂z(c^2)))
c² = c^2

∫χdV   = Integral(χ)
∫χ2dV  = Integral(χ2)
∫ε_qdV = Integral(ε_q)
∫c²dV  = Integral(c²)

outputs = (; χ, ∫χdV,
           χ2, ∫χ2dV,
           ε_q, ∫ε_qdV,
           c, c², ∫c²dV)

dt = 1e-4
simulation.output_writers[:tracer] = NetCDFOutputWriter(model, outputs;
                                                        filename = "tracer_diff.nc",
                                                        schedule = TimeInterval(dt),
                                                        overwrite_existing = true)
run!(simulation)

using NCDatasets, GLMakie

xc, yc, zc = nodes(c)

ds = NCDataset(simulation.output_writers[:tracer].filepath, "r")

lines(ds["time"], ds["∫χdV"], label="∫χdV (new conservative form)")
lines!(ds["time"], ds["∫χ2dV"], label="∫χ2dV (old non-conservative form)")

∂t∫c²dV = -diff(ds["∫c²dV"]) / dt
lines!(ds["time"][2:end], ∂t∫c²dV, label="∂(∫c²dV)/∂t")
axislegend()

@show ds["∫χdV"][2:end] ./ ∂t∫c²dV

close(ds)

With N=16 (i.e. 16x16x16 grid), we get these results:

Indicating that the manual, non-conservative calculation of the tracer dissipation rate is very close to correct (based on the rate of decrease of $\int c^2 dV$), but the new calculation proposed by this PR is way off. In fact, it's a factor of almost exactly 64 off:

julia> @show mean(ds["∫χdV"][2:end] ./ ∂t∫c²dV)
mean((ds["∫χdV"])[2:end] ./ ∂t∫c²dV) = 63.703221099778375
63.703221099778375

When I run the same example with N=4, the results are similar but instead of the results being off by a factor of 64, we get a factor of 16:

julia> @show mean(ds["∫χdV"][2:end] ./ ∂t∫c²dV)
mean((ds["∫χdV"])[2:end] ./ ∂t∫c²dV) = 15.949624383944931
15.949624383944931

I'll sit down at some point soon and work out the units of the functions diffusive_flux_z(), etc., but I'm pretty sure we need to multiply each flux by its areas here:

https://github.com/tomchor/Oceanostics.jl/blob/d1983465bf8a6688f2ecc6adade895e9e191eee2/src/FlowDiagnostics.jl#L362-L367

[tomchor]

Indeed after multiplying the fluxes by the appropriate area the lines all match

What a difference from the non-conservative form.

[tomchor]

A note is that this formulation can also work for anisotropic diffusivities, but I think at the moment the diffusive_flux_x() functions aren't defined for tuples since when I use a tuple closure they fall back to here and return zero: https://github.com/CliMA/Oceananigans.jl/blob/2cd91ab4f42e528981422171d6ab2f323a1ce044/src/TurbulenceClosures/abstract_scalar_diffusivity_closure.jl#L158-L169

[tomchor]

@glwagner I think this is ready to merge. Let me know if you have any comments/suggestions otherwise I'll merge this soon.