AL01 Namespace Reference

Functions
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	geom_length (amrex::Real l_g_raw, amrex::Real l_g_max)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	smooth_Rt (amrex::Real Rt, amrex::Real Rt_crit, amrex::Real Rt_min)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	turb_length (amrex::Real l_g, amrex::Real N2, amrex::Real tke, amrex::Real Cmu0_pow3, amrex::Real inv_Cb_sq, amrex::Real Rt_crit, amrex::Real Rt_min, amrex::Real eps)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	unstable_length_bound (amrex::Real l_g, amrex::Real Cmu0_pow3, amrex::Real inv_Cb_sq, amrex::Real Rt_min)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	dissipation (amrex::Real rho, amrex::Real Cmu0_pow3, amrex::Real tke, amrex::Real length)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	richardson (amrex::Real length, amrex::Real N2, amrex::Real tke, amrex::Real Cmu0_pow3)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	cmu (amrex::Real Rt, amrex::Real Cmu0)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	cmu_prime (amrex::Real Rt, amrex::Real Cmu0)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	cmu_numerator_root (amrex::Real Cmu0)
AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real	realizable_nut (amrex::Real nut, amrex::Real tke, amrex::Real S, amrex::Real Cr, amrex::Real &ratio)

Variables
constexpr amrex::Real	Rt_min_lower_bound = amrex::Real(-3.5986889)

Function Documentation

cmu()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::cmu	(	amrex::Real	Rt,
		amrex::Real	Cmu0
	)

Momentum stability function (AL01 Eq. 31).

{
    return (Cmu0 + amrex::Real(0.108)*Rt)
         / (one + amrex::Real(0.308)*Rt + amrex::Real(0.00837)*Rt*Rt);
}

Referenced by ComputeTurbulentViscosityRANS().

Here is the caller graph for this function:

cmu_numerator_root()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::cmu_numerator_root

(

amrex::Real

Cmu0

)

Zero of the numerator of cmu, -Cmu0/0.108 (-5.15 for the default Cmu0). Below it cmu is negative even above the pole, so Rt_min must also lie above it; this is the tighter of the two limits when Cmu0 < 0.3887.

{
    return -Cmu0 / amrex::Real(0.108);
}

Referenced by TurbChoice::init_params().

Here is the caller graph for this function:

cmu_prime()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::cmu_prime	(	amrex::Real	Rt,
		amrex::Real	Cmu0
	)

Scalar stability function (AL01 Eq. 32).

{
    return Cmu0 / (1 + amrex::Real(0.277)*Rt);
}

Referenced by ComputeTurbulentViscosityRANS().

Here is the caller graph for this function:

dissipation()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::dissipation	(	amrex::Real	rho,
		amrex::Real	Cmu0_pow3,
		amrex::Real	tke,
		amrex::Real	length
	)

Dissipation rate per unit volume, rho Cmu0^3 k^1.5 / L (AL01 Eq. 19).

{
    return rho * Cmu0_pow3 * std::pow(tke,amrex::Real(1.5)) / length;
}

Referenced by ComputeTurbulentViscosityRANS().

Here is the caller graph for this function:

geom_length()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::geom_length	(	amrex::Real	l_g_raw,
		amrex::Real	l_g_max
	)

Geometric length scale with the harmonic cap (AL01 Eq. 22 plus limiter).

Parameters

l_g_raw	kappa times (wall distance + z0)
l_g_max	upper bound, about kappa times 0.1 zi

{
    return l_g_max * l_g_raw / (l_g_max + l_g_raw);
}

Referenced by ComputeTurbulentViscosityRANS().

Here is the caller graph for this function:

realizable_nut()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::realizable_nut	(	amrex::Real	nut,
		amrex::Real	tke,
		amrex::Real	S,
		amrex::Real	Cr,
		amrex::Real &	ratio
	)

Realizability cap on the eddy viscosity (erf.rans_realizability = C_r > 0):

nu_t <- min(nu_t, C_r k / S),      S = sqrt(2 S_ij S_ij).

This is the form of the SST limiter. The closure gives nu_t S / k = sqrt(cmu Cmu0^3 P / eps) (AL01 Eqs. 11 and 19), which in a neutral equilibrium layer (cmu = Cmu0, P = eps) is Cmu0^2 = 0.3094, so C_r = Cmu0^2 leaves neutral equilibrium shear untouched and acts where P / eps > Cmu0 / cmu, as in the shear layer that leaves a crest. The cap on nu_t carries over to nu_t' (the scalar diffusivity keeps the ratio cmu'/cmu of the uncapped closure).

Parameters

nut	uncapped eddy viscosity cmu sqrt(k) L
tke	turbulent kinetic energy k (already floored)
S	strain-rate magnitude sqrt(2 S_ij S_ij)
Cr	cap constant C_r > 0 (the caller skips the cap, and the strain read, when erf.rans_realizability is 0)
ratio	out: nu_t / (C_r k / S), the plotfile diagnostic cap_ratio; above 1 the cap is active

Returns

the capped eddy viscosity

S is floored at machine epsilon before the division so the cap evaluates (to a value far above nu_t) without a divide-by-zero trap in cells of uniform flow; the compiler may evaluate both sides of a select.

{
    constexpr amrex::Real eps = std::numeric_limits<amrex::Real>::epsilon();
    AMREX_ASSERT(Cr > amrex::Real(0));
    const amrex::Real cap = Cr * tke / amrex::max(S, eps);
    ratio = nut / cap;
    return amrex::min(nut, cap);
}

Referenced by ComputeTurbulentViscosityRANS().

Here is the caller graph for this function:

richardson()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::richardson	(	amrex::Real	length,
		amrex::Real	N2,
		amrex::Real	tke,
		amrex::Real	Cmu0_pow3
	)

Turbulent Richardson number, Eq. 29 combined with Eq. 19: Rt = k^2 N^2 / eps^2 = L^2 N^2 / (k Cmu0^6).

{
    return length*length * N2 / (tke * Cmu0_pow3 * Cmu0_pow3);
}

Referenced by ComputeDiffusivityYSU(), and ComputeTurbulentViscosityRANS().

Here is the caller graph for this function:

smooth_Rt()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::smooth_Rt	(	amrex::Real	Rt,
		amrex::Real	Rt_crit,
		amrex::Real	Rt_min
	)

Burchard & Petersen smoothing of Rt below Rt_crit, which maps (-inf, Rt_crit) onto (Rt_min, Rt_crit) so the stability functions stay finite in strongly unstable air: Rt_s = Rt - (Rt - Rt_crit)^2 / (Rt + Rt_min - 2 Rt_crit). Written with x = Rt - Rt_crit and a = Rt_min - Rt_crit as Rt_s = Rt_crit + a x / (x + a), which is the same function without the cancellation of two O(|Rt|) terms that returned garbage for |Rt| beyond about 1e15 (tiny k under a strong unstable N^2). Limits: Rt_crit as x -> 0, Rt_min as x -> -inf.

{
    if (Rt >= Rt_crit) { return Rt; }
    const amrex::Real x = Rt - Rt_crit;      // < 0
    const amrex::Real a = Rt_min - Rt_crit;  // < 0 by input validation
    return Rt_crit + a * x / (x + a);
}

Referenced by ComputeTurbulentViscosityRANS(), and turb_length().

Here is the caller graph for this function:

turb_length()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::turb_length	(	amrex::Real	l_g,
		amrex::Real	N2,
		amrex::Real	tke,
		amrex::Real	Cmu0_pow3,
		amrex::Real	inv_Cb_sq,
		amrex::Real	Rt_crit,
		amrex::Real	Rt_min,
		amrex::Real	eps
	)

Turbulent length scale: neutral (l_g), stable (Eq. 26) or unstable (Eq. 28).

Unstable: Eq. 28 is Eq. 26 rewritten with the dissipation of Eq. 19 on the right-hand side. Iterating it (dissipation from the new length, length from the new Rt) is the fixed-point map of Eq. 26 with N^2 < 0, which has no fixed point once l_g^2 |N^2| exceeds Cb^2 k (AL01, p. 78), so the length grows without bound in strongly convective, weakly turbulent air. Here Eq. 28 is evaluated once with the dissipation of the geometric length and the Burchard & Petersen smoothed Rt, as in the Kynema KLAxell implementation, which bounds the unstable length by l_g sqrt(1 + Cmu0^6 |Rt_min| / Cb^2), about 1.31 l_g for the defaults.

Parameters

l_g	capped geometric length
N2	Brunt-Vaisala frequency squared
tke	turbulent kinetic energy (floored, > 0)
Cmu0_pow3	Cmu0 cubed
inv_Cb_sq	one over Cb squared
Rt_crit,Rt_min	smoothing parameters
eps	threshold below which N2 counts as neutral

{
    amrex::Real length;
    if (std::abs(N2) <= eps) {
        length = l_g;
    } else if (N2 > eps) {
        // Stable (AL01, Eqn. 26)
        length = std::sqrt(one /
                (one / (l_g * l_g) + inv_Cb_sq * N2 / tke));
    } else {
        // Unstable (AL01, Eqn. 28) with Rt from the geometric length
        // (Eq. 29 with Eq. 19: Rt = l_g^2 N^2 / (k Cmu0^6)), smoothed
        const amrex::Real Rt = smooth_Rt(l_g*l_g * N2 / (tke * Cmu0_pow3 * Cmu0_pow3),
                                         Rt_crit, Rt_min);
        length = l_g * std::sqrt(one - Cmu0_pow3*Cmu0_pow3 * inv_Cb_sq * Rt);
    }
    return length;
}

Referenced by ComputeTurbulentViscosityRANS().

Here is the call graph for this function:

Here is the caller graph for this function:

unstable_length_bound()

AMREX_GPU_HOST_DEVICE AMREX_FORCE_INLINE amrex::Real AL01::unstable_length_bound	(	amrex::Real	l_g,
		amrex::Real	Cmu0_pow3,
		amrex::Real	inv_Cb_sq,
		amrex::Real	Rt_min
	)

Upper bound of the unstable length for given smoothing parameters.

{
    return l_g * std::sqrt(one - Cmu0_pow3*Cmu0_pow3 * inv_Cb_sq * Rt_min);
}

Variable Documentation

Rt_min_lower_bound

constexpr amrex::Real AL01::Rt_min_lower_bound = amrex::Real(-3.5986889)

constexpr

Lowest admissible Rt_min (exclusive). The pole of cmu is the larger root of 1 + 0.308 Rt + 0.00837 Rt^2, (-0.308 + sqrt(0.308^2 - 4 * 0.00837)) / (2 * 0.00837) = -3.5986889850, and the pole of cmu' (-1/0.277 = -3.6101) lies below it. Just below the pole of cmu its denominator is negative, so cmu and the eddy viscosity are negative (down to the numerator's zero, -Cmu0/0.108). The literal is the root rounded toward zero, which keeps it above the pole in both double and single precision. Values just above it are accepted but give a very large cmu (about 980 at -3.598).

Referenced by TurbChoice::init_params().