171 {
174 V piNeqNormSqr { };
175 MomentaF().computePiNeqNormSqr(cell, piNeqNormSqr);
176 const V rho =
MomentaF().computeRho(cell);
177 const V omega = parameters.template get<descriptors::OMEGA>();
178 const V smagorinsky = parameters.template get<collision::LES::Smagorinsky>();
180 V conSmagoR[DESCRIPTOR::q];
182 V tauMol = V{1} / omega;
183 V
cs2 = V{1} / descriptors::invCs2<V,DESCRIPTOR>();
185 V Phi = (-0.5*(-rho*tauMol*
cs2+
util::sqrt(rho*rho*tauMol*tauMol*cs2*cs2+V{2}*(smagorinsky*smagorinsky)*rho*piNeqNorm))/(smagorinsky*smagorinsky*rho*piNeqNorm));
186 for (int n=0; n < util::TensorVal<DESCRIPTOR>::n; ++n) {
187 S[n] = Phi*pi[n];
188 }
189
190 V SNormSqr = S[0]*S[0] + 2.0*S[1]*S[1] + S[2]*S[2];
192 SNormSqr += S[2]*S[2] + S[3]*S[3] + 2.0*S[4]*S[4] + S[5]*S[5];
193 }
195 V preFactor = smagorinsky*smagorinsky
196 * descriptors::invCs2<V,DESCRIPTOR>()*descriptors::invCs2<V,DESCRIPTOR>()
198
199
200 {
202 V
t = descriptors::t<V,DESCRIPTOR>(q);
203
204 H[0] = descriptors::c<DESCRIPTOR>(q,0)*descriptors::c<DESCRIPTOR>(q,0)-
cs2;
205 H[1] = descriptors::c<DESCRIPTOR>(q,0)*descriptors::c<DESCRIPTOR>(q,1);
206 H[2] = descriptors::c<DESCRIPTOR>(q,1)*descriptors::c<DESCRIPTOR>(q,1)-
cs2;
208 H[2] = descriptors::c<DESCRIPTOR>(q,0)*descriptors::c<DESCRIPTOR>(q,2);
209 H[3] = descriptors::c<DESCRIPTOR>(q,1)*descriptors::c<DESCRIPTOR>(q,1)-
cs2;
210 H[4] = descriptors::c<DESCRIPTOR>(q,1)*descriptors::c<DESCRIPTOR>(q,2);
211 H[5] = descriptors::c<DESCRIPTOR>(q,2)*descriptors::c<DESCRIPTOR>(q,2)-
cs2;
212 }
213
214 V contractHS = H[0]*S[0] + 2.0*H[1]*S[1] + H[2]*S[2];
216 contractHS += H[2]*S[2] + H[3]*S[3] + 2.0*H[4]*S[4] + H[5]*S[5];
217 }
218
219 conSmagoR[
q] =
t*preFactor*SNorm*contractHS;
220 }
221 return conSmagoR[0];
222 }
platform_constant Fraction cs2
platform_constant Fraction t[Q]
constexpr int q() any_platform
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
typename MOMENTA::template type< DESCRIPTOR > MomentaF
static constexpr int n
result stored in n