amroc/lbm/LBMD3Q19Thermal_BC.h

Go to the documentation of this file.

// -*- C++ -*-

// Copyright (C) 2012 Oak Ridge National Laboratory
// Ralf Deiterding, deiterdingr@ornl.gov

#ifndef LBM_D3Q19_H
#define LBM_D3Q19_H

#include "LBMBase.h"
#include <cfloat>

#define LB_FACTOR 1.0e5

template <class DataType>
class LBMD3Q19Thermal : public LBMBase<Vector<DataType,26>, Vector<DataType,5>, 3> {
  typedef LBMBase<Vector<DataType,26>, Vector<DataType,5>, 3> base;
public:
  typedef typename base::vec_grid_data_type vec_grid_data_type;
  typedef typename base::grid_data_type grid_data_type;
  typedef typename base::MicroType MicroType;
  typedef typename base::MacroType MacroType;
  typedef GridData<MacroType,3> macro_grid_data_type;
  typedef typename base::SideName SideName;
  typedef typename base::point_type point_type;

  enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
  enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall, ExtrapolationEquilibrium, SlipWallTemperature, NoSlipWallTemperature, BounceBackDirichlet };
  enum GFMPredefined { GFMExtrapolation, GFMSlipWallTemp, GFMNoSlipWallTemp, GFMWallLaw, GFMBounceBack };
  enum TurbulenceModel { laminar, LES_Smagorinsky, LES_dynamic };

  LBMD3Q19Thermal() : base(), P0(1.), R0(1.), U0(1.), S0(1.), Temp0(1.), prp(0.), rhop(0.), Tpmin(0.), Tpmax(0.), Tref(0.5), csp(0.), cs2p(0.), nup (0.), nutp (0.), gp(1.), Wp(1.), Rp(1.)  {
    cs2  = DataType(1.)/DataType(3.);
    cs22 = DataType(2.)*cs2;
    cssq = DataType(2.)/DataType(9.);
    ds2  = DataType(1.)/DataType(2.);
    dcs2 = ds2/cs2;
    method[0] = 2;
    turbulence = laminar;

    mdx[0]=mdx[3]=mdx[4]=mdx[5]=mdx[6]=mdx[11]=mdx[12]=mdx[13]=mdx[14]=0;
    mdx[1]=mdx[7]=mdx[9]=mdx[15]=mdx[17]=1;
    mdx[2]=mdx[8]=mdx[10]=mdx[16]=mdx[18]=-1;

    mdy[0]=mdy[1]=mdy[2]=mdy[5]=mdy[6]=mdy[15]=mdy[16]=mdy[17]=mdy[18]=0;
    mdy[3]=mdy[7]=mdy[10]=mdy[11]=mdy[13]=1;
    mdy[4]=mdy[8]=mdy[9]=mdy[12]=mdy[14]=-1;

    mdz[0]=mdz[1]=mdz[2]=mdz[3]=mdz[4]=mdz[7]=mdz[8]=mdz[9]=mdz[10]=0;
    mdz[5]=mdz[11]=mdz[14]=mdz[15]=mdz[18]=1;
    mdz[6]=mdz[12]=mdz[13]=mdz[16]=mdz[17]=-1;

    mdx[19]=mdx[22]=mdx[23]=mdx[24]=mdx[25]=0;
    mdx[20]=1;
    mdx[21]=-1;
   
    mdy[19]=mdy[20]=mdy[21]=mdy[24]=mdy[25]=0;
    mdy[22]=1;
    mdy[23]=-1;

    mdz[19]=mdz[20]=mdz[21]=mdz[22]=mdz[23]=0;
    mdz[24]=1;
    mdz[25]=-1;

  }

  virtual ~LBMD3Q19Thermal() {}

  virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
    base::LocCtrl = Ctrl.getSubDevice(prefix+"D3Q19");
    RegisterAt(base::LocCtrl,"Method(1)",method[0]);
    RegisterAt(base::LocCtrl,"Turbulence",turbulence);
    RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
    RegisterAt(base::LocCtrl,"Gas_rho",rhop);
    RegisterAt(base::LocCtrl,"Gas_Cs",csp);
    RegisterAt(base::LocCtrl,"Gas_pr",prp);
    RegisterAt(base::LocCtrl,"Gas_nu",nup);
    RegisterAt(base::LocCtrl,"Gas_nuT",nutp); 
    RegisterAt(base::LocCtrl,"Gas_Tmin",Tpmin); 
    RegisterAt(base::LocCtrl,"Gas_Tmax",Tpmax); 
    RegisterAt(base::LocCtrl,"Gas_gamma",gp);
    RegisterAt(base::LocCtrl,"Gas_W",Wp);
    RegisterAt(base::LocCtrl,"Ideal_R",Rp);
    RegisterAt(base::LocCtrl,"SpeedUp",S0);
  }

  virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
    base::SetupData(gh,ghosts);
    if (rhop>0.) R0 = rhop;
    if (prp>0.) P0 = prp;
    if (Tpmin!=0. || Tpmax!=0.0) SetTemperatureScale(Tpmin,Tpmax);
    if (csp>0.) {
      cs2p = csp*csp;
      SetTimeScale(LatticeSpeedOfSound()/csp*L0());
    }
    std::cout << "D3Q19Thermal: Gas_rho=" << rhop << "   Gas_Cs=" << csp
              << "   Gas_nu=" << nup << std::endl;
    std::cout << "D3Q19Thermal: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0
              << "   R0=" << R0 << "   Temp0=" << Temp0 << "   S0=" << S0 << "   P0=" << P0 << std::endl;
    if ( TurbulenceType() ==  laminar ) {
      std::cout << "LBM Setup() Laminar" << std::endl;
    }
    if ( TurbulenceType() ==  LES_Smagorinsky ) {
      std::cout << "LBM Setup() LES_Smagorinsky : Smagorinsky Constant = " << SmagorinskyConstant()
                << std::endl;

    }
    if ( TurbulenceType() ==  LES_dynamic ) {
      std::cout << "LBM Setup() LES_dynamic " << std::endl;
    }
  }

  inline virtual MacroType MacroVariables(const MicroType &f) const {
    MacroType q;

    DataType part1,part2;

    q(1)=f(1)-f(2)+f(7)-f(8)+f(9)-f(10)+f(15)-f(16)+f(17)-f(18);
    q(2)=f(3)-f(4)+f(7)-f(8)-f(9)+f(10)+f(11)-f(12)+f(13)-f(14);
    q(3)=f(5)-f(6)+f(11)-f(12)-f(13)+f(14)+f(15)-f(16)-f(17)+f(18);

    part1 = f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8)+f(9)+f(10)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18);
    part2 = ds2*(std::pow(q(1),2)+std::pow(q(2),2)+std::pow(q(3),2));
   
    q(0)= (DataType(1.)/DataType(2.))*(part1-part2);
    q(4)= f(20)+f(21)+f(22)+f(23)+f(24)+f(25);
    return q;
  }

  inline virtual MicroType Equilibrium(const MacroType &q) const {
    MicroType feq;

    DataType p  = q(0);
    DataType u  = q(1);
    DataType v  = q(2);
    DataType w  = q(3);
    DataType T  = q(4);

    DataType ri, rt0, rt1, rt2, pre0, pre1, pre2;
    
    ri = DataType(0.);
    rt0 = DataType(1.) / DataType(3.);
    rt1 = DataType(1.) / DataType(18.);
    rt2 = DataType(1.) / DataType(36.);

    pre0 =  -DataType(2.)*p;
    pre1 =  (DataType(1.)/DataType(6.))*p;
    pre2 =  (DataType(1.)/DataType(12.))*p;


    DataType usq = u * u;
    DataType vsq = v * v;
    DataType wsq = w * w;
    DataType sumsq = (usq + vsq + wsq) / cs22;
    DataType sumsqxy = (usq + vsq) / cs22;
    DataType sumsqyz = (vsq + wsq) / cs22;
    DataType sumsqxz = (usq + wsq) / cs22;
    DataType sumsq2xy = sumsqxy * (DataType(1.) - cs2) / cs2;
    DataType sumsq2yz = sumsqyz * (DataType(1.) - cs2) / cs2;
    DataType sumsq2xz = sumsqxz * (DataType(1.) - cs2) / cs2;
    DataType u2 = usq / cssq;
    DataType v2 = vsq / cssq;
    DataType w2 = wsq / cssq;
    DataType u22 = usq / cs22;
    DataType v22 = vsq / cs22;
    DataType w22 = wsq / cs22;
    DataType ui = u / cs2;
    DataType vi = v / cs2;
    DataType wi = w / cs2;
    DataType uv = ui * vi;
    DataType uw = ui * wi;
    DataType vw = vi * wi;

    DataType Temp = T / DataType(6.);
    DataType uT = DataType(3.) * u;
    DataType vT = DataType(3.) * v;
    DataType wT = DataType(3.) * w;

    feq(0) = pre0 + rt0 * (ri - sumsq);
    //x
    feq(1) = pre1 + rt1 * (ri - sumsq + u2 + ui);
    feq(2) = pre1 + rt1 * (ri - sumsq + u2 - ui);
    //y
    feq(3) = pre1 + rt1 * (ri - sumsq + v2 + vi);
    feq(4) = pre1 + rt1 * (ri - sumsq + v2 - vi);
    //z 
    feq(5) = pre1 + rt1 * (ri - sumsq + w2 + wi);
    feq(6) = pre1 + rt1 * (ri - sumsq + w2 - wi);
    //x-y
    feq(7) = pre2 + rt2 * (ri + sumsq2xy -w22 + ui + vi + uv);
    feq(8) = pre2 + rt2 * (ri + sumsq2xy -w22 - ui - vi + uv);
    feq(9) = pre2 + rt2 * (ri + sumsq2xy -w22 + ui - vi - uv);
    feq(10)= pre2 + rt2 * (ri + sumsq2xy -w22 - ui + vi - uv);
    //y-z
    feq(11)= pre2 + rt2 * (ri + sumsq2yz -u22 + vi + wi + vw);
    feq(12)= pre2 + rt2 * (ri + sumsq2yz -u22 - vi - wi + vw);
    feq(13)= pre2 + rt2 * (ri + sumsq2yz -u22 + vi - wi - vw);
    feq(14)= pre2 + rt2 * (ri + sumsq2yz -u22 - vi + wi - vw);
    //x-z
    feq(15)= pre2 + rt2 * (ri + sumsq2xz -v22 + ui + wi + uw);
    feq(16)= pre2 + rt2 * (ri + sumsq2xz -v22 - ui - wi + uw);
    feq(17)= pre2 + rt2 * (ri + sumsq2xz -v22 + ui - wi - uw);
    feq(18)= pre2 + rt2 * (ri + sumsq2xz -v22 - ui + wi - uw);

    //D3Q6 
    feq(19) = DataType(0.); //Temp; 
    feq(20) = Temp * (DataType(1.) + uT);
    feq(21) = Temp * (DataType(1.) - uT);
    feq(22) = Temp * (DataType(1.) + vT);
    feq(23) = Temp * (DataType(1.) - vT);
    feq(24) = Temp * (DataType(1.) + wT);
    feq(25) = Temp * (DataType(1.) - wT);

    return feq;
  }

  inline virtual void Collision(MicroType &f, const DataType dt) const {
    MacroType q = MacroVariables(f);
    MicroType feq = Equilibrium(q);

    DataType force = DataType(0.5)*dt/S0*gbeta*(q(4)-Tref);                         

    DataType omega = Omega(dt), omegaT = OmegaT(dt);
    
    for(int i=0;i<=25;i++){ 
            if(i<19){
                f(i)=(DataType(1.)-omega)*f(i) + omega*feq(i);
                if(i==3) f(3) += force; 
                if(i==4) f(4) -= force; 
            }
            else f(i)=(DataType(1.)-omegaT)*f(i) + omegaT*feq(i);
    }
  }     

   virtual int IncomingIndices(const int side, int indices[]) const {
    switch (side) {
    case base::Left:
      indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;
      indices[5] = 20;  
      break;
    case base::Right:
      indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;   
      indices[5] = 21;
      break;
    case base::Bottom:
      indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;  
      indices[5] = 22;
      break;
    case base::Top:
      indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;  
      indices[5] = 23;
      break;
    case base::Front:
      indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;  
      indices[5] = 24;
      break;
    case base::Back:
      indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;  
      indices[5] = 25;
      break;
     }
    return 6;
  }

  virtual int OutgoingIndices(const int side, int indices[]) const {
    switch (side) {
    case base::Left:
      indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;   
      indices[5] = 21;
      break;
    case base::Right:
      indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;  
      indices[5] = 20;
      break;
    case base::Bottom:
      indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;  
      indices[5] = 23;
      break;
    case base::Top:
      indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;  
      indices[5] = 22;
      break;
    case base::Front:
      indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;  
      indices[5] = 25;
      break;
    case base::Back:
      indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;  
      indices[5] = 24;
      break;
     }
    return 6;
  }

  // Revert just one layer of cells at the boundary 
  virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
                             const int side) const {
    int Nx = fvec.extents(0), Ny = fvec.extents(1);
    int b = fvec.bottom();
    assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) 
            && fvec.stepsize(2)==bb.stepsize(2));
    int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
    int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
    int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
    int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
    int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
    int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
    MicroType *f = (MicroType *)fvec.databuffer();

#ifdef DEBUG_PRINT_FIXUP
    ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
      << fvec.bbox() << " using " << bb << "\n").flush();
#endif
    
    switch (side) {
    case base::Left:
      for (register int k=mzs; k<=mze; k++) 
        for (register int j=mys; j<=mye; j++) {
          f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](2);
          f[b+base::idx(mxs,j,k,Nx,Ny)](21) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](21);
          f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](8);
          f[b+base::idx(mxs,j,k,Nx,Ny)](10)= f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](10);
          f[b+base::idx(mxs,j,k,Nx,Ny)](16)= f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](16);
          f[b+base::idx(mxs,j,k,Nx,Ny)](18)= f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](18);
        }
      break;
    case base::Right:
      for (register int k=mzs; k<=mze; k++) 
        for (register int j=mys; j<=mye; j++) {
          f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](1);
          f[b+base::idx(mxe,j,k,Nx,Ny)](20) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](20);
          f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](7);
          f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](9);
          f[b+base::idx(mxe,j,k,Nx,Ny)](15)= f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](15);
          f[b+base::idx(mxe,j,k,Nx,Ny)](17)= f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](17);
        }
      break;
    case base::Bottom:
      for (register int k=mzs; k<=mze; k++) 
        for (register int i=mxs; i<=mxe; i++) {
          f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](4);
          f[b+base::idx(i,mys,k,Nx,Ny)](23) = f[b+base::idx(i,mys-1,k,Nx,Ny)](23);
          f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](8);
          f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](9);
          f[b+base::idx(i,mys,k,Nx,Ny)](12)= f[b+base::idx(i,mys-1,k-1,Nx,Ny)](12);
          f[b+base::idx(i,mys,k,Nx,Ny)](14)= f[b+base::idx(i,mys-1,k+1,Nx,Ny)](14);
        }
      break;
    case base::Top:
      for (register int k=mzs; k<=mze; k++) 
        for (register int i=mxs; i<=mxe; i++) {
          f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](3);
          f[b+base::idx(i,mye,k,Nx,Ny)](22) = f[b+base::idx(i,mye+1,k,Nx,Ny)](22);
          f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](7);
          f[b+base::idx(i,mye,k,Nx,Ny)](10)= f[b+base::idx(i-1,mye+1,k,Nx,Ny)](10);
          f[b+base::idx(i,mye,k,Nx,Ny)](11)= f[b+base::idx(i,mye+1,k+1,Nx,Ny)](11);
          f[b+base::idx(i,mye,k,Nx,Ny)](13)= f[b+base::idx(i,mye+1,k-1,Nx,Ny)](13);
        }
      break;
    case base::Front:
      for (register int j=mys; j<=mye; j++) 
        for (register int i=mxs; i<=mxe; i++) {
          f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](6);
          f[b+base::idx(i,j,mzs,Nx,Ny)](25) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](25);
          f[b+base::idx(i,j,mzs,Nx,Ny)](12)= f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](12);
          f[b+base::idx(i,j,mzs,Nx,Ny)](13)= f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](13);
          f[b+base::idx(i,j,mzs,Nx,Ny)](16)= f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](16);
          f[b+base::idx(i,j,mzs,Nx,Ny)](17)= f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](17);
        }
      break;
    case base::Back:
      for (register int j=mys; j<=mye; j++) 
        for (register int i=mxs; i<=mxe; i++) {
          f[b+base::idx(i,j,mzs,Nx,Ny)](5) = f[b+base::idx(i,j,mzs+1,Nx,Ny)](5);
          f[b+base::idx(i,j,mzs,Nx,Ny)](24) = f[b+base::idx(i,j,mzs+1,Nx,Ny)](24);
          f[b+base::idx(i,j,mzs,Nx,Ny)](11)= f[b+base::idx(i,j+1,mzs+1,Nx,Ny)](11);
          f[b+base::idx(i,j,mzs,Nx,Ny)](14)= f[b+base::idx(i,j-1,mzs+1,Nx,Ny)](14);
          f[b+base::idx(i,j,mzs,Nx,Ny)](15)= f[b+base::idx(i+1,j,mzs+1,Nx,Ny)](15);
          f[b+base::idx(i,j,mzs,Nx,Ny)](18)= f[b+base::idx(i-1,j,mzs+1,Nx,Ny)](18);
        }
      break;
    }
  }   

  virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec, 
                         const BBox &bb, const double& dt) const {
    int Nx = fvec.extents(0), Ny = fvec.extents(1);
    int b = fvec.bottom();
    MicroType *f = (MicroType *)fvec.databuffer();
    MicroType *of = (MicroType *)ovec.databuffer();

#ifdef DEBUG_PRINT_FIXUP
    ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox() 
      << " using " << bb << "\n").flush();
#endif
    
    assert (fvec.extents(0)==ovec.extents(0));
    assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) && fvec.stepsize(2)==bb.stepsize(2) &&
            fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1) && fvec.stepsize(2)==ovec.stepsize(2));
    int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
    int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
    int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
    int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
    int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
    int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));

    for (register int k=mzs; k<=mze; k++) 
      for (register int j=mys; j<=mye; j++) 
        for (register int i=mxs; i<=mxe; i++) {
          for (register int n=0; n<base::NMicroVar(); n++) 
            f[b+base::idx(i,j,k,Nx,Ny)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n],Nx,Ny)](n);
          Collision(f[b+base::idx(i,j,k,Nx,Ny)],dt);
        }
  }  

  virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
                      const double& t, const double& dt, const int& mpass) const {
    // Make sure that the physical length and times are scaling consistently
    // into lattice quantities
    DCoords dx = base::GH().worldStep(fvec.stepsize());
    
    //std::cout << "STEP   dt=" << dt << "   TO=" << T0() << std::endl;
    if ( (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR)
         || (std::fabs(dx(1)/L0()-dx(2)/L0()) > DBL_EPSILON*LB_FACTOR) ) {
      std::cerr << "LBMD3Q19 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
                << " and dx(2)=" << dx(2)/L0()
                << " to be equal." << std::endl << std::flush;
      std::exit(-1);
    }
    if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
      std::cerr << "LBMD3Q19 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
                << " to be equal."
                << "   dt=" << dt << "   TO=" << T0()
                << std::endl << std::flush;
      std::exit(-1);
    }

    int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
    MicroType *f = (MicroType *)fvec.databuffer();

    // Streaming - update all cells
    for (register int k=Nz-1; k>=1; k--)
      for (register int j=Ny-1; j>=1; j--)
        for (register int i=0; i<=Nx-2; i++) {
          //x-y
          f[base::idx(i,j,k,Nx,Ny)](3) = f[base::idx(i,j-1,k,Nx,Ny)](3);
          f[base::idx(i,j,k,Nx,Ny)](22) = f[base::idx(i,j-1,k,Nx,Ny)](22);
          f[base::idx(i,j,k,Nx,Ny)](10) = f[base::idx(i+1,j-1,k,Nx,Ny)](10);
          //x-z
          f[base::idx(i,j,k,Nx,Ny)](5) = f[base::idx(i,j,k-1,Nx,Ny)](5);
          f[base::idx(i,j,k,Nx,Ny)](24) = f[base::idx(i,j,k-1,Nx,Ny)](24);
          f[base::idx(i,j,k,Nx,Ny)](18) = f[base::idx(i+1,j,k-1,Nx,Ny)](18);
        }

    for (register int k=0; k<=Nz-2; k++)
      for (register int j=Ny-1; j>=1; j--)
        for (register int i=0; i<=Nx-2; i++) {
          //y-z
          f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);//%%%%% -j,+k
        }
    for (register int k=Nz-1; k>=1; k--)
      for (register int j=Ny-1; j>=1; j--)
        for (register int i=Nx-1; i>=1; i--) {
          //x-y
          f[base::idx(i,j,k,Nx,Ny)](1) = f[base::idx(i-1,j,k,Nx,Ny)](1);
          f[base::idx(i,j,k,Nx,Ny)](20) = f[base::idx(i-1,j,k,Nx,Ny)](20);
          f[base::idx(i,j,k,Nx,Ny)](7) = f[base::idx(i-1,j-1,k,Nx,Ny)](7);
          //y-z
          f[base::idx(i,j,k,Nx,Ny)](11) = f[base::idx(i,j-1,k-1,Nx,Ny)](11);
          //x-z
          f[base::idx(i,j,k,Nx,Ny)](15) = f[base::idx(i-1,j,k-1,Nx,Ny)](15);
        }

    for (register int k=0; k<=Nz-2; k++)
      for (register int j=0; j<=Ny-2; j++)
        for (register int i=Nx-1; i>=1; i--) {
          //x-y
          f[base::idx(i,j,k,Nx,Ny)](4) = f[base::idx(i,j+1,k,Nx,Ny)](4);
          f[base::idx(i,j,k,Nx,Ny)](23) = f[base::idx(i,j+1,k,Nx,Ny)](23);
          f[base::idx(i,j,k,Nx,Ny)](9) = f[base::idx(i-1,j+1,k,Nx,Ny)](9);
          //x-z
          f[base::idx(i,j,k,Nx,Ny)](6) = f[base::idx(i,j,k+1,Nx,Ny)](6);
          f[base::idx(i,j,k,Nx,Ny)](25) = f[base::idx(i,j,k+1,Nx,Ny)](25);
          f[base::idx(i,j,k,Nx,Ny)](17) = f[base::idx(i-1,j,k+1,Nx,Ny)](17);
        }
    for (register int k=Nz-1; k>=1; k--)
      for (register int j=0; j<=Ny-2; j++)
        for (register int i=0; i<=Nx-2; i++) {
          //y-z
          f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);//%%%%% +j,-k
        }
    for (register int k=0; k<=Nz-2; k++)
      for (register int j=0; j<=Ny-2; j++)
        for (register int i=0; i<=Nx-2; i++) {
          //x-y
          f[base::idx(i,j,k,Nx,Ny)](2) = f[base::idx(i+1,j,k,Nx,Ny)](2);
          f[base::idx(i,j,k,Nx,Ny)](21) = f[base::idx(i+1,j,k,Nx,Ny)](21);
          f[base::idx(i,j,k,Nx,Ny)](8) = f[base::idx(i+1,j+1,k,Nx,Ny)](8);
          //y-z
          f[base::idx(i,j,k,Nx,Ny)](12) = f[base::idx(i,j+1,k+1,Nx,Ny)](12);
          //x-z
          f[base::idx(i,j,k,Nx,Ny)](16) = f[base::idx(i+1,j,k+1,Nx,Ny)](16);
        }

     // Collision - only in the interior region
    int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
    int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
    for (register int k=mzs; k<=mze; k++)
      for (register int j=mys; j<=mye; j++)
        for (register int i=mxs; i<=mxe; i++)
          Collision(f[base::idx(i,j,k,Nx,Ny)],dt);

    return DataType(1.);
  }

  virtual void ICStandard(vec_grid_data_type &fvec, const int type,
                          DataType* aux=0, const int naux=0, const int scaling=0) const {
    int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
    int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
    int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
    MicroType *f = (MicroType *)fvec.databuffer();

    //std::cout << "IC Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << "   type=" << type << std::endl;

    if (type==ConstantMicro) {
      //std::cout << "ConstantMicro start\n";
      MicroType g;
      for (register int i=0; i<base::NMicroVar(); i++)
        g(i) = aux[i];

      for (register int k=mzs; k<=mze; k++)
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++)
            f[base::idx(i,j,k,Nx,Ny)] = g;
      //std::cout << "ConstantMicro end\n";
    }
    
    else if (type==ConstantMacro) {
      //std::cout << "ConstantMacro start\n";
      MacroType q;
      if (scaling==base::Physical) {
        q(0) = aux[0]/P0;
        q(1) = aux[1]/U0;
        q(2) = aux[2]/U0;
        q(3) = aux[3]/U0;
        q(4) = aux[4]/Temp0;
      }
      else
        for (register int i=0; i<base::NMacroVar(); i++)
          q(i) = aux[i];
      q(1) *= S0; q(2) *= S0; q(3) *= S0;
      
      for (register int k=mzs; k<=mze; k++)
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++)
            f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
      //std::cout << "ConstantMacro end\n";
    }
    
    else if (type==GasAtRest) {
      //std::cout << "GasAtRest start\n";
      MacroType q;
      q(0) = GasDensity()/R0;
      q(1) = DataType(0.);
      q(2) = DataType(0.);
      q(3) = DataType(0.);
      q(4) = DataType(0.);
      for (register int k=mzs; k<=mze; k++)
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++)
            f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
      //std::cout << "GasAtRest end\n";
    }
  }

  // Set just the one layer of ghost cells at the boundary
  virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
                          DataType* aux=0, const int naux=0, const int scaling=0) const {
    int Nx = fvec.extents(0), Ny = fvec.extents(1);
    int b = fvec.bottom();
    assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) 
            && fvec.stepsize(2)==bb.stepsize(2));
    int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
    int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
    int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
    int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
    int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
    int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
    MicroType *f = (MicroType *)fvec.databuffer();
    
    //std::cout << "BC Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << std::endl;
    if (type==Symmetry) {
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
            f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
            f[b+base::idx(mxe,j,k,Nx,Ny)](20) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](21);
            if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);

            if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
            if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
          }
        break;
      case base::Right:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
            f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
            f[b+base::idx(mxs,j,k,Nx,Ny)](21) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](20);
            if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);

            if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
            if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
          }
        break;
      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
            f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
            f[b+base::idx(i,mye,k,Nx,Ny)](22) = f[b+base::idx(i,mye+1,k,Nx,Ny)](23);
            if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
            if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
          }
        break;
      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
            f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
            f[b+base::idx(i,mys,k,Nx,Ny)](23) = f[b+base::idx(i,mys-1,k,Nx,Ny)](22);
            if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
            if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
            if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
          }
        break;
      case base::Front:
        //std::cout << "SYMMETRY Front Nz=" << Nz << std::endl;
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
            f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
            f[b+base::idx(i,j,mze,Nx,Ny)](24) = f[b+base::idx(i,j,mze+1,Nx,Ny)](25);
            if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
            if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
            if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
          }
        break;
      case base::Back:
        //std::cout << "SYMMETRY Back Nz=" << Nz << std::endl;
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
            f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
            f[b+base::idx(i,j,mzs,Nx,Ny)](25) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](24);
            if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
            if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
            if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
          }
        break;
      }
    }
    
    else if (type==SlipWall) {
      switch (side) {
      case base::Left:
        //std::cout << "S Left Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
            f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
            f[b+base::idx(mxe,j,k,Nx,Ny)](20) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](21);
            if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
            if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
          }
        break;
      case base::Right:
        //std::cout << "S Right Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
            f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
            f[b+base::idx(mxs,j,k,Nx,Ny)](21) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](20);
            if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
            if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
            if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
          }
        break;
      case base::Bottom:
        //std::cout << "S Bottom Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
            f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
            f[b+base::idx(i,mye,k,Nx,Ny)](22) = f[b+base::idx(i,mye+1,k,Nx,Ny)](23);
            if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
            if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
          }
        break;
      case base::Top:
        //std::cout << "S Top Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
            f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
            f[b+base::idx(i,mys,k,Nx,Ny)](23) = f[b+base::idx(i,mys-1,k,Nx,Ny)](22);
            if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
            if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
            if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
          }
        break;
      case base::Front:
        //std::cout << "S Front Nz=" << Nz << std::endl;
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
            f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
            f[b+base::idx(i,j,mze,Nx,Ny)](24) = f[b+base::idx(i,j,mze+1,Nx,Ny)](25);
            if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
            if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
            if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
          }
        break;
      case base::Back:
        //std::cout << "S Back Nz=" << Nz << std::endl;
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
            f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
            f[b+base::idx(i,j,mzs,Nx,Ny)](25) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](24);
            if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
            if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
            if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
          }
        break;
      }
    }

    else if (type==NoSlipWall) {
      switch (side) {
      case base::Left:
        //std::cout << "NS Left Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
            f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
            f[b+base::idx(mxe,j,k,Nx,Ny)](20) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](21);
            if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
            if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
          }
        break;
      case base::Right:
        //std::cout << "NS Right Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
            f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
            f[b+base::idx(mxs,j,k,Nx,Ny)](21) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](20);
            if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
            if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
            if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
          }
        break;
      case base::Bottom:
        //std::cout << "NS Bottom Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
            f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
            f[b+base::idx(i,mye,k,Nx,Ny)](22) = f[b+base::idx(i,mye+1,k,Nx,Ny)](23);
            if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
            if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
          }
        break;
      case base::Top:
        //std::cout << "NS Top Nz=" << Nz << std::endl;
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
            f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
            f[b+base::idx(i,mys,k,Nx,Ny)](23) = f[b+base::idx(i,mys-1,k,Nx,Ny)](22);
            if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
            if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
            if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
          }
        break;
      case base::Front:
        //std::cout << "NS Front Nz=" << Nz << std::endl;
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
            f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
            f[b+base::idx(i,j,mze,Nx,Ny)](24) = f[b+base::idx(i,j,mze+1,Nx,Ny)](25);
            if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
            if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
            if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
          }
        break;
      case base::Back:
        //std::cout << "NS Back Nz=" << Nz << std::endl;
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
            f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
            f[b+base::idx(i,j,mzs,Nx,Ny)](25) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](24);
            if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
            if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
            if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
          }
        break;
      }
    }

    else if (type==Inlet) {
      if (aux==0 || naux<=0) return;
      DataType a0=S0*aux[0], a1=0., a2=0.;
      if (naux>1) {
        a1 = S0*aux[1];
        a2 = S0*aux[2];
      }
      if (scaling==base::Physical) {
        a0 /= U0;
        a1 /= U0;
        a2 /= U0;
      }
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
            q(1) = a0;
            if (naux>1) {
              q(2) = a1;
              q(3) = a2;
            }
            f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Right:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
            q(1) = a0;
            if (naux>1) {
              q(2) = a1;
              q(3) = a2;
            }
            f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
            if (naux==1)
              q(2) = a0;
            else {
              q(1) = a0;
              q(2) = a1;
              q(3) = a2;
            }
            f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
            if (naux==1)
              q(2) = a0;
            else {
              q(1) = a0;
              q(2) = a1;
              q(3) = a2;
            }
            f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Front:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
            if (naux==1)
              q(3) = a0;
            else {
              q(1) = a1;
              q(2) = a2;
              q(3) = a0;
            }
            f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Back:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
            if (naux==1)
              q(3) = a0;
            else {
              q(1) = a1;
              q(2) = a2;
              q(3) = a0;
            }
            f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
          }
        break;
      }
    }

    else if (type==Outlet) {
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
            if (q(0)>0) f[b+base::idx(mxe,j,k,Nx,Ny)] = f[b+base::idx(mxe+1,j,k,Nx,Ny)]/q(0);
            //f[b+base::idx(mxe,j,k,Nx,Ny)] = f[b+base::idx(mxe+1,j,k,Nx,Ny)];
          }
        break;
      case base::Right:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
            if (q(0)>0) f[b+base::idx(mxs,j,k,Nx,Ny)] = f[b+base::idx(mxs-1,j,k,Nx,Ny)]/q(0);
            //f[b+base::idx(mxs,j,k,Nx,Ny)] = f[b+base::idx(mxs-1,j,k,Nx,Ny)];
          }
        break;
      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
            if (q(0)>0) f[b+base::idx(i,mye,k,Nx,Ny)] = f[b+base::idx(i,mye+1,k,Nx,Ny)]/q(0);
            //f[b+base::idx(i,mye,k,Nx,Ny)] = f[b+base::idx(i,mye+1,k,Nx,Ny)];
          }
        break;
      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
            if (q(0)>0) f[b+base::idx(i,mys,k,Nx,Ny)] = f[b+base::idx(i,mys-1,k,Nx,Ny)]/q(0);
            //f[b+base::idx(i,mys,k,Nx,Ny)] = f[b+base::idx(i,mys-1,k,Nx,Ny)];
          }
        break;
      case base::Front:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
            if (q(0)>0) f[b+base::idx(i,j,mze,Nx,Ny)] = f[b+base::idx(i,j,mze+1,Nx,Ny)]/q(0);
            //f[b+base::idx(i,j,mze,Nx,Ny)] = f[b+base::idx(i,j,mze+1,Nx,Ny)];
          }
        break;
      case base::Back:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
            if (q(0)>0) f[b+base::idx(i,j,mzs,Nx,Ny)] = f[b+base::idx(i,j,mzs-1,Nx,Ny)]/q(0);
            //f[b+base::idx(i,j,mzs,Nx,Ny)] = f[b+base::idx(i,j,mzs-1,Nx,Ny)];
          }
        break;
      }
    }
    
    // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)
    else if (type==Pressure) {
      if (aux==0 || naux<=0) return;
      DataType a0=aux[0];
      if (scaling==base::Physical)
        a0 /= R0;
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
            q(0) = a0;
            f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Right:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
            q(0) = a0;
            f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
            q(0) = a0;
            f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
            q(0) = a0;
            f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Front:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
            q(0) = a0;
            f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
          }
        break;
      case base::Back:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
            q(0) = a0;
            f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
          }
        break;
      }
    }
    
    else if (type==SlidingWall) {
      if (aux==0 || naux<=0) return;
      if (aux==0 || naux<=0) return;
      DataType a0=S0*aux[0];
      if (scaling==base::Physical)
        a0 /= U0;
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
            DataType rd1 = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10), rd2 = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
            DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
            DataType pw = DataType(1.)-qw;
            if (j<mye) f[b+base::idx(mxe,j+1,k,Nx,Ny)](9) = pw*rd1+qw*rd2;
            f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
            if (j>mys) f[b+base::idx(mxe,j-1,k,Nx,Ny)](7) = qw*rd1+pw*rd2;

            if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
            if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
          }
        break;
      case base::Right:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
            DataType rd1 = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7), rd2 = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
            DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
            DataType pw = DataType(1.)-qw;
            if (j<mye) f[b+base::idx(mxs,j+1,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
            f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
            if (j>mys) f[b+base::idx(mxs,j-1,k,Nx,Ny)](10) = qw*rd1+pw*rd2;

            if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
            if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
          }
        break;
      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
            DataType rd1 = f[b+base::idx(i,mye+1,k,Nx,Ny)](9), rd2 = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
            DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
            DataType pw = DataType(1.)-qw;
            if (i<mxe) f[b+base::idx(i+1,mye,k,Nx,Ny)](10) = pw*rd1+qw*rd2;
            f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
            if (i>mxs) f[b+base::idx(i-1,mye,k,Nx,Ny)](7) = qw*rd1+pw*rd2;

            if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
            if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
          }
        break;
      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
            DataType rd1 = f[b+base::idx(i,mys-1,k,Nx,Ny)](7), rd2 = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
            DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
            DataType pw = DataType(1.)-qw;
            if (i<mxe) f[b+base::idx(i+1,mys,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
            f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
            if (i>mxs) f[b+base::idx(i-1,mys,k,Nx,Ny)](9) = qw*rd1+pw*rd2;

            if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
            if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
          }
        break;
      case base::Front:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
            DataType rd1 = f[b+base::idx(i,j,mze+1,Nx,Ny)](7), rd2 = f[b+base::idx(i,j,mze+1,Nx,Ny)](10);
            DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
            DataType pw = DataType(1.)-qw;
            if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = pw*rd1+qw*rd2;
            f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
            if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = qw*rd1+pw*rd2;

            if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
            if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
          }
        break;
      case base::Back:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
            DataType rd1 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](7), rd2 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](10);
            DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
            DataType pw = DataType(1.)-qw;
            if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = pw*rd1+qw*rd2;
            f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
            if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = qw*rd1+pw*rd2;

            if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
            if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
          }
        break;
      }
    }

    // (Z.Guo, B. Shi and C.Zheng - A coupled lattice BGK model for the boussinesq equations)     
    else if (type==ExtrapolationEquilibrium) {
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++) 
          for (register int j=mys; j<=mye; j++) {
          MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
          MacroType qBC = q;
          if (aux==0 || naux<=0) return;
          for(int i=0; i<=4; i++){
                if(aux[i] != 0){
                        if(i==0) qBC(i) = aux[i+5]/P0;
                        else if(i>0 && i<4) qBC(i) = aux[i+5]/(U0/S0);
                        else qBC(i)=(aux[i+5]-Tpmin)/Temp0;
                }
          }

          f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(qBC)+f[b+base::idx(mxe+1,j,k,Nx,Ny)]-Equilibrium(q);
        }
        break;

     case base::Right:
        for (register int k=mzs; k<=mze; k++) 
          for (register int j=mys; j<=mye; j++) {
          MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
          MacroType qBC = q;
          if (aux==0 || naux<=0) return;
          for(int i=0; i<=4; i++){
                if(aux[i] != 0){
                        if(i==0) qBC(i) = aux[i+5]/P0;
                        else if(i>0 && i<4) qBC(i) = aux[i+5]/(U0/S0);
                        else qBC(i)=(aux[i+5]-Tpmin)/Temp0;
                }
          }

          f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(qBC)+f[b+base::idx(mxs-1,j,k,Nx,Ny)]-Equilibrium(q);
        }
        break;

      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
            MacroType qBC = q;
            if (aux==0 || naux<=0) return;
          for(int j=0; j<=4; j++){
                if(aux[j] != 0){
                        if(j==0) qBC(j) = aux[j+5]/P0;
                        else if (j>0 && j<4) qBC(j) = aux[j+5]/(U0/S0);
                        else qBC(j)=(aux[j+5]-Tpmin)/Temp0;
                }
          }

          f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(qBC)+f[b+base::idx(i,mye+1,k,Nx,Ny)]-Equilibrium(q);
         }
        break;

      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
            MacroType qBC = q;
            if (aux==0 || naux<=0) return;
          for(int j=0; j<=4; j++){
                if(aux[j] != 0){
                        if(j==0) qBC(j) = aux[j+5]/P0;
                        else if (j>0 && j<4) qBC(j) = aux[j+5]/(U0/S0);
                        else qBC(j)=(aux[j+5]-Tpmin)/Temp0;
                }
          }

          f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(qBC)+f[b+base::idx(i,mys-1,k,Nx,Ny)]-Equilibrium(q);
          }
        break;

      case base::Front:
        for (register int j=mys; j<=mye; j++) 
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
            MacroType qBC = q;
          if (aux==0 || naux<=0) return;
          for(int k=0; k<=4; k++){
                if(aux[k] != 0){
                        if(k==0) qBC(k) = aux[k+5]/P0;
                        else if (k>0 && k<4) qBC(k) = aux[k+5]/(U0/S0);
                        else qBC(k)=(aux[k+5]-Tpmin)/Temp0;
                }
          }

          f[b+base::idx(i,j,mze,Nx,Ny)] = Equilibrium(qBC)+f[b+base::idx(i,j,mze+1,Nx,Ny)]-Equilibrium(q);

        }
        break;
      case base::Back:
        for (register int j=mys; j<=mye; j++) 
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
            MacroType qBC = q;
          if (aux==0 || naux<=0) return;
          for(int k=0; k<=4; k++){
                if(aux[k] != 0){
                        if(k==0) qBC(k) = aux[k+5]/P0;
                        else if (k>0 && k<4) qBC(k) = aux[k+5]/(U0/S0);
                        else qBC(k)=(aux[k+5]-Tpmin)/Temp0;
                }
          }

          f[b+base::idx(i,j,mzs,Nx,Ny)] = Equilibrium(qBC)+f[b+base::idx(i,j,mzs-1,Nx,Ny)]-Equilibrium(q);        
        }
        break;
      }
    }

   else if (type==SlipWallTemperature) {
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
            f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
            f[b+base::idx(mxe,j,k,Nx,Ny)](20) = gBC(20);
            if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
            if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
          }
        break;
      case base::Right:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
            f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
            f[b+base::idx(mxs,j,k,Nx,Ny)](21) = gBC(21);
            if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
            if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
            if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
          }
        break;
      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
            f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
            f[b+base::idx(i,mye,k,Nx,Ny)](22) = gBC(22);
            if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
            if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
          }
        break;
      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
            f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
            f[b+base::idx(i,mys,k,Nx,Ny)](23) = gBC(23);
            if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
            if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
            if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
          }
        break;
      case base::Front:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
            f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
            f[b+base::idx(i,j,mze,Nx,Ny)](24) = gBC(24);
            if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
            if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
            if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
          }
        break;
      case base::Back:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
            f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
            f[b+base::idx(i,j,mzs,Nx,Ny)](25) = gBC(25);
            if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
            if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
            if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
          }
        break;
      }
    }

    else if (type==NoSlipWallTemperature) {
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
            f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
            f[b+base::idx(mxe,j,k,Nx,Ny)](20) = gBC(20);
            if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
            if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
          }
        break;
      case base::Right:
        for (register int k=mzs; k<=mze; k++)
          for (register int j=mys; j<=mye; j++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
            f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
            f[b+base::idx(mxs,j,k,Nx,Ny)](21) = gBC(21);
            if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
            if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
            if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
          }
        break;
      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
            f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
            f[b+base::idx(i,mye,k,Nx,Ny)](22) = gBC(22);
            if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
            if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
            if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
          }
        break;
      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
            f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
            f[b+base::idx(i,mys,k,Nx,Ny)](23) = gBC(23);
            if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
            if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
            if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
          }
        break;
      case base::Front:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
            f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
            f[b+base::idx(i,j,mze,Nx,Ny)](24) = gBC(24);
            if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
            if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
            if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
          }
        break;
      case base::Back:
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType qBC;
            qBC(4) = (aux[0]-Tpmin)/Temp0;
            MicroType gBC = Equilibrium(qBC);
            if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
            f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
            f[b+base::idx(i,j,mzs,Nx,Ny)](25) = gBC(25);
            if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
            if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
            if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
          }
        break;
      }
    }

    else if (type==BounceBackDirichlet) {
      switch (side) {
      case base::Left:
        for (register int k=mzs; k<=mze; k++) 
          for (register int j=mys; j<=mye; j++) {
          MacroType qBC;
          if (aux==0 || naux<=0) return;
          
          qBC(1) = aux[0]/(U0/S0);
          qBC(2) = aux[1]/(U0/S0);
          qBC(3) = aux[2]/(U0/S0);
          qBC(4) = (aux[3]-Tpmin)/Temp0;

          DataType Force = DataType(1.)/DataType(12.);
          DataType ForceT = (DataType(1.)/DataType(3.))*qBC(4);

          if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10)+Force*(qBC(1)-qBC(2));
          f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2)+Force/DataType(2.)*qBC(1);
          f[b+base::idx(mxe,j,k,Nx,Ny)](20) = -f[b+base::idx(mxe+1,j,k,Nx,Ny)](21)+ForceT;
          if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8)+Force*(qBC(1)+qBC(2));
          if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18)+Force*(qBC(1)-qBC(3));
          if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16)+Force*(qBC(1)+qBC(3));        

        }
        break;

     case base::Right:
        for (register int k=mzs; k<=mze; k++) 
          for (register int j=mys; j<=mye; j++) {
          MacroType qBC;
          if (aux==0 || naux<=0) return;
          
          qBC(1) = aux[0]/(U0/S0);
          qBC(2) = aux[1]/(U0/S0);
          qBC(3) = aux[2]/(U0/S0);
          qBC(4) = (aux[3]-Tpmin)/Temp0;

          DataType Force = DataType(1.)/DataType(12.);
          DataType ForceT = (DataType(1.)/DataType(3.))*qBC(4);

          if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7)+Force*(-qBC(1)-qBC(2));
          f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1)+Force/DataType(2.)*(-qBC(1));
          f[b+base::idx(mxs,j,k,Nx,Ny)](21) = -f[b+base::idx(mxs-1,j,k,Nx,Ny)](20)+ForceT;
          if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9)+Force*(-qBC(1)+qBC(2));
          if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15)+Force*(-qBC(1)-qBC(3));
          if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17)+Force*(-qBC(1)+qBC(3));

        }
        break;

      case base::Bottom:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
          MacroType qBC;
          if (aux==0 || naux<=0) return;
          
          qBC(1) = aux[0]/(U0/S0);
          qBC(2) = aux[1]/(U0/S0);
          qBC(3) = aux[2]/(U0/S0);
          qBC(4) = (aux[3]-Tpmin)/Temp0;

          DataType Force = DataType(1.)/DataType(12.);
          DataType ForceT = (DataType(1.)/DataType(3.))*qBC(4);

          if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9)+Force*(-qBC(1)+qBC(2));
          f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4)+Force/DataType(2.)*(qBC(2));
          f[b+base::idx(i,mye,k,Nx,Ny)](22) = -f[b+base::idx(i,mye+1,k,Nx,Ny)](23)+ForceT;
          if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8)+Force*(qBC(1)+qBC(2));
          if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14)+Force*(qBC(2)-qBC(3));
          if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12)+Force*(qBC(2)+qBC(3));


         }
        break;

      case base::Top:
        for (register int k=mzs; k<=mze; k++)
          for (register int i=mxs; i<=mxe; i++) {
          MacroType qBC;
          if (aux==0 || naux<=0) return;
          
          qBC(1) = aux[0]/(U0/S0);
          qBC(2) = aux[1]/(U0/S0);
          qBC(3) = aux[2]/(U0/S0);
          qBC(4) = (aux[3]-Tpmin)/Temp0;

          DataType Force = DataType(1.)/DataType(12.);
          DataType ForceT = (DataType(1.)/DataType(3.))*qBC(4);

          if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7)+Force*(-qBC(1)-qBC(2));
          f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3)+Force/DataType(2.)*(-qBC(2));
          f[b+base::idx(i,mys,k,Nx,Ny)](23) = -f[b+base::idx(i,mys-1,k,Nx,Ny)](22)+ForceT;
          if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10)+Force*(qBC(1)-qBC(2));
          if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11)+Force*(-qBC(2)-qBC(3));
          if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13)+Force*(-qBC(2)+qBC(3));


          }
        break;

      case base::Front:
        for (register int j=mys; j<=mye; j++) 
          for (register int i=mxs; i<=mxe; i++) {
          MacroType qBC;
          if (aux==0 || naux<=0) return;
          
          qBC(1) = aux[0]/(U0/S0);
          qBC(2) = aux[1]/(U0/S0);
          qBC(3) = aux[2]/(U0/S0);
          qBC(4) = (aux[3]-Tpmin)/Temp0;

          DataType Force = DataType(1.)/DataType(12.);
          DataType ForceT = (DataType(1.)/DataType(3.))*qBC(4);

          if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17)+Force*(-qBC(1)+qBC(3));
          f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6)+Force/DataType(2.)*qBC(3);
          f[b+base::idx(i,j,mze,Nx,Ny)](24) = -f[b+base::idx(i,j,mze+1,Nx,Ny)](25)+ForceT;
          if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16)+Force*(qBC(1)+qBC(3));
          if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13)+Force*(-qBC(2)+qBC(3));
          if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12)+Force*(qBC(2)+qBC(3));

        }
        break;

      case base::Back:
        for (register int j=mys; j<=mye; j++) 
          for (register int i=mxs; i<=mxe; i++) {
          MacroType qBC;
          if (aux==0 || naux<=0) return;
          
          qBC(1) = aux[0]/(U0/S0);
          qBC(2) = aux[1]/(U0/S0);
          qBC(3) = aux[2]/(U0/S0);
          qBC(4) = (aux[3]-Tpmin)/Temp0;

          DataType Force = DataType(1.)/DataType(12.);
          DataType ForceT = (DataType(1.)/DataType(3.))*qBC(4);

          if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15)+Force*(-qBC(1)-qBC(3));
          f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5)+Force/DataType(2.)*(-qBC(3));
          f[b+base::idx(i,j,mzs,Nx,Ny)](25) = -f[b+base::idx(i,j,mzs-1,Nx,Ny)](24)+ForceT;
          if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18)+Force*(qBC(1)-qBC(3));
          if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11)+Force*(-qBC(2)-qBC(3));
          if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14)+Force*(qBC(2)-qBC(3));
    
        }
        break;
      }
    }
  }

  virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
                             const MicroType* f, const point_type* xc, const DataType* distance,
                             const point_type* normal, DataType* aux=0, const int naux=0,
                             const int scaling=0) const {
    if (type==GFMNoSlipWallTemp || type==GFMSlipWallTemp) {
      DataType fac = S0;
      if (scaling==base::Physical)
        fac /= U0;
      for (int n=0; n<nc; n++) {
        MacroType q=MacroVariables(f[n]);
        DataType u=-q(1);
        DataType v=-q(2);
        DataType w=-q(3);
        DataType T= q(4);

        // Add boundary velocities if available
        if (naux>=4) {
          u += fac*aux[naux*n];
          v += fac*aux[naux*n+1];
          w += fac*aux[naux*n+2];
          T = aux[naux*n+3];
        }

        if (type==GFMNoSlipWallTemp) {
          // Prescribe entire velocity vector in ghost cells
          q(1) += DataType(2.)*u;
          q(2) += DataType(2.)*v;
          q(3) += DataType(2.)*w;
          q(4)  = T;
        }
        else if (type==GFMSlipWallTemp) {
          // Prescribe only normal velocity vector in ghost cells
          DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1)+w*normal[n](2));
          q(1) += vl*normal[n](0);
          q(2) += vl*normal[n](1);
          q(3) += vl*normal[n](2);
          q(4)  = T;
        }
        fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
      }
    }

    else if (type==GFMExtrapolation) {
      for (int n=0; n<nc; n++) {
        MacroType q=MacroVariables(f[n]);
        DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1)+q(3)*normal[n](2);
        q(1) = vl*normal[n](0);
        q(2) = vl*normal[n](1);
        q(3) = vl*normal[n](2);
        fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
      }
    }

    else if (type==GFMWallLaw ) {
      DCoords dx = base::GH().worldStep(fvec.stepsize());
      DataType fac = S0;
      if (scaling==base::Physical)
        fac /= U0;

      // DataType BLH = 10.;
      DataType yp_ = 0.,  yp = 0., yMax = 50.;
      DataType yPlus[3];
      int lc_low[3], lc_up[3], yPlus_lc[3];

      // int mxe = fvec.upper(0)/fvec.stepsize(0);
      // int mye = fvec.upper(1)/fvec.stepsize(1);
      // int mze = fvec.upper(2)/fvec.stepsize(2);

      DataType shear_vel = 1., karman = 0.41, y0 = 1., law = 1.;

      DataType shearStress;
      DataType u, v, w, u_, v_, w_, vel0;
      // DataType vel1;
      MacroType q, q_;
      int kMax = 0;

      for (int n=0; n<nc; n++) {

        q=MacroVariables(f[n]);
        q_=MacroVariables(f[n]);

        yPlus[0] = (xc[n](0) - DataType(2.)*dx(0)*normal[n](0));
        yPlus[1] = (xc[n](1) - DataType(2.)*dx(1)*normal[n](1));
        yPlus[2] = (xc[n](2) - DataType(2.)*dx(2)*normal[n](2));

        yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
        yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
        yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);

        lc_low[0] = fvec.lower(0);   lc_low[1] = fvec.lower(1);   lc_low[2] = fvec.lower(2);
        lc_up[0] = fvec.upper(0);   lc_up[1] = fvec.upper(1);   lc_up[2] = fvec.upper(2);

        //std::cout << "xc " << xc[n] << " yPlus " << yPlus[0] << " " << yPlus[1] << " " << yPlus[2] << std::endl;

        if ( yPlus_lc[0] >  lc_low[0] ) {
          if ( yPlus_lc[1] >  lc_low[1] ) {
            if ( yPlus_lc[2] >  lc_low[2] ) {

              if ( yPlus_lc[0] < lc_up[0] ) {
                if ( yPlus_lc[1] <  lc_up[1] ) {
                  if ( yPlus_lc[2] <  lc_up[2] ) {
                    q=MacroVariables( fvec( base::GH().localCoords((const DataType *) yPlus) ) );
                  }
                }
              }

            }
          }
        }

        u_=-q(1);
        v_=-q(2);
        w_=-q(3);
        // vel1 = sqrt(u_*u_+v_*v_+w_*w_);

        q=MacroVariables(f[n]);
        u=-q(1);
        v=-q(2);
        w=-q(3);
        vel0 = sqrt(u*u+v*v+w*w);

        u_+=q(1);
        v_+=q(2);
        w_+=q(3);
        vel0 = ( u_*normal[n](0) + v_*normal[n](1) + w_*normal[n](2) );

        u_ =  u_ - vel0 * normal[n](0);
        v_ =  v_ - vel0 * normal[n](1);
        w_ =  w_ - vel0 * normal[n](2);
        vel0 = sqrt(u_*u_ + v_*v_ + w_*w_)*U0/S0;

        shearStress = (DataType(0.5)*vel0/dx(0) * q(0)*R0*GasViscosity() );

        // Add boundary velocities if available
        if (naux>=3) {
          u += fac*aux[naux*n];
          v += fac*aux[naux*n+1];
          w += fac*aux[naux*n+2];
        }
        q(1) += DataType(2.)*u;
        q(2) += DataType(2.)*v;
        q(3) += DataType(2.)*w;
        fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);

        shear_vel = std::sqrt( shearStress/(q(0)*R0) );

        y0 = karman/30.;

        //std::cout << n << " vel0 " << vel0 << " vel1 " << vel1  << " shear_vel = " << shear_vel << std::endl;

        if ( shear_vel > 1.e-5 ) {
          kMax = (int) std::ceil( (yMax*GasViscosity())/(dx(0)*shear_vel) );

          for (int k=1; k < kMax; k++) {

            yPlus[0] = (xc[n](0) - k*dx(0)*normal[n](0));
            yPlus[1] = (xc[n](1) - k*dx(1)*normal[n](1));
            yPlus[2] = (xc[n](2) - k*dx(2)*normal[n](2));

            yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
            yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
            yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);

            if ( yPlus_lc[0] >  lc_low[0] ) {
              if ( yPlus_lc[1] >  lc_low[1] ) {
                if ( yPlus_lc[2] >  lc_low[2] ) {

                  if ( yPlus_lc[0] < lc_up[0] ) {
                    if ( yPlus_lc[1] <  lc_up[1] ) {
                      if ( yPlus_lc[2] <  lc_up[2] ) {
                        yp = k*dx(0);
                        yp_ = yp*shear_vel/GasViscosity();

                        if ( yp_ <= 5.0 ) {
                          law = yp_;
                        }
                        else if (yp_ > 5 && yp_ < 30.0 ) {
                          //law = (shear_vel/karman*std::log( yp/y0 ) + yp_)*DataType(0.5);
                          law = std::min( yp_ , shear_vel/karman*std::log( yp/y0 ) );

                        }
                        else if (yp_ >= 30.0 ) {
                          law = shear_vel/karman*std::log( yp/y0 );
                        }
                        //std::cout << "k = " << k << " yp = " << yp << " shear_vel = " << shear_vel << " yp_ = " << yp_ << " law = " << law << std::endl;
                        q_=MacroVariables(f[n]);
                        u = law*(-q_(1) + fac*aux[naux*n]);
                        v = law*(-q_(2) + fac*aux[naux*n+1]);
                        w = law*(-q_(3) + fac*aux[naux*n+2]);

                        q_(1) += DataType(2.)*u;
                        q_(2) += DataType(2.)*v;
                        q_(3) += DataType(2.)*w;
                        fvec( base::GH().localCoords((const DataType *) yPlus) ) = Equilibrium(q_);
                      }
                    }
                  }


                }
              }
            }

          }

        }

      }
      //std::cout << "====== WallLaw Done\n";
    }


    else if (type==GFMBounceBack ) {

     DCoords dx = base::GH().worldStep(fvec.stepsize());
     DataType fac = S0;
     if(scaling==base::Physical) fac/=U0;

      DataType dt_temp = dx(0)*S0/U0;
   
      //xc : Cell midpoints of the ghost cells
      //nc : Number of ghost cells
      //normal : normalized normal vectors, pointed inside
      //distance: Level-Set distance between ghost cell midpoint and boundary
       
      DataType distB[3], distBScalar, cosa, leng, d_dist, qNorm, mdx2, mdy2, mdz2, zaehler, nenner, lengXi;
        
      for (int n=0; n<nc; n++){
        distB[0] = (double &)distance[n]*normal[n](0)*dx(0);
        distB[1] = (double &)distance[n]*normal[n](1)*dx(1);
        distB[2] = (double &)distance[n]*normal[n](2)*dx(2);

        
        //std::ofstream deb("debugOut.txt", std::ios::app);


        //alpha: Angle between discrete velocity direction (Ghostcell) and inner normal direction (Boundary)
        //leng: Distance between ghost cell midpoint and curved boundary along discrete velocity direction
        //d_dist: Distance between fluid-cellmidpoint and boundary along discrete velocity direction
        //qNorm: Normalized length from fluid-cellmidpoint to boundary along discrete velocity direction
 
        for(int ind=0; ind<base::NMicroVar(); ind++){
                
                mdx2 = mdx[ind]*mdx[ind];
                mdy2 = mdy[ind]*mdy[ind];
                mdz2 = mdz[ind]*mdz[ind];
                zaehler = (distB[0]*mdx[ind]+distB[1]*mdy[ind]+distB[2]*mdz[ind]);
                distBScalar = std::sqrt(std::pow(distB[0],2)+std::pow(distB[1],2)+std::pow(distB[2],2));
                nenner = std::sqrt(mdx2+mdy2+mdz2)*distBScalar;
                lengXi = std::sqrt(mdx2*dx(0)*dx(0)+mdy2*dx(1)*dx(1)+mdz2*dx(2)*dx(2));

                if(nenner!=0) cosa = zaehler/nenner;
                else cosa = 0;
                if(cosa!=0) {
                        leng = distBScalar*(DataType(1.)/cosa);
                        d_dist = lengXi-leng;
                        qNorm = d_dist/lengXi;
                }
                else {
                        leng = 0;
                        d_dist = 0;
                        qNorm = 0;
                }
                
                if(qNorm>0 && qNorm < std::sqrt(2.)){
                //The positive lengths are important only       
                        
                        //Two cases: qNorm greater equal or lower than 0.5
                        int idxA, idxB;
                        if(ind<19){
                                if(ind==0)      idxA=0;
                                else if(ind % 2 == 0 && ind <19) idxA = ind-1;
                                else idxA = ind+1;
                        }
                        else{
                                if(ind==19) idxA=19;
                                else if(ind % 2 == 0 && ind >19) idxA = ind+1;
                                else idxA = ind-1;
                        }

                        //std::cout << "ind = " << ind <<  ", idxA = " << idxA << ", idxB = " << idxB << std::endl;
                        idxB=idxA;
                        if(qNorm>=DataType(0.5)) idxB=ind;

                        //Calculate Coordinates of needed fluidcell                     
                        DataType cn[3], cnA[3], cnB[3];
                        cn[0] = xc[n](0);
                        cn[1] = xc[n](1);
                        cn[2] = xc[n](2);

                        //cnA = cn; //Coords of fluidcell
                        cnA[0] = cn[0] + mdx[ind]*dx(0);
                        cnA[1] = cn[1] + mdy[ind]*dx(1);
                        cnA[2] = cn[2] + mdz[ind]*dx(2);

                        cnB[0] = cnA[0]; //Coords of 2nd fluidcell
                        cnB[1] = cnA[1]; 
                        cnB[2] = cnA[2]; 

                        if(qNorm<DataType(0.5) && qNorm > DataType(0.)){
                                cnB[0] += mdx[ind]*dx(0);
                                cnB[1] += mdy[ind]*dx(1);
                                cnB[2] += mdz[ind]*dx(2);
                        }

                        //Getting the corresponding partial probability distribution functions
                        DataType distFuncA, distFuncB;

                        distFuncA = fvec(base::GH().localCoords((const DataType*) cnA))(idxA);
                        distFuncB = fvec(base::GH().localCoords((const DataType*) cnB))(idxB);
        
                        DataType intp1,intp2,intp,q2;
                        q2 = DataType(2.)*qNorm;

                        //Using formular from Eitel-Amor for velocity
                        if(ind<19){     
                                if(qNorm<DataType(0.5) && qNorm >= DataType(0.)){
                                        intp1=q2*distFuncA;
                                        intp2=(DataType(1.)-q2)*distFuncB;
                                }
                                else{
                                        if(q2!=0){
                                                intp1 = (DataType(1.)/q2)*distFuncA;
                                                intp2 = ((q2-DataType(1.))/q2)*distFuncB;
                                        }
                                        else{
                                        intp1 = DataType(0.);
                                        intp2 = intp1;
                                        }
                                }
                        }
                        //Using formular from L. Li et al. for temperature distributions
                        else{           
                                if(qNorm<=DataType(0.5) && qNorm >= DataType(0.)){
                                        intp1 = -q2*distFuncA;
                                        intp2 = (q2-DataType(1.))*distFuncB;
                                }
                                else{
                                        if(q2!=0){
                                                intp1 = -(DataType(1.)/q2)*distFuncA;
                                                intp2 = (DataType(1.)-(DataType(1.)/q2))*distFuncB;
                                        }
                                        else{
                                                intp1 = DataType(0.);
                                                intp2 = intp1;
                                        }
                                }
                        }

                        intp = intp1+intp2;
                        
                        fvec(base::GH().localCoords((const DataType *) cn))(ind) = intp;

                        //ADDITION PART FOR EXTERNAL FORCE LIKE VELOCITY OR TEMPERATURE
                        //From: Momentum transfer of a Boltzmann-lattice fluid with boundaries - M. Bouzidi et al. (2001) Physics of Fluids

                        if(naux>=3 && ind<19){
                                DataType u = fac*aux[naux*n];
                                DataType v = fac*aux[naux*n+1];
                                DataType w = fac*aux[naux*n+2];
        
                                DataType weight;
                                weight = DataType(1.)/DataType(12.);
                                if(ind==0) weight *= 0;
                                else if(ind>=1 && ind<=6) weight *= DataType(2.);
        
                                DataType velForce;
                                        
                                if(qNorm<DataType(0.5)) velForce = DataType(2.)*weight;
                                else velForce = (DataType(1.)/qNorm)*weight;
                                        
                                velForce *= ((mdx[ind]*u)+(mdy[ind]*v)+(mdz[ind]*w));
                                fvec(base::GH().localCoords((const DataType*) cn))(ind) +=velForce;
                        }

                        if(naux>=4 && ind>=19){

                                DataType T = (aux[naux*n+3]-Tpmin)/(Tpmax-Tpmin);
                                DataType TForce = DataType(1.)/DataType(3.);

                                if(qNorm>DataType(0.5)) TForce /= q2;

                                fvec(base::GH().localCoords((const DataType*) cn))(ind) +=TForce*T;
                        
                        }       
                }       
        }
      }

    }//End GFMBounceBack

  }

  virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
                      const int skip_ghosts=1) const {
    int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
    int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
    int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
    MicroType *f = (MicroType *)fvec.databuffer();
    DataType *work = (DataType *)workvec.databuffer();
    DCoords dx = base::GH().worldStep(fvec.stepsize());
    DataType dt_temp = dx(0)*S0/U0;

    if (cnt>=1 && cnt<=10) { 
      for (register int k=mzs; k<=mze; k++)
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            MacroType q=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
            switch(cnt) {
            case 1:
              work[base::idx(i,j,k,Nx,Ny)]=0;
              break;
            case 2:
              work[base::idx(i,j,k,Nx,Ny)]=q(1)*U0/S0;
              break;
            case 3:
              work[base::idx(i,j,k,Nx,Ny)]=q(2)*U0/S0;
              break;
            case 4:
              work[base::idx(i,j,k,Nx,Ny)]=q(3)*U0/S0;
              break;
            case 6://Temperature
              work[base::idx(i,j,k,Nx,Ny)]=q(4)*(Tpmax-Tpmin)+Tpmin;
              break;
            case 7://Pressure
              work[base::idx(i,j,k,Nx,Ny)]=q(0)*P0; 
              break;
            case 9: {
              MicroType feq = Equilibrium(q);     
              work[base::idx(i,j,k,Nx,Ny)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt_temp),
                                                        GasSpeedofSound(),dt_temp );
              break;
            }
            case 10:
              work[base::idx(i,j,k,Nx,Ny)]=std::sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3))*U0/S0;
              break;
            }     
          }
    }
    else if (cnt>=11 && cnt<=14) {
      macro_grid_data_type qgrid(fvec.bbox());
      MacroType *q = (MacroType *)qgrid.databuffer();
      for (register int k=mzs; k<=mze; k++)
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) {
            q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
            q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
            q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
            q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
          }
      
      if (skip_ghosts==0) {
        switch(cnt) {
        case 11:
          mys+=1; mye-=1; mzs+=1; mze-=1;
          break;
        case 12:
          mxs+=1; mxe-=1; mzs+=1; mze-=1;
          break;
        case 13:
          mxs+=1; mxe-=1; mys+=1; mye-=1; 
          break;
        case 14:
          mxs+=1; mxe-=1; mys+=1; mye-=1; mzs+=1; mze-=1;
          break;
        }
      }

      for (register int k=mzs; k<=mze; k++)
        for (register int j=mys; j<=mye; j++)
          for (register int i=mxs; i<=mxe; i++) 
            switch(cnt) {
            case 11: 
              work[base::idx(i,j,k,Nx,Ny)]=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(2.*dx(1))-
                (q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(2.*dx(2));
              break;
            case 12: 
              work[base::idx(i,j,k,Nx,Ny)]=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(2.*dx(2))-
                (q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(2.*dx(0));
              break;
            case 13: 
              work[base::idx(i,j,k,Nx,Ny)]=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(2.*dx(0))-
                (q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(2.*dx(1));
              break;
            case 14:
              DataType ox = (q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(2.*dx(1))-
                (q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(2.*dx(2));
              DataType oy = (q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(2.*dx(2))-
                (q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(2.*dx(0));
              DataType oz = (q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(2.*dx(0))-
                (q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(2.*dx(1));
              work[base::idx(i,j,k,Nx,Ny)] = std::sqrt(ox*ox+oy*oy+oz*oz);
              break;
            }
    }
  }

  virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
                     const int skip_ghosts=1) const {
    int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
    int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
    int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
    MicroType *f = (MicroType *)fvec.databuffer();
    DataType *work = (DataType *)workvec.databuffer();
    DCoords dx = base::GH().worldStep(fvec.stepsize());

    for (register int k=mzs; k<=mze; k++)
      for (register int j=mys; j<=mye; j++)
        for (register int i=mxs; i<=mxe; i++) {
          switch(cnt) {
          case 2:
            f[base::idx(i,j,k,Nx,Ny)](1)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
            break;
          case 3:
            f[base::idx(i,j,k,Nx,Ny)](2)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
            break;
          case 4:
            f[base::idx(i,j,k,Nx,Ny)](3)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
            f[base::idx(i,j,k,Nx,Ny)] = Equilibrium((const MacroType &)f[base::idx(i,j,k,Nx,Ny)]);
            break;
          case 6://Temperature
            f[base::idx(i,j,k,Nx,Ny)](3)=(work[base::idx(i,j,k,Nx,Ny)]-Tpmin)/(Tpmax-Tpmin);
            break;
          case 7://Pressure
            f[base::idx(i,j,k,Nx,Ny)](0)=work[base::idx(i,j,k,Nx,Ny)]/P0;
            break;
          }
        }
  }

  virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
                    const int verbose) const {
    int Nx = fvec.extents(0), Ny = fvec.extents(1);
    int b = fvec.bottom();
    assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
    BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
    int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
    int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
    int mzs = std::max(bb.lower(2)/bb.stepsize(2),bbmax.lower(2)/bbmax.stepsize(2));
    int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
    int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
    int mze = std::min(bb.upper(2)/bb.stepsize(2),bbmax.upper(2)/bbmax.stepsize(2));
    MicroType *f = (MicroType *)fvec.databuffer();

    //std::cout << "CHECK Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << std::endl;
    int result = 1;
    for (register int k=mzs; k<=mze; k++)
      for (register int j=mys; j<=mye; j++)
        for (register int i=mxs; i<=mxe; i++)
          for (register int l=0; l<base::NMicroVar(); l++)
            if (!(f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l)>-1.e37)) {
              result = 0;
              if (verbose)
                std::cerr << "D3Q19-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << "," << k-mdz[l] << ")(" << l << ")="
                          << f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l) << " " << fvec.bbox() << std::endl;
            }
    return result;
  }

  virtual int NMethodOrder() const { return 2; }

  inline const DataType& L0() const { return base::LengthScale(); }
  inline const DataType& T0() const { return base::TimeScale(); }
  inline void SetDensityScale(const DataType r0) { R0 = r0; }
  inline void SetPressureScale(const DataType p0) { P0 = p0; }
  inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
  inline void SetSpeedUp(const DataType s0) { S0 = s0; }
  inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }

  inline void SetTemperatureScale(const DataType tpmin, const DataType tpmax) { 
    Tpmin = tpmin; Tpmax = tpmax; 
    Temp0 = tpmax-tpmin;
    Tref = DataType(0.5)*(Tpmax+Tpmin)/Temp0;
    //Tref = DataType(0.5);
  }

  inline const DataType& DensityScale() const { return R0; }
  inline const DataType& PressureScale() const { return P0; }
  inline const DataType VelocityScale() const { return U0/S0; }
  inline const DataType& TemperatureScale() const { return Temp0; }
  // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
  inline const DataType& SpeedUp() const { return S0; }

  // Lattice quantities for the operator
  inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
  inline DataType LatticeViscosityT(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*ds2); }
  //inline DataType LatticeViscosityT(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
  inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
  inline DataType LatticeSpeedOfSoundT() const { return (std::sqrt(ds2)); }
  //inline DataType LatticeSpeedOfSoundT() const { return (std::sqrt(cs2)); }

  // Physical quantities for the operator
  inline void SetGas(DataType pr, DataType rho, DataType nu, DataType cs) {
    rhop = rho; nup = nu; csp = cs; cs2p = cs*cs; prp = pr;
    SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
    SetPressureScale(prp);
  }
  inline void SetThermalGas(const DataType tmin, const DataType tmax, const DataType diff, 
                            const DataType gp, const DataType betap) {    
    nutp = diff;
    SetTemperatureScale(tmin,tmax);
    //DataType H = width;

    DataType g = gp/(L0()/T0()/T0()); //g in Lattice Coordinates
    DataType beta = betap*Temp0; //beta in Lattice Coordinates
    gbeta = g*beta;

    //Force input is Rayleigh-Number => Calculate gbeta directly in Lattice
    //gbeta = force*LatticeViscosity(Omega(T0()))*LatticeViscosityT(OmegaT(T0()))/(H/L0())/(H/L0())/(H/L0());
  }

  inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
  inline virtual const DataType OmegaT(const DataType dt) const { return (dcs2*cs2p*dt/S0/(nutp*S0+dcs2*cs2p*dt/S0*DataType(0.5))); }
  //inline virtual const DataType OmegaT(const DataType dt) const { return (cs2p*dt/S0/(nutp*S0+cs2p*dt/S0*DataType(0.5))); }

  inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q,
                                              const DataType dt) const {
    DataType M2 = 0.0;
    //x-y
    M2 += (f(7) - feq(7))*(f(7) - feq(7));
    M2 += (f(8) - feq(8))*(f(8) - feq(8));
    M2 += (f(9) - feq(9))*(f(9) - feq(9));
    M2 += (f(10) - feq(10))*(f(10) - feq(10));
    //y-z
    M2 += (f(11) - feq(11))*(f(11) - feq(11));
    M2 += (f(12) - feq(12))*(f(12) - feq(12));
    M2 += (f(13) - feq(13))*(f(13) - feq(13));
    M2 += (f(14) - feq(14))*(f(14) - feq(14));
    //x-z
    M2 += (f(15) - feq(15))*(f(15) - feq(15));
    M2 += (f(16) - feq(16))*(f(16) - feq(16));
    M2 += (f(17) - feq(17))*(f(17) - feq(17));
    M2 += (f(18) - feq(18))*(f(18) - feq(18));
    M2 = std::sqrt(M2);

    DataType tau = DataType(1.0)/Omega(dt);
    DataType om_turb =  (DataType(2.0)/( tau
                                         + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));

    return ( om_turb );
  }
  inline const DataType Omega_LES_dynamic(const MicroType &f, const MicroType &feq, const MacroType& q,
                                          const DataType dt) const {
    DataType M2 = 0.0;
    //x-y
    M2 += (f(7) - feq(7))*(f(7) - feq(7));
    M2 += (f(8) - feq(8))*(f(8) - feq(8));
    M2 += (f(9) - feq(9))*(f(9) - feq(9));
    M2 += (f(10) - feq(10))*(f(10) - feq(10));
    //y-z
    M2 += (f(11) - feq(11))*(f(11) - feq(11));
    M2 += (f(12) - feq(12))*(f(12) - feq(12));
    M2 += (f(13) - feq(13))*(f(13) - feq(13));
    M2 += (f(14) - feq(14))*(f(14) - feq(14));
    //x-z
    M2 += (f(15) - feq(15))*(f(15) - feq(15));
    M2 += (f(16) - feq(16))*(f(16) - feq(16));
    M2 += (f(17) - feq(17))*(f(17) - feq(17));
    M2 += (f(18) - feq(18))*(f(18) - feq(18));
    M2 = std::sqrt(M2);

    DataType tau = DataType(1.0)/Omega(dt);
    DataType vel = sqrt( q(1)*q(1)+q(2)*q(2)+q(3)*q(3) );
    DataType Cs_SmagDYN = DataType(0.5)*L0()*L0()*( (vel*vel/csp*(M2-vel)) / ((M2-vel)*(M2-vel))  );
    DataType om_turb =  (DataType(2.0)/( tau
                                         + sqrt(tau*tau + (DataType(18.)*Cs_SmagDYN*M2)/(R0*q(0)*cs2/S0/S0)) ));

    return ( om_turb );
  }
  inline const int TurbulenceType() const { return turbulence; }
  inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
  inline const DataType& GasDensity() const { return rhop; }
  inline const DataType& GasPressure() const { return prp; }
  inline const DataType& GasViscosity() const { return nup; }
  inline const DataType& GasViscosityT() const { return nutp; }
  inline const DataType& GasSpeedofSound() const { return csp; }
  inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
  { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
  inline const DataType GasViscosityT(const DataType omegat, const DataType cs, const DataType dt) const
  { return ((DataType(1.)/omegat-DataType(0.5))*dcs2*cs*cs*dt/S0/S0); }
  // Quantities for output only
  inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }

protected:
  DataType cs2, cs22, cssq, ds2, dcs2;
  DataType P0, R0, U0, S0, Temp0, prp, rhop, Tpmin, Tpmax, Tref, csp, cs2p, nup, nutp, gbeta, gp, Wp, Rp;
  DataType Cs_Smagorinsky, turbulence;
  int method[1], mdx[26], mdy[26], mdz[26];
};


#endif