Loading [MathJax]/extensions/tex2jax.js
MAIA bb96820c
Multiphysics at AIA
•All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Pages
lbfunctions.h
Go to the documentation of this file.
1// Copyright (C) 2024 The m-AIA AUTHORS
2//
3// This file is part of m-AIA (https://git.rwth-aachen.de/aia/m-AIA/m-AIA)
4//
5// SPDX-License-Identifier: LGPL-3.0-only
6
7#ifndef LBFUNCTIONS_H
8#define LBFUNCTIONS_H
9#include <array>
10#include <numeric>
11#include "lbconstants.h"
12#include "lblatticedescriptor.h"
13
14// Within this namespace inline functions are provided for several LB
15// calculations given as templates
16// Definitions must be independent from anything besides lbconstants.h.
17namespace lbfunc {
18using namespace lbconstants;
19
28template <MInt nDim, MInt nDist>
29inline void calcEqDistsIncompressible(const MFloat& p_rho, const MFloat& p_squaredVelocity, MFloat const* const p_u,
30 MFloat* const eqDist) {
31 using Ld = LbLatticeDescriptor<nDim, nDist>;
32 const MFloat squaredVelocityB2 = p_squaredVelocity * F1B2;
33
34 const MFloat rhomulCSsq = p_rho * CSsq;
35
36 const MFloat preTerm = rhomulCSsq - squaredVelocityB2;
37
38 MFloat b[2 * nDim];
39 for(MInt n = 0; n < nDim; n++) {
40 b[2 * n] = -p_u[n];
41 b[2 * n + 1] = p_u[n];
42 }
43
44 // Calculation of distributions for directions with only one component
45 constexpr MFloat lb_tp_coef1 = Ld::tp(1) * F1BCSsq;
46 for(MInt j = 0; j < Ld::distFld(0); j++) {
47 eqDist[j] = lb_tp_coef1 * (preTerm + b[j] + b[j] * b[j] * F1B2mulF1BCSsq);
48 }
49
50 // Calculation of distributions for directions with two components
51 constexpr MFloat lb_tp_coef2 = Ld::tp(2) * F1BCSsq;
52 for(MInt j = 0; j < Ld::distFld(1); j++) {
53 const MFloat tmp = (b[Ld::mFld1(2 * j)] + b[Ld::mFld1(2 * j + 1)]);
54 eqDist[Ld::distFld(0) + j] = lb_tp_coef2 * (preTerm + tmp + tmp * tmp * F1B2mulF1BCSsq);
55 }
56
57 // Calculation of distributions for directions with three components
58 constexpr MFloat lb_tp_coef3 = Ld::tp(3) * F1BCSsq;
59 for(MInt j = 0; j < Ld::distFld(2); j++) {
60 const MFloat tmp = (b[Ld::mFld2(3 * j)] + b[Ld::mFld2(3 * j + 1)] + b[Ld::mFld2(3 * j + 2)]);
61 eqDist[Ld::distFld(0) + Ld::distFld(1) + j] = lb_tp_coef3 * (preTerm + tmp + tmp * tmp * F1B2mulF1BCSsq);
62 }
63
64 // Calculation of distribution for rest particle distribution (center)
65 eqDist[Ld::lastId()] = Ld::tp(0) * (p_rho - F1B2mulF1BCSsq * p_squaredVelocity);
66}
67
76template <MInt nDim, MInt nDist>
77inline void calcEqDistsCompressible(const MFloat& p_rho, const MFloat& p_squaredVelocity, MFloat const* const p_u,
78 MFloat* const eqDist) {
79 using Ld = LbLatticeDescriptor<nDim, nDist>;
80 const MFloat squaredVelocityB2 = p_squaredVelocity * F1B2;
81
82 MFloat b[2 * nDim];
83 for(MInt n = 0; n < nDim; n++) {
84 b[2 * n] = -p_u[n];
85 b[2 * n + 1] = p_u[n];
86 }
87
88 // Calculation of distributions for directions with only one component
89 constexpr MFloat lb_tp_coef1 = Ld::tp(1) * F1BCSsq;
90 for(MInt j = 0; j < Ld::distFld(0); j++) {
91 eqDist[j] = lb_tp_coef1 * (CSsq + b[j] + b[j] * b[j] * F1B2mulF1BCSsq - squaredVelocityB2) * p_rho;
92 }
93
94 // Calculation of distributions for directions with two components
95 constexpr MFloat lb_tp_coef2 = Ld::tp(2) * F1BCSsq;
96 for(MInt j = 0; j < Ld::distFld(1); j++) {
97 const MFloat tmp = (b[Ld::mFld1(2 * j)] + b[Ld::mFld1(2 * j + 1)]);
98 eqDist[Ld::distFld(0) + j] = lb_tp_coef2 * (CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2) * p_rho;
99 }
100
101 // Calculation of distributions for directions with three components
102 constexpr MFloat lb_tp_coef3 = Ld::tp(3) * F1BCSsq;
103 for(MInt j = 0; j < Ld::distFld(2); j++) {
104 const MFloat tmp = (b[Ld::mFld2(3 * j)] + b[Ld::mFld2(3 * j + 1)] + b[Ld::mFld2(3 * j + 2)]);
105 eqDist[Ld::distFld(0) + Ld::distFld(1) + j] =
106 lb_tp_coef3 * (CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2) * p_rho;
107 }
108
109 // Calculation of distribution for rest particle distribution (center)
110 eqDist[Ld::lastId()] = Ld::tp(0) * (1.0 - F1B2mulF1BCSsq * p_squaredVelocity) * p_rho;
111}
112
121template <MInt nDim, MInt nDist, MBool compressible = false>
122inline void calcEqDists(const MFloat& p_rho, const MFloat& p_squaredVelocity, MFloat const* const p_u,
123 MFloat* const eqDist) {
124 if constexpr(compressible)
125 calcEqDistsCompressible<nDim, nDist>(p_rho, p_squaredVelocity, p_u, eqDist);
126 else
127 calcEqDistsIncompressible<nDim, nDist>(p_rho, p_squaredVelocity, p_u, eqDist);
128}
129
137template <MInt nDim, MInt nDist, MBool compressible = false>
138inline void calcEqDists(const MFloat& p_rho, MFloat const* const p_u, MFloat* const eqDist) {
139 const MFloat p_squaredVelocity = std::inner_product(&p_u[0], &p_u[nDim], &p_u[0], .0);
140
141 calcEqDists<nDim, nDist, compressible>(p_rho, p_squaredVelocity, p_u, eqDist);
142}
143
144#ifdef WAR_NVHPC_PSTL
151template <MInt nDim, MInt nDist>
152inline void calcEqDistsIncompressible(const MFloat& p_rho, const MFloat& p_squaredVelocity, MFloat const* const p_u,
153 MFloat* const eqDist, const MInt* myMFld1, const MInt* myMFld2,
154 const MFloat* myTp, const MInt* myDistFld) {
155 using Ld = LbLatticeDescriptor<nDim, nDist>;
156 const MFloat squaredVelocityB2 = p_squaredVelocity * F1B2;
157
158 const MFloat rhomulCSsq = p_rho * CSsq;
159
160 MFloat b[2 * nDim];
161 for(MInt n = 0; n < nDim; n++) {
162 b[2 * n] = -p_u[n];
163 b[2 * n + 1] = p_u[n];
164 }
165
166 // Calculation of distributions for directions with only one component
167 const MFloat lb_tp_coef1 = myTp[1] * F1BCSsq;
168 for(MInt j = 0; j < myDistFld[0]; j++) {
169 eqDist[j] = lb_tp_coef1 * (rhomulCSsq + b[j] + b[j] * b[j] * F1B2mulF1BCSsq - squaredVelocityB2);
170 }
171
172 // Calculation of distributions for directions with two components
173 const MFloat lb_tp_coef2 = myTp[2] * F1BCSsq;
174 for(MInt j = 0; j < myDistFld[1]; j++) {
175 const MFloat tmp = (b[myMFld1[2 * j]] + b[myMFld1[2 * j + 1]]);
176 eqDist[myDistFld[0] + j] = lb_tp_coef2 * (rhomulCSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2);
177 }
178
179 // Calculation of distributions for directions with three components
180 const MFloat lb_tp_coef3 = myTp[3] * F1BCSsq;
181 for(MInt j = 0; j < myDistFld[2]; j++) {
182 const MFloat tmp = (b[myMFld2[3 * j]] + b[myMFld2[3 * j + 1]] + b[myMFld2[3 * j + 2]]);
183 eqDist[myDistFld[0] + myDistFld[1] + j] =
184 lb_tp_coef3 * (rhomulCSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2);
185 }
186
187 // Calculation of distribution for rest particle distribution (center)
188 eqDist[Ld::lastId()] = myTp[0] * (p_rho - F1B2mulF1BCSsq * p_squaredVelocity);
189}
190
197template <MInt nDim, MInt nDist>
198inline void calcEqDistsCompressible(const MFloat& p_rho, const MFloat& p_squaredVelocity, MFloat const* const p_u,
199 MFloat* const eqDist, const MInt* myMFld1, const MInt* myMFld2, const MFloat* myTp,
200 const MInt* myDistFld) {
201 using Ld = LbLatticeDescriptor<nDim, nDist>;
202 const MFloat squaredVelocityB2 = p_squaredVelocity * F1B2;
203
204 MFloat b[2 * nDim];
205 for(MInt n = 0; n < nDim; n++) {
206 b[2 * n] = -p_u[n];
207 b[2 * n + 1] = p_u[n];
208 }
209
210 // Calculation of distributions for directions with only one component
211 const MFloat lb_tp_coef1 = myTp[1] * F1BCSsq;
212 for(MInt j = 0; j < myDistFld[0]; j++) {
213 eqDist[j] = lb_tp_coef1 * (CSsq + b[j] + b[j] * b[j] * F1B2mulF1BCSsq - squaredVelocityB2) * p_rho;
214 }
215
216 // Calculation of distributions for directions with two components
217 const MFloat lb_tp_coef2 = myTp[2] * F1BCSsq;
218 for(MInt j = 0; j < myDistFld[1]; j++) {
219 const MFloat tmp = (b[myMFld1[2 * j]] + b[myMFld1[2 * j + 1]]);
220 eqDist[myDistFld[0] + j] = lb_tp_coef2 * (CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2) * p_rho;
221 }
222
223 // Calculation of distributions for directions with three components
224 const MFloat lb_tp_coef3 = myTp[3] * F1BCSsq;
225 for(MInt j = 0; j < myDistFld[2]; j++) {
226 const MFloat tmp = (b[myMFld2[3 * j]] + b[myMFld2[3 * j + 1]] + b[myMFld2[3 * j + 2]]);
227 eqDist[myDistFld[0] + myDistFld[1] + j] =
228 lb_tp_coef3 * (CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2) * p_rho;
229 }
230
231 // Calculation of distribution for rest particle distribution (center)
232 eqDist[Ld::lastId()] = myTp[0] * (1.0 - F1B2mulF1BCSsq * p_squaredVelocity) * p_rho;
233}
234
243template <MInt nDim, MInt nDist, MBool compressible = false>
244inline void calcEqDists(const MFloat& p_rho, const MFloat& p_squaredVelocity, MFloat const* const p_u,
245 MFloat* const eqDist, const MInt* myMFld1, const MInt* myMFld2, const MFloat* myTp,
246 const MInt* myDistFld) {
247 if constexpr(compressible)
248 calcEqDistsCompressible<nDim, nDist>(p_rho, p_squaredVelocity, p_u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
249 else
250 calcEqDistsIncompressible<nDim, nDist>(p_rho, p_squaredVelocity, p_u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
251}
252
260template <MInt nDim, MInt nDist, MBool compressible = false>
261inline void calcEqDists(const MFloat& p_rho, MFloat const* const p_u, MFloat* const eqDist, const MInt* myMFld1,
262 const MInt* myMFld2, const MFloat* myTp, const MInt* myDistFld) {
263 const MFloat p_squaredVelocity = std::inner_product(&p_u[0], &p_u[nDim], &p_u[0], .0);
264
265 calcEqDists<nDim, nDist, compressible>(p_rho, p_squaredVelocity, p_u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
266}
267#endif
268
281template <MInt nDim, MInt nDist, MBool compressible = false>
282inline void calcNonEqDists(const MFloat& rho, const MFloat* const u, const MFloat* const dist,
283 MFloat* const nonEqDist) {
284 // Calculate equlibirum distribution
285 std::array<MFloat, nDist> eqDist{};
286 calcEqDists<nDim, nDist, compressible>(rho, u, eqDist.data());
287
288 // Calculate non-equilibium part
289 for(MInt j = 0; j < nDist; j++) {
290 nonEqDist[j] = dist[j] - eqDist[j];
291 }
292}
293
304template <MInt nDim, MInt nDist, MBool compressible = false>
305inline void calcNonEqDists(const MFloat* const dist, MFloat* const nonEqDist) {
306 // Calculate macroscopic variables
307 MFloat rho{};
308 std::array<MFloat, nDim> u{};
309 calcMacroscopicVars(dist, rho, u.data());
310
311 calcNonEqDists<nDim, nDist, compressible>(rho, u.data(), dist, nonEqDist);
312}
313
324template <MInt nDim, MInt nDist>
325inline void calcMomentumFluxFromNonEq(const MFloat* const nonEqDist, MFloat* const momentumFlux) {
326 using Ld = LbLatticeDescriptor<nDim, nDist>;
327
328 for(MInt d = 0; d < nDim * nDim; d++) {
329 momentumFlux[d] = F0;
330 }
331 for(MInt j = 0; j < nDist; j++) {
332 for(MInt k = 0; k < nDim; k++) {
333 for(MInt l = 0; l < nDim; l++) {
334 momentumFlux[k * nDim + l] += nonEqDist[j] * Ld::ppdfDir(j, k) * Ld::ppdfDir(j, l);
335 }
336 }
337 }
338}
339
352template <MInt nDim, MInt nDist, MBool compressible = false>
353inline void calcMomentumFlux(const MFloat& rho, const MFloat* const u, const MFloat* const dist,
354 MFloat* const momentumFlux) {
355 std::array<MFloat, nDist> nonEqDist{};
356 calcNonEqDists<nDim, nDist, compressible>(rho, u, dist, nonEqDist.data());
357
358 calcMomentumFluxFromNonEq<nDim, nDist>(nonEqDist.data(), momentumFlux);
359}
360
371template <MInt nDim, MInt nDist, MBool compressible = false>
372inline void calcMomentumFlux(const MFloat* const dist, MFloat* const momentumFlux) {
373 std::array<MFloat, nDist> nonEqDist{};
374 calcNonEqDists<nDim, nDist, compressible>(dist, nonEqDist.data());
375
376 calcMomentumFluxFromNonEq<nDim, nDist>(nonEqDist.data(), momentumFlux);
377}
378
392template <MInt nDim, MInt nDist>
393inline void calcEqDistsThermal(const MFloat& T, const MFloat squaredVelocity, const MFloat* const u,
394 MFloat* const eqDist) {
395 calcEqDists<nDim, nDist, true>(T, squaredVelocity, u, eqDist);
396}
397
409template <MInt nDim, MInt nDist>
410inline void calcEqDistsThermal(const MFloat& T, const MFloat* const u, MFloat* const eqDist) {
411 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
412
413 calcEqDists<nDim, nDist, true>(T, squaredVelocity, u, eqDist);
414}
415
416#ifdef WAR_NVHPC_PSTL
430template <MInt nDim, MInt nDist>
431inline void calcEqDistsThermal(const MFloat& T, const MFloat squaredVelocity, const MFloat* const u,
432 MFloat* const eqDist, const MInt* myMFld1, const MInt* myMFld2, const MFloat* myTp,
433 const MInt* myDistFld) {
434 calcEqDists<nDim, nDist, true>(T, squaredVelocity, u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
435}
436
448template <MInt nDim, MInt nDist>
449inline void calcEqDistsThermal(const MFloat& T, const MFloat* const u, MFloat* const eqDist, const MInt* myMFld1,
450 const MInt* myMFld2, const MFloat* myTp, const MInt* myDistFld) {
451 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
452
453 calcEqDists<nDim, nDist, true>(T, squaredVelocity, u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
454}
455#endif
456
457template <MInt nDim, MInt nDist>
458inline void calcEqDistsInnerEnergy(const MFloat& T, const MFloat& rho, const MFloat squaredVelocity,
459 const MFloat* const u, MFloat* const eqDist) {
460 using Ld = LbLatticeDescriptor<nDim, nDist>;
461
462 const MFloat F1BD = F1 / nDim;
463 const MFloat squaredVelocityB2 = 0.5 * squaredVelocity;
464 const MFloat innerEnergy = rho * T * F1B2 * (MFloat)nDim;
465
466 std::array<MFloat, 2 * nDim> b{};
467 for(MInt d = 0; d < nDim; d++) {
468 b[2 * d] = -u[d];
469 b[2 * d + 1] = u[d];
470 }
471
472 // Calculate all equilibrium distributions
473 // Calculation of eq. distributions for directions with only one component
474 for(MInt j = 0; j < Ld::distFld(0); j++) {
475 const MFloat l_innerterm_ie =
476 F1BCSsq * F1BD * CSsq + (F1BCSsq * F1BD - F2 * F1BD) * b[j] + b[j] * b[j] * F1B2 * F1BCSsq - squaredVelocityB2;
477 eqDist[j] = Ld::tp(1) * F1BCSsq * innerEnergy * l_innerterm_ie;
478 }
479
480 // Calculation of eq. distributions for directions with two components
481 for(MInt j = 0; j < Ld::distFld(1); j++) {
482 const MFloat tmp = (b[Ld::mFld1(2 * j)] + b[Ld::mFld1(2 * j + 1)]);
483 const MFloat l_innerterm_ie = F2 * F1BCSsq * F1BD * CSsq + (F2 * F1BCSsq * F1BD - F2 * F1BD) * tmp
484 + tmp * tmp * F1B2 * F1BCSsq - squaredVelocityB2;
485 eqDist[Ld::distFld(0) + j] = Ld::tp(2) * F1BCSsq * innerEnergy * l_innerterm_ie;
486 }
487
488 // Calculation of eq. distributions for directions with three components
489 for(MInt j = 0; j < Ld::distFld(2); j++) {
490 const MFloat tmp = (b[Ld::mFld2(3 * j)] + b[Ld::mFld2(3 * j + 1)] + b[Ld::mFld2(3 * j + 2)]);
491 const MFloat l_innerterm_ie = F3 * F1BCSsq * F1BD * CSsq + (F3 * F1BCSsq * F1BD - F2 * F1BD) * tmp
492 + tmp * tmp * F1B2 * F1BCSsq - squaredVelocityB2;
493 eqDist[Ld::distFld(0) + Ld::distFld(1) + j] = Ld::tp(3) * F1BCSsq * innerEnergy * l_innerterm_ie;
494 }
495
496 // Rest distribution
497 eqDist[Ld::lastId()] = -Ld::tp(0) * innerEnergy * F1B2 * F1BCSsq * squaredVelocity;
498}
499
500template <MInt nDim, MInt nDist>
501inline void calcEqDistsInnerEnergy(const MFloat& T, const MFloat& rho, const MFloat* const u, MFloat* const eqDist) {
502 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
503
504 calcEqDistsInnerEnergy<nDim, nDist>(T, rho, squaredVelocity, u, eqDist);
505}
506
507#ifdef WAR_NVHPC_PSTL
508template <MInt nDim, MInt nDist>
509inline void calcEqDistsInnerEnergy(const MFloat& T, const MFloat& rho, const MFloat squaredVelocity,
510 const MFloat* const u, MFloat* const eqDist, const MInt* myMFld1,
511 const MInt* myMFld2, const MFloat* myTp, const MInt* myDistFld) {
512 using Ld = LbLatticeDescriptor<nDim, nDist>;
513
514 const MFloat F1BD = F1 / nDim;
515 const MFloat squaredVelocityB2 = 0.5 * squaredVelocity;
516 const MFloat innerEnergy = rho * T * F1B2 * nDim;
517
518 std::array<MFloat, 2 * nDim> b{};
519 for(MInt d = 0; d < nDim; d++) {
520 b[2 * d] = -u[d];
521 b[2 * d + 1] = u[d];
522 }
523
524 // Calculate all equilibrium distributions
525 // Calculation of eq. distributions for directions with only one component
526 const MFloat lb_tp_coef1 = myTp[1] * F1BCSsq;
527 for(MInt j = 0; j < myDistFld[0]; j++) {
528 const MFloat l_innerterm_ie =
529 F1BCSsq * F1BD * CSsq + (F1BCSsq * F1BD - F2 * F1BD) * b[j] + b[j] * b[j] * F1B2 * F1BCSsq - squaredVelocityB2;
530 eqDist[j] = lb_tp_coef1 * innerEnergy * l_innerterm_ie;
531 }
532
533 // Calculation of eq. distributions for directions with two components
534 const MFloat lb_tp_coef2 = myTp[2] * F1BCSsq;
535 for(MInt j = 0; j < myDistFld[1]; j++) {
536 const MFloat tmp = (b[myMFld1[2 * j]] + b[myMFld1[2 * j + 1]]);
537 const MFloat l_innerterm_ie = F2 * F1BCSsq * F1BD * CSsq + (F2 * F1BCSsq * F1BD - F2 * F1BD) * tmp
538 + tmp * tmp * F1B2 * F1BCSsq - squaredVelocityB2;
539 eqDist[myDistFld[0] + j] = lb_tp_coef2 * innerEnergy * l_innerterm_ie;
540 }
541
542 // Calculation of eq. distributions for directions with three components
543 const MFloat lb_tp_coef3 = myTp[3] * F1BCSsq;
544 for(MInt j = 0; j < myDistFld[2]; j++) {
545 const MFloat tmp = (b[myMFld2[3 * j]] + b[myMFld2[3 * j + 1]] + b[myMFld2[3 * j + 2]]);
546 const MFloat l_innerterm_ie = F3 * F1BCSsq * F1BD * CSsq + (F3 * F1BCSsq * F1BD - F2 * F1BD) * tmp
547 + tmp * tmp * F1B2 * F1BCSsq - squaredVelocityB2;
548 eqDist[myDistFld[0] + myDistFld[1] + j] = lb_tp_coef3 * innerEnergy * l_innerterm_ie;
549 }
550
551 // Rest distribution
552 eqDist[Ld::lastId()] = -myTp[0] * innerEnergy * F1B2 * F1BCSsq * squaredVelocity;
553}
554
555template <MInt nDim, MInt nDist>
556inline void calcEqDistsInnerEnergy(const MFloat& T, const MFloat& rho, const MFloat* const u, MFloat* const eqDist,
557 const MInt* myMFld1, const MInt* myMFld2, const MFloat* myTp,
558 const MInt* myDistFld) {
559 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
560
561 calcEqDistsInnerEnergy<nDim, nDist>(T, rho, squaredVelocity, u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
562}
563#endif
564
565template <MInt nDim, MInt nDist>
566inline void calcEqDistsTotalEnergy(const MFloat& T, const MFloat& rho, const MFloat squaredVelocity,
567 const MFloat* const u, MFloat* const eqDistT) {
568 using Ld = LbLatticeDescriptor<nDim, nDist>;
569
570 const MFloat squaredVelocityB2 = 0.5 * squaredVelocity;
571 const MFloat p0 = rho * CSsq;
572 const MFloat totalEnergy = T * nDim * F1B2 + squaredVelocityB2;
573
574 std::array<MFloat, 2 * nDim> b{};
575 for(MInt d = 0; d < nDim; d++) {
576 b[d * 2] = -u[d];
577 b[d * 2 + 1] = u[d];
578 }
579
580 // Calculate all equilibrium distributions
581 // Calculation of eq. distributions for directions with only one component
582 constexpr MFloat lb_tp_coef1 = Ld::tp(1) * F1BCSsq;
583 for(MInt j = 0; j < Ld::distFld(0); j++) {
584 const MFloat l_innerterm = CSsq + b[j] + b[j] * b[j] * F1B2mulF1BCSsq - squaredVelocityB2;
585 const MFloat eq_dist = lb_tp_coef1 * rho * l_innerterm;
586 const MFloat l_innerterm_te = b[j] + b[j] * b[j] * F1BCSsq - squaredVelocity + F1B2 * (F1 - nDim * CSsq);
587 eqDistT[j] = lb_tp_coef1 * p0 * l_innerterm_te + totalEnergy * eq_dist;
588 }
589
590 // Calculation of eq. distributions for directions with two components
591 constexpr MFloat lb_tp_coef2 = Ld::tp(2) * F1BCSsq;
592 for(MInt j = 0; j < Ld::distFld(1); j++) {
593 const MFloat tmp = (b[Ld::mFld1(2 * j)] + b[Ld::mFld1(2 * j + 1)]);
594 const MFloat l_innerterm = CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2;
595 const MFloat eq_dist = lb_tp_coef2 * rho * l_innerterm;
596 const MFloat l_innerterm_te = tmp + tmp * tmp * F1BCSsq - squaredVelocity + F1B2 * (F2 - nDim * CSsq);
597 eqDistT[Ld::distFld(0) + j] = lb_tp_coef2 * p0 * l_innerterm_te + totalEnergy * eq_dist;
598 }
599
600 // Calculation of eq. distributions for directions with three components
601 constexpr MFloat lb_tp_coef3 = Ld::tp(3) * F1BCSsq;
602 for(MInt j = 0; j < Ld::distFld(2); j++) {
603 const MFloat tmp = (b[Ld::mFld2(3 * j)] + b[Ld::mFld2(3 * j + 1)] + b[Ld::mFld2(3 * j + 2)]);
604 const MFloat l_innerterm = CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2;
605 const MFloat eq_dist = lb_tp_coef3 * rho * l_innerterm;
606 const MFloat l_innerterm_te = tmp + tmp * tmp * F1BCSsq - squaredVelocity + F1B2 * (F3 - nDim * CSsq);
607 eqDistT[Ld::distFld(0) + Ld::distFld(1) + j] = lb_tp_coef3 * p0 * l_innerterm_te + totalEnergy * eq_dist;
608 }
609
610 // Rest distribution
611 constexpr MFloat lb_tp_coef0 = Ld::tp(0);
612 const MFloat l_innerterm = F1 - F1B2mulF1BCSsq * squaredVelocity;
613 const MFloat eq_dist = lb_tp_coef0 * rho * l_innerterm;
614 const MFloat l_innerterm_te = -F1BCSsq * squaredVelocity + F1B2mulF1BCSsq * (F0 - nDim * CSsq);
615 eqDistT[Ld::lastId()] = lb_tp_coef0 * p0 * l_innerterm_te + totalEnergy * eq_dist;
616}
617
618template <MInt nDim, MInt nDist>
619inline void calcEqDistsTotalEnergy(const MFloat& T, const MFloat& rho, const MFloat* const u, MFloat* const eqDist) {
620 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
621
622 calcEqDistsTotalEnergy<nDim, nDist>(T, rho, squaredVelocity, u, eqDist);
623}
624
625#ifdef WAR_NVHPC_PSTL
626template <MInt nDim, MInt nDist>
627inline void calcEqDistsTotalEnergy(const MFloat& T, const MFloat& rho, const MFloat squaredVelocity,
628 const MFloat* const u, MFloat* const eqDistT, const MInt* myMFld1,
629 const MInt* myMFld2, const MFloat* myTp, const MInt* myDistFld) {
630 using Ld = LbLatticeDescriptor<nDim, nDist>;
631
632 const MFloat squaredVelocityB2 = 0.5 * squaredVelocity;
633 const MFloat p0 = rho * CSsq;
634 const MFloat totalEnergy = T * nDim * F1B2 + squaredVelocityB2;
635
636 std::array<MFloat, 2 * nDim> b{};
637 for(MInt d = 0; d < nDim; d++) {
638 b[2 * d] = -u[d];
639 b[2 * d + 1] = u[d];
640 }
641
642 // Calculate all equilibrium distributions
643 // Calculation of eq. distributions for directions with only one component
644 const MFloat lb_tp_coef1 = myTp[1] * F1BCSsq;
645 for(MInt j = 0; j < myDistFld[0]; j++) {
646 const MFloat l_innerterm = CSsq + b[j] + b[j] * b[j] * F1B2mulF1BCSsq - squaredVelocityB2;
647 const MFloat eq_dist = lb_tp_coef1 * rho * l_innerterm;
648 const MFloat l_innerterm_te = b[j] + b[j] * b[j] * F1BCSsq - squaredVelocity + F1B2 * (F1 - nDim * CSsq);
649 eqDistT[j] = lb_tp_coef1 * p0 * l_innerterm_te + totalEnergy * eq_dist;
650 }
651
652 // Calculation of eq. distributions for directions with two components
653 const MFloat lb_tp_coef2 = myTp[2] * F1BCSsq;
654 for(MInt j = 0; j < myDistFld[1]; j++) {
655 const MFloat tmp = (b[myMFld1[2 * j]] + b[myMFld1[2 * j + 1]]);
656 const MFloat l_innerterm = CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2;
657 const MFloat eq_dist = lb_tp_coef2 * rho * l_innerterm;
658 const MFloat l_innerterm_te = tmp + tmp * tmp * F1BCSsq - squaredVelocity + F1B2 * (F2 - nDim * CSsq);
659 eqDistT[myDistFld[0] + j] = lb_tp_coef2 * p0 * l_innerterm_te + totalEnergy * eq_dist;
660 }
661
662 // Calculation of eq. distributions for directions with three components
663 const MFloat lb_tp_coef3 = myTp[3] * F1BCSsq;
664 for(MInt j = 0; j < myDistFld[2]; j++) {
665 const MFloat tmp = (b[myMFld2[3 * j]] + b[myMFld2[3 * j + 1]] + b[myMFld2[3 * j + 2]]);
666 const MFloat l_innerterm = CSsq + tmp + tmp * tmp * F1B2mulF1BCSsq - squaredVelocityB2;
667 const MFloat eq_dist = lb_tp_coef3 * rho * l_innerterm;
668 const MFloat l_innerterm_te = tmp + tmp * tmp * F1BCSsq - squaredVelocity + F1B2 * (F3 - nDim * CSsq);
669 eqDistT[myDistFld[0] + myDistFld[1] + j] = lb_tp_coef3 * p0 * l_innerterm_te + totalEnergy * eq_dist;
670 }
671
672 // Rest distribution
673 const MFloat l_innerterm = F1 - F1B2mulF1BCSsq * squaredVelocity;
674 const MFloat eq_dist = myTp[0] * rho * l_innerterm;
675 const MFloat l_innerterm_te = -F1BCSsq * squaredVelocity + F1B2mulF1BCSsq * (F0 - nDim * CSsq);
676 eqDistT[Ld::lastId()] = myTp[0] * p0 * l_innerterm_te + totalEnergy * eq_dist;
677}
678
679template <MInt nDim, MInt nDist>
680inline void calcEqDistsTotalEnergy(const MFloat& T, const MFloat& rho, const MFloat* const u, MFloat* const eqDist,
681 const MInt* myMFld1, const MInt* myMFld2, const MFloat* myTp,
682 const MInt* myDistFld) {
683 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
684
685 calcEqDistsTotalEnergy<nDim, nDist>(T, rho, squaredVelocity, u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
686}
687#endif
688
702template <MInt nDim, MInt nDist>
703inline void calcEqDistsTransport(const MFloat& C, const MFloat squaredVelocity, const MFloat* const u,
704 MFloat* const eqDist) {
705 calcEqDists<nDim, nDist, true>(C, squaredVelocity, u, eqDist);
706}
707
719template <MInt nDim, MInt nDist>
720inline void calcEqDistsTransport(const MFloat& C, const MFloat* const u, MFloat* const eqDist) {
721 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
722
723 calcEqDists<nDim, nDist, true>(C, squaredVelocity, u, eqDist);
724}
725
726#ifdef WAR_NVHPC_PSTL
740template <MInt nDim, MInt nDist>
741inline void calcEqDistsTransport(const MFloat& C, const MFloat squaredVelocity, const MFloat* const u,
742 MFloat* const eqDist, const MInt* myMFld1, const MInt* myMFld2, const MFloat* myTp,
743 const MInt* myDistFld) {
744 calcEqDists<nDim, nDist, true>(C, squaredVelocity, u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
745}
746
758template <MInt nDim, MInt nDist>
759inline void calcEqDistsTransport(const MFloat& C, const MFloat* const u, MFloat* const eqDist, const MInt* myMFld1,
760 const MInt* myMFld2, const MFloat* myTp, const MInt* myDistFld) {
761 const MFloat squaredVelocity = std::inner_product(&u[0], &u[nDim], &u[0], .0);
762
763 calcEqDists<nDim, nDist, true>(C, squaredVelocity, u, eqDist, myMFld1, myMFld2, myTp, myDistFld);
764}
765#endif
766
774template <MInt nDim, MInt nDist, MBool compressible = false>
775inline void calcMacroVars(MFloat const* const p_dist, MFloat& p_rho, MFloat* const p_u) {
776 using Ld = LbLatticeDescriptor<nDim, nDist>;
777 // density
778 p_rho = 0.0;
779 for(MInt j = 0; j < nDist; j++) {
780 p_rho += p_dist[j];
781 }
782 // velocities [rho*u]
783 for(MInt i = 0; i < nDim; i++) {
784 p_u[i] = 0.0;
785 for(MInt j = 0; j < Ld::dxQyFld(); j++) {
786 p_u[i] += p_dist[Ld::pFld(i, j)];
787 p_u[i] -= p_dist[Ld::nFld(i, j)];
788 }
789 if constexpr(compressible) p_u[i] /= p_rho;
790 }
791}
792
793#ifdef WAR_NVHPC_PSTL
801template <MInt nDim, MInt nDist, MBool compressible = false>
802inline void calcMacroVars(MFloat const* const p_dist, MFloat& p_rho, MFloat* const p_u, MInt* const myPFld,
803 MInt* const myNFld, const MInt fldlen) {
804 using Ld = LbLatticeDescriptor<nDim, nDist>;
805 // density
806 p_rho = 0.0;
807 for(MInt j = 0; j < nDist; j++) {
808 p_rho += p_dist[j];
809 }
810 // velocities [rho*u]
811 for(MInt i = 0; i < nDim; i++) {
812 p_u[i] = 0.0;
813 for(MInt j = 0; j < fldlen; j++) {
814 p_u[i] += p_dist[myPFld[i * fldlen + j]];
815 p_u[i] -= p_dist[myNFld[i * fldlen + j]];
816 }
817 if constexpr(compressible) p_u[i] /= p_rho;
818 }
819}
820#endif
821} // namespace lbfunc
822
823#endif
MInt
int32_t MInt
Definition: maiatypes.h:62
MFloat
double MFloat
Definition: maiatypes.h:52
lbfunc::calcMacroVars
void calcMacroVars(MFloat const *const p_dist, MFloat &p_rho, MFloat *const p_u)
Calculate macroscopic variables from distribution.
Definition: lbfunctions.h:775
lbfunc::calcEqDistsIncompressible
void calcEqDistsIncompressible(const MFloat &p_rho, const MFloat &p_squaredVelocity, MFloat const *const p_u, MFloat *const eqDist)
Calculate equilibrium function (incompressible version). Squared velocity has already been computed.
Definition: lbfunctions.h:29
lbfunc::calcMomentumFlux
void calcMomentumFlux(const MFloat &rho, const MFloat *const u, const MFloat *const dist, MFloat *const momentumFlux)
Calculate the momentum flux for given macroscopic variables and distribution.
Definition: lbfunctions.h:353
lbfunc::calcEqDistsTotalEnergy
void calcEqDistsTotalEnergy(const MFloat &T, const MFloat &rho, const MFloat squaredVelocity, const MFloat *const u, MFloat *const eqDistT)
Definition: lbfunctions.h:566
lbfunc::calcEqDistsThermal
void calcEqDistsThermal(const MFloat &T, const MFloat squaredVelocity, const MFloat *const u, MFloat *const eqDist)
Calculate thermal equilibrium distribution.
Definition: lbfunctions.h:393
lbfunc::calcMomentumFluxFromNonEq
void calcMomentumFluxFromNonEq(const MFloat *const nonEqDist, MFloat *const momentumFlux)
Calculate the momentum flux for a given non-equilibirum distribution.
Definition: lbfunctions.h:325
lbfunc::calcEqDists
void calcEqDists(const MFloat &p_rho, const MFloat &p_squaredVelocity, MFloat const *const p_u, MFloat *const eqDist)
Calculate equilibrium function. Squared velocity has already been computed.
Definition: lbfunctions.h:122
lbfunc::calcNonEqDists
void calcNonEqDists(const MFloat &rho, const MFloat *const u, const MFloat *const dist, MFloat *const nonEqDist)
Calculate the non-equilibrium part for given macroscopic variables and distribution.
Definition: lbfunctions.h:282
lbfunc::calcEqDistsInnerEnergy
void calcEqDistsInnerEnergy(const MFloat &T, const MFloat &rho, const MFloat squaredVelocity, const MFloat *const u, MFloat *const eqDist)
Definition: lbfunctions.h:458
lbfunc::calcEqDistsTransport
void calcEqDistsTransport(const MFloat &C, const MFloat squaredVelocity, const MFloat *const u, MFloat *const eqDist)
Calculate transport equilibrium distribution.
Definition: lbfunctions.h:703
lbfunc::calcEqDistsCompressible
void calcEqDistsCompressible(const MFloat &p_rho, const MFloat &p_squaredVelocity, MFloat const *const p_u, MFloat *const eqDist)
Calculate equilibrium function (compressible version). Squared velocity has already been computed.
Definition: lbfunctions.h:77
dist
MFloat dist(const Point< DIM > &p, const Point< DIM > &q)
Definition: pointbox.h:54
LbLatticeDescriptor
LB lattice descriptor for arrays depending on D and Q.
Definition: lblatticedescriptor.h:297