54 V fEq[DESCRIPTOR::q] { };
56 V rho, u[DESCRIPTOR::d];
57 MomentaF().computeRhoU(cell, rho, u);
59 const V omega =
parameters.template get<descriptors::OMEGA>();
62 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
63 df[iPop] = cell[iPop] - fEq[iPop];
64 fEq[iPop] += descriptors::t<V,DESCRIPTOR>(iPop);
67 V M200 = 0 , M011 = 0;
70 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
71 V c0 = descriptors::c<DESCRIPTOR::d, DESCRIPTOR::q>(iPop, 0);
72 V c1 = descriptors::c<DESCRIPTOR::d, DESCRIPTOR::q>(iPop, 1);
73 V c2 = descriptors::c<DESCRIPTOR::d, DESCRIPTOR::q>(iPop, 2);
75 M200 += c0 * c0 * df[iPop];
76 M020 += c1 * c1 * df[iPop];
77 M002 += c2 * c2 * df[iPop];
79 M110 += c0 * c1 * df[iPop];
80 M011 += c1 * c2 * df[iPop];
81 M101 += c0 * c2 * df[iPop];
93 V ds[DESCRIPTOR::q] = {
95 Nxz/3 - Nyz/6, Nyz/3 - Nxz/6, - Nxz/6 - Nyz/6,
97 -Pxz/4, Pyz/4, -Pyz/4,
98 Nxz/3 - Nyz/6, Nyz/3 - Nxz/6, - Nxz/6 - Nyz/6,
100 -Pxz/4, Pyz/4, -Pyz/4,
105 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
106 dh[iPop] = df[iPop] - ds[iPop];
112 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
113 dsdh += ds[iPop] * dh[iPop] / fEq[iPop];
114 dhdh += dh[iPop] * dh[iPop] / fEq[iPop];
117 V gamma = dhdh < 1.0e-10 ? 2.0: 1.0/beta - (2 - 1./beta) * dsdh/dhdh;
120 cell.template setField<descriptors::GAMMA>(gamma);
121 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
122 cell[iPop] -= beta * (2*ds[iPop] + gamma*dh[iPop]);