33template <
typename DESCRIPTOR>
36 template <
typename RHO,
typename U,
typename V=RHO>
40 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
41 c_u += descriptors::c<DESCRIPTOR>(iPop,iD)*u[iD];
44 * descriptors::t<V,DESCRIPTOR>(iPop)
45 * ( V{1} + c_u * descriptors::invCs2<V,DESCRIPTOR>() )
46 - descriptors::t<V,DESCRIPTOR>(iPop);
50 template <
typename RHO,
typename U,
typename USQR,
typename V=RHO>
54 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
55 c_u += descriptors::c<DESCRIPTOR>(iPop,iD)*u[iD];
58 * descriptors::t<V,DESCRIPTOR>(iPop)
60 + descriptors::invCs2<V,DESCRIPTOR>() * c_u
61 + descriptors::invCs2<V,DESCRIPTOR>() * descriptors::invCs2<V,DESCRIPTOR>() * V{0.5} * c_u * c_u
62 - descriptors::invCs2<V,DESCRIPTOR>() * V{0.5} * uSqr)
63 - descriptors::t<V,DESCRIPTOR>(iPop);
66 template <
typename RHO,
typename U,
typename V=RHO>
69 const V uSqr = util::normSqr<U,DESCRIPTOR::d>(u);
74 template <
typename RHO,
typename U,
typename V=RHO>
79 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
80 c_u += descriptors::c<DESCRIPTOR>(iPop,iD)*u[iD];
82 return descriptors::t<V,DESCRIPTOR>(iPop) * (rho + c_u)
83 - descriptors::t<V,DESCRIPTOR>(iPop);
86 template <
typename J,
typename JSQR,
typename PRESSURE,
typename V=PRESSURE>
90 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
91 c_j += descriptors::c<DESCRIPTOR>(iPop,iD)*j[iD];
93 return descriptors::t<V,DESCRIPTOR>(iPop)
94 * ( descriptors::invCs2<V,DESCRIPTOR>() * pressure
95 + descriptors::invCs2<V,DESCRIPTOR>() * c_j
96 + descriptors::invCs2<V,DESCRIPTOR>() * descriptors::invCs2<V,DESCRIPTOR>()/V{2} * c_j * c_j
97 - descriptors::invCs2<V,DESCRIPTOR>()/V{2} * jSqr )
98 - descriptors::t<V,DESCRIPTOR>(iPop);
100 template <
typename J,
typename PRESSURE,
typename V=PRESSURE>
103 const V jSqr = util::normSqr<J,DESCRIPTOR::d>(j);
109 template <
typename V>
111 const V Jgrad[DESCRIPTOR::d * DESCRIPTOR::d],
114 using L = DESCRIPTOR;
117 for (
int iAlpha=0; iAlpha < L::d; ++iAlpha) {
118 for (
int iBeta=0; iBeta < L::d; ++iBeta) {
119 V toAdd = descriptors::c<L>(iPop,iAlpha) * descriptors::c<L>(iPop,iBeta);
120 if (iAlpha == iBeta) {
121 toAdd -= 1./descriptors::invCs2<V,L>();
126 toAdd *= Jgrad[iJgrad++];
130 fNeq *= - descriptors::t<V,L>(iPop) * descriptors::invCs2<V,L>() / omega;
142 template <
typename V,
typename PI>
145 using L = DESCRIPTOR;
150 for (
int iAlpha=0; iAlpha < L::d; ++iAlpha) {
151 for (
int iBeta=iAlpha; iBeta < L::d; ++iBeta) {
152 V toAdd = descriptors::c<L>(iPop,iAlpha)*descriptors::c<L>(iPop,iBeta);
153 if (iAlpha == iBeta) {
154 toAdd -= V{1}/descriptors::invCs2<V,L>();
163 fNeq *= descriptors::t<V,L>(iPop) * descriptors::invCs2<V,L>() * descriptors::invCs2<V,L>() / V{2};
167 template <
typename V>
171 for (
int iD = 0; iD < DESCRIPTOR::d; ++iD) {
172 fNeq += descriptors::c<DESCRIPTOR>(iPop,iD) * jNeq[iD];
174 fNeq *= descriptors::t<V,DESCRIPTOR>(iPop) * descriptors::invCs2<V,DESCRIPTOR>();
181template <
typename DESCRIPTOR>
184 template <
typename CELL,
typename V=
typename CELL::value_t>
188 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
196 template <
typename CELL,
typename J,
typename V=
typename CELL::value_t>
199 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
202 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
203 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
204 j[iD] += cell[iPop]*descriptors::c<DESCRIPTOR>(iPop,iD);
210 template <
typename CELL,
typename RHO,
typename J,
typename V=
typename CELL::value_t>
218 template <
typename CELL,
typename RHO,
typename U,
typename V=
typename CELL::value_t>
222 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
228 template <
typename CELL,
typename RHO,
typename U,
typename PI,
typename V=
typename CELL::value_t>
232 for (
int iAlpha=0; iAlpha < DESCRIPTOR::d; ++iAlpha) {
233 for (
int iBeta=iAlpha; iBeta < DESCRIPTOR::d; ++iBeta) {
235 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
236 pi[iPi] += descriptors::c<DESCRIPTOR>(iPop,iAlpha)*
237 descriptors::c<DESCRIPTOR>(iPop,iBeta) * cell[iPop];
240 pi[iPi] -= rho*u[iAlpha]*u[iBeta];
242 pi[iPi] -= V{1} / descriptors::invCs2<V,DESCRIPTOR>()*(rho-V{1});
250 template <
typename CELL,
typename RHO,
typename U,
typename PI,
typename V=
typename CELL::value_t>
257 template <
typename CELL,
typename FEQ,
typename V=
typename CELL::value_t>
261 V u[DESCRIPTOR::d] {};
263 const V uSqr = util::normSqr<V,DESCRIPTOR::d>(u);
264 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
270 template <
typename CELL,
typename FNEQ,
typename RHO,
typename U,
typename V=
typename CELL::value_t>
273 const V uSqr = util::normSqr<U,DESCRIPTOR::d>(u);
274 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
279 template <
typename CELL,
typename FNEQ,
typename V=
typename CELL::value_t>
283 V u[DESCRIPTOR::d] {};
289 template <
typename CELL,
typename RHO,
typename VELOCITY,
typename OMEGA,
typename V=
typename CELL::value_t>
292 const V uSqr = util::normSqr<VELOCITY,DESCRIPTOR::d>(u);
293 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
294 cell[iPop] *= V{1} - omega;
301 template <
typename CELL,
typename RHO,
typename VELOCITY,
typename OMEGA,
typename V=
typename CELL::value_t>
304 const V uSqr = util::normSqr<VELOCITY,DESCRIPTOR::d>(u);
305 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
306 cell[iPop] *= V{1} - omega;
313 template <
typename CELL,
typename PRESSURE,
typename J,
typename OMEGA,
typename V=
typename CELL::value_t>
316 const V jSqr = util::normSqr<J,DESCRIPTOR::d>(j);
317 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
318 cell[iPop] *= V{1} - omega;
325 template <
typename CELL,
typename RHO,
typename U,
typename RATIORHO,
typename OMEGA,
typename V=
typename CELL::value_t>
327 const RATIORHO& ratioRho,
const OMEGA& omega)
any_platform
329 const V uSqr = util::normSqr<U,DESCRIPTOR::d>(u);
330 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
333 ratioRho*(feq+descriptors::t<V,DESCRIPTOR>(iPop))-descriptors::t<V,DESCRIPTOR>(iPop) +
334 (V{1}-omega)*(cell[iPop]-feq);
340 template <
typename CELL,
typename RHO,
typename U,
typename OMEGA,
typename V=
typename CELL::value_t>
343 const V uSqr = util::normSqr<U,DESCRIPTOR::d>(u);
346 for (
int iD = 0; iD < DESCRIPTOR::d; ++iD ) {
350 V fEq[DESCRIPTOR::q];
351 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop ) {
353 for (
int iD = 0; iD < DESCRIPTOR::d; ++iD ) {
354 j1[iD] += descriptors::c<DESCRIPTOR>(iPop,iD) * ( cell[iPop] - fEq[iPop] );
359 for (
int iPop = 0; iPop < DESCRIPTOR::q; ++iPop ) {
361 for (
int iD = 0; iD < DESCRIPTOR::d; ++iD ) {
362 fNeq += descriptors::c<DESCRIPTOR>(iPop,iD) * j1[iD];
364 fNeq *= descriptors::t<V,DESCRIPTOR>(iPop) * descriptors::invCs2<V,DESCRIPTOR>();
365 cell[iPop] = fEq[iPop] + ( V{1} - omega ) * fNeq;
371 template <
typename CELL,
typename RHO,
typename U,
typename PI,
typename OMEGA,
typename V=
typename CELL::value_t>
374 const V uSqr = util::normSqr<U,DESCRIPTOR::d>(u);
377 for (
int iPop=1; iPop <= DESCRIPTOR::q/2; ++iPop) {
383 cell[iPop+DESCRIPTOR::q/2] += fNeq;
388 template <
typename CELL,
typename NEWRHO,
typename NEWU,
typename V=
typename CELL::value_t>
391 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
396 template <
typename CELL,
typename OLDRHO,
typename OLDU,
typename NEWRHO,
typename NEWU,
typename V=
typename CELL::value_t>
398 const OLDRHO& oldRho,
const OLDU& oldU,
401 const V oldUSqr = util::normSqr<OLDU,DESCRIPTOR::d>(oldU);
402 const V newUSqr = util::normSqr<NEWU,DESCRIPTOR::d>(newU);
403 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
409 template <
typename CELL,
typename RHO,
typename U,
typename PI,
typename V=
typename CELL::value_t>
415 const V uSqr = util::normSqr<U,DESCRIPTOR::d>(u);
416 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
423 template <
typename CELL,
typename FORCE,
typename V=
typename CELL::value_t>
431 for (
int Alpha=0; Alpha<DESCRIPTOR::d; ++Alpha) {
432 for (
int Beta=Alpha; Beta<DESCRIPTOR::d; ++Beta) {
433 ForceTensor[iPi] = rho/2.*(force[Alpha]*u[Beta] + u[Alpha]*force[Beta]);
438 for (
int iPi=0; iPi < util::TensorVal<DESCRIPTOR >::n; ++iPi) {
439 pi[iPi] += ForceTensor[iPi];
441 V PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2];
443 PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] +pi[5]*pi[5];
449 template <
typename CELL,
typename V=
typename CELL::value_t>
454 V PiNeqNormSqr = pi[0]*pi[0] + 2.*pi[1]*pi[1] + pi[2]*pi[2];
456 PiNeqNormSqr += pi[2]*pi[2] + pi[3]*pi[3] + 2.*pi[4]*pi[4] +pi[5]*pi[5];
462 template <
typename CELL,
typename RHO,
typename U,
typename OMEGA,
typename FORCE,
typename V=
typename CELL::value_t>
465 for (
int iPop=0; iPop < DESCRIPTOR::q; ++iPop) {
467 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
468 c_u += descriptors::c<DESCRIPTOR>(iPop,iD)*u[iD];
470 c_u *= descriptors::invCs2<V,DESCRIPTOR>() * descriptors::invCs2<V,DESCRIPTOR>();
472 for (
int iD=0; iD < DESCRIPTOR::d; ++iD) {
474 ( (descriptors::c<DESCRIPTOR>(iPop,iD) - u[iD]) * descriptors::invCs2<V,DESCRIPTOR>()
475 + c_u * descriptors::c<DESCRIPTOR>(iPop,iD)
479 forceTerm *= descriptors::t<V,DESCRIPTOR>(iPop);
480 forceTerm *= V{1} - omega * V{0.5};
482 cell[iPop] += forceTerm;
Top level namespace for all of OpenLB.
static V fromPiToFneq(int iPop, const PI &pi) any_platform
Compute off-equilibrium part of the f's from the stress tensor Pi.
static V secondOrder(int iPop, const RHO &rho, const U &u) any_platform
Computation of equilibrium distribution, second order in u.
static V incompressible(int iPop, const J &j, const JSQR &jSqr, const PRESSURE &pressure) any_platform
static V fromJneqToFneq(int iPop, const V jNeq[DESCRIPTOR::d]) any_platform
static V P1(int iPop, const RHO &rho, const U &u) any_platform
static V secondOrder(int iPop, const RHO &rho, const U &u, const USQR &uSqr) any_platform
Computation of equilibrium distribution, second order in u.
static V fromJgradToFneq(int iPop, const V Jgrad[DESCRIPTOR::d *DESCRIPTOR::d], V omega) any_platform
compute off-equilibrium part of the populations from gradient of the flux for asymmetric regularizati...
static V firstOrder(int iPop, const RHO &rho, const U &u) any_platform
Computation of equilibrium distribution, first order in u.
static V incompressible(int iPop, const J &j, const PRESSURE &pressure) any_platform
Collection of common computations for LBM.
static V adeBgkCollision(CELL &cell, const RHO &rho, const VELOCITY &u, const OMEGA &omega) any_platform
Advection diffusion BGK collision step.
static void defineEqFirstOrder(CELL &cell, const NEWRHO &newRho, const NEWU &newU) any_platform
static V constRhoBgkCollision(CELL &cell, const RHO &rho, const U &u, const RATIORHO &ratioRho, const OMEGA &omega) any_platform
BGK collision step with density correction.
static void computeFeq(CELL &cell, FEQ &fEq) any_platform
static V computePiNeqNormSqr(CELL &cell) any_platform
Computes squared norm of non-equilibrium part of 2nd momentum for standard (non-forced) dynamics.
static void computeJ(CELL &cell, J &j) any_platform
Computation of momentum.
static void computeFneq(CELL &cell, FNEQ &fNeq, const RHO &rho, const U &u) any_platform
Computation of non-equilibrium distribution.
static void computeAllMomenta(CELL &cell, RHO &rho, U &u, PI &pi) any_platform
Computation of all hydrodynamic variables.
static V rlbCollision(CELL &cell, const RHO &rho, const U &u, const PI &pi, const OMEGA &omega) any_platform
Renormalized DESCRIPTOR Boltzmann collision operator, fIn --> fOut.
static void computeFneq(CELL &cell, FNEQ &fNeq) any_platform
static V computePiNeqNormSqr(CELL &cell, const FORCE &force) any_platform
Computes squared norm of non-equilibrium part of 2nd momentum for forced dynamics.
static V rlbCollision(CELL &cell, const RHO &rho, const U &u, const OMEGA &omega) any_platform
RLB advection diffusion collision step.
static void computeRhoJ(CELL &cell, RHO &rho, J &j) any_platform
Computation of hydrodynamic variables.
static V computeRho(CELL &cell) any_platform
Computation of density.
static void addExternalForce(CELL &cell, const RHO &rho, const U &u, const OMEGA &omega, const FORCE &force) any_platform
Add a force term after BGK collision.
static V bgkCollision(CELL &cell, const RHO &rho, const VELOCITY &u, const OMEGA &omega) any_platform
BGK collision step.
static void computeStress(CELL &cell, const RHO &rho, const U &u, PI &pi) any_platform
Computation of stress tensor.
static void defineNEq(CELL &cell, const OLDRHO &oldRho, const OLDU &oldU, const NEWRHO &newRho, const NEWU &newU) any_platform
static V incBgkCollision(CELL &cell, const PRESSURE &pressure, const J &j, const OMEGA &omega) any_platform
Incompressible BGK collision step.
static void defineNEqFromPi(CELL &cell, const RHO &rho, const U &u, const PI &pi) any_platform
static void computeRhoU(CELL &cell, RHO &rho, U &u) any_platform
Computation of hydrodynamic variables.
Compute number of elements of a symmetric d-dimensional tensor.
Set of functions commonly used in LB computations – header file.