📄 lbhelpers3d.h
字号:
grid[iX][iY][iZ][iPop+9] = grid[nX][nY][nZ][iPop]; grid[nX][nY][nZ][iPop] = fTmp; } static void swapAndStream3D(Cell<T,descriptors::ForcedD3Q19Descriptor> ***grid, int iX, int iY, int iZ) { T fTmp; swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY, iZ, 1, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX, iY-1, iZ, 2, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX, iY , iZ-1, 3, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY-1, iZ, 4, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY+1, iZ, 5, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY , iZ-1, 6, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY , iZ+1, 7, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX , iY-1, iZ-1, 8, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX , iY-1, iZ+1, 9, fTmp); }};// Efficient specialization for D3Q15 latticetemplate<typename T>struct lbDynamicsHelpers<T, descriptors::D3Q15DescriptorBase<T> > { static T equilibrium( int iPop, T rho, const T u[3], const T uSqr ) { typedef descriptors::D3Q15DescriptorBase<T> L; T c_u = L::c[iPop][0]*u[0] + L::c[iPop][1]*u[1] + L::c[iPop][2]*u[2]; return rho * L::t[iPop] * ( 1. + 3.*c_u + 4.5*c_u*c_u - 1.5*uSqr ) - L::t[iPop]; } static T incEquilibrium(int iPop, const T j[3], const T jSqr, const T pressure) { typedef descriptors::D3Q15DescriptorBase<T> L; T c_j = L::c[iPop][0]*j[0] + L::c[iPop][1]*j[1] + L::c[iPop][2]*j[2]; return L::t[iPop] * ( 3.*pressure + 3.*c_j + 4.5*c_j*c_j - 1.5*jSqr ) - L::t[iPop]; } static void computeFneq(T const* cell, T fNeq[15], T rho, const T u[3] ) { const T uSqr = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; for (int iPop=0; iPop < 15; ++iPop) { fNeq[iPop] = cell[iPop] - equilibrium(iPop, rho, u, uSqr); } } static T bgkCollision(T* cell, T rho, const T u[3], T omega) { const T uSqr = u[0]*u[0] + u[1]*u[1] + u[2]*u[2]; for (int iPop=0; iPop < 15; ++iPop) { cell[iPop] *= (T)1-omega; cell[iPop] += omega * lbDynamicsHelpers<T,descriptors::D3Q15DescriptorBase<T> >::equilibrium(iPop, rho, u, uSqr); } return uSqr; } static T incBgkCollision(T* cell, T pressure, const T j[3], T omega) { const T jSqr = util::normSqr<T,descriptors::D3Q15DescriptorBase<T>::d>(j); for (int iPop=0; iPop < descriptors::D3Q15DescriptorBase<T>::q; ++iPop) { cell[iPop] *= (T)1-omega; cell[iPop] += omega * lbDynamicsHelpers<T,descriptors::D3Q15DescriptorBase<T> >::incEquilibrium ( iPop, j, jSqr, pressure ); } return jSqr; } static T constRhoBgkCollision(T* cell, T rho, const T u[3], T ratioRho, T omega) { const T uSqr = util::normSqr<T,descriptors::D3Q15DescriptorBase<T>::d>(u); for (int iPop=0; iPop < descriptors::D3Q15DescriptorBase<T>::q; ++iPop) { T feq = lbHelpers<T,descriptors::D3Q15DescriptorBase>:: equilibrium(iPop, rho, u, uSqr ); cell[iPop] = ratioRho*(feq+descriptors::D3Q15DescriptorBase<T>::t[iPop]) -descriptors::D3Q15DescriptorBase<T>::t[iPop] + ((T)1-omega)*(cell[iPop]-feq); } return uSqr; } static void partial_rho(T const* cell, T& surfX_M1, T& surfX_0, T& surfX_P1, T& surfY_M1, T& surfY_P1, T& surfZ_M1, T& surfZ_P1 ) { surfX_M1 = cell[1] + cell[4] + cell[5] + cell[6] + cell[7]; surfX_0 = cell[0] + cell[2] + cell[3] + cell[9] + cell[10]; surfX_P1 = cell[8] + cell[11] + cell[12] + cell[13] + cell[14]; surfY_M1 = cell[2] + cell[4] + cell[5] + cell[13] + cell[14]; surfY_P1 = cell[6] + cell[7] + cell[9] + cell[11] + cell[12]; surfZ_M1 = cell[3] + cell[4] + cell[6] + cell[12] + cell[14]; surfZ_P1 = cell[5] + cell[7] + cell[10] + cell[11] + cell[13]; } static T computeRho(T const* cell) { T rho = cell[0] + cell[1] + cell[2] + cell[3] + cell[4] + cell[5] + cell[6] + cell[7] + cell[8] + cell[9] + cell[10] + cell[11] + cell[12] + cell[13] + cell[14] + (T)1; return rho; } static void computeRhoU(T const* cell, T& rho, T u[3]) { T surfX_M1, surfX_0, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1; partial_rho(cell, surfX_M1, surfX_0, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1); rho = surfX_M1 + surfX_0 + surfX_P1 + (T)1; u[0] = ( surfX_P1 - surfX_M1 ) / rho; u[1] = ( surfY_P1 - surfY_M1 ) / rho; u[2] = ( surfZ_P1 - surfZ_M1 ) / rho; } static void computeJ(T const* cell, T j[3]) { T surfX_M1, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1; surfX_M1 = cell[1] + cell[4] + cell[5] + cell[6] + cell[7]; surfX_P1 = cell[8] + cell[11] + cell[12] + cell[13] + cell[14]; surfY_M1 = cell[2] + cell[4] + cell[5] + cell[13] + cell[14]; surfY_P1 = cell[6] + cell[7] + cell[9] + cell[11] + cell[12]; surfZ_M1 = cell[3] + cell[4] + cell[6] + cell[12] + cell[14]; surfZ_P1 = cell[5] + cell[7] + cell[10] + cell[11] + cell[13]; j[0] = ( surfX_P1 - surfX_M1 ); j[1] = ( surfY_P1 - surfY_M1 ); j[2] = ( surfZ_P1 - surfZ_M1 ); } static void computeStress(T const* cell, T rho, const T u[3], T pi[6]) { typedef descriptors::D3Q15DescriptorBase<T> L; using namespace util::tensorIndices3D; T surfX_M1, surfX_0, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1; partial_rho(cell, surfX_M1, surfX_0, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1); pi[xx] = surfX_P1+surfX_M1 - 1./L::invCs2*(rho-(T)1) - rho*u[0]*u[0]; pi[yy] = surfY_P1+surfY_M1 - 1./L::invCs2*(rho-(T)1) - rho*u[1]*u[1]; pi[zz] = surfZ_P1+surfZ_M1 - 1./L::invCs2*(rho-(T)1) - rho*u[2]*u[2]; pi[xy] = cell[4] + cell[5] - cell[6] - cell[7] + cell[11] + cell[12] - cell[13] - cell[14] - rho*u[0]*u[1]; pi[xz] = cell[4] - cell[5] + cell[6] - cell[7] + cell[11] - cell[12] + cell[13] - cell[14] - rho*u[0]*u[2]; pi[yz] = cell[4] - cell[5] - cell[6] + cell[7] + cell[11] - cell[12] - cell[13] + cell[14] - rho*u[1]*u[2]; } static void computeAllMomenta(T const* cell, T& rho, T u[3], T pi[6]) { typedef descriptors::D3Q15DescriptorBase<T> L; using namespace util::tensorIndices3D; T surfX_M1, surfX_0, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1; partial_rho(cell, surfX_M1, surfX_0, surfX_P1, surfY_M1, surfY_P1, surfZ_M1, surfZ_P1); rho = surfX_M1 + surfX_0 + surfX_P1 + (T)1; T rhoU0 = ( surfX_P1 - surfX_M1 ) / rho; T rhoU1 = ( surfY_P1 - surfY_M1 ) / rho; T rhoU2 = ( surfZ_P1 - surfZ_M1 ) / rho; u[0] = rhoU0 / rho; u[1] = rhoU1 / rho; u[2] = rhoU2 / rho; pi[xx] = surfX_P1+surfX_M1 - 1./L::invCs2*(rho-(T)1) - rhoU0*u[0]; pi[yy] = surfY_P1+surfY_M1 - 1./L::invCs2*(rho-(T)1) - rhoU1*u[1]; pi[zz] = surfZ_P1+surfZ_M1 - 1./L::invCs2*(rho-(T)1) - rhoU2*u[2]; pi[xy] = cell[4] + cell[5] - cell[6] - cell[7] + cell[11] + cell[12] - cell[13] - cell[14] - rhoU0*u[1]; pi[xz] = cell[4] - cell[5] + cell[6] - cell[7] + cell[11] - cell[12] + cell[13] - cell[14] - rhoU0*u[2]; pi[yz] = cell[4] - cell[5] - cell[6] + cell[7] + cell[11] - cell[12] - cell[13] + cell[14] - rhoU1*u[2]; } static void modifyVelocity(T* cell, const T newU[3]) { T rho, oldU[3]; computeRhoU(cell, rho, oldU); const T oldUSqr = util::normSqr<T,3>(oldU); const T newUSqr = util::normSqr<T,3>(newU); for (int iPop=0; iPop<15; ++iPop) { cell[iPop] = cell[iPop] - equilibrium(iPop, rho, oldU, oldUSqr) + equilibrium(iPop, rho, newU, newUSqr); } }}; //struct lbDynamicsHelpers<D3Q15DescriptorBase>// Efficient specialization for D3Q15 lattice and for forced D3Q15 lattice// (operations applying to the whole lattice)template<typename T>struct lbLatticeHelpers<T, descriptors::D3Q15Descriptor> { static void swapAndStreamCell ( Cell<T,descriptors::D3Q15Descriptor> ***grid, int iX, int iY, int iZ, int nX, int nY, int nZ, int iPop, T& fTmp ) { fTmp = grid[iX][iY][iZ][iPop]; grid[iX][iY][iZ][iPop] = grid[iX][iY][iZ][iPop+7]; grid[iX][iY][iZ][iPop+7] = grid[nX][nY][nZ][iPop]; grid[nX][nY][nZ][iPop] = fTmp; } static void swapAndStream3D(Cell<T,descriptors::D3Q15Descriptor> ***grid, int iX, int iY, int iZ) { T fTmp; swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY, iZ, 1, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX, iY-1, iZ, 2, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX, iY , iZ-1, 3, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY-1, iZ-1, 4, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY-1, iZ+1, 5, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY+1, iZ-1, 6, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY+1, iZ+1, 7, fTmp); }};template<typename T>struct lbLatticeHelpers<T, descriptors::ForcedD3Q15Descriptor> { static void swapAndStreamCell ( Cell<T,descriptors::D3Q15Descriptor> ***grid, int iX, int iY, int iZ, int nX, int nY, int nZ, int iPop, T& fTmp ) { fTmp = grid[iX][iY][iZ][iPop]; grid[iX][iY][iZ][iPop] = grid[iX][iY][iZ][iPop+7]; grid[iX][iY][iZ][iPop+7] = grid[nX][nY][nZ][iPop]; grid[nX][nY][nZ][iPop] = fTmp; } static void swapAndStream3D(Cell<T,descriptors::D3Q15Descriptor> ***grid, int iX, int iY, int iZ) { T fTmp; swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY, iZ, 1, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX, iY-1, iZ, 2, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX, iY , iZ-1, 3, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY-1, iZ-1, 4, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY-1, iZ+1, 5, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY+1, iZ-1, 6, fTmp); swapAndStreamCell(grid, iX, iY, iZ, iX-1, iY+1, iZ+1, 7, fTmp); }};} // namespace olb#endif
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -