blocklattice3d.hh

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

void BlockLattice3D<T,Lattice>::collideAndStream(int x0, int x1, int y0, int y1, int z0, int z1) {
    OLB_PRECONDITION(x0>=0 && x1<nx);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=0 && y1<ny);
    OLB_PRECONDITION(y1>=y0);
    OLB_PRECONDITION(z0>=0 && z1<nz);
    OLB_PRECONDITION(z1>=z0);

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

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

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

/** Post-processing steps are called at the end of this method.
 * \sa collideAndStream(int,int,int,int,int,int) */
template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::collideAndStream(bool periodic) {
    collideAndStream(0, nx-1, 0, ny-1, 0, nz-1);

    if (periodic) {
        makePeriodic();
    }

    postProcess();
}

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

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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::releaseMemory() {
    delete [] rawData;
    for (int iX=0; iX<nx; ++iX) {
      delete [] grid[iX];
    }
    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 BlockLattice3D<T,Lattice>::boundaryStream (
    int lim_x0, int lim_x1, int lim_y0, int lim_y1, int lim_z0, int lim_z1,
    int x0, int x1, int y0, int y1, int z0, int z1 )
{
    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(lim_z0>=0 && lim_z1<nz);
    OLB_PRECONDITION(lim_z1>=lim_z0);

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

    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 iZ=z0; iZ<=z1; ++iZ) {
                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];
                    int nextZ = iZ + Lattice<T>::c[iPop][2];
                    if ( nextX>=lim_x0 && nextX<=lim_x1 &&
                         nextY>=lim_y0 && nextY<=lim_y1 &&
                         nextZ>=lim_z0 && nextZ<=lim_z1 )
                    {
                        std::swap(grid[iX][iY][iZ][iPop+Lattice<T>::q/2],
                                  grid[nextX][nextY][nextZ][iPop]);
                    }
                }
            }
        }
    }
}

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

    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 iZ=z0; iZ<=z1; ++iZ) {
                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];
                    int nextZ = iZ + Lattice<T>::c[iPop][2];
                    std::swap(grid[iX][iY][iZ][iPop+Lattice<T>::q/2],
                              grid[nextX][nextY][nextZ][iPop]);
                }
            }
        }
    }
}

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

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