Smagorinsky LES adds each velocity gradient to itself, so mu_t uses |grad u| not |S|

[WeiqunZhang]

Severity: high/medium · Category: physics-numerics · Fix order: 8 of 21 — fix this 8th.

_{Filenames are numbered in reverse fix order: 001 = fix last, 021 = fix first. This file is 014.}

Location: Source/NS_LES.cpp:128

Based on commit 9bf664bf (line numbers refer to that tree).

All five findings are the same single-line defect in the Smagorinsky branch of calc_mut_LES, reported independently; one patch closes them.

The defect

What is wrong — Smagorinsky model adds a velocity-gradient component to itself instead of its transpose: 'src(i,j,k,i_symij) + src(i,j,k,i_symij)' computes 2*dui/dxj, not the strain rate dui/dxj + duj/dxi, so mu_t is based on |grad u| (strain plus rotation) rather than |S|. _{Reported independently at this same line by F048, F049, F051, F057; each reviewer's own wording and evidence is under Verification evidence below.}

Why it matters

Any run with ns.do_LES=1 and LES_model=Smagorinsky (2D or 3D). AMReX compVelGrad stores du_m/dx_n at comp SPACEDIMn+m, so the code yields mu_t=(Csdx)^2sqrt(2(S:S+W:W)). In solid-body rotation (u=-Omegay, v=Omegax) true |S|=0 but the code produces mu_t=(Csdx)^22Omega, spuriously damping vortices; mu_t is overpredicted in every rotational flow.

How to reach it

• ns.do_LES=1 with default ns.LES_model (Smagorinsky) and nonzero ns.vel_visc_coef; e.g. Tutorials/HotSpot/inputs.3d.LES_hotspot with its line 160 (Sigma override) removed. No assert/abort gates it.

Suggested fix

Replace the self-add with the transpose partner. compVelGrad's documented layout (AMReX_MLTensorOp_grad.cpp: component SPACEDIM*n+m holds du_m/dx_n) makes the transpose of flat index c equal (c%SPACEDIM)*SPACEDIM + c/SPACEDIM, so symij becomes g_mn + g_nm = 2*S_mn; the existing 0.5* and sqrt then give exactly sqrt(2 S_ij S_ij), i.e. the standard mu_t = (Cs*Delta)^2 |S|. Only that one term changes, in 2D and 3D.

This is not a Fortran-to-C++ regression: the self-add appears in the introducing commit 9831fa0 ("more validation needed now"), so there is no earlier correct kernel to restore. Leave the Sigma branch untouched, it legitimately uses the raw gradient tensor.

Maintainer calls: whether smago_Cs_cst = 0.18 needs recalibration against a correct |S|, and whether the filter width should stay dx[idim]. Separately, NavierStokes.cpp:2150 adds this kinematic mu_t to a dynamic viscosity with no rho factor. No Exec input sets do_LES, so no baselines move; a solid-body-rotation check (mu_t == 0) would lock the fix in.

For Source/NS_LES.cpp:128 (F003):

--- a/Source/NS_LES.cpp
+++ b/Source/NS_LES.cpp
@@ -125,6 +125,8 @@
             Real smag = 0;
             for (int i_symij = 0; i_symij < dim_fluxes; ++i_symij)
             {
-              Real symij = src(i,j,k,i_symij) + src(i,j,k,i_symij);
+              // compVelGrad stores du_m/dx_n in component AMREX_SPACEDIM*n+m
+              int i_symji = (i_symij%AMREX_SPACEDIM)*AMREX_SPACEDIM + i_symij/AMREX_SPACEDIM;
+              Real symij = src(i,j,k,i_symij) + src(i,j,k,i_symji);
               smag += symij * symij;
             }

_{Pairs each velocity-gradient component with its transpose so symij = 2S_mn, giving smag = 2S:S and mu_t = (Csdx)^2sqrt(2 S_ij S_ij). Relies on MLTensorOp/MLEBTensorOp::compVelGrad ordering documented in AMReX_MLTensorOp_grad.cpp (comp = deriv_dir*SPACEDIM + vel_comp). Rotation contribution now cancels; verify Smagorinsky Cs calibration is for this |S| definition.}

_{Diff(s) are against 9bf664bf, written from the current source and verified only with git apply --check — never compiled, never run, never applied to the tree. Treat them as precise intent, not tested patches.}

Verification evidence

F003 — confirmed (two independent verifier lenses)

Lens 1 (refutation attempt): Line 128: 'Real symij = src(i,j,k,i_symij) + src(i,j,k,i_symij);'. AMReX_MLTensorOp_grad.cpp:164-167 documents compVelGrad output as raw dU/dx, dV/dx, dU/dy,... so self-add gives 2g_c, yielding sqrt(2(S:S+W:W)) after the 0.5/sqrt scaffolding, exactly as claimed. Solid-body rotation gives 2Omega where |S|=0. Reachable via ns.do_LES=1 (NavierStokes.cpp:2139). git log -S shows never fixed since introduction (9831fa0).

Lens 2 (reachability/intent): Source/NS_LES.cpp:128 symij = src(i,j,k,i_symij) + src(i,j,k,i_symij); AMReX MLTensorOp_grad.cpp documents comp layout "dU/dx, dV/dx, dU/dy, dV/dy...", so the transpose is a different comp. mu_LES feeds diffuse_velocity via getViscosity (NavierStokes.cpp:2139-2151, 1024-1031), never overwritten. Introduced verbatim in 9831fa0 ("more validation needed now") — no deliberate-choice evidence.

F048 — confirmed (one verifier lens)

Reported as: Smagorinsky |S| computation adds each velocity-gradient component to itself ('src(i,j,k,i_symij) + src(i,j,k,i_symij)') instead of adding its transpose, so the eddy viscosity is based on |grad u| (including the antisymmetric rotation part) rather than the strain-rate magnitude |S|.

Failure scenario: ns.do_LES=1 with LES_model=Smagorinsky (the default model): for solid-body rotation dU/dy=-Omega, dV/dx=Omega the true S_ij=0 but the code yields sqrt(smag)=2*Omega, so vortex cores receive large spurious SGS viscosity; every LES-Smagorinsky result is quantitatively wrong wherever the flow rotates.

Lens 1 (refutation attempt): Duplicate of F003; same line 128 self-add confirmed in current source. For u=(-Omy,Omx): g:g=2Om^2, code computes sqrt(0.54g:g)=2Om while true |S|=0 — the claimed sqrt(smag)=2*Omega is exact. LES_model defaults to "Smagorinsky" (NavierStokesBase.cpp:144), path gated only by ns.do_LES=1.

F049 — confirmed (one verifier lens)

Reported as: Smagorinsky model builds the strain-rate magnitude from 'src(i,j,k,i_symij) + src(i,j,k,i_symij)' (2x the same raw velocity-gradient component) instead of adding the transpose component, so mu_t uses the full gradient tensor, not the symmetric strain rate.

Failure scenario: Run with ns.do_LES=1, ns.LES_model=Smagorinsky. MLTensorOp::compVelGrad returns raw dU_i/dx_j (AMReX comment confirms layout). For rotation-dominated flow (e.g., solid-body rotation u=(-y,x)) true S_ij=0 so mu_t must vanish, but the code yields mu_t ~ (Csdx)^2|grad u| > 0; eddy viscosity is wrong wherever vorticity is nonzero.

Lens 1 (refutation attempt): Duplicate of F003. AMReX comment in MLTensorOp_grad.cpp:164-167 confirms raw-gradient layout as the claim states. Code yields mu_t=(Csdx)^2sqrt(2)*|grad u|_F, nonzero wherever vorticity is nonzero, instead of |S|-based mu_t. Not fixed in current source.

F051 — confirmed (one verifier lens)

Reported as: The Smagorinsky LES model computes 'symij' as src(i,j,k,i_symij) + src(i,j,k,i_symij) — the gradient component added to ITSELF instead of its transpose — so the strain-rate tensor S_ij = g_ij + g_ji is never formed; the model acts on |grad U| including the antisymmetric (rotation) part.

Failure scenario: Any run with ns.do_LES=1 and LES_model=Smagorinsky (the default model): in solid-body rotation (u=y,-x) where S=0, the code produces nonzero eddy viscosity mu_t = (Csdx)^22|omega| instead of 0; in simple shear mu_t is overpredicted by sqrt(2). All Smagorinsky LES results are systematically wrong in rotational flow.

Lens 1 (refutation attempt): Duplicate of F003; self-add at line 128 confirmed. Simple shear dU/dy=gamma: code gives sqrt(2)*gamma vs correct |S|=gamma — the claimed sqrt(2) overprediction is exact. Solid-body u=(y,-x): code gives 2 = 2|Omega|, matching the claim with omega read as angular velocity. Core defect (transpose never added, S_ij never formed) is correct.

F057 — confirmed (one verifier lens)

Reported as: The Smagorinsky kernel forms the strain tensor as src(n)+src(n) (the same velocity-gradient component added to itself) instead of adding the transpose component g_ji, so it computes |2*gradU| rather than the symmetrized rate-of-strain magnitude.

Failure scenario: Any run with ns.do_LES=1 (LES_model defaults to Smagorinsky). MLTensorOp::compVelGrad returns the raw gradient (dU/dx,dV/dx,...,order documented in AMReX_MLTensorOp_grad.cpp). For solid-body rotation (u=-Omy, v=Omx), S_ij=0 so mu_t should be 0, but the code produces mu_t=(Csdx)^22*Om; simple shear is overestimated by sqrt(2). Silently wrong eddy viscosity in every LES run.

Lens 1 (refutation attempt): Duplicate of F003; line 128 unchanged in current source. Layout claim verified against AMReX_MLTensorOp_grad.cpp comment. For u=(-Omy,Omx) code produces mu_t=(Csdx)^22*Om, exact as claimed (precise magnitude is sqrt(2)|gradU|_F rather than |2 gradU|, a minor wording slip that does not affect the defect). Reachable via ns.do_LES=1 with default LES_model.

---

_{Based on commit 9bf664bf, which is also the tree the audit verified against. From an automated audit of Source/, Tutorials/ and Util/. Audit finding ids: F003, F048, F049, F051, F057. Reviewer unit(s): Diffusion-2 + LES, theme:eb-smallcell, theme:ghost-fillpatch, theme:gpu-capture, theme:state-comp-indexing. Nothing here was compiled or run — the failure scenarios are code reasoning, so the reaching configuration above is the cheapest way to confirm or refute it.}