59 {
60
61
62 for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) {
63 V cell_iPop = cell[iPop];
64 cell[iPop] = cell[descriptors::opposite<DESCRIPTOR>(iPop)];
65 cell[descriptors::opposite<DESCRIPTOR>(iPop)] = cell_iPop;
66 }
67
68
69 const auto pop_f = cell.template getField<descriptors::F>();
70 const auto dJdF = cell.template getFieldPointer<descriptors::DJDF>();
71 const V
d = cell.template getField<descriptors::POROSITY>();
72 const V omega =
parameters.template get<descriptors::OMEGA>();
73
74
75 V rho_f { };
76 V u_f[DESCRIPTOR::d] { };
77
78 rho_f = V();
79 for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
80 u_f[iD] = V();
81 }
82 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
83 rho_f += pop_f[iPop];
84 for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
85 u_f[iD] += pop_f[iPop]*descriptors::c<DESCRIPTOR>(iPop,iD);
86 }
87 }
88 rho_f += (V)1;
89 for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
90 u_f[iD] /= rho_f;
91 }
92
93
94 V pheq[DESCRIPTOR::q] { };
95 for (int i=0; i < DESCRIPTOR::q; ++i) {
96 pheq[i] = V{0};
97 for (int j=0; j < DESCRIPTOR::q; ++j) {
99 + descriptors::t<V,DESCRIPTOR>(j);
100 V dot_ij = V();
101 for (int iD=0; iD < DESCRIPTOR::d; ++iD) {
102 dot_ij += ( descriptors::c<DESCRIPTOR>(j,iD) -
d*u_f[iD] )
103 * ( descriptors::c<DESCRIPTOR>(i,iD) - u_f[iD] );
104 }
105 pheq[i] += cell[j]*feq_j*( V{1} + descriptors::invCs2<V,DESCRIPTOR>()*
d*dot_ij );
106 }
107 pheq[i] /= rho_f;
108 }
109
110
111 for (int i=0; i < DESCRIPTOR::q; ++i) {
112 cell[i] = cell[i] - omega*( cell[i] - pheq[i] ) + dJdF[i];
113 }
114
115
116 V rho_phi, u_phi[DESCRIPTOR::d];
117 MomentaF().computeRhoU( cell, rho_phi, u_phi );
118 V uSqr_phi = util::normSqr<V,DESCRIPTOR::d>( u_phi );
119
120
121
122 for (int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) {
123 V cell_iPop = cell[iPop];
124 cell[iPop] = cell[descriptors::opposite<DESCRIPTOR>(iPop)];
125 cell[descriptors::opposite<DESCRIPTOR>(iPop)] = cell_iPop;
126 }
127 return {rho_phi, uSqr_phi};
128 };
constexpr int d() any_platform
typename MOMENTA::template type< DESCRIPTOR > MomentaF
typename meta::list< descriptors::OMEGA > parameters
static V secondOrder(int iPop, const RHO &rho, const U &u, const USQR &uSqr) any_platform
Computation of equilibrium distribution, second order in u.