lbhelpers.h

open lattice boltzmann project www.openlb.org

        T rho = T();
        for (int iPop=0; iPop < Descriptor::q; ++iPop) {
            rho += cell[iPop];
        }
        rho += (T)1;
        return rho;
    }

    /// Computation of momentum
    static void computeJ(T const* cell, T j[Descriptor::d]) {
        for (int iD=0; iD < Descriptor::d; ++iD) {
            j[iD] = T();
        }
        for (int iPop=0; iPop < Descriptor::q; ++iPop) {
            for (int iD=0; iD < Descriptor::d; ++iD) {
                j[iD] += cell[iPop]*Descriptor::c[iPop][iD];
            }
        }
    }

    /// Computation of hydrodynamic variables
    static void computeRhoU(T const* cell, T& rho, T u[Descriptor::d]) {
        rho = T();
        for (int iD=0; iD < Descriptor::d; ++iD) {
            u[iD] = T();
        }
        for (int iPop=0; iPop < Descriptor::q; ++iPop) {
            rho += cell[iPop];
            for (int iD=0; iD < Descriptor::d; ++iD) {
                u[iD] += cell[iPop]*Descriptor::c[iPop][iD];
            }
        }
        rho += (T)1;
        for (int iD=0; iD < Descriptor::d; ++iD) {
            u[iD] /= rho;
        }
    }

    /// Computation of stress tensor
    static void computeStress(T const* cell, T rho, const T u[Descriptor::d],
                              T pi[util::TensorVal<Descriptor>::n] )
    {
        int iPi = 0;
        for (int iAlpha=0; iAlpha < Descriptor::d; ++iAlpha) {
            for (int iBeta=iAlpha; iBeta < Descriptor::d; ++iBeta) {
                pi[iPi] = T();
                for (int iPop=0; iPop < Descriptor::q; ++iPop) {
                    pi[iPi] += Descriptor::c[iPop][iAlpha]*
                               Descriptor::c[iPop][iBeta] * cell[iPop];
                }
                // stripe off equilibrium contribution
                pi[iPi] -= rho*u[iAlpha]*u[iBeta];
                if (iAlpha==iBeta) {
                    pi[iPi] -= 1./Descriptor::invCs2*(rho-(T)1);
                }
                ++iPi;
            }
        }
    }
    
    /// Computation of all hydrodynamic variables
    static void computeAllMomenta(T const* cell, T& rho, T u[Descriptor::d],
                                  T pi[util::TensorVal<Descriptor>::n] )
    {
        computeRhoU(cell, rho, u);
        computeStress(cell, rho, u, pi);
    }

    static void modifyVelocity(T* cell, const T newU[Descriptor::d]) {
        T rho, oldU[Descriptor::d];
        computeRhoU(cell, rho, oldU);
        const T oldUSqr = util::normSqr<T,Descriptor::d>(oldU);
        const T newUSqr = util::normSqr<T,Descriptor::d>(newU);
        for (int iPop=0; iPop<Descriptor::q; ++iPop) {
            cell[iPop] = cell[iPop]
                             - equilibrium(iPop, rho, oldU, oldUSqr)
                             + equilibrium(iPop, rho, newU, newUSqr);
        }
    }

};  // struct lbDynamicsHelpers

/// Helper functions for dynamics that access external field
template<typename T, template<typename U> class Lattice>
struct lbExternalHelpers {
    /// Add a force term after BGK collision
    static void addExternalForce(Cell<T,Lattice>& cell, const T u[Lattice<T>::d], T omega, T amplitude) {
        static const int forceBeginsAt = Lattice<T>::ExternalField::forceBeginsAt;
        T* force = cell.getExternal(forceBeginsAt);
        for (int iPop=0; iPop < Lattice<T>::q; ++iPop) {
            T c_u = T();
            for (int iD=0; iD < Lattice<T>::d; ++iD) {
               c_u += Lattice<T>::c[iPop][iD]*u[iD];
            }
            c_u *= Lattice<T>::invCs2 * Lattice<T>::invCs2;
            T forceTerm = T();
            for (int iD=0; iD < Lattice<T>::d; ++iD) {
                forceTerm +=
                    (   (Lattice<T>::c[iPop][iD]-u[iD]) * Lattice<T>::invCs2
                      + c_u * Lattice<T>::c[iPop][iD]
                    )
                    * force[iD];
            }
            forceTerm *= Lattice<T>::t[iPop];
            forceTerm *= 1-omega/(T)2;
            forceTerm *= amplitude;
            cell[iPop] += forceTerm;
        }
    }
};  // struct externalFieldHelpers

/// Helper functions with full-lattice access
template<typename T, template<typename U> class Lattice>
struct lbLatticeHelpers {
    /// Swap ("bounce-back") values of a cell (2D), and apply streaming step
    static void swapAndStream2D(Cell<T,Lattice> **grid, int iX, int iY)
    {
        const int half = Lattice<T>::q/2;
        for (int iPop=1; iPop<=half; ++iPop) {
            int nextX = iX + Lattice<T>::c[iPop][0];
            int nextY = iY + Lattice<T>::c[iPop][1];
            T fTmp                   = grid[iX][iY][iPop];
            grid[iX][iY][iPop]       = grid[iX][iY][iPop+half];
            grid[iX][iY][iPop+half]  = grid[nextX][nextY][iPop];
            grid[nextX][nextY][iPop] = fTmp;
         }
    }

    /// Swap ("bounce-back") values of a cell (3D), and apply streaming step
    static void swapAndStream3D(Cell<T,Lattice> ***grid,
                                int iX, int iY, int iZ)
    {
        const int half = Lattice<T>::q/2;
        for (int iPop=1; iPop<=half; ++iPop) {
            int nextX = iX + Lattice<T>::c[iPop][0];
            int nextY = iY + Lattice<T>::c[iPop][1];
            int nextZ = iZ + Lattice<T>::c[iPop][2];
            T fTmp                          = grid[iX][iY][iZ][iPop];
            grid[iX][iY][iZ][iPop]          = grid[iX][iY][iZ][iPop+half];
            grid[iX][iY][iZ][iPop+half]     = grid[nextX][nextY][nextZ][iPop];
            grid[nextX][nextY][nextZ][iPop] = fTmp;
         }
    }
};

/// All boundary helper functions are inside this structure
template<typename T, template<typename U> class Lattice, int direction, int orientation>
struct BoundaryHelpers {
    static void computeStress (
         Cell<T,Lattice> const& cell, T rho, const T u[Lattice<T>::d],
         T pi[util::TensorVal<Lattice<T> >::n] )
    {
        typedef Lattice<T> L;
        const T uSqr = util::normSqr<T,L::d>(u);

        std::vector<int> const& onWallIndices = util::subIndex<L, direction, 0>();
        std::vector<int> const& normalIndices = util::subIndex<L, direction, orientation>();

        T fNeq[Lattice<T>::q];
        for (unsigned fIndex=0; fIndex<onWallIndices.size(); ++fIndex) {
            int iPop = onWallIndices[fIndex];
            fNeq[iPop] =
                cell[iPop] -
                lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr);
        }
        for (unsigned fIndex=0; fIndex<normalIndices.size(); ++fIndex) {
            int iPop = normalIndices[fIndex];
            if (iPop == 0) {
                fNeq[iPop] = T();  // fNeq[0] will not be used anyway
            }
            else {
                fNeq[iPop] =
                    cell[iPop] -
                    lbHelpers<T,Lattice>::equilibrium(iPop, rho, u, uSqr);
            }
        }

        int iPi = 0;
        for (int iAlpha=0; iAlpha<L::d; ++iAlpha) {
            for (int iBeta=iAlpha; iBeta<L::d; ++iBeta) {
                pi[iPi] = T();
                for (unsigned fIndex=0; fIndex<onWallIndices.size(); ++fIndex)
                {
                    const int iPop = onWallIndices[fIndex];
                    pi[iPi] +=
                        L::c[iPop][iAlpha]*L::c[iPop][iBeta]*fNeq[iPop];
                }
                for (unsigned fIndex=0; fIndex<normalIndices.size(); ++fIndex)
                {
                    const int iPop = normalIndices[fIndex];
                    pi[iPi] += (T)2 * L::c[iPop][iAlpha]*L::c[iPop][iBeta]*
                               fNeq[iPop];
                }
                ++iPi;
            }
        }
    }

};  // struct boundaryHelpers

}  // namespace olb

// The specialized code is directly included. That is because we never want
// it to be precompiled so that in both the precompiled and the
// "include-everything" version, the compiler can apply all the
// optimizations it wants.
#include "lbHelpers2D.h"
#include "lbHelpers3D.h"

#endif