multiblocklattice2d.hh

open lattice boltzmann project www.openlb.org

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2007, 2008 Jonas Latt
 *  Address: Rue General Dufour 24,  1211 Geneva 4, Switzerland 
 *  E-mail: [e-mail]
 *
 *  This program is free software; you can redistribute it and/or
 *  modify it under the terms of the GNU General Public License
 *  as published by the Free Software Foundation; either version 2
 *  of the License, or (at your option) any later version.
 *
 *  This program is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public 
 *  License along with this program; if not, write to the Free 
 *  Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
 *  Boston, MA  02110-1301, USA.
*/

/** \file
 * A 2D multiblock lattice -- generic implementation.
 */
#ifndef MULTI_BLOCK_LATTICE_2D_HH
#define MULTI_BLOCK_LATTICE_2D_HH

#include "multiBlockLattice2D.h"
#include "multiDataAnalysis2D.h"
#include <limits>


namespace olb {

////////////////////// Class MultiBlockLattice2D /////////////////////////

template<typename T, template<typename U> class Lattice>
MultiBlockLattice2D<T,Lattice>::MultiBlockLattice2D(MultiDataDistribution2D const& dataDistribution_)
    : locatedBlock(0),
      statisticsOn(true),
      periodicCommunicationOn(true),
      serializer(0), unSerializer(0),
      serializerPolicy(*this), unSerializerPolicy(*this),
      dataAnalysis(0)
{
#ifdef PARALLEL_MODE_MPI
    multiBlockHandler = new ParallelMultiBlockHandler2D<T,Lattice>(dataDistribution_);
#else
    multiBlockHandler = new SerialMultiBlockHandler2D<T,Lattice>(dataDistribution_);
#endif
    statistics = new LatticeStatistics<T>;
    allocateBlocks();
    eliminateStatisticsInEnvelope();
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock<getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        reductor.startNewSubscription();
        blockLattices[iBlock] -> subscribeReductions(reductor);
    }
    dataAnalysis = new MultiDataAnalysis2D<T,Lattice>(*this);
}

template<typename T, template<typename U> class Lattice>
MultiBlockLattice2D<T,Lattice>::~MultiBlockLattice2D() {
    delete serializer;
    delete unSerializer;
    delete statistics;
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock<getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        delete blockLattices[iBlock];
    }
    delete multiBlockHandler;
    delete dataAnalysis;
}

template<typename T, template<typename U> class Lattice>
MultiBlockLattice2D<T,Lattice>::MultiBlockLattice2D(MultiBlockLattice2D<T,Lattice> const& rhs)
    : locatedBlock(rhs.locatedBlock),
      statisticsOn(rhs.statisticsOn),
      periodicCommunicationOn(rhs.periodicCommunicationOn),
      serializer(0), unSerializer(0),
      serializerPolicy(*this), unSerializerPolicy(*this),
      dataAnalysis(0)
{
    statistics = new LatticeStatistics<T>;
#ifdef PARALLEL_MODE_MPI
    multiBlockHandler = new ParallelMultiBlockHandler2D<T,Lattice> (
            rhs.multiBlockHandler->getMultiDataDistribution() );
#else
    multiBlockHandler = new SerialMultiBlockHandler2D<T,Lattice> (
            rhs.multiBlockHandler->getMultiDataDistribution() );
#endif
    allocateBlocks();
    eliminateStatisticsInEnvelope();
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock<getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        *(blockLattices[iBlock]) = *(rhs.blockLattices[iBlock]);
    }
    for (int rBlock=0; rBlock<getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        reductor.startNewSubscription();
        blockLattices[iBlock] -> subscribeReductions(reductor);
    }
    dataAnalysis = new MultiDataAnalysis2D<T,Lattice>(*this);
}

template<typename T, template<typename U> class Lattice>
MultiBlockLattice2D<T,Lattice>& MultiBlockLattice2D<T,Lattice>::operator= (
        MultiBlockLattice2D<T,Lattice> const& rhs )
{
    MultiBlockLattice2D<T,Lattice>(rhs).swap(*this);
    return *this;
}

template<typename T, template<typename U> class Lattice>
void MultiBlockLattice2D<T,Lattice>::swap(MultiBlockLattice2D<T,Lattice>& rhs) {
    std::swap(locatedBlock, rhs.locatedBlock);
    std::swap(multiBlockHandler, rhs.multiBlockHandler);
    blockLattices.swap(rhs.blockLattices);
    std::swap(statistics, rhs.statistics);
    std::swap(statisticsOn, rhs.statisticsOn);
    std::swap(periodicCommunicationOn, rhs.periodicCommunicationOn);
    std::swap(reductor, rhs.reductor);
    std::swap(serializer, rhs.serializer);
    std::swap(unSerializer, rhs.unSerializer);
    std::swap(dataAnalysis, rhs.dataAnalysis);
}

template<typename T, template<typename U> class Lattice>
Cell<T,Lattice>& MultiBlockLattice2D<T,Lattice>::get(int iX, int iY) {
    std::vector<int> foundId;
    std::vector<int> foundX, foundY;
    bool hasBulkCell;
    locatedBlock = multiBlockHandler -> locateLocally(iX, iY, foundId, foundX, foundY, hasBulkCell, locatedBlock);
    returnCells.clear();
    for (unsigned iBlock=0; iBlock<foundId.size(); ++iBlock) {
        int foundBlock = foundId[iBlock];
        returnCells.push_back ( &blockLattices[foundBlock] -> get ( foundX[iBlock], foundY[iBlock] ) );
    }
    return multiBlockHandler -> getDistributedCell(returnCells, hasBulkCell);
}

template<typename T, template<typename U> class Lattice>
Cell<T,Lattice> const& MultiBlockLattice2D<T,Lattice>::get(int iX, int iY) const {
    std::vector<int> foundId;
    std::vector<int> foundX, foundY;
    bool hasBulkCell;
    locatedBlock = multiBlockHandler -> locateLocally(iX, iY, foundId, foundX, foundY, hasBulkCell, locatedBlock);
    constReturnCells.clear();
    for (unsigned iBlock=0; iBlock<foundId.size(); ++iBlock) {
        int foundBlock = foundId[iBlock];
        constReturnCells.push_back ( &blockLattices[foundBlock] -> get ( foundX[iBlock], foundY[iBlock] ) );
    }
    return multiBlockHandler -> getDistributedCell(constReturnCells, hasBulkCell);
}

template<typename T, template<typename U> class Lattice>
void MultiBlockLattice2D<T,Lattice>::initialize() {
    // Invoke postProcessMultiBlock(), which fills the envelope of
    //   each sub-block. This needs to be done in the first place,
    //   because the method initialize() of each sub-block may want
    //   to access the envelope. An additional postProcessMultiBlock()
    //   is invoked in the end to copy the result of initialize() to
    //   all processors
    postProcessMultiBlock();
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock < getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        blockLattices[iBlock] -> initialize();
    }
    postProcessMultiBlock();
}

template<typename T, template<typename U> class Lattice>
void MultiBlockLattice2D<T,Lattice>::defineDynamics (
        int x0_, int x1_, int y0_, int y1_, Dynamics<T,Lattice>* dynamics )
{
    BlockCoordinates2D domain(x0_, x1_, y0_, y1_), inters;
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock < getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        BlockParameters2D const& params = getParameters(iBlock);
        if (util::intersect(domain, params.getNonPeriodicEnvelope(), inters ) ) {
            inters = params.toLocal(inters);
            blockLattices[iBlock] -> defineDynamics (
                    inters.x0, inters.x1, inters.y0, inters.y1, dynamics );
        }
    }
    std::vector<int> const& periodicOverlapWithMe
        = multiBlockHandler->getRelevantIndexes().getPeriodicOverlapWithMe();
    for (unsigned rOverlap=0; rOverlap < periodicOverlapWithMe.size(); ++rOverlap) {
        int iOverlap = periodicOverlapWithMe[rOverlap];
        Overlap2D const& overlap = getPeriodicOverlap(iOverlap);
        int overlapId = overlap.getOverlapId();
        if (util::intersect(domain, overlap.getOriginalCoordinates(), inters ) ) {
            inters = inters.shift( -overlap.getShiftX(), -overlap.getShiftY() );
            inters = getParameters(overlapId).toLocal(inters);
            blockLattices[overlapId] -> defineDynamics (
                    inters.x0, inters.x1, inters.y0, inters.y1, dynamics );
        }
    }
}

template<typename T, template<typename U> class Lattice>
void MultiBlockLattice2D<T,Lattice>::defineDynamics (
        int iX, int iY, Dynamics<T,Lattice>* dynamics )
{
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock < getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        BlockParameters2D const& params = getParameters(iBlock);
        if (util::contained(iX, iY, params.getNonPeriodicEnvelope()) ) {
            int localX = params.toLocalX(iX);
            int localY = params.toLocalY(iY);
            blockLattices[iBlock] -> defineDynamics (localX,localY, dynamics );
        }
    }
    std::vector<int> const& periodicOverlapWithMe
        = multiBlockHandler->getRelevantIndexes().getPeriodicOverlapWithMe();
    for (unsigned rOverlap=0; rOverlap < periodicOverlapWithMe.size(); ++rOverlap) {
        int iOverlap = periodicOverlapWithMe[rOverlap];
        Overlap2D const& overlap = getPeriodicOverlap(iOverlap);
        int overlapId = overlap.getOverlapId();
        BlockParameters2D const& params = getParameters(overlapId);
        if (util::contained(iX,iY, overlap.getOriginalCoordinates()) ) {
            int localX = params.toLocalX(iX-overlap.getShiftX());
            int localY = params.toLocalY(iY-overlap.getShiftY());
            blockLattices[overlapId] -> defineDynamics (localX, localY, dynamics);
        }
    }
}

template<typename T, template<typename U> class Lattice>
void MultiBlockLattice2D<T,Lattice>::specifyStatisticsStatus (
        int x0_, int x1_, int y0_, int y1_, bool status )
{
    BlockCoordinates2D domain(x0_, x1_, y0_, y1_), inters;
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock < getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        BlockParameters2D const& params = getParameters(iBlock);
        if (util::intersect(domain, params.getBulk(), inters ) ) {
            inters = params.toLocal(inters);
            blockLattices[iBlock] -> specifyStatisticsStatus (
                    inters.x0, inters.x1, inters.y0, inters.y1, status );
        }
    }
}

template<typename T, template<typename U> class Lattice>
void MultiBlockLattice2D<T,Lattice>::collide(int x0_, int x1_, int y0_, int y1_)
{
    BlockCoordinates2D domain(x0_, x1_, y0_, y1_), inters;
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock < getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        BlockParameters2D const& params = getParameters(iBlock);
        if (util::intersect(domain, params.getEnvelope(), inters ) ) {
            inters = params.toLocal(inters);
            blockLattices[iBlock] -> collide(inters.x0, inters.x1, inters.y0, inters.y1);
        }
    }
}

template<typename T, template<typename U> class Lattice>
void MultiBlockLattice2D<T,Lattice>::collide() {
    std::vector<int> const& relevantBlocks = getRelevantBlocks();
    for (int rBlock=0; rBlock < getNumRelevantBlocks(); ++rBlock) {
        int iBlock = relevantBlocks[rBlock];
        blockLattices[iBlock] -> collide();