476 {
477 using DESCRIPTOR = typename CELLS::template value_t<names::NavierStokes>::descriptor_t;
479 cells.template get<names::NavierStokes>().computeAllMomenta(rho, u, pi);
480
481 T PiNeqNormSqr = pi[0]*pi[0] + 2.0*pi[1]*pi[1] + pi[2]*pi[2];
483 PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2*pi[4]*pi[4] +pi[5]*pi[5];
484 }
486
487 if(PiNeqNorm != T{0.}){
488 auto omega = cells.template get<names::NavierStokes>().template getFieldPointer<descriptors::OMEGA>();
489 PiNeqNorm /= (-rho/omega[0]/descriptors::invCs2<T,DESCRIPTOR>());
490 T angle =
parameters.template get<FRICTION_ANGLE>();
491 T ps = ( T{0.5}*omega[0] - T{1} ) * T{1./3.} * ( pi[0] + pi[3] + pi[5] );
492 T granularTau = ps/rho *
util::sin(T{3.14}*angle/T{180}) / PiNeqNorm * descriptors::invCs2<T,DESCRIPTOR>() + T{0.5};
493 omega[0] = T{1} / granularTau;
494 }
495
496
497 auto force = cells.template get<names::NavierStokes>().template getFieldPointer<descriptors::FORCE>();
498 auto forcePrefactor =
parameters.template get<FORCE_PREFACTOR>();
499 T densityDifference = (rho -
parameters.template get<RHO0>()) / rho;
500 for (unsigned iD = 0; iD < DESCRIPTOR::d; ++iD) {
501 force[iD] = forcePrefactor[iD] * densityDifference;
502 }
503 }
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
ADf< T, DIM > sin(const ADf< T, DIM > &a)
meta::list< FRICTION_ANGLE, RHO0, FORCE_PREFACTOR > parameters
static constexpr int n
result stored in n