Shear used in EDMF sources from precomputed args becomes negative

[ilopezgp]

Description

See draft PR #2148 to reproduce this behavior. An erroneous behavior is encountered when drawing the squared shear from the precomputed arguments in the computation of sources. It is observed that the shear is negative when unpacked, but not when computed.

The ekman_layer.jl simulation in this case has the following characteristics:

• Make BCs consistent with ICs to ensure there are no oscillations coming from the ground. This was tested offline without other contributions to the TKE equation.
• EDMF is tested in decoupled mode, with the grid-mean being evolved by a constant kinematic viscosity closure. For the results in the grid mean variables, refer to Grid imprinting in specific variables (per unit mass) due to filtering, and energy-momentum leakage #2116 .
• Make mixing length constant and equal to 5 m.
• Strip the TKE equations of all sources except for the shear production. It is observed that the TKE goes negative, which is inconsistent with the fact that the shear production is always positive.

[yairchn]

adding this print in edmf_kernels.jl line 1171 - in "flux(::en_ρatke, ::Diffusion, atmos, args)"
""""
Shear² = diffusive.turbconv.S²
if Shear² < 0
@print("Shear² in flux(en_ρatke) = ", Shear², "\n")
end
"""
this print still gets negative values:

Shear² in flux(en_ρatke) = -0.12189892688029821
Shear² in flux(en_ρatke) = -0.12189892688020655
Shear² in flux(en_ρatke) = -0.12189892688031359
Shear² in flux(en_ρatke) = -0.12189892688012727
Shear² in flux(en_ρatke) = -0.12189892688031001
Shear² in flux(en_ρatke) = -0.12189892688062727
Shear² in flux(en_ρatke) = -0.1218989268805853
Shear² in flux(en_ρatke) = -0.12189892688023629
Shear² in flux(en_ρatke) = -0.12189892688034511
Shear² in flux(en_ρatke) = -0.1218989268809321
Shear² in flux(en_ρatke) = -0.12189892688045434
Shear² in flux(en_ρatke) = -0.12189892688071481
Shear² in flux(en_ρatke) = -0.12189892688107767

[charleskawczynski]

I can confirm that this is indeed strange.

[charleskawczynski]

From digging around for the day, I've narrowed it down to here. cc @mwarusz. Going to continue looking tomorrow. Thanks for setting up the PR @ilopezgp, it's been helpful so far.

[mwarusz]

It makes sense that this is happening after the gradient interface contribution is added. I am not sure that it is worth investigating any further. This computation needs to happen after the full gradient is available inside source!. Otherwise it is wrong because compute_gradient_flux! is assumed linear. Furthermore, I am not aware of any work using nonlinear auxiliary LDG variables and trying to compute the squared shear directly (and making sure that it stays positive) is bound to cause problems. The solution is to calculate and store velocity gradients and then compute the squared shear from them.

[yairchn]

Thanks @mwarusz we will work towards it (tore velocity gradients and then compute the squared shear from them).
It seems to me that a warning in the form of a comment should be left in the compute_gradient_flux as this might not be trivial to every user / developer.

[mwarusz]

This is mentioned in compute_gradient_flux! docstring

ClimateMachine.jl/src/BalanceLaws/interface.jl

Lines 205 to 206 in bb51cbb

	[`BalanceLaw`](@ref) subtype `BL`. This should be a linear function of
	`∇transformstate`

but I sort of wonder if renaming this function to something like linear_gradient_transform would help prevent this kind of issues.

[ilopezgp]

This issue has been resolved by #2158 and #2160 .