blocklattice2d.hh

open lattice boltzmann project www.openlb.org

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2006-2008 Jonas Latt
 *  OMP parallel code by Mathias Krause, Copyright (C) 2007
 *  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
 * The dynamics of a 2D block lattice -- generic implementation.
 */
#ifndef BLOCK_LATTICE_2D_HH
#define BLOCK_LATTICE_2D_HH

#include <algorithm>
#include "blockLattice2D.h"
#include "dynamics.h"
#include "cell.h"
#include "lbHelpers.h"
#include "util.h"
#include "ompManager.h"
#include "loadBalancer.h"
#include "dataAnalysis2D.h"

namespace olb {

////////////////////// Class BlockLattice2D /////////////////////////

/** \param nx_ lattice width (first index)
 *  \param ny_ lattice height (second index)
 */
template<typename T, template<typename U> class Lattice>
BlockLattice2D<T,Lattice>::BlockLattice2D(int nx_, int ny_)
        : nx(nx_), ny(ny_),
          serializer(0), unSerializer(0),
          dataAnalysis( new DataAnalysis2D<T,Lattice>(*this) )
{
    allocateMemory();
    resetPostProcessors();
    #ifdef PARALLEL_MODE_OMP
        statistics = new LatticeStatistics<T>* [3*omp.get_size()];
        #pragma omp parallel
        {   statistics[omp.get_rank() + omp.get_size()]
                                       = new LatticeStatistics<T>;
            statistics[omp.get_rank()] = new LatticeStatistics<T>;
            statistics[omp.get_rank() + 2*omp.get_size()]
                                       = new LatticeStatistics<T>;
        }
    #else
        statistics = new LatticeStatistics<T>;
    #endif
}

/** During destruction, the memory for the lattice and the contained
 * cells is released. However, the dynamics objects pointed to by
 * the cells must be deleted manually by the user.
 */
template<typename T, template<typename U> class Lattice>
BlockLattice2D<T,Lattice>::~BlockLattice2D()
{
    releaseMemory();
    clearPostProcessors();
    clearLatticeCouplings();
    #ifdef PARALLEL_MODE_OMP
        #pragma omp parallel
        {
            delete statistics[omp.get_rank()];
        }
        delete statistics;
    #else
        delete statistics;
    #endif
    delete serializer;
    delete unSerializer;
    delete dataAnalysis;
}

/** The whole data of the lattice is duplicated. This includes
 * both particle distribution function and external fields.
 * \warning The dynamics objects and postProcessors are not copied
 * \param rhs the lattice to be duplicated
 */
template<typename T, template<typename U> class Lattice>
BlockLattice2D<T,Lattice>::BlockLattice2D(BlockLattice2D<T,Lattice> const& rhs)
    : serializer(0), unSerializer(0),
      dataAnalysis( new DataAnalysis2D<T,Lattice>(*this) )
{
    nx = rhs.nx;
    ny = rhs.ny;
    allocateMemory();
    resetPostProcessors();
    for (int iX=0; iX<nx; ++iX) {
        for (int iY=0; iY<ny; ++iY) {
            grid[iX][iY] = rhs.grid[iX][iY];
        }
    }
    #ifdef PARALLEL_MODE_OMP
        statistics = new LatticeStatistics<T>* [3*omp.get_size()];
        #pragma omp parallel
        {   statistics[omp.get_rank() + omp.get_size()]
                                       = new LatticeStatistics<T>;
            statistics[omp.get_rank()] = new LatticeStatistics<T> (**rhs.statistics);
            statistics[omp.get_rank() + 2*omp.get_size()]
                                       = new LatticeStatistics<T>;
        }
    #else
        statistics = new LatticeStatistics<T> (*rhs.statistics);
    #endif
}

/** The current lattice is deallocated, then the lattice from the rhs
 * is duplicated. This includes both particle distribution function
 * and external fields. 
 * \warning The dynamics objects and postProcessors are not copied
 * \param rhs the lattice to be duplicated
 */
template<typename T, template<typename U> class Lattice>
BlockLattice2D<T,Lattice>& BlockLattice2D<T,Lattice>::operator= (
        BlockLattice2D<T,Lattice> const& rhs )
{
    BlockLattice2D<T,Lattice> tmp(rhs);
    swap(tmp);
    return *this;
}

/** The swap is efficient, in the sense that only pointers to the 
 * lattice are copied, and not the lattice itself.
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::swap(BlockLattice2D& rhs) {
    std::swap(nx, rhs.nx);
    std::swap(ny, rhs.ny);
    std::swap(rawData, rhs.rawData);
    std::swap(grid, rhs.grid);
    postProcessors.swap(rhs.postProcessors);
    std::swap(statistics, rhs.statistics);
    std::swap(serializer, rhs.serializer);
    std::swap(unSerializer, rhs.unSerializer);
    std::swap(dataAnalysis, rhs.dataAnalysis);
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::initialize() {
    postProcess();
}

/** The dynamics object is not duplicated: all cells of the rectangular
 * domain point to the same dynamics.
 *
 * The dynamics object is not owned by the BlockLattice2D object, its
 * memory management is in charge of the user.
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::defineDynamics (
        int x0, int x1, int y0, int y1, Dynamics<T,Lattice>* dynamics )
{
    OLB_PRECONDITION(x0>=0 && x1<nx);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=0 && y1<ny);
    OLB_PRECONDITION(y1>=y0);

    for (int iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            grid[iX][iY].defineDynamics(dynamics);
        }
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::defineDynamics (
        int iX, int iY, Dynamics<T,Lattice>* dynamics )
{
    OLB_PRECONDITION(iX>=0 && iX<nx);
    OLB_PRECONDITION(iY>=0 && iY<ny);

    grid[iX][iY].defineDynamics(dynamics);
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::specifyStatisticsStatus (
        int x0, int x1, int y0, int y1, bool status )
{
    OLB_PRECONDITION(x0>=0 && x1<nx);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=0 && y1<ny);
    OLB_PRECONDITION(y1>=y0);

    for (int iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            grid[iX][iY].specifyStatisticsStatus(status);
        }
    }
}

template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::collide(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);

    int iX;
    #ifdef PARALLEL_MODE_OMP
    #pragma omp parallel for schedule(dynamic,1)
    #endif
    for (iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            grid[iX][iY].collide(getStatistics());
            grid[iX][iY].revert();
        }
    }
}

/** \sa collide(int,int,int,int) */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::collide() {
    collide(0, nx-1, 0, ny-1);
}

/** 
 * A useful method for initializing the flow field to a given velocity
 * profile.
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::staticCollide(int x0, int x1, int y0, int y1,
                                              TensorFieldBase2D<T,2> const& u)
{
    OLB_PRECONDITION(x0>=0 && x1<nx);
    OLB_PRECONDITION(x1>=x0);
    OLB_PRECONDITION(y0>=0 && y1<ny);
    OLB_PRECONDITION(y1>=y0);

    int iX;

    #ifdef PARALLEL_MODE_OMP
    #pragma omp parallel for schedule(dynamic,1)
    #endif
    for (iX=x0; iX<=x1; ++iX) {
        for (int iY=y0; iY<=y1; ++iY) {
            grid[iX][iY].staticCollide(u.get(iX,iY), getStatistics());
            grid[iX][iY].revert();
        }
    }
}

/** \sa collide(int,int,int,int) */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::staticCollide(TensorFieldBase2D<T,2> const& u) {
    staticCollide(0, nx-1, 0, ny-1, u);
}

/** The distribution functions never leave the rectangular domain. On the
 * domain boundaries, the (outgoing) distribution functions that should
 * be streamed outside are simply left untouched.
 * The post-processing steps are not automatically invoked by this method,
 * as they are in the method stream(). If you want them to be executed, you
 * must explicitly call the method postProcess().
 * \sa stream()
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::stream(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);

    bulkStream(x0+1,x1-1,y0+1,y1-1);

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

/** At the end of this method, the post-processing steps are automatically
 * invoked.
 * \sa stream(int,int,int,int)
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice2D<T,Lattice>::stream(bool periodic) {
    stream(0, nx-1, 0, ny-1);

    if (periodic) {
        makePeriodic();
    }

    postProcess();
}

/** This operation is more efficient than a successive application of