48 {
49 using V = typename CELL::value_t;
50
51
52 auto j = cell.template getField<descriptors::VELOCITY>();
54 cell.template setField<descriptors::VELOCITY>(j);
55
57 const V rho0 =
parameters.template get<RHO0>();
58 const V omega =
parameters.template get<OMEGA>();
59
60 const V rho = cell.template getFieldComponent<descriptors::STATISTIC>(0) * rho0;
61
62 const V rhoTot = rho * omega;
63
64 Vector<V,DESCRIPTOR::d> uTot{};
65 auto u = cell.template getField<descriptors::VELOCITY>();
66 uTot = (u*rho0*omega) / rhoTot;
67
68
69 Vector<V,DESCRIPTOR::d> rhoContribution;
70 const V psi = POTENTIAL().compute(rho,
parameters);
71
72 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
73 const V rho_ = cell.neighbor(descriptors::c<DESCRIPTOR>(iPop))
74 .template getFieldComponent<descriptors::STATISTIC>(0)
75 * rho0;
76 const V psi_ = POTENTIAL().compute(rho_,
parameters);
77 rhoContribution += psi * psi_
78 * descriptors::c<DESCRIPTOR>(iPop)
79 * descriptors::t<V,DESCRIPTOR>(iPop);
80 }
81
82
83
84
85 auto force = cell.template getField<descriptors::EXTERNAL_FORCE>();
86 cell.template setField<descriptors::VELOCITY>(uTot);
87 cell.template setField<descriptors::FORCE>(force - g*rhoContribution/rho);
88 }
typename POTENTIAL::parameters::template include< G, RHO0, OMEGA > parameters
static void computeJ(CELL &cell, J &j) any_platform
Computation of momentum.