52 {
53
54 V fEq[DESCRIPTOR::q] { };
56 V rho, u[DESCRIPTOR::d];
57 MomentaF().computeRhoU(cell, rho, u);
58
59 const V omega =
parameters.template get<descriptors::OMEGA>();
60
61 V df[DESCRIPTOR::q];
62 for (int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
63 df[iPop] = cell[iPop] - fEq[iPop];
64 fEq[iPop] += descriptors::t<V,DESCRIPTOR>(iPop);
65 }
66
67 V M200 = 0 , M011 = 0;
68 V M002 = 0, M110 = 0;
69 V M020 = 0, M101 = 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);
74
75 M200 += c0 * c0 * df[iPop];
76 M020 += c1 * c1 * df[iPop];
77 M002 += c2 * c2 * df[iPop];
78
79 M110 += c0 * c1 * df[iPop];
80 M011 += c1 * c2 * df[iPop];
81 M101 += c0 * c2 * df[iPop];
82 }
83
84
85 V Nxz = M200 - M002;
86 V Nyz = M020 - M002;
87
88 V Pxy = M110;
89 V Pxz = M101;
90 V Pyz = M011;
91
92
93 V ds[DESCRIPTOR::q] = {
94 0,
95 Nxz/3 - Nyz/6, Nyz/3 - Nxz/6, - Nxz/6 - Nyz/6,
96 Pxy/4, -Pxy/4, Pxz/4,
97 -Pxz/4, Pyz/4, -Pyz/4,
98 Nxz/3 - Nyz/6, Nyz/3 - Nxz/6, - Nxz/6 - Nyz/6,
99 Pxy/4, -Pxy/4, Pxz/4,
100 -Pxz/4, Pyz/4, -Pyz/4,
101 };
102
103
104 V dh[DESCRIPTOR::q];
105 for (int iPop = 0; iPop < DESCRIPTOR::q; ++iPop) {
106 dh[iPop] = df[iPop] - ds[iPop];
107 }
108
109 V beta = omega / 2;
110 V dsdh(0);
111 V dhdh(0);
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];
115 }
116
117 V gamma = dhdh < 1.0e-10 ? 2.0: 1.0/beta - (2 - 1./beta) * dsdh/dhdh;
118
119
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]);
123 }
124
125 return statistic;
126 };
typename EQUILIBRIUM::template type< DESCRIPTOR, MOMENTA > EquilibriumF
typename MOMENTA::template type< DESCRIPTOR > MomentaF
typename meta::list< descriptors::OMEGA > parameters