blocklattice3d.hh

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

/*  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 3D block lattice -- generic implementation.
 */
#ifndef BLOCK_LATTICE_3D_HH
#define BLOCK_LATTICE_3D_HH

#include <algorithm>
#include "blockLattice3D.h"
#include "dynamics.h"
#include "cell.h"
#include "lbHelpers.h"
#include "util.h"
#include "ompManager.h"
#include "loadBalancer.h"
#include "dataAnalysis3D.h"


namespace olb {

////////////////////// Class BlockLattice3D /////////////////////////

/** \param nx_ lattice width (first index)
 *  \param ny_ lattice height (second index)
 *  \param nz_ lattice depth (third index)
 */
template<typename T, template<typename U> class Lattice>
BlockLattice3D<T,Lattice>::BlockLattice3D(int nx_, int ny_, int nz_)
        : nx(nx_), ny(ny_), nz(nz_),
          serializer(0), unSerializer(0),
          dataAnalysis(new DataAnalysis3D<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>
BlockLattice3D<T,Lattice>::~BlockLattice3D()
{
    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>
BlockLattice3D<T,Lattice>::BlockLattice3D(BlockLattice3D<T,Lattice> const& rhs)
    :  serializer(0), unSerializer(0),
       dataAnalysis(new DataAnalysis3D<T,Lattice>(*this) )
{
    nx = rhs.nx;
    ny = rhs.ny;
    nz = rhs.nz;
    allocateMemory();
    resetPostProcessors();
    for (int iX=0; iX<nx; ++iX) {
        for (int iY=0; iY<ny; ++iY) {
            for (int iZ=0; iZ<nz; ++iZ) {
                grid[iX][iY][iZ] = rhs.grid[iX][iY][iZ];
            }
        }
    }
    #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>
BlockLattice3D<T,Lattice>& BlockLattice3D<T,Lattice>::operator= (
        BlockLattice3D<T,Lattice> const& rhs )
{
    BlockLattice3D<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 BlockLattice3D<T,Lattice>::swap(BlockLattice3D& rhs) {
    std::swap(nx, rhs.nx);
    std::swap(ny, rhs.ny);
    std::swap(nz, rhs.nz);
    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 BlockLattice3D<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 BlockLattice3D object, its
 * memory management must be taken care of by the user.
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::defineDynamics (
        int x0, int x1, int y0, int y1, int z0, int z1,
        Dynamics<T,Lattice>* dynamics )
{
    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);

    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].defineDynamics(dynamics);
            }
        }
    }
}

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

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

template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::specifyStatisticsStatus (
        int x0, int x1, int y0, int y1, int z0, int z1, bool status)
{
    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);

    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].specifyStatisticsStatus(status);
            }
        }
    }
}

/** 
 * This method is automatically parallelized if your compiler understands
 * OpenMP
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::collide (
    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);

    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) {
            for (int iZ=z0; iZ<=z1; ++iZ) {
                grid[iX][iY][iZ].collide(getStatistics());
                grid[iX][iY][iZ].revert();
            }
        }
    }
}

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

/** 
 * A useful method for initializing the flow field to a given velocity
 * profile.
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::staticCollide (
        int x0, int x1, int y0, int y1, int z0, int z1,
        TensorFieldBase3D<T,3> const& u )
{
    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);

    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) {
            for (int iZ=z0; iZ<=z1; ++iZ) {
                grid[iX][iY][iZ].staticCollide(u.get(iX,iY,iZ), getStatistics());
                grid[iX][iY][iZ].revert();
            }
        }
    }
}

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

/** The distribution function never leave the rectangular domain. On the
 * domain boundaries, the (outgoing) distribution functions that should 
 * be streamed outside are simply left untouched.
 */
template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::stream(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);

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

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

/** Post-processing steps are called at the end of this method.
 * \sa stream(int,int,int,int,int,int) */
template<typename T, template<typename U> class Lattice>
void BlockLattice3D<T,Lattice>::stream(bool periodic) {
    stream(0, nx-1, 0, ny-1, 0, nz-1);

    if (periodic) {
        makePeriodic();
    }

    postProcess();
}

/** This operation is more efficient than a successive application of
 * collide(int,int,int,int,int,int) and stream(int,int,int,int,int,int),
 * because memory is traversed only once instead of twice.
 */
template<typename T, template<typename U> class Lattice>