178 {
179 using V = typename CELLS::template value_t<names::Component1>::value_t;
180 using DESCRIPTOR = typename CELLS::template value_t<names::Component1>::descriptor_t;
181
186 auto eps =
parameters.template get<EPSILON>();
193 auto alpha =
parameters.template get<ALPHA>();
196
197 Vector<V,10> rhoField{};
198 V rhoTot = 0;
199 Vector<V,DESCRIPTOR::d> u_sum{};
200 auto& cell1 = cells.template get<names::Component1>();
201 rhoField[0] = cell1.template getFieldComponent<descriptors::STATISTIC>(0);
202 Vector<V,DESCRIPTOR::d> u1{};
204 u_sum += u1;
205
206 auto& cell2 = cells.template get<names::Component2>();
207 rhoField[1] = cell2.template getFieldComponent<descriptors::STATISTIC>(0);
208 Vector<V,DESCRIPTOR::d> u2{};
210 u_sum += u2;
211
212#ifdef THIRD_COMPONENT
213 auto& cell3 = cells.template get<names::Component3>();
214 rhoField[2] = cell3.template getFieldComponent<descriptors::STATISTIC>(0);
215 Vector<V,DESCRIPTOR::d> u3{};
217 u_sum += u3;
218#endif
219#ifdef FOURTH_COMPONENT
220 auto& cell4 = cells.template get<names::Component4>();
221 rhoField[3] = cell4.template getFieldComponent<descriptors::STATISTIC>(0);
222 Vector<V,DESCRIPTOR::d> u4{};
224 u_sum += u4;
225#endif
226 for (unsigned n = 0; n < N_COMPONENTS; ++n) {
227 rhoTot += rhoField[n];
228 }
229
230
231 Vector<V,DESCRIPTOR::d> totalForce{};
232 Vector<V,DESCRIPTOR::d> A_ij{};
233
234 V p = POTENTIAL().computeP(rhoField, t, a, b, T_c, m, alpha, g_I, g_II, M);
236
237
238 for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
239 Vector<V,10> rhoField_i{};
240 auto nextCell1 = cell1.neighbor(descriptors::c<DESCRIPTOR>(iPop));
241 rhoField_i[0] = nextCell1.template getFieldComponent<descriptors::STATISTIC>(0);
242 auto nextCell2 = cell2.neighbor(descriptors::c<DESCRIPTOR>(iPop));
243 rhoField_i[1] = nextCell2.template getFieldComponent<descriptors::STATISTIC>(0);
244#ifdef THIRD_COMPONENT
245 auto nextCell3 = cell3.neighbor(descriptors::c<DESCRIPTOR>(iPop));
246 rhoField_i[2] = nextCell3.template getFieldComponent<descriptors::STATISTIC>(0);
247#endif
248#ifdef FOURTH_COMPONENT
249 auto nextCell4 = cell4.neighbor(descriptors::c<DESCRIPTOR>(iPop));
250 rhoField_i[3] = nextCell4.template getFieldComponent<descriptors::STATISTIC>(0);
251#endif
252 V rhoTot_i = 0;
253 for (unsigned n = 0; n < N_COMPONENTS; ++n) {
254 rhoTot_i += rhoField_i[n];
255 }
256 V p_i = POTENTIAL().computeP(rhoField_i, t, a, b, T_c, m, alpha, g_I, g_II, M);
257
258 V phi_i = 6.*(k*p_i - rhoTot_i/3.)/g;
259
260 V psi_i = 0;
261 if (phi_i >= 0){
263 }
264 else {
266 }
267 totalForce += -g * psi * psi_i * descriptors::c<DESCRIPTOR>(iPop) * descriptors::t<V,DESCRIPTOR>(iPop);
268 for (int jPop = 0; jPop < DESCRIPTOR::q; ++jPop) {
269 Vector<V,10> rhoField_j{};
270 auto nextCell1 = cell1.neighbor(descriptors::c<DESCRIPTOR>(jPop));
271 rhoField_j[0] = nextCell1.template getFieldComponent<descriptors::STATISTIC>(0);
272 auto nextCell2 = cell2.neighbor(descriptors::c<DESCRIPTOR>(jPop));
273 rhoField_j[1] = nextCell2.template getFieldComponent<descriptors::STATISTIC>(0);
274#ifdef THIRD_COMPONENT
275 auto nextCell3 = cell3.neighbor(descriptors::c<DESCRIPTOR>(jPop));
276 rhoField_j[2] = nextCell3.template getFieldComponent<descriptors::STATISTIC>(0);
277#endif
278#ifdef FOURTH_COMPONENT
279 auto nextCell4 = cell4.neighbor(descriptors::c<DESCRIPTOR>(jPop));
280 rhoField_j[3] = nextCell4.template getFieldComponent<descriptors::STATISTIC>(0);
281#endif
282 V rhoTot_j = 0;
283 for (unsigned n = 0; n < N_COMPONENTS; ++n) {
284 rhoTot_j += rhoField_j[n];
285 }
286 V p_j = POTENTIAL().computeP(rhoField_j, t, a, b, T_c, m, alpha, g_I, g_II, M);
287
288 V phi_j = 6.*(k*p_j - (rhoTot_j)/3.)/g;
289
290 V psi_j = 0;
291 if (phi_j >= 0){
293 }
294 else {
296 }
297 Vector<V,DESCRIPTOR::d> c1{};
298 Vector<V,DESCRIPTOR::d> c2{};
299 for (
int m = 0;
m < DESCRIPTOR::d; ++
m) {
300 Vector<V,DESCRIPTOR::d> C1{};
301 V C2 = 0.;
302 for (int n = 0; n < DESCRIPTOR::d; ++n) {
303 C1[n] += descriptors::c<DESCRIPTOR>(jPop)[
m] * descriptors::c<DESCRIPTOR>(jPop)[n];
304 C2 += descriptors::c<DESCRIPTOR>(jPop)[n] * descriptors::c<DESCRIPTOR>(jPop)[n] - 1./descriptors::invCs2<V,DESCRIPTOR>();
305 if (n==m){
306 C1[n] -= 1./descriptors::invCs2<V,DESCRIPTOR>();
307 }
308 }
309 for (int n = 0; n < DESCRIPTOR::d; ++n) {
310 c1[
m] += descriptors::c<DESCRIPTOR>(iPop)[n] * C1[n];
311 }
312 c2[
m] += descriptors::c<DESCRIPTOR>(iPop)[
m] * C2;
313 }
314 A_ij = 0.5*g*descriptors::t<V,DESCRIPTOR>(iPop)*descriptors::t<V,DESCRIPTOR>(jPop)*descriptors::invCs2<V,DESCRIPTOR>()*
315 ((3./2.*eps-(sig-1))*c1 + (sig-1)*c2);
316 totalForce += A_ij * psi_i * psi_j;
317 }
318
319 }
320
321
322 auto externalForce1 = cell1.template getField<descriptors::EXTERNAL_FORCE>();
323 Vector<V,DESCRIPTOR::d> uTot = (u_sum + totalForce/2.) / rhoTot + externalForce1/2.;
324 cell1.template setField<descriptors::FORCE>(externalForce1 + chi[0]*totalForce/rhoField[0]);
325 cell1.template setField<descriptors::VELOCITY>(uTot);
326
327 auto externalForce2 = cell2.template getField<descriptors::EXTERNAL_FORCE>();
328 cell2.template setField<descriptors::FORCE>(externalForce2 + chi[1]*totalForce/rhoField[1]);
329 cell2.template setField<descriptors::VELOCITY>(uTot);
330
331#ifdef THIRD_COMPONENT
332 auto externalForce3 = cell3.template getField<descriptors::EXTERNAL_FORCE>();
333 cell3.template setField<descriptors::FORCE>(externalForce3 + chi[2]*totalForce/rhoField[2]);
334 cell3.template setField<descriptors::VELOCITY>(uTot);
335#endif
336#ifdef FOURTH_COMPONENT
337 auto externalForce4 = cell4.template getField<descriptors::EXTERNAL_FORCE>();
338 cell4.template setField<descriptors::FORCE>(externalForce4 + chi[3]*totalForce/rhoField[3]);
339 cell4.template setField<descriptors::VELOCITY>(uTot);
340#endif
341
342 cell1.template setField<descriptors::SCALAR>(p);
343 cell2.template setField<descriptors::SCALAR>(psi);
344 }
platform_constant Fraction t[Q]
platform_constant Fraction M[Q][Q]
constexpr T m(unsigned iPop, unsigned jPop, tag::MRT)
cpu::simd::Pack< T > sqrt(cpu::simd::Pack< T > value)
ADf< T, DIM > abs(const ADf< T, DIM > &a)
meta::list< CHI, TEMPERATURE, G, SIGMA, EPSILON, K, A, B, MM, TCRIT, DEVI, ALPHA, GI, GII > parameters
static void computeJ(CELL &cell, J &j) any_platform
Computation of momentum.