186 {
187 V rho, u[DESCRIPTOR::d];
188 const V epsilon = 1. - cell.template getField<descriptors::POROSITY>();
189 const V omega =
parameters.template get<descriptors::OMEGA>();
190 const V paramA = V{1.} / omega - V{0.5};
191 MomentaF().computeRhoU(cell, rho, u);
192
193 auto u_s = cell.template getField<descriptors::VELOCITY_SOLID>();
194 if (epsilon < 1e-5) {
196 }
197 else {
198
199 V paramB = (epsilon * paramA) / ((1. - epsilon) + paramA);
200
201 V paramC = (1. - paramB);
202 V omega_s[DESCRIPTOR::q];
203 V cell_tmp[DESCRIPTOR::q];
204 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
205 cell_tmp[iPop] = cell[iPop];
207
208 omega_s[iPop] = (
EquilibriumF().compute(iPop, rho, u_s) - cell[iPop])
209 + (V{1} - omega) * (cell[iPop] -
EquilibriumF().compute(iPop, rho, u));
210
211
212
213
214 }
216 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
217 cell[iPop] = cell_tmp[iPop] + paramC * (cell[iPop] - cell_tmp[iPop]);
218 cell[iPop] += paramB * omega_s[iPop];
219 }
220 for (int iVel=0; iVel<DESCRIPTOR::d; ++iVel) {
221 u[iVel] = paramC * u[iVel] + paramB * u_s[iVel];
222 }
223 }
224
225 return {rho, util::normSqr<V,DESCRIPTOR::d>(u)};
226 }
typename MOMENTA::template type< DESCRIPTOR > MomentaF
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
typename COLLISION::template type< DESCRIPTOR, MOMENTA, EQUILIBRIUM > CollisionO
typename meta::list< descriptors::OMEGA > parameters