fix: normalize closed-form Smagorinsky strain magnitude

docs/development/les_normalization_candidate.md

@@ -0,0 +1,70 @@
+# Minimal Smagorinsky normalization integration candidate
+
+This local candidate changes only the two closed-form D2Q9/D3Q19 coefficients
+from 18 to 18 sqrt(2), plus their convention comments. The intended convention
+is |S| = sqrt(2 S:S), with self-consistent effective tau. No Cs default,
+iterative formula/cap/iteration count, collision dispatch, AMR guard, inlet,
+channel-parity rule, or unrelated source file changes.
+
+The resulting turbulence_smag.cuh SHA256 is
+`a3db46c2e7d27530b2d1436e65cdbed0b72cc7e047bfbeef2a640877271b51aa`,
+exactly the independently reviewed temporary correction used in the paired
+scalar and end-to-end flow experiments. Its uncorrected baseline remains at
+validation checkpoint `9805383f9c50d867e8555471baf5dbda8333a652`.
+
+## Evidence and limits
+
+The separately audited scalar-helper GPU diagnostic on the cd3 source campaign
+covered 48 cases per variant. All 32 paired iterative results were bit-identical;
+maximum patched closed-form error was approximately 1.19e-7, versus 0.06454 for
+the baseline against the same standard convention. This is scalar evidence,
+not full-collision or turbulence qualification.
+
+The independently reviewed wall-bounded, x-periodic SRT/shifted-FP32/double-buffer
+shear experiment distinguished baseline from corrected source under identical
+predeclared criteria. Baseline Cs=1 SGS decrement errors at the last two
+checkpoints were approximately 28.50% and 28.05%; corrected errors were 0.0628%
+and 0.00612%. Corrected maximum Cs=1 profile error was approximately 0.11344%,
+and the Cs=0 observations were identical between variants. Cs=1 is an amplified
+normalization diagnostic, not a recommended physical coefficient. These data
+qualify this bounded normalization comparison, not general turbulence accuracy.
+
+The current direct-header host regression uses the existing independent
+Hermite tensor construction and quadratic bisection oracle. It covers 120
+parameter subcases across D2/D3, Cs=0/.16, rho=.7/1/1.3, tau0=.56/.8,
+equilibrium, both signs of shear and trace-free extension, including a 3D axis.
+It also checks converged iterative agreement, Cs=0 omega guard bounds, and
+explicitly proves that removing sqrt(2) in temporary headers fails the oracle.
+No stored candidate header or full solver fixture is added.
+
+MRT full-trace/bulk coupling, forcing/Guo correction, other streaming/storage
+paths, three-dimensional turbulent flows, and AMR LES remain uncertified.
+Full-stress closed form and deviatoric/capped iteration intentionally need not
+agree on isotropic or cap-active stresses. Existing guards remain unchanged.
+
+## Run and retain provenance
+
+From this candidate checkout:
+
+    PYTHONPATH=src python tests/unit/test_smagorinsky_normalization.py -v
+    PYTHONPATH=src python tests/unit/test_les_flow_runner.py -v
+
+The flow runner already supports `--variant corrected` against the exact
+reviewed helper hash. Its source-pin tests now exercise that current variant;
+they still verify that the archived baseline variant rejects the corrected
+source and that unrelated source tampering is rejected.
+
+Historical Re100, Jet, matched-DFG and paired-scalar manifests/runners are not
+repinned or weakened. Their baseline-source guards are expected to reject this
+changed production tree. Run those archived source-specific checks on the clean
+9805383 checkpoint; their old GPU evidence does not become evidence for this
+candidate. The new direct-current-header test is the regression for corrected
+production. Do not describe the historical full suite as passing on this tree.
+
+Before publication or main integration, obtain final independent diff review
+and rerun the applicable scalar/flow GPU commands on the exact final candidate
+commit after the executor reconnects. Retain new checkout/import/source-hash
+receipts and compare against the frozen baseline using separate processes.
+A new commit receipt is required even though its corrected header matches the
+previously tested temporary correction. No GPU rerun, publication, or main
+merge is claimed by this candidate preparation.

src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh

@@ -5,7 +5,9 @@
 //   nu_t = (C_s * Delta)^2 * |S|,  Delta = 1 (lattice unit)
 //   omega_eff = 1 / (3*(nu + nu_t) + 0.5)
 //
-// |S| is estimated from non-equilibrium second-order moments using tau0.
+// |S| = sqrt(2 S_ij S_ij), with Delta = 1 in lattice units.
+// The stress estimate uses self-consistent tau_eff; tau0 is molecular.
+// Closed form uses full stress. Iteration uses deviatoric stress and a cap.
 // ============================================================================
 
 #ifndef CELERIS_OPERATORS_TURBULENCE_SMAG_CUH
@@ -65,12 +67,12 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
 
 #if LES_CLOSED_FORM
     // Closed-form Smagorinsky (no iteration):
-    //   ω = 2 / (τ₀ + √(τ₀² + 18·Cs²·√Q / ρ))
+    //   ω = 2 / (τ₀ + √(τ₀² + 18·√2·Cs²·√Q / ρ))
     // where Q = Πneq_ij · Πneq_ij (Frobenius norm squared)
     const float Q = pixx*pixx + piyy*piyy + 2.0f*pixy*pixy;
     const float sqrtQ = sqrtf(fmaxf(Q, 0.0f));
     const float smag_const_sq = LES_CS * LES_CS;
-    const float discriminant = tau0*tau0 + 18.0f * smag_const_sq * sqrtQ / rho_safe;
+    const float discriminant = tau0*tau0 + 18.0f * sqrtf(2.0f) * smag_const_sq * sqrtQ / rho_safe;
     const float tau_eff = 0.5f * (tau0 + sqrtf(fmaxf(discriminant, tau0*tau0)));
     return clamp_omega(1.0f / fmaxf(tau_eff, 0.500001f));
 #else
@@ -140,7 +142,7 @@ __device__ __forceinline__ float compute_omega_smag(const float* __restrict__ f,
                   + 2.0f*(pixy*pixy + pixz*pixz + piyz*piyz);
     const float sqrtQ = sqrtf(fmaxf(Q, 0.0f));
     const float smag_const_sq = LES_CS * LES_CS;
-    const float discriminant = tau0*tau0 + 18.0f * smag_const_sq * sqrtQ / rho_safe;
+    const float discriminant = tau0*tau0 + 18.0f * sqrtf(2.0f) * smag_const_sq * sqrtQ / rho_safe;
     const float tau_eff = 0.5f * (tau0 + sqrtf(fmaxf(discriminant, tau0*tau0)));
     return clamp_omega(1.0f / fmaxf(tau_eff, 0.500001f));
 #else

tests/unit/test_les_flow_runner.py

@@ -41,9 +41,9 @@ class GuardTests(unittest.TestCase):
         with self.assertRaisesRegex(RuntimeError, "velocity/profile"):
             decode_macros(invalid, w, cx, cy)
 
-    def test_baseline_pins(self):
+    def test_current_corrected_pins(self):
         self.assertEqual(
-            source_guard(ROOT, "baseline")["base_commit"],
+            source_guard(ROOT, "corrected")["base_commit"],
             "cd3a657f9d3833e812c3cca173525f165b48710a",
         )
 
@@ -55,20 +55,7 @@ class GuardTests(unittest.TestCase):
                 dest = root / name
                 dest.parent.mkdir(parents=True, exist_ok=True)
                 shutil.copyfile(ROOT / name, dest)
-            # Rebuild the reviewed candidate from actual pinned headers, not fixtures.
-            candidate = root / "header_candidate"
-            subprocess.run(
-                [
-                    sys.executable,
-                    str(HERE / "les/prepare.py"),
-                    "--repo",
-                    str(ROOT),
-                    "--output-dir",
-                    str(candidate),
-                ],
-                check=True,
-            )
-            shutil.copyfile(candidate / "patched" / HEADER, root / HEADER)
+            # Current production is the corrected variant; archived baseline pins stay frozen.
             with patch(
                 "run_uniform_les_flow.subprocess.check_output",
                 return_value="isolated-fixture\n",

tests/unit/test_smagorinsky_normalization.py

@@ -0,0 +1,99 @@
+"""Direct-current-header regression; historical candidate manifests stay frozen."""
+
+import importlib.util
+from pathlib import Path
+import shutil
+import sys
+from tempfile import TemporaryDirectory
+import unittest
+
+ROOT = Path(__file__).resolve().parents[2]
+HELPERS = ROOT / "tests/validation/les"
+sys.path.insert(0, str(HELPERS))
+spec = importlib.util.spec_from_file_location(
+    "smagorinsky_oracle", HELPERS / "test_candidate.py"
+)
+oracle = importlib.util.module_from_spec(spec)
+spec.loader.exec_module(oracle)
+HEADER = Path("src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh")
+DESCRIPTORS = Path("src/CelerisLab/lbm/kernels/core/descriptors.cuh")
+
+
+class SmagorinskyNormalization(unittest.TestCase):
+    def test_current_header_matches_independent_tensor_oracle(self):
+        with TemporaryDirectory() as directory:
+            build_dir = Path(directory)
+            for dim in (2, 3):
+                for cs in (0.0, 0.16):
+                    # Compile directly against this checkout, not saved candidate headers.
+                    closed = oracle.build(
+                        ROOT.parent, ROOT.name, dim, cs, 1, 3, build_dir
+                    )
+                    converged = oracle.build(
+                        ROOT.parent, ROOT.name, dim, cs, 0, 200, build_dir
+                    )
+                    for rho in (0.7, 1.0, 1.3):
+                        for tau0 in (0.56, 0.8):
+                            for kind, magnitude in (
+                                ("equilibrium", 0.0),
+                                ("shear", -0.02),
+                                ("shear", 0.02),
+                                ("extension", -0.02),
+                                ("extension", 0.02),
+                            ):
+                                stress = oracle.tensor(
+                                    dim, kind, magnitude * rho, 0, dim - 1
+                                )
+                                f, eq = oracle.state(dim, rho, stress)
+                                expected = 1 / oracle.standard_tau(
+                                    tau0, cs, oracle.norm(stress), rho
+                                )
+                                with self.subTest(
+                                    dim=dim,
+                                    cs=cs,
+                                    rho=rho,
+                                    tau=tau0,
+                                    kind=kind,
+                                    sign=magnitude,
+                                ):
+                                    self.assertAlmostEqual(
+                                        closed(f, eq, rho, tau0), expected, delta=3e-6
+                                    )
+                                    self.assertAlmostEqual(
+                                        converged(f, eq, rho, tau0),
+                                        expected,
+                                        delta=3e-6,
+                                    )
+                    if cs == 0:
+                        f, eq = oracle.state(
+                            dim, 1.0, oracle.tensor(dim, "equilibrium", 0.0)
+                        )
+                        self.assertAlmostEqual(
+                            closed(f, eq, 1.0, 0.50001), 1.99, delta=3e-6
+                        )
+                        self.assertAlmostEqual(
+                            closed(f, eq, 1.0, 200.0), 0.01, delta=3e-6
+                        )
+
+    def test_removing_sqrt_two_is_detected(self):
+        with TemporaryDirectory() as directory:
+            root = Path(directory)
+            for relative in (HEADER, DESCRIPTORS):
+                target = root / "old_convention" / relative
+                target.parent.mkdir(parents=True, exist_ok=True)
+                shutil.copyfile(ROOT / relative, target)
+            target = root / "old_convention" / HEADER
+            text = target.read_text()
+            factor = "18.0f * sqrtf(2.0f) * smag_const_sq"
+            self.assertEqual(text.count(factor), 2)
+            target.write_text(text.replace(factor, "18.0f * smag_const_sq"))
+            for dim in (2, 3):
+                old = oracle.build(root, "old_convention", dim, 0.16, 1, 3, root)
+                stress = oracle.tensor(dim, "shear", 0.02)
+                f, eq = oracle.state(dim, 1.0, stress)
+                expected = 1 / oracle.standard_tau(0.8, 0.16, oracle.norm(stress), 1.0)
+                self.assertGreater(abs(old(f, eq, 1.0, 0.8) - expected), 1e-3)
+
+
+if __name__ == "__main__":
+    unittest.main()