diff --git a/docs/development/les_normalization_candidate.md b/docs/development/les_normalization_candidate.md new file mode 100644 index 0000000..412956e --- /dev/null +++ b/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. diff --git a/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh b/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh index 0ff99e8..76e2747 100644 --- a/src/CelerisLab/lbm/kernels/operators/turbulence_smag.cuh +++ b/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 diff --git a/tests/unit/test_les_flow_runner.py b/tests/unit/test_les_flow_runner.py index c8584ea..98f4d7f 100644 --- a/tests/unit/test_les_flow_runner.py +++ b/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", diff --git a/tests/unit/test_smagorinsky_normalization.py b/tests/unit/test_smagorinsky_normalization.py new file mode 100644 index 0000000..300856c --- /dev/null +++ b/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()