multiGridLattice3D.hh

Go to the documentation of this file.

/* This file is part of the Palabos library.
 *
 * Copyright (C) 2011 FlowKit Sarl
 * Avenue de Chailly 23
 * 1012 Lausanne, Switzerland
 * E-mail contact: contact@flowkit.com
 *
 * The most recent release of Palabos can be downloaded at 
 * <http://www.palabos.org/>
 *
 * The library Palabos is free software: you can redistribute it and/or
 * modify it under the terms of the GNU Affero General Public License as
 * published by the Free Software Foundation, either version 3 of the
 * License, or (at your option) any later version.
 *
 * The library 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 Affero General Public License for more details.
 *
 * You should have received a copy of the GNU Affero General Public License
 * along with this program.  If not, see <http://www.gnu.org/licenses/>.
*/

/* Main author: Daniel Lagrava
 **/

#ifndef MULTI_GRID_LATTICE_3D_HH
#define MULTI_GRID_LATTICE_3D_HH

#include "multiGrid/multiGridLattice3D.h"
#include "multiGrid/multiGridParameterManager.h"
#include "io/parallelIO.h"
#include "multiGrid/multiGridGenerator3D.h"


namespace plb {

/* ****** MultiGridLattice3D ****** */

template <typename T, template <typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>::MultiGridLattice3D( MultiGridManagement3D management,
                                                      Dynamics<T,Descriptor>* backgroundDynamics,
                                                      plint behaviorLevel )
       : MultiGrid3D (
             management,
             behaviorLevel )
{
    std::vector<Dynamics<T,Descriptor>*> dynamicsVector(this->getNumLevels());
    for (plint iLevel=0; iLevel<this->getNumLevels(); ++iLevel){
        Dynamics<T,Descriptor>* cloned = backgroundDynamics->clone();
        cloned->rescale(behaviorLevel-iLevel,behaviorLevel-iLevel);
        dynamicsVector[iLevel] = cloned;
    }
    delete backgroundDynamics;
    // generate the lattices that correspond to the blocks contained inside the management object
    lattices = generateLattices( this->getMultiGridManagement(), dynamicsVector,
                                 defaultMultiGridPolicy3D().getBlockCommunicator(management.getNumLevels()),
                                 defaultMultiGridPolicy3D().getCombinedStatistics(management.getNumLevels()) );
                                 
    // initialize the internalStatSubscription (multiGrid) once all the lattices have been created
    this->internalStatSubscription().initialize();    
}

template <typename T, template <typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>::MultiGridLattice3D (
        MultiGridManagement3D management,
        std::vector<BlockCommunicator3D*> communicators_,
        std::vector<CombinedStatistics*> combinedStatistics_, 
        Dynamics<T,Descriptor>* backgroundDynamics,
        plint behaviorLevel )
    : MultiGrid3D( management,
                   behaviorLevel )
{
    
    std::vector<Dynamics<T,Descriptor>*> dynamicsVector(this->getNumLevels());
    for (plint iLevel=0; iLevel<this->getNumLevels(); ++iLevel){
        Dynamics<T,Descriptor>* cloned = backgroundDynamics->clone();
        cloned->rescale(behaviorLevel-iLevel,behaviorLevel-iLevel);
        dynamicsVector[iLevel] = cloned;
    }
    delete backgroundDynamics;
    // generate the lattices that correspond to the blocks contained inside the management object
    lattices = generateLattices( this->getMultiGridManagement(), dynamicsVector,
                                 communicators_, combinedStatistics_ );
                                 
    // initialize the internalStatSubscription (multiGrid) once all the lattices have been created
    this->internalStatSubscription().initialize();
}


template <typename T, template <typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>::~MultiGridLattice3D(){
    for (plint i = 0; i < this->getNumLevels(); ++i){
        delete lattices[i];
    }
}

template <typename T, template <typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>::MultiGridLattice3D(MultiGridLattice3D<T,Descriptor> const& rhs)
  : MultiGrid3D(rhs)
{
    lattices.resize(this->getNumLevels());
    for (plint iLattice=0; iLattice<this->getNumLevels(); ++iLattice) {
        lattices[iLattice] = new MultiBlockLattice3D<T,Descriptor>(rhs.lattices[iLattice]);
    }
}

template <typename T, template <typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>::MultiGridLattice3D( MultiGridLattice3D<T,Descriptor> const& rhs, 
                                                      Box3D subDomain, bool crop )
  : MultiGrid3D(rhs, subDomain, crop)
{
    std::vector<Dynamics<T,Descriptor>*> backgroundDynamics;
    for (plint iDyn=0; iDyn<this->getNumLevels(); ++iDyn) {
      backgroundDynamics[iDyn] = new NoDynamics<T,Descriptor>();
    }
    // generate the lattices that correspond to the blocks contained inside the management object
    lattices = generateLattices( this->getMultiGridManagement(), backgroundDynamics,
                                 this->getCommunicators(), this->getCombinedStatistics() );
    // initialize the internalStatSubscription (multiGrid) once all the lattices have been created
    this->internalStatSubscription().initialize();
}


template <typename T, template <typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>::MultiGridLattice3D(MultiGrid3D const& rhs)
  : MultiGrid3D(rhs)
{
    std::vector<Dynamics<T,Descriptor>*> backgroundDynamics;
    for (plint iDyn=0; iDyn<this->getNumLevels(); ++iDyn) {
      backgroundDynamics[iDyn] = new NoDynamics<T,Descriptor>();
    }
    // generate the lattices that correspond to the blocks contained inside the management object
    lattices = generateLattices( this->getMultiGridManagement(), backgroundDynamics,
                                 this->getCommunicators(), this->getCombinedStatistics() );
    // initialize the internalStatSubscription (multiGrid) once all the lattices have been created
    this->internalStatSubscription().initialize();
}

template <typename T, template <typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>::MultiGridLattice3D( MultiGrid3D const& rhs, 
                                                      Box3D subDomain, bool crop )
  : MultiGrid3D(rhs, subDomain, crop)
{
    std::vector<Dynamics<T,Descriptor>*> backgroundDynamics;
    for (plint iDyn=0; iDyn<this->getNumLevels(); ++iDyn) {
      backgroundDynamics[iDyn] = new NoDynamics<T,Descriptor>();
    }
    // generate the lattices that correspond to the blocks contained inside the management object
    lattices = generateLattices( this->getMultiGridManagement(), backgroundDynamics,
                                 this->getCommunicators(), this->getCombinedStatistics() );
    // initialize the internalStatSubscription (multiGrid) once all the lattices have been created
    this->internalStatSubscription().initialize();
}

template<typename T, template<typename U> class Descriptor>
MultiGridLattice3D<T,Descriptor>& MultiGridLattice3D<T,Descriptor>::operator= (
        MultiGridLattice3D<T,Descriptor> const& rhs )
{
    MultiGridLattice3D<T,Descriptor> tmp(rhs);
    swap(tmp);
    return *this;
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::initialize()
{
    this->createInterfaces();
    
    // Execute this data processor once, because this allocates memory on the fine grid dynamcis
    //   for time t0 (after we have values for t1)
    lattices[0]->executeInternalProcessors();
    std::vector<std::vector<Box3D> > fineGridInterfaces = this->getMultiGridManagement().getFineInterface();
    plint timeForImmediateExecution = 0;
    plint numTimeSteps=2;
    // fineGrid interfaces start at level 1
    for (pluint iLevel=1; iLevel<fineGridInterfaces.size(); ++iLevel){
        std::vector<Box3D> levelInterfaces=fineGridInterfaces[iLevel];
        for (pluint iInterf=0; iInterf<levelInterfaces.size(); ++iInterf){
            applyProcessingFunctional (
                new Copy_t1_to_t0_3D<T, Descriptor>(numTimeSteps, timeForImmediateExecution),
                levelInterfaces[iInterf].multiply(2),
                *lattices[iLevel] );
                
             Box3D box = levelInterfaces[iInterf].multiply(2);
        }
        
        lattices[iLevel]->executeInternalProcessors();
    }
    
    // toggle all stats on
    for (int iLevel = 0; iLevel < (plint)lattices.size(); ++iLevel){
        lattices[iLevel]->initialize();
        lattices[iLevel]->toggleInternalStatistics(true);
    }
    
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::createInterfaces(){
    plb::createInterfaces( lattices, this->getMultiGridManagement());
};

/* *** MultiGrid3D methods *** */
template <typename T, template <typename U> class Descriptor>
MultiBlockLattice3D<T,Descriptor>& MultiGridLattice3D<T,Descriptor>::getComponent(plint iBlock){
    PLB_PRECONDITION( iBlock >= 0 && iBlock < (plint) lattices.size() );
    return *lattices[iBlock];
}

template <typename T, template <typename U> class Descriptor>
const MultiBlockLattice3D<T,Descriptor>& MultiGridLattice3D<T,Descriptor>::getComponent(plint iBlock) const
{
    PLB_PRECONDITION( iBlock >= 0 && iBlock < (plint) lattices.size() );
    return *lattices[iBlock];
}

/* **** BlockLatticeBase3D methods **** */

template <typename T, template <typename U> class Descriptor>
Cell<T,Descriptor>& MultiGridLattice3D<T,Descriptor>::get(plint iX, plint iY, plint iZ){
    return lattices[this->getBehaviorLevel()]->get(iX,iY,iZ);
}

template <typename T, template <typename U> class Descriptor>
Cell<T,Descriptor> const& MultiGridLattice3D<T,Descriptor>::get(plint iX, plint iY, plint iZ) const{
    return lattices[this->getBehaviorLevel()]->get(iX,iY,iZ);
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::specifyStatisticsStatus(Box3D domain, bool status){
    for (plint iLevel = 0; iLevel < (plint)lattices.size(); ++iLevel){
        Box3D rescaledDomain = this->getScaleManager().scaleBox(domain,iLevel);
        lattices[iLevel]->specifyStatisticsStatus(rescaledDomain,status);
    }
}


template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::collide(Box3D domain){
    plint levelNumber = (plint)lattices.size();
    PLB_PRECONDITION( levelNumber >= 0 && levelNumber <= (plint)lattices.size() );
    for (plint iLevel = 0; iLevel < levelNumber; ++iLevel){
        // each level will iterate 2 times the iterations of level - 1
        plint fineGridIt = util::roundToInt(util::twoToThePower(iLevel));
        for (plint iterations = 0; iterations < fineGridIt; ++iterations){
            lattices[iLevel]->collide(domain);
        }
    }
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::collide() {
    plint levelNumber = (plint)lattices.size();
    PLB_PRECONDITION( levelNumber >= 0 && levelNumber <= (plint)lattices.size() );
    for (plint iLevel = 0; iLevel < levelNumber; ++iLevel){
        // each level will iterate 2 times the iterations of level - 1
        plint fineGridIt = util::roundToInt(util::twoToThePower(iLevel));
        for (plint iterations = 0; iterations < fineGridIt; ++iterations){
            lattices[iLevel]->collide();
        }
    }
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::stream(Box3D domain) {
    plint levelNumber = (plint)lattices.size();
    PLB_PRECONDITION( levelNumber >= 0 && levelNumber <= (plint)lattices.size() );
    for (plint iLevel = 0; iLevel < levelNumber; ++iLevel){
        // each level will iterate 2 times the iterations of level - 1
        plint fineGridIt = util::roundToInt(util::twoToThePower(iLevel));
        for (plint iterations = 0; iterations < fineGridIt; ++iterations){
            lattices[iLevel]->stream(domain);
        }
    }
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::stream() {
    plint levelNumber = (plint)lattices.size();
    PLB_PRECONDITION( levelNumber >= 0 && levelNumber <= (plint)lattices.size() );
    for (plint iLevel = 0; iLevel < levelNumber; ++iLevel){
        // each level will iterate 2 times the iterations of level - 1
        plint fineGridIt = util::roundToInt(util::twoToThePower(iLevel));
        for (plint iterations = 0; iterations < fineGridIt; ++iterations){
                lattices[iLevel]->stream();
        }
        // update the values for parallelism
        lattices[iLevel]->getBlockCommunicator().duplicateOverlaps(*lattices[iLevel], modif::staticVariables);
    }
    
    this->evaluateStatistics();
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::collideAndStream(Box3D domain){
    plint levelNumber = (plint)lattices.size();
    PLB_PRECONDITION( levelNumber >= 0 && levelNumber <= (plint)lattices.size() );
    for (plint iLevel = 0; iLevel < levelNumber; ++iLevel){
        // each level will iterate 2 times the iterations of level - 1
        plint fineGridIt = util::roundToInt(util::twoToThePower(iLevel));
        for (plint iterations = 0; iterations < fineGridIt; ++iterations){
            lattices[iLevel]->collideAndStream(domain);
        }
    }
}


template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::iterateMultiGrid(plint level){
    PLB_PRECONDITION( level>=0 && level<(plint)lattices.size() );
    lattices[level]->collideAndStream();
    if ((pluint)level<lattices.size()-1) {
        iterateMultiGrid(level+1);
        iterateMultiGrid(level+1);
        // Overlaps must be duplicated on the coarse lattice, because the
        //   fine->coarse copy acts on bulk nodes only.
        lattices[level]->getBlockCommunicator().duplicateOverlaps(*lattices[level],modif::staticVariables);
    }
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::collideAndStream(){
    iterateMultiGrid(0);
    this->evaluateStatistics();
}

template <typename T, template <typename U> class Descriptor>
void MultiGridLattice3D<T,Descriptor>::incrementTime()
{}

template <typename T, template <typename U> class Descriptor>
TimeCounter& MultiGridLattice3D<T,Descriptor>::getTimeCounter(){
    return lattices[this->getBehaviorLevel()]->getTimeCounter();
}

template <typename T, template <typename U> class Descriptor>
TimeCounter const& MultiGridLattice3D<T,Descriptor>::getTimeCounter() const{
    return lattices[this->getBehaviorLevel()]->getTimeCounter();
}

template <typename T, template <typename U> class Descriptor>
int  MultiGridLattice3D<T,Descriptor>::getBlockId () const{
    return lattices[this->getReferenceLevel()]->getStaticId();
}



template<typename T, template<typename U> class Descriptor>
double getStoredAverageDensity(MultiGridLattice3D<T,Descriptor> const& multiGrid){
    return Descriptor<T>::fullRho (
                 multiGrid.getInternalStatistics().getAverage(
                    LatticeStatistics::avRhoBar ) );
}

template<typename T, template<typename U> class Descriptor>
double getStoredAverageEnergy(MultiGridLattice3D<T,Descriptor> const& multiGrid){
    return 0.5 * multiGrid.getInternalStatistics().getAverage(
                                            LatticeStatistics::avUSqr );
}

template<typename T, template<typename U> class Descriptor>
double getStoredMaxVelocity(MultiGridLattice3D<T,Descriptor> const& multiGrid){
    return std::sqrt( multiGrid.getInternalStatistics().getMax (
                               LatticeStatistics::maxUSqr ) );
}


}// namespace plb

#endif  // MULTI_GRID_LATTICE_3D_HH