One-equation k RANS: wall condition, closure fixes, terrain wall distance, Canonical_RANS regression cases

[hgopalan]

One-equation k RANS: a minimal working model with regression cases

This PR makes the Axell & Liungman (2001) one-equation RANS closure (erf.rans_type = kEqn) a working, tested option on flat and terrain-fitted meshes, and adds Exec/CanonicalTests/Canonical_RANS: five cases laid out like Canonical_LES, each with a Python check script that compares numbers against stated targets, nine CTest entries (label rans), seven gtests, and a theory page. The development record with per-phase findings is Exec/CanonicalTests/Canonical_RANS/PLAN.md; the numbers are in RESULTS.md.

Every behaviour change is opt-in except two bug fixes noted below; the existing kEqn deck (Askervein) only gains erf.dirichlet_k = true.

Fixes

defect	before	after
Wall value of k (dirichlet_k) was written into the first cell once per step, then the ext_dir BC applied a one-sided flux against a zero ghost cell and drained a fifth of it every step	KE(0)/u*^2 = 2.60 against the AL01 value 3.23	3.232; the value is held through every RK stage and in the implicit solve, wall flux of k is zero
Burchard & Petersen smoothing of Rt cancelled catastrophically	-2 at Rt = -1e16, +1.4e14 at -1e30	exact rearrangement, holds Rt_min
Unstable length iterated the fixed-point map of AL01 Eq. 26, which has no fixed point in strong convection	length unbounded	Eq. 28 evaluated once with the smoothed Rt, bounded at about 1.31 l_g (matches the Kynema KLAxell implementation)
Poisson wall distance took the low-face fluxes as the cell gradient, half a cell off centre	every distance short by z dz/(2H), 0.8 % on 64 cells	cell-centred metric gradient: exact to 1e-6 m on a flat fitted mesh

New options (all default off)

erf.tke_floor, erf.implicit_tke_dissipation, erf.rans_consistent_diffusivities, erf.rans_lscale_from_pblh + erf.rans_lscale_min, erf.wall_dist_type = poisson | terrain_height; input validation for the closure constants (the stability functions have poles near Rt = -3.6); a start-up warning when kEqn runs under a surface layer without dirichlet_k; plot variables Rt, cmu, cmu_prime.

The terrain_height wall distance (height above the local surface on its normal, no linear solve) is exact to 1e-10 on a flat fitted mesh and closer to the true distance than the Poisson solve on both hills (mean 0.02 % and 0.01 % against 1.0 % and 0.3 %); the hill decks use it, and _Poisson CTest variants keep the Poisson path covered.

Cases and results

case	run	key results
Neutral_ABL_Flat	12 h Ekman layer	log law within 7 %, KE(0)/u*^2 = 3.232, 1 and 2 ranks agree to 2e-15
Stable_ABL_Flat	GABLS1, 9 h	u* 0.244, jet 1.23 Ug at 154 m, depth 134 m (LES ensemble: 0.26-0.30, 150-200 m)
Convective_ABL_Flat	0.24 K m/s, 4 h	column heat budget closes to 0.02 %, inversion at 1020 m, AL01 buoyancy term at the wall (3.82)
Neutral_Hill_2D	Witch of Agnesi ridge, 6 h	crest speed-up 0.55/0.45/0.39 vs Jackson-Hunt 2h/L = 0.4, upstream log law within 4 %
Neutral_Hill_3D	radial hill, 4 h	crest speed-up 0.41/0.28/0.22 vs 1.6 h/L = 0.32, upstream log law within 10 %

Restart is bit-exact on the flat and both terrain decks.

Findings outside this PR's scope

• Under the anelastic integrator every vertical diffusion is explicit (the implicit column solve is switched off for anelastic, PR Generalize vert_implicit_fac to be per level #3329), so RANS runs are bound by dz^2/(2K): the convective deck aborts at dt = 5 s with K reaching 40 m2/s and runs at 2 s. The column tridiagonal works with the RANS diffusivities (compressible, explicit vs implicit at the same dt: 3e-6 m/s in wind). Enabling it under anelastic is planned separately (PLAN.md phases 9-10).
• On a 3D terrain-fitted mesh with dz different from dx the initial field's divergence is a deterministic 1.788e139 before the first projection and the run aborts, independent of the closure, the boundaries, the terrain source, stretching, the box layout and the projection solver; the Askervein deck is unaffected at dx/dz of 0.5, 1 and 2. A separate uninitialised read in the w boundary fill trips amrex.fpe_trap_invalid under amrex.init_snan in the first advance. Reproducer in PLAN.md phase 6. The 3D hill deck therefore uses unit aspect ratio.
• A local-K closure keeps a superadiabatic lapse of about -2 K/km through a convective mixed layer (1.1 K spread where LES gives under 0.3 K); the check bounds it at 2 K.

Limitations

Dry buoyancy only; no hybrid RANS-LES; the Poisson wall distance is not available with EB or thin bodies; the wall distance is not recomputed after a regrid.

Verification

ctest -L rans (9 entries) and erf_unit_tests --gtest_filter=RANSClosure* (7 tests) pass; the Release build with all warnings on is clean in the changed files; single-precision build passes; codespell clean; sphinx-build reports no error on the touched pages.

🤖 Generated with Claude Code