77 T rho_old, rho_new, u_old[DESCRIPTOR::d], u_new[DESCRIPTOR::d];;
79 T fNeq[DESCRIPTOR::q];
81 T uSqr_old = util::normSqr<T,DESCRIPTOR::d>(u_old);
83 T tau_new = _converter.getLatticeRelaxationTime();
85 T physDeltaT_new = _converter.getPhysDeltaT();
86 T physDeltaT_old = physDeltaT_new * (_tau_old-0.5) / (tau_new-0.5);
87 T neqScaling = ((physDeltaT_new*tau_new)/(physDeltaT_old*_tau_old));
89 for (
int iDim=0; iDim < DESCRIPTOR::d; iDim++) {
90 u_new[iDim] = u_old[iDim] * physDeltaT_new / physDeltaT_old ;
93 rho_new = (rho_old-1.0) * (physDeltaT_new / physDeltaT_old) * (physDeltaT_new / physDeltaT_old) + 1.0;
95 T uSqr_new = util::normSqr<T,DESCRIPTOR::d>(u_new);
98 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
99 fNeq[iPop] = cell[iPop] - dynamics -> computeEquilibrium(iPop,rho_old,u_old,uSqr_old);
101 output[iPop] = (dynamics -> computeEquilibrium(iPop,rho_old,u_old,uSqr_old)+descriptors::t<T,DESCRIPTOR>(iPop)
102 + fNeq[iPop]*neqScaling)
103 * ((dynamics -> computeEquilibrium(iPop,rho_new,u_new,uSqr_new)+descriptors::t<T,DESCRIPTOR>(iPop))/(dynamics -> computeEquilibrium(iPop,rho_old,u_old,uSqr_old)+descriptors::t<T,DESCRIPTOR>(iPop))) - descriptors::t<T,DESCRIPTOR>(iPop);