⭐ 欢迎来到虫虫下载站! | 📦 资源下载 📁 资源专辑 ℹ️ 关于我们
⭐ 虫虫下载站

📄 lbhelpers3d.h

📁 open lattice boltzmann project www.openlb.org
💻 H
📖 第 1 页 / 共 2 页
字号:
        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 + -