fix: normalize closed-form Smagorinsky strain magnitude
This commit is contained in:
@@ -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.
|
||||
@@ -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
|
||||
|
||||
@@ -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",
|
||||
|
||||
@@ -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()
|
||||
Reference in New Issue
Block a user