blocklattice3d.hh

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

template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::makePeriodic() {
    int maxX = getNx()-1;
    int maxY = getNy()-1;
    int maxZ = getNz()-1;
    periodicSurface(0,      0,      0,    maxY,   0,    maxZ);
    periodicSurface(maxX,   maxX,   0,    maxY,   0,    maxZ);
    periodicSurface(1,      maxX-1, 0,    0,      0,    maxZ);
    periodicSurface(1,      maxX-1, maxY, maxY,   0,    maxZ);
    periodicSurface(1,      maxX-1, 1,    maxY-1, 0,    0);
    periodicSurface(1,      maxX-1, 1,    maxY-1, maxZ, maxZ);
}

template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::periodicSurface (
        int x0, int x1, int y0, int y1, int z0, int z1)
{
    for (int 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<0 || nextX>=getNx() ||
                         nextY<0 || nextY>=getNy() ||
                         nextZ<0 || nextZ>=getNz() )
                    {
                        nextX = (nextX+getNx())%getNx();
                        nextY = (nextY+getNy())%getNy();
                        nextZ = (nextZ+getNz())%getNz();
                        std::swap (
                            grid[iX][iY][iZ]         [iPop+Lattice<T>::q/2],
                            grid[nextX][nextY][nextZ][iPop] );
                    }
                }
            }
        }
    }
}

template<typename T, template<typename U> class Lattice>
DataAnalysisBase3D<T,Lattice> const& BlockLattice3D<T,Lattice>::getDataAnalysis() const {
    dataAnalysis -> reset();
    return *dataAnalysis;
}

template<typename T, template<typename U> class Lattice>
DataSerializer<T> const& BlockLattice3D<T,Lattice>::getSerializer(IndexOrdering::OrderingT ordering) const {
    delete serializer;
    serializer = new BlockLatticeSerializer3D<T,Lattice>(*this, ordering);
    return *serializer;
}

template<typename T, template<typename U> class Lattice>
DataUnSerializer<T>& BlockLattice3D<T,Lattice>::getUnSerializer(IndexOrdering::OrderingT ordering) {
    delete unSerializer;
    unSerializer = new BlockLatticeUnSerializer3D<T,Lattice>(*this, ordering);
    return *unSerializer;
}

template<typename T, template<typename U> class Lattice>
DataSerializer<T> const& BlockLattice3D<T,Lattice>::getSubSerializer (
            int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
            IndexOrdering::OrderingT ordering ) const
{
    delete serializer;
    serializer = new BlockLatticeSerializer3D<T,Lattice> (
            *this, x0_, x1_, y0_, y1_, z0_, z1_, ordering );
    return *serializer;
}

template<typename T, template<typename U> class Lattice>
DataUnSerializer<T>& BlockLattice3D<T,Lattice>::getSubUnSerializer (
            int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
            IndexOrdering::OrderingT ordering )
{
    delete unSerializer;
    unSerializer = new BlockLatticeUnSerializer3D<T,Lattice> (
            *this, x0_, x1_, y0_, y1_, z0_, z1_, ordering );
    return *unSerializer;
}

template<typename T, template<typename U> class Lattice>
MultiDataDistribution3D BlockLattice3D<T,Lattice>::getDataDistribution() const {
    return MultiDataDistribution3D(getNx(), getNy(), getNz());
}

template<typename T, template<typename U> class Lattice>
SpatiallyExtendedObject3D* BlockLattice3D<T,Lattice>::getComponent(int iBlock) {
    OLB_PRECONDITION( iBlock==0 );
    return this;
}

template<typename T, template<typename U> class Lattice>
SpatiallyExtendedObject3D const* BlockLattice3D<T,Lattice>::getComponent(int iBlock) const {
    OLB_PRECONDITION( iBlock==0 );
    return this;
}

template<typename T, template<typename U> class Lattice>
multiPhysics::MultiPhysicsId BlockLattice3D<T,Lattice>::getMultiPhysicsId() const {
    return multiPhysics::getMultiPhysicsBlockId<T,Lattice>();
}



////////// class BlockLatticeSerializer3D ////////////////////////////

template<typename T, template<typename U> class Lattice>
BlockLatticeSerializer3D<T,Lattice>::BlockLatticeSerializer3D (
        BlockLattice3D<T,Lattice> const& blockLattice_, IndexOrdering::OrderingT ordering_ )
    : blockLattice(blockLattice_), ordering(ordering_),
      x0(0), x1(blockLattice.getNx()-1), y0(0), y1(blockLattice.getNy()-1), z0(0), z1(blockLattice.getNz()-1),
      iX(x0), iY(y0), iZ(z0),
      sizeOfCell(Lattice<T>::q + Lattice<T>::ExternalField::numScalars)
{ }

template<typename T, template<typename U> class Lattice>
BlockLatticeSerializer3D<T,Lattice>::BlockLatticeSerializer3D (
        BlockLattice3D<T,Lattice> const& blockLattice_,
        int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
        IndexOrdering::OrderingT ordering_ )
    : blockLattice(blockLattice_), ordering(ordering_),
      x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
      iX(x0), iY(y0), iZ(z0),
      sizeOfCell(Lattice<T>::q + Lattice<T>::ExternalField::numScalars)
{ }

template<typename T, template<typename U> class Lattice>
size_t BlockLatticeSerializer3D<T,Lattice>::getSize() const {
    return (size_t)(x1-x0+1) * (size_t)(y1-y0+1) * (size_t)(z1-z0+1) * (size_t)sizeOfCell;
}

template<typename T, template<typename U> class Lattice>
const T* BlockLatticeSerializer3D<T,Lattice>::getNextDataBuffer(size_t& bufferSize) const {
    OLB_PRECONDITION( !isEmpty() );
    if (ordering==IndexOrdering::forward || ordering==IndexOrdering::memorySaving) {
        bufferSize = (size_t)(z1-z0+1)*(size_t)sizeOfCell;
        buffer.resize(bufferSize);
        for (iZ=z0; iZ<=z1; ++iZ) {
            blockLattice.get(iX,iY,iZ).serialize(&buffer[(size_t)(iZ-z0)*(size_t)sizeOfCell]);
        }
        ++iY;
        if (iY > y1) {
            iY = 0;
            ++iX;
        }
    }
    else {
        bufferSize = (size_t)(x1-x0+1)*(size_t)sizeOfCell;
        buffer.resize(bufferSize);
        for (iX=x0; iX<=x1; ++iX) {
            blockLattice.get(iX,iY,iZ).serialize(&buffer[(size_t)(iX-x0)*(size_t)sizeOfCell]);
        }
        ++iY;
        if (iY > y1) {
            iY = 0;
            ++iZ;
        }
    }
    return &buffer[0];
}

template<typename T, template<typename U> class Lattice>
bool BlockLatticeSerializer3D<T,Lattice>::isEmpty() const {
    if (ordering==IndexOrdering::forward || ordering==IndexOrdering::memorySaving) {
        return iX > x1;
    }
    else {
        return iZ > z1;
    }
}



////////// class BlockLatticeUnSerializer3D ////////////////////////////

template<typename T, template<typename U> class Lattice>
BlockLatticeUnSerializer3D<T,Lattice>::BlockLatticeUnSerializer3D (
        BlockLattice3D<T,Lattice>& blockLattice_, IndexOrdering::OrderingT ordering_ )
    : blockLattice(blockLattice_), ordering(ordering_),
      x0(0), x1(blockLattice.getNx()-1), y0(0), y1(blockLattice.getNy()-1), z0(0), z1(blockLattice.getNz()-1),
      iX(x0), iY(y0), iZ(z0),
      sizeOfCell(Lattice<T>::q + Lattice<T>::ExternalField::numScalars)
{ }

template<typename T, template<typename U> class Lattice>
BlockLatticeUnSerializer3D<T,Lattice>::BlockLatticeUnSerializer3D (
        BlockLattice3D<T,Lattice>& blockLattice_,
        int x0_, int x1_, int y0_, int y1_, int z0_, int z1_,
        IndexOrdering::OrderingT ordering_ )
    : blockLattice(blockLattice_), ordering(ordering_),
      x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_),
      iX(x0), iY(y0), iZ(z0),
      sizeOfCell(Lattice<T>::q + Lattice<T>::ExternalField::numScalars)
{ }

template<typename T, template<typename U> class Lattice>
size_t BlockLatticeUnSerializer3D<T,Lattice>::getSize() const {
    return (size_t)(x1-x0+1) * (size_t)(y1-y0+1) * (size_t)(z1-z0+1) * (size_t)sizeOfCell;
}

template<typename T, template<typename U> class Lattice>
T* BlockLatticeUnSerializer3D<T,Lattice>::getNextDataBuffer(size_t& bufferSize) {
    OLB_PRECONDITION( !isFull() );
    if (ordering==IndexOrdering::forward || ordering==IndexOrdering::memorySaving) {
        bufferSize = (size_t)(z1-z0+1)*(size_t)sizeOfCell;
    }
    else {
        bufferSize = (size_t)(x1-x0+1)*(size_t)sizeOfCell;
    }
    buffer.resize(bufferSize);
    return &buffer[0];
}

template<typename T, template<typename U> class Lattice>
void BlockLatticeUnSerializer3D<T,Lattice>::commitData() {
    OLB_PRECONDITION( !isFull() );
    if (ordering==IndexOrdering::forward || ordering==IndexOrdering::memorySaving) {
        for (iZ=z0; iZ<=z1; ++iZ) {
            blockLattice.get(iX,iY,iZ).unSerialize(&buffer[(iZ-z0)*sizeOfCell]);
        }
        ++iY;
        if (iY > y1) {
            iY = 0;
            ++iX;
        }
    }
    else {
        for (iX=x0; iX<=x1; ++iX) {
            blockLattice.get(iX,iY,iZ).unSerialize(&buffer[(iX-x0)*sizeOfCell]);
        }
        ++iY;
        if (iY > y1) {
            iY = 0;
            ++iZ;
        }
    }
}

template<typename T, template<typename U> class Lattice>
bool BlockLatticeUnSerializer3D<T,Lattice>::isFull() const {
    if (ordering==IndexOrdering::forward || ordering==IndexOrdering::memorySaving) {
        return iX > x1;
    }
    else {
        return iZ > z1;
    }
}



//// OpenMP implementation of the method bulkCollideAndStream,
//   by Mathias Krause                                         ////
#ifdef PARALLEL_MODE_OMP
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);

    if (omp.get_size() <= x1-x0+1) {
        #pragma omp parallel
            {
            loadBalancer loadbalance(omp.get_rank(), omp.get_size(), x1-x0+1, x0);
            int iX, iY, iZ, iPop;

            iX=loadbalance.get_firstGlobNum();
            for (int iY=y0; iY<=y1; ++iY) {
                for (int iZ=z0; iZ<=z1; ++iZ) {
                    grid[iX][iY][iZ].collide(getStatistics());
                    grid[iX][iY][iZ].revert();
                    }
            }

            for (iX=loadbalance.get_firstGlobNum()+1; iX<=loadbalance.get_lastGlobNum(); ++iX) {
                for (iY=y0; iY<=y1; ++iY) {
                    for (iZ=z0; iZ<=z1; ++iZ) {
                        grid[iX][iY][iZ].collide(getStatistics());
                        /** The method beneath doesnt work with Intel
                         *  compiler 9.1044 and 9.1046 for Itanium prozessors
                         *    lbHelpers<T,Lattice>::swapAndStream3D(grid, iX, iY, iZ);
                         *  Therefore we use:
                         */
                        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;
                        }
                    }
                }
            }

            #pragma omp barrier
            iX=loadbalance.get_firstGlobNum();
            for (iY=y0; iY<=y1; ++iY) {
                for (iZ=z0; iZ<=z1; ++iZ) {
                    for (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]);
                    }
                }
            }
        }
    }
    else {
        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 // defined PARALLEL_MODE_OMP

}  // namespace olb

#endif