boundarypostprocessors3d.hh

open lattice boltzmann project www.openlb.org

/*  This file is part of the OpenLB library
 *
 *  Copyright (C) 2006, 2007 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.
*/

#ifndef BOUNDARY_POST_PROCESSORS_3D_HH
#define BOUNDARY_POST_PROCESSORS_3D_HH

#include "boundaryPostProcessors3D.h"
#include "finiteDifference3D.h"
#include "blockLattice3D.h"
#include "firstOrderLbHelpers.h"
#include "util.h"

namespace olb {

////////  PlaneFdBoundaryProcessor3D ///////////////////////////////////

template<typename T, template<typename U> class Lattice, int direction, int orientation>
PlaneFdBoundaryProcessor3D<T,Lattice,direction,orientation>::
    PlaneFdBoundaryProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
    : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_)
{
    OLB_PRECONDITION(x0==x1 || y0==y1 || z0==z1);
}

template<typename T, template<typename U> class Lattice, int direction, int orientation>
void PlaneFdBoundaryProcessor3D<T,Lattice,direction,orientation>::
    processSubDomain(BlockLattice3D<T,Lattice>& blockLattice,
                     int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
{
    using namespace olb::util::tensorIndices3D;

    int newX0, newX1, newY0, newY1, newZ0, newZ1;
    if ( util::intersect (
                x0, x1, y0, y1, z0, z1,
                x0_, x1_, y0_, y1_, z0_, z1_,
                newX0, newX1, newY0, newY1, newZ0, newZ1 ) )
    {
        int iX;

        #ifdef PARALLEL_MODE_OMP
        #pragma omp parallel for
        #endif
        for (iX=newX0; iX<=newX1; ++iX) {
            T dx_u[Lattice<T>::d], dy_u[Lattice<T>::d], dz_u[Lattice<T>::d];
            for (int iY=newY0; iY<=newY1; ++iY) {
                for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
                    Cell<T,Lattice>& cell = blockLattice.get(iX,iY,iZ);
                    Dynamics<T,Lattice>* dynamics = cell.getDynamics();
                    T rho, u[Lattice<T>:: d];
                    cell.computeRhoU(rho,u);

                    interpolateGradients<0> ( blockLattice, dx_u, iX, iY, iZ );
                    interpolateGradients<1> ( blockLattice, dy_u, iX, iY, iZ );
                    interpolateGradients<2> ( blockLattice, dz_u, iX, iY, iZ );
                    T dx_ux = dx_u[0];
                    T dy_ux = dy_u[0];
                    T dz_ux = dz_u[0];
                    T dx_uy = dx_u[1];
                    T dy_uy = dy_u[1];
                    T dz_uy = dz_u[1];
                    T dx_uz = dx_u[2];
                    T dy_uz = dy_u[2];
                    T dz_uz = dz_u[2];
                    T omega = cell.getDynamics()->getOmega();
                    T sToPi = - rho / Lattice<T>::invCs2 / omega;
                    T pi[util::TensorVal<Lattice<T> >::n];
                    pi[xx] = (T)2 * dx_ux * sToPi;
                    pi[yy] = (T)2 * dy_uy * sToPi;
                    pi[zz] = (T)2 * dz_uz * sToPi;
                    pi[xy] = (dx_uy + dy_ux) * sToPi;
                    pi[xz] = (dx_uz + dz_ux) * sToPi;
                    pi[yz] = (dy_uz + dz_uy) * sToPi;

                    // Computation of the particle distribution functions
                    // according to the regularized formula
                    T uSqr = util::normSqr<T,Lattice<T>::d>(u);
                    for (int iPop = 0; iPop < Lattice<T>::q; ++iPop)
                        cell[iPop] = dynamics -> computeEquilibrium(iPop,rho,u,uSqr) +
                                     firstOrderLbHelpers<T,Lattice>::fromPiToFneq(iPop, pi);
                }
            }
        }
    }
}

template<typename T, template<typename U> class Lattice, int direction, int orientation>
void PlaneFdBoundaryProcessor3D<T,Lattice,direction,orientation>::
    process(BlockLattice3D<T,Lattice>& blockLattice)
{
    processSubDomain(blockLattice, x0, x1, y0, y1, z0, z1);
}


template<typename T, template<typename U> class Lattice, int direction, int orientation>
template<int deriveDirection>
void PlaneFdBoundaryProcessor3D<T,Lattice,direction,orientation>::
    interpolateGradients(BlockLattice3D<T,Lattice> const& blockLattice, T velDeriv[Lattice<T>::d],
                         int iX, int iY, int iZ) const
{
    fd::DirectedGradients3D<T, Lattice, direction, orientation, deriveDirection, direction==deriveDirection>::
        interpolateVector(velDeriv, blockLattice, iX, iY, iZ);
}


////////  PlaneFdBoundaryProcessorGenerator3D ///////////////////////////////

template<typename T, template<typename U> class Lattice, int direction, int orientation>
PlaneFdBoundaryProcessorGenerator3D<T,Lattice,direction,orientation>::
    PlaneFdBoundaryProcessorGenerator3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
    : PostProcessorGenerator3D<T,Lattice>(x0_, x1_, y0_, y1_, z0_, z1_)
{ }

template<typename T, template<typename U> class Lattice, int direction, int orientation>
PostProcessor3D<T,Lattice>* PlaneFdBoundaryProcessorGenerator3D<T,Lattice,direction,orientation>::
        generate() const
{
    return new PlaneFdBoundaryProcessor3D<T,Lattice, direction,orientation>
                   ( this->x0, this->x1, this->y0, this->y1, this->z0, this->z1 );
}

template<typename T, template<typename U> class Lattice, int direction, int orientation>
PostProcessorGenerator3D<T,Lattice>*
    PlaneFdBoundaryProcessorGenerator3D<T,Lattice,direction,orientation>::clone() const
{
    return new PlaneFdBoundaryProcessorGenerator3D<T,Lattice,direction,orientation>
                   (this->x0, this->x1, this->y0, this->y1, this->z0, this->z1);
}
 

////////  OuterVelocityEdgeProcessor3D ///////////////////////////////////

template<typename T, template<typename U> class Lattice, int plane, int normal1, int normal2>
OuterVelocityEdgeProcessor3D<T,Lattice, plane,normal1,normal2>::
    OuterVelocityEdgeProcessor3D(int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
    : x0(x0_), x1(x1_), y0(y0_), y1(y1_), z0(z0_), z1(z1_)
{
    OLB_PRECONDITION (
            (plane==2 && x0==x1 && y0==y1) ||
            (plane==1 && x0==x1 && z0==z1) ||
            (plane==0 && y0==y1 && z0==z1)     );

}

template<typename T, template<typename U> class Lattice, int plane, int normal1, int normal2>
void OuterVelocityEdgeProcessor3D<T,Lattice, plane,normal1,normal2>::
    processSubDomain(BlockLattice3D<T,Lattice>& blockLattice,
                     int x0_, int x1_, int y0_, int y1_, int z0_, int z1_)
{
    using namespace olb::util::tensorIndices3D;

    int newX0, newX1, newY0, newY1, newZ0, newZ1;
    if ( util::intersect ( x0, x1, y0, y1, z0, z1,
                           x0_, x1_, y0_, y1_, z0_, z1_,
                           newX0, newX1, newY0, newY1, newZ0, newZ1 ) )
    {
        int iX;

        #ifdef PARALLEL_MODE_OMP
        #pragma omp parallel for
        #endif
        for (iX=newX0; iX<=newX1; ++iX) {
            for (int iY=newY0; iY<=newY1; ++iY) {
                for (int iZ=newZ0; iZ<=newZ1; ++iZ) {
                    Cell<T,Lattice>& cell = blockLattice.get(iX,iY,iZ);
                    Dynamics<T,Lattice>* dynamics = cell.getDynamics();

                    T rho10 = getNeighborRho(iX,iY,iZ,1,0, blockLattice);
                    T rho01 = getNeighborRho(iX,iY,iZ,0,1, blockLattice);
                    T rho20 = getNeighborRho(iX,iY,iZ,2,0, blockLattice);
                    T rho02 = getNeighborRho(iX,iY,iZ,0,2, blockLattice);
                    T rho = (T)2/(T)3*(rho01+rho10)-(T)1/(T)6*(rho02+rho20);