blocklattice2d.hh

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

template<typename T, template<typename U> class Lattice>
template<int normalX, int normalY>
void BlockLattice2D<T,Lattice>::periodicEdge(int from, int to) {
    int x0 = (normalX==1) * (getNx()-1);
    int y0 = (normalY==1) * (getNy()-1);

    int xStep = (normalX==0);
    int yStep = (normalY==0);

    int x = x0 + xStep*from;
    int y = y0 + yStep*from;
    for (; from <= to; ++from) {
        for (int iPop=1; iPop<=Lattice<T>::q/2; ++iPop) {
            int nextX = x + Lattice<T>::c[iPop][0];
            int nextY = y + Lattice<T>::c[iPop][1];
            if ( nextX<0 || nextX>=getNx() ||
                 nextY<0 || nextY>=getNy() )
            {
                nextX = (nextX+getNx())%getNx();
                nextY = (nextY+getNy())%getNy();
                std::swap (
                    grid[x][y]        [iPop+Lattice<T>::q/2],
                    grid[nextX][nextY][iPop] );
            }
        }
        x += xStep;
        y += yStep;
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::makePeriodic() {
    periodicEdge<-1, 0>(0, getNy()-1);
    periodicEdge< 1, 0>(0, getNy()-1);
    periodicEdge< 0,-1>(1, getNx()-2);
    periodicEdge< 0, 1>(1, getNx()-2);
}

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

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

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

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

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

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

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

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

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



////////// class BlockLatticeSerializer2D ////////////////////////////

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

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

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

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

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


////////// class BlockLatticeUnSerializer2D ////////////////////////////

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

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

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

template<typename T, template<typename U> class Lattice>
T* BlockLatticeUnSerializer2D<T,Lattice>::getNextDataBuffer(size_t& bufferSize) {
    OLB_PRECONDITION( !isFull() );
    if (ordering==IndexOrdering::forward || ordering==IndexOrdering::memorySaving) {
        bufferSize = (size_t)(y1-y0+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 BlockLatticeUnSerializer2D<T,Lattice>::commitData() {
    OLB_PRECONDITION( !isFull() );
    if (ordering==IndexOrdering::forward || ordering==IndexOrdering::memorySaving) {
        for (iY=y0; iY<=y1; ++iY) {
            blockLattice.get(iX,iY).unSerialize(&buffer[(size_t)(iY-y0)*(size_t)sizeOfCell]);
        }
        ++iX;
    }
    else {
        for (iX=x0; iX<=x1; ++iX) {
            blockLattice.get(iX,iY).unSerialize(&buffer[(size_t)(iX-x0)*(size_t)sizeOfCell]);
        }
        ++iY;
    }
}

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


//// OpenMP implementation of the method bulkCollideAndStream,
//   by Mathias Krause                                         ////

#ifdef PARALLEL_MODE_OMP
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);

    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, iPop;

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

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

            #pragma omp barrier

            iX=loadbalance.get_firstGlobNum();
            for (iY=y0; iY<=y1; ++iY) {
                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];
                    std::swap(grid[iX][iY][iPop+Lattice<T>::q/2],
                        grid[nextX][nextY][iPop]);
                }
            }
        }
    }
    else {
        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 // defined PARALLEL_MODE_OMP


}  // namespace olb

#endif