142 {
143 V drho = rho - 1;
144 V inverse_rho = 1. / rho;
145 V uSqr = util::normSqr<V,DESCRIPTOR::d>(u);
146
147
148
149 V moments[DESCRIPTOR::q];
151 auto [mbbb, mabb, mbab, mbba, maab, macb, maba, mabc, mbaa, mbac, maaa, maac, maca, macc, mcbb, mbcb, mbbc, mccb, mcab, mcbc, mcba, mbcc, mbca, mccc, mcca, mcac, mcaa] = moments;
152
153 const V omega2 = 1;
154 const V omega3 = 1;
155 const V omega4 = 1;
156 const V omega5 = 1;
157 const V omega6 = 1;
158 const V omega7 = 1;
159 const V omega10 = 1;
160
161
162 V CUMcbb = mcbb - ((mcaa + 1. / 3) * mabb + 2 * mbba * mbab) * inverse_rho;
163 V CUMbcb = mbcb - ((maca + 1. / 3) * mbab + 2 * mbba * mabb) * inverse_rho;
164 V CUMbbc = mbbc - ((maac + 1. / 3) * mbba + 2 * mbab * mabb) * inverse_rho;
165
166 V CUMcca = mcca - (((mcaa * maca + 2 * mbba * mbba) + 1. / 3 * (mcaa + maca)) * inverse_rho - 1. / 9 * (drho * inverse_rho));
167 V CUMcac = mcac - (((mcaa * maac + 2 * mbab * mbab) + 1. / 3 * (mcaa + maac)) * inverse_rho - 1. / 9 * (drho * inverse_rho));
168 V CUMacc = macc - (((maac * maca + 2 * mabb * mabb) + 1. / 3 * (maac + maca)) * inverse_rho - 1. / 9 * (drho * inverse_rho));
169
170
171 V CUMbcc = mbcc - ((maac * mbca + maca * mbac + 4 * mabb * mbbb + 2 * (mbab * macb + mbba * mabc)) + 1. / 3 * (mbca + mbac)) * inverse_rho;
172 V CUMcbc = mcbc - ((maac * mcba + mcaa * mabc + 4 * mbab * mbbb + 2 * (mabb * mcab + mbba * mbac)) + 1. / 3 * (mcba + mabc)) * inverse_rho;
173 V CUMccb = mccb - ((mcaa * macb + maca * mcab + 4 * mbba * mbbb + 2 * (mbab * mbca + mabb * mcba)) + 1. / 3 * (macb + mcab)) * inverse_rho;
174
175
176 V CUMccc = mccc + ((-4 * mbbb * mbbb - (mcaa * macc + maca * mcac + maac * mcca) - 4 * (mabb * mcbb + mbab * mbcb + mbba * mbbc)
177 - 2 * (mbca * mbac + mcba * mabc + mcab * macb)) * inverse_rho
178 + (4 * (mbab * mbab * maca + mabb * mabb * mcaa + mbba * mbba * maac)
179 + 2 * (mcaa * maca * maac) + 16 * mbba * mbab * mabb) * inverse_rho * inverse_rho
180 - 1. / 3 * (macc + mcac + mcca) * inverse_rho
181 - 1. / 9 * (mcaa + maca + maac) * inverse_rho
182 + (2 * (mbab * mbab + mabb * mabb + mbba * mbba)
183 + (maac * maca + maac * mcaa + maca * mcaa)
184 + 1. / 3 * (maac + maca + mcaa)) * inverse_rho * inverse_rho * 2. / 3
185 + 1. / 27 * ((drho * drho - drho) * inverse_rho * inverse_rho));
186
187
188 V mxxPyyPzz = mcaa + maca + maac;
189 V mxxMyy = mcaa - maca;
190 V mxxMzz = mcaa - maac;
191
192 V mxxyPyzz = mcba + mabc;
193 V mxxyMyzz = mcba - mabc;
194
195 V mxxzPyyz = mcab + macb;
196 V mxxzMyyz = mcab - macb;
197
198 V mxyyPxzz = mbca + mbac;
199 V mxyyMxzz = mbca - mbac;
200
201
202 mxxPyyPzz += omega2 * (maaa - mxxPyyPzz);
203 mxxMyy += -(-omega) * (-mxxMyy);
204 mxxMzz += -(-omega) * (-mxxMzz);
205
206 mabb += omega * (-mabb);
207 mbab += omega * (-mbab);
208 mbba += omega * (-mbba);
209
210
211
212 mbbb += omega4 * (-mbbb);
213 mxxyPyzz += omega3 * (-mxxyPyzz);
214 mxxyMyzz += omega4 * (-mxxyMyzz);
215 mxxzPyyz += omega3 * (-mxxzPyyz);
216 mxxzMyyz += omega4 * (-mxxzMyyz);
217 mxyyPxzz += omega3 * (-mxyyPxzz);
218 mxyyMxzz += omega4 * (-mxyyMxzz);
219
220
221 mcaa = 1. / 3 * (mxxMyy + mxxMzz + mxxPyyPzz);
222 maca = 1. / 3 * (-2 * mxxMyy + mxxMzz + mxxPyyPzz);
223 maac = 1. / 3 * (mxxMyy - 2 * mxxMzz + mxxPyyPzz);
224
225 mcba = (mxxyMyzz + mxxyPyzz) * 0.5;
226 mabc = (-mxxyMyzz + mxxyPyzz) * 0.5;
227 mcab = (mxxzMyyz + mxxzPyyz) * 0.5;
228 macb = (-mxxzMyyz + mxxzPyyz) * 0.5;
229 mbca = (mxyyMxzz + mxyyPxzz) * 0.5;
230 mbac = (-mxyyMxzz + mxyyPxzz) * 0.5;
231
232 CUMacc = (1 - omega6) * (CUMacc);
233 CUMcac = (1 - omega6) * (CUMcac);
234 CUMcca = (1 - omega6) * (CUMcca);
235 CUMbbc = (1 - omega6) * (CUMbbc);
236 CUMbcb = (1 - omega6) * (CUMbcb);
237 CUMcbb = (1 - omega6) * (CUMcbb);
238
239 CUMbcc += omega7 * (-CUMbcc);
240 CUMcbc += omega7 * (-CUMcbc);
241 CUMccb += omega7 * (-CUMccb);
242
243 CUMccc += omega10 * (-CUMccc);
244
245
246 mcbb = CUMcbb + 1. / 3 * ((3 * mcaa + 1) * mabb + 6 * mbba * mbab) * inverse_rho;
247 mbcb = CUMbcb + 1. / 3 * ((3 * maca + 1) * mbab + 6 * mbba * mabb) * inverse_rho;
248 mbbc = CUMbbc + 1. / 3 * ((3 * maac + 1) * mbba + 6 * mbab * mabb) * inverse_rho;
249
250 mcca = CUMcca + (((mcaa * maca + 2 * mbba * mbba) * 9 + 3 * (mcaa + maca)) * inverse_rho - (drho * inverse_rho)) * 1. / 9;
251 mcac = CUMcac + (((mcaa * maac + 2 * mbab * mbab) * 9 + 3 * (mcaa + maac)) * inverse_rho - (drho * inverse_rho)) * 1. / 9;
252 macc = CUMacc + (((maac * maca + 2 * mabb * mabb) * 9 + 3 * (maac + maca)) * inverse_rho - (drho * inverse_rho)) * 1. / 9;
253
254 mbcc = CUMbcc + 1. / 3 * (3 * (maac * mbca + maca * mbac + 4 * mabb * mbbb + 2 * (mbab * macb + mbba * mabc)) + (mbca + mbac)) * inverse_rho;
255 mcbc = CUMcbc + 1. / 3 * (3 * (maac * mcba + mcaa * mabc + 4 * mbab * mbbb + 2 * (mabb * mcab + mbba * mbac)) + (mcba + mabc)) * inverse_rho;
256 mccb = CUMccb + 1. / 3 * (3 * (mcaa * macb + maca * mcab + 4 * mbba * mbbb + 2 * (mbab * mbca + mabb * mcba)) + (macb + mcab)) * inverse_rho;
257
258 mccc = CUMccc - ((-4 * mbbb * mbbb - (mcaa * macc + maca * mcac + maac * mcca)
259 - 4 * (mabb * mcbb + mbab * mbcb + mbba * mbbc)
260 - 2 * (mbca * mbac + mcba * mabc + mcab * macb)) * inverse_rho
261 + (4 * (mbab * mbab * maca + mabb * mabb * mcaa + mbba * mbba * maac)
262 + 2 * (mcaa * maca * maac) + 16 * mbba * mbab * mabb) * inverse_rho * inverse_rho
263 - 1. / 9 * (macc + mcac + mcca) * inverse_rho
264 - 1. / 9 * (mcaa + maca + maac) * inverse_rho
265 + (2 * (mbab * mbab + mabb * mabb + mbba * mbba) + (maac * maca + maac * mcaa + maca * mcaa)
266 + 1. / 3 * (maac + maca + mcaa)) * inverse_rho * inverse_rho * 2. / 3
267 + 1. / 27 * ((drho * drho - drho) * inverse_rho * inverse_rho));
268
269
270 mbaa = -mbaa;
271 maba = -maba;
272 maab = -maab;
273
274
276
277 return uSqr;
278 }
static void computePopulations(MOMENTA &momenta, CELL &cell, U &u) any_platform
Backwards Chimera Transformation to compute population distribution from the central moments FIXME: S...
static void computeMomenta(MOMENTA &momenta, CELL &cell, U &u) any_platform
Uses the Chimera Transformation to compute central moments from the population distribution.