49 int x0_,
int x1_,
int y0_,
int y1_,
int z0_,
int z1_)
57 for (
int iX=x0; iX<=x1; ++iX) {
58 for (
int iY=y0; iY<=y1; ++iY) {
59 for (
int iZ=z0; iZ<=z1; ++iZ) {
63 T dx_U[DESCRIPTOR::d], dy_U[DESCRIPTOR::d], dz_U[DESCRIPTOR::d];
64 interpolateGradients<0>(blockLattice, dx_U, iX, iY, iZ);
65 interpolateGradients<1>(blockLattice, dy_U, iX, iY, iZ);
66 interpolateGradients<2>(blockLattice, dz_U, iX, iY, iZ);
68 T rhoGradU[L::d][L::d];
69 rhoGradU[x][x] = rho *dx_U[x];
70 rhoGradU[x][y] = rho *dx_U[y];
71 rhoGradU[x][z] = rho *dx_U[z];
72 rhoGradU[y][x] = rho *dy_U[x];
73 rhoGradU[y][y] = rho *dy_U[y];
74 rhoGradU[y][z] = rho *dy_U[z];
75 rhoGradU[z][x] = rho *dz_U[x];
76 rhoGradU[z][y] = rho *dz_U[y];
77 rhoGradU[z][z] = rho *dz_U[z];
79 T omega = blockLattice.
getDynamics(iX, iY, iZ) -> getOmega();
80 T sToPi = - (T)1 / descriptors::invCs2<T,DESCRIPTOR>() / omega;
83 pi[xx] = (T)2 * rhoGradU[x][x] * sToPi;
84 pi[yy] = (T)2 * rhoGradU[y][y] * sToPi;
85 pi[zz] = (T)2 * rhoGradU[z][z] * sToPi;
86 pi[xy] = (rhoGradU[x][y] + rhoGradU[y][x]) * sToPi;
87 pi[xz] = (rhoGradU[x][z] + rhoGradU[z][x]) * sToPi;
88 pi[yz] = (rhoGradU[y][z] + rhoGradU[z][y]) * sToPi;
93 T uSqr = util::normSqr<T,DESCRIPTOR::d>(u);
97 T dx_rho, dy_rho, dz_rho;
98 interpolateGradients<0>(blockLattice, dx_rho, iX, iY, iZ);
99 interpolateGradients<1>(blockLattice, dy_rho, iX, iY, iZ);
100 interpolateGradients<2>(blockLattice, dz_rho, iX, iY, iZ);
101 for (
int iPop = 0; iPop < L::q; ++iPop) {
103 for (
int iAlpha=0; iAlpha < L::d; ++iAlpha) {
104 cGradRhoUU += descriptors::c<L>(iPop,iAlpha) * (
105 dx_rho*u[iAlpha]*u[x] +
106 dx_U[iAlpha]*rho*u[x] +
107 dx_U[x]*rho*u[iAlpha] +
108 dy_rho*u[iAlpha]*u[y] +
109 dy_U[iAlpha]*rho*u[y] +
110 dy_U[y]*rho*u[iAlpha] +
111 dz_rho*u[iAlpha]*u[z] +
112 dz_U[iAlpha]*rho*u[z] +
113 dz_U[z]*rho*u[iAlpha]);
118 T cDivRhoUU[L::d][L::d];
119 for (
int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
120 for (
int iBeta = 0; iBeta < L::d; ++iBeta) {
121 cDivRhoUU[iAlpha][iBeta] = descriptors::c<L>(iPop,x)*
122 (dx_rho*u[iAlpha]*u[iBeta] +
123 dx_U[iAlpha]*rho*u[iBeta] +
124 dx_U[iBeta]*rho*u[iAlpha])
125 +descriptors::c<L>(iPop,y)*
126 (dy_rho*u[iAlpha]*u[iBeta] +
127 dy_U[iAlpha]*rho*u[iBeta] +
128 dy_U[iBeta]*rho*u[iAlpha])
129 +descriptors::c<L>(iPop,z)*
130 (dz_rho*u[iAlpha]*u[iBeta] +
131 dz_U[iAlpha]*rho*u[iBeta] +
132 dz_U[iBeta]*rho*u[iAlpha]);
141 for (
int iAlpha = 0; iAlpha < L::d; ++iAlpha) {
142 for (
int iBeta = 0; iBeta < L::d; ++iBeta) {
143 int ci_ci = descriptors::c<L>(iPop,iAlpha)*descriptors::c<L>(iPop,iBeta);
144 qCdivRhoUU += ci_ci * cDivRhoUU[iAlpha][iBeta];
145 qRhoGradU += ci_ci * rhoGradU[iAlpha][iBeta];
146 if (iAlpha == iBeta) {
147 qCdivRhoUU -= cDivRhoUU[iAlpha][iBeta]/descriptors::invCs2<T,L>();
148 qRhoGradU -= rhoGradU[iAlpha][iBeta]/descriptors::invCs2<T,L>();
155 cell[iPop] = lbH::equilibrium(iPop,rho,u,uSqr)
156 - descriptors::t<T,L>(iPop) * descriptors::invCs2<T,L>() / omega
157 * (qRhoGradU - cGradRhoUU + 0.5*descriptors::invCs2<T,L>()*qCdivRhoUU);
This class computes the finite difference approximation to LB boundary conditions on a plane wall in ...
void processSubDomain(BlockLattice< T, DESCRIPTOR > &blockLattice, int x0_, int x1_, int y0_, int y1_, int z0_, int z1_) override
Execute post-processing step on a sublattice.
void process(BlockLattice< T, DESCRIPTOR > &blockLattice) override
Execute post-processing step.
std::string & getName()
read and write access to name