blocklattice2d.hh

来自「open lattice boltzmann project www.open」· HH 代码 · 共 915 行 · 第 1/3 页

 * collide(int,int,int,int) and stream(int,int,int,int), because memory
 * is traversed only once instead of twice.
 * The post-processing steps are not automatically invoked by this method,
 * as they are in the method stream(). If you want them to be executed, you
 * must explicitly call the method postProcess().
 * \sa collideAndStream()
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::collideAndStream(int x0, int x1, int y0, int y1)
{
    OLB_PRECONDITION(x0>=0 && x1<nx);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=0 && y1<ny);
    OLB_PRECONDITION(y1>=y0);

    collide(x0,x0, y0,y1);
    collide(x1,x1, y0,y1);
    collide(x0+1,x1-1, y0,y0);
    collide(x0+1,x1-1, y1,y1);

    bulkCollideAndStream(x0+1,x1-1,y0+1,y1-1);

    boundaryStream(x0,x1,y0,y1, x0,x0,y0,y1);
    boundaryStream(x0,x1,y0,y1, x1,x1,y0,y1);
    boundaryStream(x0,x1,y0,y1, x0+1,x1-1,y0,y0);
    boundaryStream(x0,x1,y0,y1, x0+1,x1-1,y1,y1);
}

/** At the end of this method, the post-processing steps are automatically
 * invoked.
 * \sa collideAndStream(int,int,int,int) */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::collideAndStream(bool periodic) {
    collideAndStream(0, nx-1, 0, ny-1);
    
    if (periodic) {
        makePeriodic();
    }

    postProcess();
}

template<typename T, template<typename U> class Lattice>
T BlockLattice2D<T,Lattice>::computeAverageDensity (
        int x0, int x1, int y0, int y1) const
{
    T sumRho = T();
    for (int iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            T rho, u[Lattice<T>::d];
            get(iX,iY).computeRhoU(rho, u);
            sumRho += rho;
        }
    }
    return sumRho / (T)(x1-x0+1) / (T)(y1-y0+1);
}

template<typename T, template<typename U> class Lattice>
T BlockLattice2D<T,Lattice>::computeAverageDensity() const {
    return computeAverageDensity(0, nx-1, 0, ny-1);
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::stripeOffDensityOffset (
        int x0, int x1, int y0, int y1, T offset )
{
    for (int iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            for (int iPop=0; iPop<Lattice<T>::q; ++iPop) {
                get(iX,iY)[iPop] -= Lattice<T>::t[iPop] * offset;
            }
        }
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::stripeOffDensityOffset(T offset) {
    stripeOffDensityOffset(0, nx-1, 0, ny-1, offset);
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::forAll (
        int x0, int x1, int y0, int y1, WriteCellFunctional<T,Lattice> const& application )
{
    for (int iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            application.apply( get(iX,iY) );
        }
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::forAll(WriteCellFunctional<T,Lattice> const& application) {
    forAll(0, nx-1, 0, ny-1, application);
}


template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::addPostProcessor (
    PostProcessorGenerator2D<T,Lattice> const& ppGen )
{
    postProcessors.push_back(ppGen.generate());
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::resetPostProcessors() {
    clearPostProcessors();
    StatPPGenerator2D<T,Lattice> statPPGenerator;
    addPostProcessor(statPPGenerator);
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::clearPostProcessors() {
    typename PostProcVector::iterator ppIt = postProcessors.begin();
    for (; ppIt != postProcessors.end(); ++ppIt) {
        delete *ppIt;
    }
    postProcessors.clear();
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::postProcess() {
    for (unsigned iPr=0; iPr<postProcessors.size(); ++iPr) {
        postProcessors[iPr] -> process(*this);
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::postProcess(int x0_, int x1_, int y0_, int y1_)
{
    for (unsigned iPr=0; iPr<postProcessors.size(); ++iPr) {
        postProcessors[iPr] -> processSubDomain(*this, x0_, x1_, y0_, y1_);
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::addLatticeCoupling (
    LatticeCouplingGenerator2D<T,Lattice> const& lcGen,
    std::vector<SpatiallyExtendedObject2D*> partners )
{
    latticeCouplings.push_back(lcGen.generate(partners));
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::executeCoupling() {
    for (unsigned iPr=0; iPr<latticeCouplings.size(); ++iPr) {
        latticeCouplings[iPr] -> process(*this);
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::executeCoupling(int x0_, int x1_, int y0_, int y1_) {
    for (unsigned iPr=0; iPr<latticeCouplings.size(); ++iPr) {
        latticeCouplings[iPr] -> processSubDomain(*this, x0_, x1_, y0_, y1_);
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::clearLatticeCouplings() {
    typename PostProcVector::iterator ppIt = latticeCouplings.begin();
    for (; ppIt != latticeCouplings.end(); ++ppIt) {
        delete *ppIt;
    }
    latticeCouplings.clear();
}


template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::subscribeReductions(Reductor<T>& reductor) {
    for (unsigned iPr=0; iPr<postProcessors.size(); ++iPr) {
        postProcessors[iPr] -> subscribeReductions(*this, &reductor);
    }
}

template<typename T, template<typename U> class Lattice>
LatticeStatistics<T>& BlockLattice2D<T,Lattice>::getStatistics() {
    #ifdef PARALLEL_MODE_OMP
        return *statistics[omp.get_rank()];
    #else
        return *statistics;
    #endif
}

template<typename T, template<typename U> class Lattice>
LatticeStatistics<T> const& BlockLattice2D<T,Lattice>::getStatistics() const {
    #ifdef PARALLEL_MODE_OMP
        return *statistics[omp.get_rank()];
    #else
        return *statistics;
    #endif
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::allocateMemory() {
    // The conversions to size_t ensure 64-bit compatibility. Note that
    //   nx and ny are of type int, which might by 32-bit types, even on
    //   64-bit platforms. Therefore, nx*ny may lead to a type overflow.
    rawData = new Cell<T,Lattice> [(size_t)nx*(size_t)ny];
    grid    = new Cell<T,Lattice>* [(size_t)nx];
    for (int iX=0; iX<nx; ++iX) {
        grid[iX] = rawData + (size_t)iX*(size_t)ny;
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::releaseMemory() {
    delete [] rawData;
    delete [] grid;
}

/** This method is slower than bulkStream(int,int,int,int), because it must
 * be verified which distribution functions are to be kept from leaving
 * the domain.
 * \sa stream(int,int,int,int)
 * \sa stream()
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::boundaryStream (
    int lim_x0, int lim_x1, int lim_y0, int lim_y1,
    int x0, int x1, int y0, int y1 )
{
    OLB_PRECONDITION(lim_x0>=0 && lim_x1<nx);
    OLB_PRECONDITION(lim_x1>=lim_x0);
    OLB_PRECONDITION(lim_y0>=0 && lim_y1<ny);
    OLB_PRECONDITION(lim_y1>=lim_y0);

    OLB_PRECONDITION(x0>=lim_x0 && x1<=lim_x1);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=lim_y0 && y1<=lim_y1);
    OLB_PRECONDITION(y1>=y0);

    int iX;

    #ifdef PARALLEL_MODE_OMP
    #pragma omp parallel for
    #endif
    for (iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            for (int iPop=1; iPop<=Lattice<T>::q/2; ++iPop) {
                int nextX = iX + Lattice<T>::c[iPop][0];
                int nextY = iY + Lattice<T>::c[iPop][1];
                if (nextX>=lim_x0 && nextX<=lim_x1 && nextY>=lim_y0 && nextY<=lim_y1) {
                    std::swap(grid[iX][iY][iPop+Lattice<T>::q/2],
                              grid[nextX][nextY][iPop]);
                }
            }
        }
    }
}

/** This method is faster than boundaryStream(int,int,int,int), but it
 * is erroneous when applied to boundary cells.
 * \sa stream(int,int,int,int)
 * \sa stream()
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::bulkStream (
    int x0, int x1, int y0, int y1 )
{
    OLB_PRECONDITION(x0>=0 && x1<nx);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=0 && y1<ny);
    OLB_PRECONDITION(y1>=y0);

    int iX;
    #ifdef PARALLEL_MODE_OMP
    #pragma omp parallel for
    #endif
    for (iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            for (int iPop=1; iPop<=Lattice<T>::q/2; ++iPop) {
                int nextX = iX + Lattice<T>::c[iPop][0];
                int nextY = iY + Lattice<T>::c[iPop][1];
                std::swap(grid[iX][iY][iPop+Lattice<T>::q/2],
                          grid[nextX][nextY][iPop]);
            }
        }
    }
}

#ifndef PARALLEL_MODE_OMP  // OpenMP parallel version is and the
                           // end of this file
/** This method is fast, but it is erroneous when applied to boundary
 * cells.
 * \sa collideAndStream(int,int,int,int)
 * \sa collideAndStream()
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::bulkCollideAndStream (
    int x0, int x1, int y0, int y1 )
{
    OLB_PRECONDITION(x0>=0 && x1<nx);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=0 && y1<ny);
    OLB_PRECONDITION(y1>=y0);

    for (int iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            grid[iX][iY].collide(getStatistics());
            lbHelpers<T,Lattice>::swapAndStream2D(grid, iX, iY);
        }
    }
}
#endif // not defined PARALLEL_MODE_OMP