amroc/lbm/LBMD3Q19_Test.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2012 Oak Ridge National Laboratory
00004 // Ralf Deiterding, deiterdingr@ornl.gov
00005 
00006 #ifndef LBM_D3Q19_H
00007 #define LBM_D3Q19_H
00008 
00016 #include "LBMBase.h"
00017 #include <cfloat>
00018 
00019 #define LB_FACTOR 1.0e5
00020 
00050 template <class DataType>
00051 class LBMD3Q19 : public LBMBase<Vector<DataType,19>, Vector<DataType,4>, 3> {
00052   typedef LBMBase<Vector<DataType,19>, Vector<DataType,4>, 3> base;
00053 public:
00054   typedef typename base::vec_grid_data_type vec_grid_data_type;
00055   typedef typename base::grid_data_type grid_data_type;
00056   typedef typename base::MicroType MicroType;
00057   typedef typename base::MacroType MacroType;
00058   typedef GridData<MacroType,3> macro_grid_data_type;
00059   typedef typename base::SideName SideName;
00060   typedef typename base::point_type point_type;
00061   typedef Vector<DataType,6> TensorType;
00062 
00063   enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
00064   enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall };
00065   enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall, GFMWallLaw, GFMBounceBack };
00066   enum TurbulenceModel { laminar, LES_Smagorinsky, LES_dynamic };
00067 
00068   LBMD3Q19() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.)  {
00069     cs2  = DataType(1.)/DataType(3.);
00070     cs22 = DataType(2.)*cs2;
00071     cssq = DataType(2.)/DataType(9.);
00072     method[0] = 2; method[1] = 0; method[2] = 0;
00073     turbulence = laminar;
00074     mdx[0]=mdx[3]=mdx[4]=mdx[5]=mdx[6]=mdx[11]=mdx[12]=mdx[13]=mdx[14]=0;
00075     mdx[1]=mdx[7]=mdx[9]=mdx[15]=mdx[17]=1;
00076     mdx[2]=mdx[8]=mdx[10]=mdx[16]=mdx[18]=-1;
00077 
00078     mdy[0]=mdy[1]=mdy[2]=mdy[5]=mdy[6]=mdy[15]=mdy[16]=mdy[17]=mdy[18]=0;
00079     mdy[3]=mdy[7]=mdy[10]=mdy[11]=mdy[13]=1;
00080     mdy[4]=mdy[8]=mdy[9]=mdy[12]=mdy[14]=-1;
00081 
00082     mdz[0]=mdz[1]=mdz[2]=mdz[3]=mdz[4]=mdz[7]=mdz[8]=mdz[9]=mdz[10]=0;
00083     mdz[5]=mdz[11]=mdz[14]=mdz[15]=mdz[18]=1;
00084     mdz[6]=mdz[12]=mdz[13]=mdz[16]=mdz[17]=-1;
00085 
00086   }
00087 
00088   virtual ~LBMD3Q19() {}
00089 
00090   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00091     base::LocCtrl = Ctrl.getSubDevice(prefix+"D3Q19");
00092     RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00093     RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00094     RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00095     RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00096     RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00097     RegisterAt(base::LocCtrl,"Gas_nu",nup);
00098     RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00099     RegisterAt(base::LocCtrl,"Gas_W",Wp);
00100     RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00101     RegisterAt(base::LocCtrl,"SpeedUp",S0);
00102   }
00103 
00104   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00105     base::SetupData(gh,ghosts);
00106     if (rhop>0.) R0 = rhop;
00107     if (csp>0.) {
00108       cs2p = csp*csp;
00109       SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00110     }
00111     std::cout << "D3Q19: Gas_rho=" << rhop << "   Gas_Cs=" << csp
00112               << "   Gas_nu=" << nup << std::endl;
00113     std::cout << "D3Q19: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0
00114               << "   R0=" << R0 << std::endl;
00115     if ( TurbulenceType() ==  laminar ) {
00116       std::cout << "LBM Setup() Laminar" << std::endl;
00117     }
00118     if ( TurbulenceType() ==  LES_Smagorinsky ) {
00119       std::cout << "LBM Setup() LES_Smagorinsky : Smagorinsky Constant = " << SmagorinskyConstant()
00120                 << std::endl;
00121 
00122     }
00123     if ( TurbulenceType() ==  LES_dynamic ) {
00124       std::cout << "LBM Setup() LES_dynamic " << std::endl;
00125     }
00126     WriteInit();
00127   }
00128 
00129   virtual void WriteInit() const {
00130     int me = MY_PROC; 
00131     if (me == VizServer) {
00132       std::ofstream ofs("chem.dat", std::ios::out);
00133       ofs << "RU         " << Rp << std::endl;
00134       ofs << "PA         " << BasePressure() << std::endl;
00135       ofs << "Sp(01)     Vicosity" << std::endl;
00136       ofs << "W(01)      " << Wp << std::endl;
00137       ofs << "Sp(02)     |Velocity|" << std::endl;
00138       ofs << "W(02)      " << 1. << std::endl;
00139       ofs << "Sp(03)     Vorticity-x" << std::endl;
00140       ofs << "W(03)      " << 1. << std::endl;
00141       ofs << "Sp(04)     Vorticity-y" << std::endl;
00142       ofs << "W(04)      " << 1. << std::endl;
00143       ofs << "Sp(05)     Vorticity-z" << std::endl;
00144       ofs << "W(05)      " << 1. << std::endl;
00145       ofs << "Sp(06)     |Vorticity|" << std::endl;
00146       ofs << "W(06)      " << 1. << std::endl;
00147       ofs.close();
00148     }    
00149   }
00150 
00151   inline virtual MacroType MacroVariables(const MicroType &f) const {
00152     MacroType q;
00153     q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8)+f(9)
00154         +f(10)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18);
00155     q(1)=(f(1)-f(2)+f(7)-f(8)+f(9)-f(10)+f(15)-f(16)+f(17)-f(18))/q(0);
00156     q(2)=(f(3)-f(4)+f(7)-f(8)-f(9)+f(10)+f(11)-f(12)+f(13)-f(14))/q(0);
00157     q(3)=(f(5)-f(6)+f(11)-f(12)-f(13)+f(14)+f(15)-f(16)-f(17)+f(18))/q(0);
00158     return q;
00159   }
00160 
00161   inline virtual MicroType Equilibrium(const MacroType &q) const {
00162     MicroType feq;
00163 
00164     DataType rho = q(0);
00165     DataType u   = q(1);
00166     DataType v   = q(2);
00167     DataType w   = q(3);
00168 
00169     DataType ri, rt0, rt1, rt2;
00170     if (method[0]==0) {
00171       ri = DataType(1.);
00172       rt0 = R0 / DataType(3.);
00173       rt1 = R0 / DataType(18.);
00174       rt2 = R0 / DataType(36.);
00175     }
00176     if (method[0]==1) {
00177       ri = rho;
00178       rt0 = DataType(1.) / DataType(3.);
00179       rt1 = DataType(1.) / DataType(18.);
00180       rt2 = DataType(1.) / DataType(36.);
00181     }
00182     else if (method[0]==2) {
00183       ri = DataType(1.);
00184       rt0 = rho / DataType(3.);
00185       rt1 = rho / DataType(18.);
00186       rt2 = rho / DataType(36.);
00187     }
00188 
00189     DataType usq = u * u;
00190     DataType vsq = v * v;
00191     DataType wsq = w * w;
00192     DataType sumsq = (usq + vsq + wsq) / cs22;
00193     DataType sumsqxy = (usq + vsq) / cs22;
00194     DataType sumsqyz = (vsq + wsq) / cs22;
00195     DataType sumsqxz = (usq + wsq) / cs22;
00196     DataType sumsq2xy = sumsqxy * (DataType(1.) - cs2) / cs2;
00197     DataType sumsq2yz = sumsqyz * (DataType(1.) - cs2) / cs2;
00198     DataType sumsq2xz = sumsqxz * (DataType(1.) - cs2) / cs2;
00199     DataType u2 = usq / cssq;
00200     DataType v2 = vsq / cssq;
00201     DataType w2 = wsq / cssq;
00202     DataType u22 = usq / cs22;
00203     DataType v22 = vsq / cs22;
00204     DataType w22 = wsq / cs22;
00205     DataType ui = u / cs2;
00206     DataType vi = v / cs2;
00207     DataType wi = w / cs2;
00208     DataType uv = ui * vi;
00209     DataType uw = ui * wi;
00210     DataType vw = vi * wi;
00211 
00212     feq(0) = rt0 * (ri - sumsq);
00213     //x
00214     feq(1) = rt1 * (ri - sumsq + u2 + ui);
00215     feq(2) = rt1 * (ri - sumsq + u2 - ui);
00216     //y
00217     feq(3) = rt1 * (ri - sumsq + v2 + vi);
00218     feq(4) = rt1 * (ri - sumsq + v2 - vi);
00219     //z
00220     feq(5) = rt1 * (ri - sumsq + w2 + wi);
00221     feq(6) = rt1 * (ri - sumsq + w2 - wi);
00222     //x-y
00223     feq(7) = rt2 * (ri + sumsq2xy -w22 + ui + vi + uv);
00224     feq(8) = rt2 * (ri + sumsq2xy -w22 - ui - vi + uv);
00225     feq(9) = rt2 * (ri + sumsq2xy -w22 + ui - vi - uv);
00226     feq(10) = rt2 * (ri + sumsq2xy -w22 - ui + vi - uv);
00227     //y-z
00228     feq(11) = rt2 * (ri + sumsq2yz -u22 + vi + wi + vw);
00229     feq(12) = rt2 * (ri + sumsq2yz -u22 - vi - wi + vw);
00230     feq(13) = rt2 * (ri + sumsq2yz -u22 + vi - wi - vw);
00231     feq(14) = rt2 * (ri + sumsq2yz -u22 - vi + wi - vw);
00232     //x-z
00233     feq(15) = rt2 * (ri + sumsq2xz -v22 + ui + wi + uw);
00234     feq(16) = rt2 * (ri + sumsq2xz -v22 - ui - wi + uw);
00235     feq(17) = rt2 * (ri + sumsq2xz -v22 + ui - wi - uw);
00236     feq(18) = rt2 * (ri + sumsq2xz -v22 - ui + wi - uw);
00237 
00238     return feq;
00239   }
00240 
00241   inline virtual void Collision(MicroType &f, const DataType dt) const {
00242     MacroType q = MacroVariables(f);
00243     MicroType feq = Equilibrium(q);
00244 
00245     DataType omega;
00246     if (turbulence == laminar) 
00247       omega = Omega(dt);
00248     else if (turbulence == LES_Smagorinsky) 
00249       omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00250     else if (turbulence == LES_dynamic) 
00251       omega = Omega_LES_dynamic(f,feq,q,dt);
00252 
00253     f=(DataType(1.)-omega)*f + omega*feq;
00254   }     
00255 
00256   virtual int IncomingIndices(const int side, int indices[]) const {
00257     switch (side) {
00258     case base::Left:
00259       indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;  
00260       break;
00261     case base::Right:
00262       indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;   
00263       break;
00264     case base::Bottom:
00265       indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;  
00266       break;
00267     case base::Top:
00268       indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;  
00269       break;
00270     case base::Front:
00271       indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;  
00272       break;
00273     case base::Back:
00274       indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;  
00275       break;
00276      }
00277     return 5;
00278   }
00279 
00280   virtual int OutgoingIndices(const int side, int indices[]) const {
00281     switch (side) {
00282     case base::Left:
00283       indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;   
00284       break;
00285     case base::Right:
00286       indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;  
00287       break;
00288     case base::Bottom:
00289       indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;  
00290       break;
00291     case base::Top:
00292       indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;  
00293       break;
00294     case base::Front:
00295       indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;  
00296       break;
00297     case base::Back:
00298       indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;  
00299       break;
00300      }
00301     return 5;
00302   }
00303 
00304   // Revert just one layer of cells at the boundary 
00305   virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
00306                              const int side) const {
00307     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00308     int b = fvec.bottom();
00309     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) 
00310             && fvec.stepsize(2)==bb.stepsize(2));
00311     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00312     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00313     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00314     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00315     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00316     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00317     MicroType *f = (MicroType *)fvec.databuffer();
00318 
00319 #ifdef DEBUG_PRINT_FIXUP
00320     ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
00321       << fvec.bbox() << " using " << bb << "\n").flush();
00322 #endif
00323     
00324     switch (side) {
00325     case base::Left:
00326       for (register int k=mzs; k<=mze; k++) 
00327         for (register int j=mys; j<=mye; j++) {
00328           f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](2);
00329           f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](8);
00330           f[b+base::idx(mxs,j,k,Nx,Ny)](10)= f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](10);
00331           f[b+base::idx(mxs,j,k,Nx,Ny)](16)= f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](16);
00332           f[b+base::idx(mxs,j,k,Nx,Ny)](18)= f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](18);
00333         }
00334       break;
00335     case base::Right:
00336       for (register int k=mzs; k<=mze; k++) 
00337         for (register int j=mys; j<=mye; j++) {
00338           f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](1);
00339           f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](7);
00340           f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](9);
00341           f[b+base::idx(mxe,j,k,Nx,Ny)](15)= f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](15);
00342           f[b+base::idx(mxe,j,k,Nx,Ny)](17)= f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](17);
00343         }
00344       break;
00345     case base::Bottom:
00346       for (register int k=mzs; k<=mze; k++) 
00347         for (register int i=mxs; i<=mxe; i++) {
00348           f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](4);
00349           f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](8);
00350           f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](9);
00351           f[b+base::idx(i,mys,k,Nx,Ny)](12)= f[b+base::idx(i,mys-1,k-1,Nx,Ny)](12);
00352           f[b+base::idx(i,mys,k,Nx,Ny)](14)= f[b+base::idx(i,mys-1,k+1,Nx,Ny)](14);
00353         }
00354       break;
00355     case base::Top:
00356       for (register int k=mzs; k<=mze; k++) 
00357         for (register int i=mxs; i<=mxe; i++) {
00358           f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](3);
00359           f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](7);
00360           f[b+base::idx(i,mye,k,Nx,Ny)](10)= f[b+base::idx(i-1,mye+1,k,Nx,Ny)](10);
00361           f[b+base::idx(i,mye,k,Nx,Ny)](11)= f[b+base::idx(i,mye+1,k+1,Nx,Ny)](11);
00362           f[b+base::idx(i,mye,k,Nx,Ny)](13)= f[b+base::idx(i,mye+1,k-1,Nx,Ny)](13);
00363         }
00364       break;
00365     case base::Front:
00366       for (register int j=mys; j<=mye; j++) 
00367         for (register int i=mxs; i<=mxe; i++) {
00368           f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](6);
00369           f[b+base::idx(i,j,mzs,Nx,Ny)](12)= f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](12);
00370           f[b+base::idx(i,j,mzs,Nx,Ny)](13)= f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](13);
00371           f[b+base::idx(i,j,mzs,Nx,Ny)](16)= f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](16);
00372           f[b+base::idx(i,j,mzs,Nx,Ny)](17)= f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](17);
00373         }
00374       break;
00375     case base::Back:
00376       for (register int j=mys; j<=mye; j++) 
00377         for (register int i=mxs; i<=mxe; i++) {
00378           f[b+base::idx(i,j,mzs,Nx,Ny)](5) = f[b+base::idx(i,j,mzs+1,Nx,Ny)](5);
00379           f[b+base::idx(i,j,mzs,Nx,Ny)](11)= f[b+base::idx(i,j+1,mzs+1,Nx,Ny)](11);
00380           f[b+base::idx(i,j,mzs,Nx,Ny)](14)= f[b+base::idx(i,j-1,mzs+1,Nx,Ny)](14);
00381           f[b+base::idx(i,j,mzs,Nx,Ny)](15)= f[b+base::idx(i+1,j,mzs+1,Nx,Ny)](15);
00382           f[b+base::idx(i,j,mzs,Nx,Ny)](18)= f[b+base::idx(i-1,j,mzs+1,Nx,Ny)](18);
00383         }
00384       break;
00385     }
00386   }   
00387 
00388   virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec, 
00389                          const BBox &bb, const double& dt) const {
00390     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00391     int b = fvec.bottom();
00392     MicroType *f = (MicroType *)fvec.databuffer();
00393     MicroType *of = (MicroType *)ovec.databuffer();
00394 
00395 #ifdef DEBUG_PRINT_FIXUP
00396     ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox() 
00397       << " using " << bb << "\n").flush();
00398 #endif
00399     
00400     assert (fvec.extents(0)==ovec.extents(0));
00401     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) && fvec.stepsize(2)==bb.stepsize(2) &&
00402             fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1) && fvec.stepsize(2)==ovec.stepsize(2));
00403     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00404     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00405     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00406     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00407     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00408     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00409 
00410     for (register int k=mzs; k<=mze; k++) 
00411       for (register int j=mys; j<=mye; j++) 
00412         for (register int i=mxs; i<=mxe; i++) {
00413           for (register int n=0; n<base::NMicroVar(); n++) 
00414             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);
00415           Collision(f[b+base::idx(i,j,k,Nx,Ny)],dt);
00416         }
00417   }  
00418 
00419   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00420                       const double& t, const double& dt, const int& mpass) const {
00421     // Make sure that the physical length and times are scaling consistently
00422     // into lattice quantities
00423     DCoords dx = base::GH().worldStep(fvec.stepsize());
00424     
00425     //std::cout << "STEP   dt=" << dt << "   TO=" << T0() << std::endl;
00426     if ( (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR)
00427          || (std::fabs(dx(1)/L0()-dx(2)/L0()) > DBL_EPSILON*LB_FACTOR) ) {
00428       std::cerr << "LBMD3Q19 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00429                 << " and dx(2)=" << dx(2)/L0()
00430                 << " to be equal." << std::endl << std::flush;
00431       std::exit(-1);
00432     }
00433     if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
00434       std::cerr << "LBMD3Q19 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00435                 << " to be equal."
00436                 << "   dt=" << dt << "   TO=" << T0()
00437                 << std::endl << std::flush;
00438       std::exit(-1);
00439     }
00440 
00441     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00442     MicroType *f = (MicroType *)fvec.databuffer();
00443 
00444     // Streaming - update all cells
00445     for (register int k=Nz-1; k>=1; k--)
00446       for (register int j=Ny-1; j>=1; j--)
00447         for (register int i=0; i<=Nx-2; i++) {
00448           //x-y
00449           f[base::idx(i,j,k,Nx,Ny)](3) = f[base::idx(i,j-1,k,Nx,Ny)](3);
00450           f[base::idx(i,j,k,Nx,Ny)](10) = f[base::idx(i+1,j-1,k,Nx,Ny)](10);
00451           //y-z
00452           //f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);//%%%%% -j,+k
00453           //x-z
00454           f[base::idx(i,j,k,Nx,Ny)](5) = f[base::idx(i,j,k-1,Nx,Ny)](5);
00455           f[base::idx(i,j,k,Nx,Ny)](18) = f[base::idx(i+1,j,k-1,Nx,Ny)](18);
00456         }
00457 
00458     for (register int k=0; k<=Nz-2; k++)
00459       for (register int j=Ny-1; j>=1; j--)
00460         for (register int i=0; i<=Nx-2; i++) {
00461           //y-z
00462           f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);//%%%%% -j,+k
00463         }
00464     for (register int k=Nz-1; k>=1; k--)
00465       for (register int j=Ny-1; j>=1; j--)
00466         for (register int i=Nx-1; i>=1; i--) {
00467           //x-y
00468           f[base::idx(i,j,k,Nx,Ny)](1) = f[base::idx(i-1,j,k,Nx,Ny)](1);
00469           f[base::idx(i,j,k,Nx,Ny)](7) = f[base::idx(i-1,j-1,k,Nx,Ny)](7);
00470           //y-z
00471           f[base::idx(i,j,k,Nx,Ny)](11) = f[base::idx(i,j-1,k-1,Nx,Ny)](11);
00472           //x-z
00473           f[base::idx(i,j,k,Nx,Ny)](15) = f[base::idx(i-1,j,k-1,Nx,Ny)](15);
00474         }
00475 
00476     for (register int k=0; k<=Nz-2; k++)
00477       for (register int j=0; j<=Ny-2; j++)
00478         for (register int i=Nx-1; i>=1; i--) {
00479           //x-y
00480           f[base::idx(i,j,k,Nx,Ny)](4) = f[base::idx(i,j+1,k,Nx,Ny)](4);
00481           f[base::idx(i,j,k,Nx,Ny)](9) = f[base::idx(i-1,j+1,k,Nx,Ny)](9);
00482           //y-z
00483           //f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);//%%%%% +j,-k
00484           //x-z
00485           f[base::idx(i,j,k,Nx,Ny)](6) = f[base::idx(i,j,k+1,Nx,Ny)](6);
00486           f[base::idx(i,j,k,Nx,Ny)](17) = f[base::idx(i-1,j,k+1,Nx,Ny)](17);
00487         }
00488     for (register int k=Nz-1; k>=1; k--)
00489       for (register int j=0; j<=Ny-2; j++)
00490         for (register int i=0; i<=Nx-2; i++) {
00491           //y-z
00492           f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);//%%%%% +j,-k
00493         }
00494     for (register int k=0; k<=Nz-2; k++)
00495       for (register int j=0; j<=Ny-2; j++)
00496         for (register int i=0; i<=Nx-2; i++) {
00497           //x-y
00498           f[base::idx(i,j,k,Nx,Ny)](2) = f[base::idx(i+1,j,k,Nx,Ny)](2);
00499           f[base::idx(i,j,k,Nx,Ny)](8) = f[base::idx(i+1,j+1,k,Nx,Ny)](8);
00500           //y-z
00501           f[base::idx(i,j,k,Nx,Ny)](12) = f[base::idx(i,j+1,k+1,Nx,Ny)](12);
00502           //x-z
00503           f[base::idx(i,j,k,Nx,Ny)](16) = f[base::idx(i+1,j,k+1,Nx,Ny)](16);
00504         }
00505 
00506      // Collision - only in the interior region
00507     int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00508     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00509     for (register int k=mzs; k<=mze; k++)
00510       for (register int j=mys; j<=mye; j++)
00511         for (register int i=mxs; i<=mxe; i++)
00512           Collision(f[base::idx(i,j,k,Nx,Ny)],dt);
00513 
00514     return DataType(1.);
00515   }
00516 
00517   virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00518                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00519     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00520     int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00521     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00522     MicroType *f = (MicroType *)fvec.databuffer();
00523 
00524     //std::cout << "IC Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << "   type=" << type << std::endl;
00525 
00526     if (type==ConstantMicro) {
00527       //std::cout << "ConstantMicro start\n";
00528       MicroType g;
00529       for (register int i=0; i<base::NMicroVar(); i++)
00530         g(i) = aux[i];
00531 
00532       for (register int k=mzs; k<=mze; k++)
00533         for (register int j=mys; j<=mye; j++)
00534           for (register int i=mxs; i<=mxe; i++)
00535             f[base::idx(i,j,k,Nx,Ny)] = g;
00536       //std::cout << "ConstantMicro end\n";
00537     }
00538     
00539     else if (type==ConstantMacro) {
00540       //std::cout << "ConstantMacro start\n";
00541       MacroType q;
00542       if (scaling==base::Physical) {
00543         q(0) = aux[0]/R0;
00544         q(1) = aux[1]/U0;
00545         q(2) = aux[2]/U0;
00546         q(3) = aux[3]/U0;
00547       }
00548       else
00549         for (register int i=0; i<base::NMacroVar(); i++)
00550           q(i) = aux[i];
00551       q(1) *= S0; q(2) *= S0; q(3) *= S0;
00552       
00553 
00554       /*std::cout << "aux[0]=" << aux[0] << "  aux[1]=" << aux[1]
00555   << "  aux[2]=" << aux[2] << "  aux[3]=" << aux[3]
00556   << "  equlibrium=" << Equilibrium(q)
00557   << "  q=" <<  MacroVariables(Equilibrium(q)) << std::endl;*/
00558 
00559       for (register int k=mzs; k<=mze; k++)
00560         for (register int j=mys; j<=mye; j++)
00561           for (register int i=mxs; i<=mxe; i++)
00562             f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00563       //std::cout << "ConstantMacro end\n";
00564     }
00565     
00566     else if (type==GasAtRest) {
00567       //std::cout << "GasAtRest start\n";
00568       MacroType q;
00569       q(0) = GasDensity()/R0;
00570       q(1) = DataType(0.);
00571       q(2) = DataType(0.);
00572       q(3) = DataType(0.);
00573       for (register int k=mzs; k<=mze; k++)
00574         for (register int j=mys; j<=mye; j++)
00575           for (register int i=mxs; i<=mxe; i++)
00576             f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00577       //std::cout << "GasAtRest end\n";
00578     }
00579   }
00580 
00581   // Set just the one layer of ghost cells at the boundary
00582   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00583                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00584     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00585     int b = fvec.bottom();
00586     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) 
00587             && fvec.stepsize(2)==bb.stepsize(2));
00588     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00589     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00590     int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00591     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00592     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00593     int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00594     MicroType *f = (MicroType *)fvec.databuffer();
00595     
00596     //std::cout << "BC Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << std::endl;
00597     if (type==Symmetry) {
00598       switch (side) {
00599       case base::Left:
00600         for (register int k=mzs; k<=mze; k++)
00601           for (register int j=mys; j<=mye; j++) {
00602             if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
00603             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00604             if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
00605 
00606             if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
00607             if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
00608           }
00609         break;
00610       case base::Right:
00611         for (register int k=mzs; k<=mze; k++)
00612           for (register int j=mys; j<=mye; j++) {
00613             if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
00614             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00615             if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
00616 
00617             if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
00618             if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
00619           }
00620         break;
00621       case base::Bottom:
00622         for (register int k=mzs; k<=mze; k++)
00623           for (register int i=mxs; i<=mxe; i++) {
00624             if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
00625             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00626             if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
00627 
00628             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);
00629             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);
00630           }
00631         break;
00632       case base::Top:
00633         for (register int k=mzs; k<=mze; k++)
00634           for (register int i=mxs; i<=mxe; i++) {
00635             if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
00636             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00637             if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
00638 
00639             if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
00640             if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
00641           }
00642         break;
00643       case base::Front:
00644         //std::cout << "SYMMETRY Front Nz=" << Nz << std::endl;
00645         for (register int j=mys; j<=mye; j++)
00646           for (register int i=mxs; i<=mxe; i++) {
00647             if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
00648             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00649             if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
00650 
00651             if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
00652             if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
00653           }
00654         break;
00655       case base::Back:
00656         //std::cout << "SYMMETRY Back Nz=" << Nz << std::endl;
00657         for (register int j=mys; j<=mye; j++)
00658           for (register int i=mxs; i<=mxe; i++) {
00659             /*std::cout << "NS Back Nz=" << Nz << "  {"<< i << "," << j << "," << mzs << "}"
00660         << " mxs=" << mxs << " mxe=" << mxe
00661         << " mys=" << mys << " mye=" << mye
00662         << " mzs=" << mzs << " mze=" << mze << std::endl;*/
00663             if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
00664             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00665             if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
00666 
00667             if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
00668             if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
00669           }
00670         break;
00671       }
00672     }
00673     
00674     else if (type==SlipWall) {
00675       switch (side) {
00676       case base::Left:
00677         //std::cout << "S Left Nz=" << Nz << std::endl;
00678         for (register int k=mzs; k<=mze; k++)
00679           for (register int j=mys; j<=mye; j++) {
00680             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);
00681             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00682             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);
00683 
00684             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);
00685             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);
00686           }
00687         break;
00688       case base::Right:
00689         //std::cout << "S Right Nz=" << Nz << std::endl;
00690         for (register int k=mzs; k<=mze; k++)
00691           for (register int j=mys; j<=mye; j++) {
00692             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);
00693             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00694             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);
00695 
00696             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);
00697             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);
00698           }
00699         break;
00700       case base::Bottom:
00701         //std::cout << "S Bottom Nz=" << Nz << std::endl;
00702         for (register int k=mzs; k<=mze; k++)
00703           for (register int i=mxs; i<=mxe; i++) {
00704             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);
00705             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00706             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);
00707 
00708             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);
00709             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);
00710           }
00711         break;
00712       case base::Top:
00713         //std::cout << "S Top Nz=" << Nz << std::endl;
00714         for (register int k=mzs; k<=mze; k++)
00715           for (register int i=mxs; i<=mxe; i++) {
00716             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);
00717             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00718             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);
00719 
00720             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);
00721             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);
00722           }
00723         break;
00724       case base::Front:
00725         //std::cout << "S Front Nz=" << Nz << std::endl;
00726         for (register int j=mys; j<=mye; j++)
00727           for (register int i=mxs; i<=mxe; i++) {
00728             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);
00729             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00730             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);
00731 
00732             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);
00733             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);
00734           }
00735         break;
00736       case base::Back:
00737         //std::cout << "S Back Nz=" << Nz << std::endl;
00738         for (register int j=mys; j<=mye; j++)
00739           for (register int i=mxs; i<=mxe; i++) {
00740             /*std::cout << "NS Back Nz=" << Nz << "  {"<< i << "," << j << "," << mzs << "}"
00741         << " mxs=" << mxs << " mxe=" << mxe
00742         << " mys=" << mys << " mye=" << mye
00743         << " mzs=" << mzs << " mze=" << mze << std::endl;*/
00744             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);
00745             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00746             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);
00747 
00748             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);
00749             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);
00750           }
00751         break;
00752       }
00753     }
00754 
00755     else if (type==NoSlipWall) {
00756       switch (side) {
00757       case base::Left:
00758         //std::cout << "NS Left Nz=" << Nz << std::endl;
00759         for (register int k=mzs; k<=mze; k++)
00760           for (register int j=mys; j<=mye; j++) {
00761             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);
00762             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00763             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);
00764 
00765             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);
00766             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);
00767           }
00768         break;
00769       case base::Right:
00770         //std::cout << "NS Right Nz=" << Nz << std::endl;
00771         for (register int k=mzs; k<=mze; k++)
00772           for (register int j=mys; j<=mye; j++) {
00773             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);
00774             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00775             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);
00776 
00777             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);
00778             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);
00779           }
00780         break;
00781       case base::Bottom:
00782         //std::cout << "NS Bottom Nz=" << Nz << std::endl;
00783         for (register int k=mzs; k<=mze; k++)
00784           for (register int i=mxs; i<=mxe; i++) {
00785             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);
00786             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00787             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);
00788 
00789             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);
00790             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);
00791           }
00792         break;
00793       case base::Top:
00794         //std::cout << "NS Top Nz=" << Nz << std::endl;
00795         for (register int k=mzs; k<=mze; k++)
00796           for (register int i=mxs; i<=mxe; i++) {
00797             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);
00798             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00799             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);
00800 
00801             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);
00802             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);
00803           }
00804         break;
00805       case base::Front:
00806         //std::cout << "NS Front Nz=" << Nz << std::endl;
00807         for (register int j=mys; j<=mye; j++)
00808           for (register int i=mxs; i<=mxe; i++) {
00809             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);
00810             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00811             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);
00812 
00813             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);
00814             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);
00815           }
00816         break;
00817       case base::Back:
00818         //std::cout << "NS Back Nz=" << Nz << std::endl;
00819         for (register int j=mys; j<=mye; j++)
00820           for (register int i=mxs; i<=mxe; i++) {
00821             /*std::cout << "NS Back Nz=" << Nz << "  {"<< i << "," << j << "," << mzs << "}"
00822         << " mxs=" << mxs << " mxe=" << mxe
00823         << " mys=" << mys << " mye=" << mye
00824         << " mzs=" << mzs << " mze=" << mze << std::endl;*/
00825             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);
00826             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00827             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);
00828 
00829             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);
00830             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);
00831           }
00832         break;
00833       }
00834     }
00835 
00836     else if (type==Inlet) {
00837       if (aux==0 || naux<=0) return;
00838       DataType a0=S0*aux[0], a1=0., a2=0.;
00839       if (naux>1) {
00840         a1 = S0*aux[1];
00841         a2 = S0*aux[2];
00842       }
00843       if (scaling==base::Physical) {
00844         a0 /= U0;
00845         a1 /= U0;
00846         a2 /= U0;
00847       }
00848       switch (side) {
00849       case base::Left:
00850         for (register int k=mzs; k<=mze; k++)
00851           for (register int j=mys; j<=mye; j++) {
00852             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00853             q(1) = a0;
00854             if (naux>1) {
00855               q(2) = a1;
00856               q(3) = a2;
00857             }
00858             f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
00859           }
00860         break;
00861       case base::Right:
00862         for (register int k=mzs; k<=mze; k++)
00863           for (register int j=mys; j<=mye; j++) {
00864             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00865             q(1) = a0;
00866             if (naux>1) {
00867               q(2) = a1;
00868               q(3) = a2;
00869             }
00870             f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
00871           }
00872         break;
00873       case base::Bottom:
00874         for (register int k=mzs; k<=mze; k++)
00875           for (register int i=mxs; i<=mxe; i++) {
00876             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
00877             if (naux==1)
00878               q(2) = a0;
00879             else {
00880               q(1) = a0;
00881               q(2) = a1;
00882               q(3) = a2;
00883             }
00884             f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
00885           }
00886         break;
00887       case base::Top:
00888         for (register int k=mzs; k<=mze; k++)
00889           for (register int i=mxs; i<=mxe; i++) {
00890             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
00891             if (naux==1)
00892               q(2) = a0;
00893             else {
00894               q(1) = a0;
00895               q(2) = a1;
00896               q(3) = a2;
00897             }
00898             f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
00899           }
00900         break;
00901       case base::Front:
00902         for (register int j=mys; j<=mye; j++)
00903           for (register int i=mxs; i<=mxe; i++) {
00904             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
00905             if (naux==1)
00906               q(3) = a0;
00907             else {
00908               q(1) = a1;
00909               q(2) = a2;
00910               q(3) = a0;
00911             }
00912             f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
00913           }
00914         break;
00915       case base::Back:
00916         for (register int j=mys; j<=mye; j++)
00917           for (register int i=mxs; i<=mxe; i++) {
00918             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
00919             if (naux==1)
00920               q(3) = a0;
00921             else {
00922               q(1) = a1;
00923               q(2) = a2;
00924               q(3) = a0;
00925             }
00926             f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
00927           }
00928         break;
00929       }
00930     }
00931 
00932     else if (type==Outlet) {
00933       //std::cout << "Outlet" << std::endl;
00934       switch (side) {
00935       case base::Left:
00936         //std::cout << "Outlet Left" << std::endl;
00937         for (register int k=mzs; k<=mze; k++)
00938           for (register int j=mys; j<=mye; j++) {
00939             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00940             //std::cout << "Outlet Left q(0)=" << q(0) << std::endl;
00941             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);
00942           }
00943         break;
00944       case base::Right:
00945         //std::cout << "Outlet Right" << std::endl;
00946         for (register int k=mzs; k<=mze; k++)
00947           for (register int j=mys; j<=mye; j++) {
00948             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00949             //std::cout << "Outlet Right q(0)=" << q(0) << std::endl;
00950             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);
00951           }
00952         break;
00953       case base::Bottom:
00954         for (register int k=mzs; k<=mze; k++)
00955           for (register int i=mxs; i<=mxe; i++) {
00956             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
00957             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);
00958           }
00959         break;
00960       case base::Top:
00961         for (register int k=mzs; k<=mze; k++)
00962           for (register int i=mxs; i<=mxe; i++) {
00963             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
00964             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);
00965           }
00966         break;
00967       case base::Front:
00968         for (register int j=mys; j<=mye; j++)
00969           for (register int i=mxs; i<=mxe; i++) {
00970             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
00971             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);
00972           }
00973         break;
00974       case base::Back:
00975         for (register int j=mys; j<=mye; j++)
00976           for (register int i=mxs; i<=mxe; i++) {
00977             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
00978             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);
00979           }
00980         break;
00981       }
00982     }
00983     
00984     // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)
00985     else if (type==Pressure) {
00986       if (aux==0 || naux<=0) return;
00987       DataType a0=aux[0];
00988       if (scaling==base::Physical)
00989         a0 /= R0;
00990       switch (side) {
00991       case base::Left:
00992         for (register int k=mzs; k<=mze; k++)
00993           for (register int j=mys; j<=mye; j++) {
00994             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00995             q(0) = a0;
00996             f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
00997           }
00998         break;
00999       case base::Right:
01000         for (register int k=mzs; k<=mze; k++)
01001           for (register int j=mys; j<=mye; j++) {
01002             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01003             q(0) = a0;
01004             f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
01005           }
01006         break;
01007       case base::Bottom:
01008         for (register int k=mzs; k<=mze; k++)
01009           for (register int i=mxs; i<=mxe; i++) {
01010             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01011             q(0) = a0;
01012             f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
01013           }
01014         break;
01015       case base::Top:
01016         for (register int k=mzs; k<=mze; k++)
01017           for (register int i=mxs; i<=mxe; i++) {
01018             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01019             q(0) = a0;
01020             f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
01021           }
01022         break;
01023       case base::Front:
01024         for (register int j=mys; j<=mye; j++)
01025           for (register int i=mxs; i<=mxe; i++) {
01026             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01027             q(0) = a0;
01028             f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
01029           }
01030         break;
01031       case base::Back:
01032         for (register int j=mys; j<=mye; j++)
01033           for (register int i=mxs; i<=mxe; i++) {
01034             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01035             q(0) = a0;
01036             f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
01037           }
01038         break;
01039       }
01040     }
01041     
01042     else if (type==SlidingWall) {
01043       if (aux==0 || naux<=0) return;
01044       if (aux==0 || naux<=0) return;
01045       DataType a0=S0*aux[0];
01046       if (scaling==base::Physical)
01047         a0 /= U0;
01048       switch (side) {
01049       case base::Left:
01050         for (register int k=mzs; k<=mze; k++)
01051           for (register int j=mys; j<=mye; j++) {
01052             MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01053             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);
01054             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01055             DataType pw = DataType(1.)-qw;
01056             if (j<mye) f[b+base::idx(mxe,j+1,k,Nx,Ny)](9) = pw*rd1+qw*rd2;
01057             f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01058             if (j>mys) f[b+base::idx(mxe,j-1,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01059 
01060             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);
01061             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);
01062           }
01063         break;
01064       case base::Right:
01065         for (register int k=mzs; k<=mze; k++)
01066           for (register int j=mys; j<=mye; j++) {
01067             MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01068             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);
01069             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01070             DataType pw = DataType(1.)-qw;
01071             if (j<mye) f[b+base::idx(mxs,j+1,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01072             f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01073             if (j>mys) f[b+base::idx(mxs,j-1,k,Nx,Ny)](10) = qw*rd1+pw*rd2;
01074 
01075             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);
01076             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);
01077           }
01078         break;
01079       case base::Bottom:
01080         for (register int k=mzs; k<=mze; k++)
01081           for (register int i=mxs; i<=mxe; i++) {
01082             MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01083             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);
01084             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01085             DataType pw = DataType(1.)-qw;
01086             if (i<mxe) f[b+base::idx(i+1,mye,k,Nx,Ny)](10) = pw*rd1+qw*rd2;
01087             f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01088             if (i>mxs) f[b+base::idx(i-1,mye,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01089 
01090             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);
01091             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);
01092           }
01093         break;
01094       case base::Top:
01095         for (register int k=mzs; k<=mze; k++)
01096           for (register int i=mxs; i<=mxe; i++) {
01097             MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01098             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);
01099             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01100             DataType pw = DataType(1.)-qw;
01101             if (i<mxe) f[b+base::idx(i+1,mys,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01102             f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01103             if (i>mxs) f[b+base::idx(i-1,mys,k,Nx,Ny)](9) = qw*rd1+pw*rd2;
01104 
01105             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);
01106             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);
01107           }
01108         break;
01109       case base::Front:
01110         for (register int j=mys; j<=mye; j++)
01111           for (register int i=mxs; i<=mxe; i++) {
01112             MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01113             DataType rd1 = f[b+base::idx(i,j,mze+1,Nx,Ny)](17), rd2 = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
01114             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01115             DataType pw = DataType(1.)-qw;
01116             if (i<mxe) f[b+base::idx(i+1,j,mze,Nx,Ny)](18) = pw*rd1+qw*rd2;
01117             f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01118             if (i>mxs) f[b+base::idx(i-1,j,mze,Nx,Ny)](15) = qw*rd1+pw*rd2;
01119 
01120             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);
01121             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);
01122           }
01123         break;
01124       case base::Back:
01125         for (register int j=mys; j<=mye; j++)
01126           for (register int i=mxs; i<=mxe; i++) {
01127             MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01128             DataType rd1 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15), rd2 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
01129             DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01130             DataType pw = DataType(1.)-qw;
01131             if (i<mxe) f[b+base::idx(i+1,j,mzs,Nx,Ny)](16) = pw*rd1+qw*rd2;
01132             f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
01133             if (i>mxs) f[b+base::idx(i-1,j,mzs,Nx,Ny)](17) = qw*rd1+pw*rd2;
01134 
01135             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);
01136             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);
01137           }
01138         break;
01139       }
01140     }
01141   }
01142 
01143   virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
01144                              const MicroType* f, const point_type* xc, const DataType* distance,
01145                              const point_type* normal, DataType* aux=0, const int naux=0,
01146                              const int scaling=0) const {
01147     if (type==GFMNoSlipWall || type==GFMSlipWall) {
01148       DataType fac = S0;
01149       if (scaling==base::Physical)
01150         fac /= U0;
01151       for (int n=0; n<nc; n++) {
01152         MacroType q=MacroVariables(f[n]);
01153         DataType u=-q(1);
01154         DataType v=-q(2);
01155         DataType w=-q(3);
01156 
01157         // Add boundary velocities if available
01158         if (naux>=3) {
01159           u += fac*aux[naux*n];
01160           v += fac*aux[naux*n+1];
01161           w += fac*aux[naux*n+2];
01162         }
01163         if (type==GFMNoSlipWall) {
01164           // Prescribe entire velocity vector in ghost cells
01165           q(1) += DataType(2.)*u;
01166           q(2) += DataType(2.)*v;
01167           q(3) += DataType(2.)*w;
01168         }
01169         else if (type==GFMSlipWall) {
01170           // Prescribe only normal velocity vector in ghost cells
01171           DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1)+w*normal[n](2));
01172           q(1) += vl*normal[n](0);
01173           q(2) += vl*normal[n](1);
01174           q(3) += vl*normal[n](2);
01175         }
01176         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01177       }
01178     }
01179 
01180     else if (type==GFMExtrapolation) {
01181       for (int n=0; n<nc; n++) {
01182         MacroType q=MacroVariables(f[n]);
01183         DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1)+q(3)*normal[n](2);
01184         q(1) = vl*normal[n](0);
01185         q(2) = vl*normal[n](1);
01186         q(3) = vl*normal[n](2);
01187         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01188       }
01189     }
01190 
01191     else if (type==GFMWallLaw ) {
01192       DCoords dx = base::GH().worldStep(fvec.stepsize());
01193       DataType fac = S0;
01194       if (scaling==base::Physical)
01195         fac /= U0;
01196 
01197       // DataType BLH = 10.;
01198       DataType yp_ = 0.,  yp = 0., yMax = 50.;
01199       DataType yPlus[3];
01200       int lc_low[3], lc_up[3], yPlus_lc[3];
01201 
01202       // int mxe = fvec.upper(0)/fvec.stepsize(0);
01203       // int mye = fvec.upper(1)/fvec.stepsize(1);
01204       // int mze = fvec.upper(2)/fvec.stepsize(2);
01205 
01206       DataType shear_vel = 1., karman = 0.41, y0 = 1., law = 1.;
01207 
01208       DataType shearStress;
01209       DataType u, v, w, u_, v_, w_, vel0;
01210       // DataType vel1;
01211       MacroType q, q_;
01212       int kMax = 0;
01213 
01214       for (int n=0; n<nc; n++) {
01215 
01217         q=MacroVariables(f[n]);
01218         q_=MacroVariables(f[n]);
01219 
01220         yPlus[0] = (xc[n](0) - DataType(2.)*dx(0)*normal[n](0));
01221         yPlus[1] = (xc[n](1) - DataType(2.)*dx(1)*normal[n](1));
01222         yPlus[2] = (xc[n](2) - DataType(2.)*dx(2)*normal[n](2));
01223 
01224         yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
01225         yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
01226         yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
01227 
01228         lc_low[0] = fvec.lower(0);   lc_low[1] = fvec.lower(1);   lc_low[2] = fvec.lower(2);
01229         lc_up[0] = fvec.upper(0);   lc_up[1] = fvec.upper(1);   lc_up[2] = fvec.upper(2);
01230 
01231         //std::cout << "xc " << xc[n] << " yPlus " << yPlus[0] << " " << yPlus[1] << " " << yPlus[2] << std::endl;
01232 
01233         if ( yPlus_lc[0] >  lc_low[0] ) {
01234           if ( yPlus_lc[1] >  lc_low[1] ) {
01235             if ( yPlus_lc[2] >  lc_low[2] ) {
01236 
01237               if ( yPlus_lc[0] < lc_up[0] ) {
01238                 if ( yPlus_lc[1] <  lc_up[1] ) {
01239                   if ( yPlus_lc[2] <  lc_up[2] ) {
01240                     q=MacroVariables( fvec( base::GH().localCoords((const DataType *) yPlus) ) );
01241                   }
01242                 }
01243               }
01244 
01245             }
01246           }
01247         }
01248 
01249         u_=-q(1);
01250         v_=-q(2);
01251         w_=-q(3);
01252         // vel1 = sqrt(u_*u_+v_*v_+w_*w_);
01253 
01254         q=MacroVariables(f[n]);
01255         u=-q(1);
01256         v=-q(2);
01257         w=-q(3);
01258         vel0 = sqrt(u*u+v*v+w*w);
01259 
01260         u_+=q(1);
01261         v_+=q(2);
01262         w_+=q(3);
01263         vel0 = ( u_*normal[n](0) + v_*normal[n](1) + w_*normal[n](2) );
01264 
01265         u_ =  u_ - vel0 * normal[n](0);
01266         v_ =  v_ - vel0 * normal[n](1);
01267         w_ =  w_ - vel0 * normal[n](2);
01268         vel0 = sqrt(u_*u_ + v_*v_ + w_*w_)*U0/S0;
01269 
01270         shearStress = (DataType(0.5)*vel0/dx(0) * q(0)*R0*GasViscosity() );
01271 
01272         // Add boundary velocities if available
01273         if (naux>=3) {
01274           u += fac*aux[naux*n];
01275           v += fac*aux[naux*n+1];
01276           w += fac*aux[naux*n+2];
01277         }
01278         q(1) += DataType(2.)*u;
01279         q(2) += DataType(2.)*v;
01280         q(3) += DataType(2.)*w;
01281         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01282 
01284         shear_vel = std::sqrt( shearStress/(q(0)*R0) );
01285 
01287         y0 = karman/30.;
01288 
01289         //std::cout << n << " vel0 " << vel0 << " vel1 " << vel1  << " shear_vel = " << shear_vel << std::endl;
01290 
01291         if ( shear_vel > 1.e-5 ) {
01292           kMax = (int) std::ceil( (yMax*GasViscosity())/(dx(0)*shear_vel) );
01293 
01294           for (int k=1; k < kMax; k++) {
01295 
01296             yPlus[0] = (xc[n](0) - k*dx(0)*normal[n](0));
01297             yPlus[1] = (xc[n](1) - k*dx(1)*normal[n](1));
01298             yPlus[2] = (xc[n](2) - k*dx(2)*normal[n](2));
01299 
01300             yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
01301             yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
01302             yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
01303 
01304             if ( yPlus_lc[0] >  lc_low[0] ) {
01305               if ( yPlus_lc[1] >  lc_low[1] ) {
01306                 if ( yPlus_lc[2] >  lc_low[2] ) {
01307 
01308                   if ( yPlus_lc[0] < lc_up[0] ) {
01309                     if ( yPlus_lc[1] <  lc_up[1] ) {
01310                       if ( yPlus_lc[2] <  lc_up[2] ) {
01311                         yp = k*dx(0);
01312                         yp_ = yp*shear_vel/GasViscosity();
01313 
01316                         if ( yp_ <= 5.0 ) {
01317                           law = yp_;
01318                         }
01320                         else if (yp_ > 5 && yp_ < 30.0 ) {
01321                           //law = (shear_vel/karman*std::log( yp/y0 ) + yp_)*DataType(0.5);
01322                           law = std::min( yp_ , shear_vel/karman*std::log( yp/y0 ) );
01323 
01324                         }
01326                         else if (yp_ >= 30.0 ) {
01327                           law = shear_vel/karman*std::log( yp/y0 );
01328                         }
01329                         //std::cout << "k = " << k << " yp = " << yp << " shear_vel = " << shear_vel << " yp_ = " << yp_ << " law = " << law << std::endl;
01330                         q_=MacroVariables(f[n]);
01331                         u = law*(-q_(1) + fac*aux[naux*n]);
01332                         v = law*(-q_(2) + fac*aux[naux*n+1]);
01333                         w = law*(-q_(3) + fac*aux[naux*n+2]);
01334 
01335                         q_(1) += DataType(2.)*u;
01336                         q_(2) += DataType(2.)*v;
01337                         q_(3) += DataType(2.)*w;
01338                         fvec( base::GH().localCoords((const DataType *) yPlus) ) = Equilibrium(q_);
01339                       }
01340                     }
01341                   }
01342 
01343 
01344                 }
01345               }
01346             }
01347 
01348           }
01349 
01350         }
01351 
01352       }
01353       //std::cout << "====== WallLaw Done\n";
01354     }
01355 
01356 /*************************************************************************************************
01357 *************************************************************************************************/
01358 
01359     else if (type==GFMBounceBack ) {
01360 
01361      DCoords dx = base::GH().worldStep(fvec.stepsize());
01362      DataType fac = S0;
01363      if(scaling==base::Physical) fac/=U0;
01364    
01365       //xc : Cell midpoints of the ghost cells
01366       //nc : Number of ghost cells
01367       //normal : normalized normal vectors, pointed inside
01368       //distance: Level-Set distance between ghost cell midpoint and boundary
01369        
01370       DataType distB[3], distBScalar, cosa, leng, d_dist, qNorm, mdx2, mdy2, mdz2, zaehler, nenner, lengXi;
01371         
01372       for (int n=0; n<nc; n++){
01373         distB[0] = (double &)distance[n]*normal[n](0);
01374         distB[1] = (double &)distance[n]*normal[n](1);
01375         distB[2] = (double &)distance[n]*normal[n](2);
01376 
01377         //alpha: Angle between discrete velocity direction (Ghostcell) and inner normal direction (Boundary)
01378         //leng: Distance between ghost cell midpoint and curved boundary along discrete velocity direction
01379         //d_dist: Distance between fluid-cellmidpoint and boundary along discrete velocity direction
01380         //qNorm: Normalized length from fluid-cellmidpoint to boundary along discrete velocity direction
01381  
01382         for(int ind=1; ind<base::NMicroVar(); ind++){
01383                 
01384                 mdx2 = mdx[ind]*mdx[ind];
01385                 mdy2 = mdy[ind]*mdy[ind];
01386                 mdz2 = mdz[ind]*mdz[ind];
01387                 distBScalar = std::sqrt(std::pow(distB[0],2)+std::pow(distB[1],2)+std::pow(distB[2],2));
01388                 leng = distBScalar*distBScalar*std::sqrt(mdx2+mdy2+mdz2)/(distB[0]*mdx[ind]+distB[1]*mdy[ind]+distB[2]*mdz[ind]);
01389                 lengXi = std::sqrt(mdx2*dx(0)*dx(0)+mdy2*dx(1)*dx(1)+mdz2*dx(2)*dx(2));
01390                 d_dist = lengXi-leng;
01391                 qNorm = d_dist/lengXi;
01392 
01393                 if(qNorm > 0 && qNorm < 1.){ 
01394                 //The positive lengths are important only       
01395                         
01396                         //Two cases: qNorm greater equal or lower than 0.5
01397                         int idxA, idxB;
01398                         if(ind==0.)     idxA=0.;
01399                         else if(ind % 2 == 0.) idxA = ind-1.;
01400                         else idxA = ind+1.;
01401                         
01402                         idxB=idxA;
01403                         if(qNorm>=DataType(0.5)) idxB=ind;
01404 
01405                         //Coords of ghostcell           
01406                         DataType cn[3], cnA[3], cnB[3];
01407                         cn[0] = xc[n](0);
01408                         cn[1] = xc[n](1);
01409                         cn[2] = xc[n](2);
01410 
01411                         //Coords of fluidcell
01412                         cnA[0] = cn[0] + mdx[ind]*dx(0);
01413                         cnA[1] = cn[1] + mdy[ind]*dx(1);
01414                         cnA[2] = cn[2] + mdz[ind]*dx(2);
01415 
01416                         //Coords of 2nd fluidcell
01417                         cnB[0] = cnA[0]; 
01418                         cnB[1] = cnA[1]; 
01419                         cnB[2] = cnA[2]; 
01420 
01421                         if(qNorm<DataType(0.5)){
01422                                 cnB[0] = cn[0] + DataType(2.)*mdx[ind]*dx(0);
01423                                 cnB[1] = cn[1] + DataType(2.)*mdy[ind]*dx(1);
01424                                 cnB[2] = cn[2] + DataType(2.)*mdz[ind]*dx(2);   
01425                         }
01426 /*
01427                         if(ind==0){
01428                                 idxA = ind; idxB = ind;
01429                                 cnA[0] = cn[0]; cnA[1] = cn[1]; cnA[2] = cn[2];
01430                                 cnB[0] = cn[0]; cnB[1] = cn[1]; cnB[2] = cn[2];
01431                         }
01432 */
01433                         //Getting the corresponding partial probability distribution functions
01434 
01435                         DataType distFuncA, distFuncB;
01436 
01437                         //distFuncA = fvec( base::GH().localCoords((const DataType *) cnA) )(idxA);
01438                         //distFuncB = fvec( base::GH().localCoords((const DataType *) cnB) )(idxB);
01439 
01440                         int cnA_lc[3], cnB_lc[3];
01441                         cnA_lc[0] = base::GH().localCoords((const DataType *) cnA)(0);
01442                         cnA_lc[1] = base::GH().localCoords((const DataType *) cnA)(1);
01443                         cnA_lc[2] = base::GH().localCoords((const DataType *) cnA)(2);
01444                         cnB_lc[0] = base::GH().localCoords((const DataType *) cnB)(0);
01445                         cnB_lc[1] = base::GH().localCoords((const DataType *) cnB)(1);
01446                         cnB_lc[2] = base::GH().localCoords((const DataType *) cnB)(2);
01447 
01448                         int lc_low[3], lc_up[3];
01449                         lc_low[0] = fvec.lower(0);   lc_low[1] = fvec.lower(1);   lc_low[2] = fvec.lower(2);
01450                         lc_up[0] = fvec.upper(0);   lc_up[1] = fvec.upper(1);   lc_up[2] = fvec.upper(2);
01451 
01452                         if ( (cnA_lc[0] >  lc_low[0]) && (cnB_lc[0] >  lc_low[0]) ) {
01453                           if ( (cnA_lc[1] >  lc_low[1]) && (cnB_lc[1] >  lc_low[1]) ) {
01454                             if ( (cnA_lc[2] >  lc_low[2]) && (cnB_lc[2] >  lc_low[2]) ) {
01455 
01456                               if ( (cnA_lc[0] < lc_up[0]) && (cnB_lc[0] < lc_up[0]) ) {
01457                                 if ( (cnA_lc[1] <  lc_up[1]) && (cnB_lc[1] < lc_up[1]) ) {
01458                                   if ( (cnA_lc[2] <  lc_up[2]) && (cnB_lc[2] < lc_up[2]) ) {
01459 
01460                                         distFuncA = fvec( base::GH().localCoords((const DataType *) cnA) )(idxA);
01461                                         distFuncB = fvec( base::GH().localCoords((const DataType *) cnB) )(idxB);
01462 
01463                                         DataType intp1,intp2,intp,q2;
01464                                         q2 = DataType(2.)*qNorm;
01465 
01466                                         //Using formular from Eitel-Amor/Bouzidi for density distribution functions
01467                                         if(qNorm<DataType(0.5)){
01468                                                 intp1=q2*distFuncA;
01469                                                 intp2=(DataType(1.)-q2)*distFuncB;
01470                                         }
01471                                         else{
01472                                                 if(q2!=0){
01473                                                         intp1 = (DataType(1.)/q2)*distFuncA;
01474                                                         intp2 = ((q2-DataType(1.))/q2)*distFuncB;
01475                                                 }
01476                                                 else{
01477                                                 intp1 = DataType(0.);
01478                                                         
01479                                                 }
01480                                         }
01481                 
01482                                         intp = intp1+intp2;
01483                                         
01484                                         fvec( base::GH().localCoords((const DataType *) cn) )(ind) = intp;
01485 
01486                                   }
01487                                 }
01488                               }
01489                 
01490                             }
01491                           }
01492                         }
01493 
01494                         //ADDITION PART FOR EXTERNAL FORCE LIKE VELOCITY
01495                         //From: Momentum transfer of a Boltzmann-lattice fluid with boundaries - M. Bouzidi et al. (2001) Physics of Fluids
01496 
01497                         if(naux>=3){
01498                                 DataType u = fac*aux[naux*n];
01499                                 DataType v = fac*aux[naux*n+1];
01500                                 DataType w = fac*aux[naux*n+2];
01501         
01502                                 DataType weight;
01503                                 weight = DataType(1.)/DataType(12.);
01504                                 if(ind==0) weight *= 0;
01505                                 else if(ind>=1 && ind<=6) weight *= DataType(2.);
01506         
01507                                 DataType velForce;
01508                                         
01509                                 if(qNorm<DataType(0.5)) velForce = DataType(2.)*weight;
01510                                 else velForce = (DataType(1.)/qNorm)*weight;
01511                                         
01512                                 velForce *= ((mdx[ind]*u)+(mdy[ind]*v)+(mdz[ind]*w));
01513                                 fvec( base::GH().localCoords((const DataType*) cn) )(ind) +=velForce;
01514                         }
01515 
01516                 }       
01517         }
01518       }
01519 
01520     }//End GFMBounceBack
01521 
01522   }
01523 
01524   virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01525                       const int skip_ghosts=1) const {
01526     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
01527     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
01528     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
01529     MicroType *f = (MicroType *)fvec.databuffer();
01530     DataType *work = (DataType *)workvec.databuffer();
01531     DCoords dx = base::GH().worldStep(fvec.stepsize());
01532     DataType dt = dx(0)*S0/U0;
01533 
01534     if ((cnt>=1 && cnt<=10) || (cnt>=19 && cnt<=24) || (cnt>=28 && cnt<=33)) { 
01535       for (register int k=mzs; k<=mze; k++)
01536         for (register int j=mys; j<=mye; j++)
01537           for (register int i=mxs; i<=mxe; i++) {
01538             MacroType q=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
01539             switch(cnt) {
01540             case 1:
01541               if (method[0]==0) work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0+rhop;
01542               else 
01543                 work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0;
01544               break;
01545             case 2:
01546               work[base::idx(i,j,k,Nx,Ny)]=q(1)*U0/S0;
01547               break;
01548             case 3:
01549               work[base::idx(i,j,k,Nx,Ny)]=q(2)*U0/S0;
01550               break;
01551             case 4:
01552               work[base::idx(i,j,k,Nx,Ny)]=q(3)*U0/S0;
01553               break;
01554             case 6:
01555               work[base::idx(i,j,k,Nx,Ny)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure()); 
01556               break;
01557             case 7:
01558               work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure(); 
01559               break;
01560             case 9: {
01561               MicroType feq = Equilibrium(q);     
01562               work[base::idx(i,j,k,Nx,Ny)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt),
01563                                                         GasSpeedofSound(),dt);
01564               break;
01565             }
01566             case 10:
01567               work[base::idx(i,j,k,Nx,Ny)]=std::sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3))*U0/S0;
01568               break;
01569             case 19:
01570             case 20:
01571             case 21:
01572             case 22:
01573             case 23:
01574             case 24: {
01575               MicroType feq = Equilibrium(q);
01576               DataType omega;
01577               if (turbulence == laminar) 
01578                 omega = Omega(dt);
01579               else if (turbulence == LES_Smagorinsky) 
01580                 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01581               else if (turbulence == LES_dynamic) 
01582                 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01583               work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*DeviatoricStress(f[base::idx(i,j,k,Nx,Ny)],feq,omega)(cnt-19);
01584               break;
01585             }
01586             case 28: 
01587             case 29:
01588             case 30: 
01589             case 31: 
01590             case 32: 
01591             case 33: {
01592               DataType omega;
01593               if (turbulence == laminar) 
01594                 omega = Omega(dt);
01595               else if (turbulence == LES_Smagorinsky) {
01596                 MicroType feq = Equilibrium(q);         
01597                 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01598               }
01599               else if (turbulence == LES_dynamic) {
01600                 MicroType feq = Equilibrium(q);         
01601                 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01602               }
01603               work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*Stress(f[base::idx(i,j,k,Nx,Ny)],q,omega)(cnt-28);
01604               if (cnt==28 || cnt==30 || cnt==33) work[base::idx(i,j,k,Nx,Ny)]-=BasePressure();
01605               break;
01606             }
01607             }
01608           }
01609     }
01610     else if ((cnt>=11 && cnt<=14) || cnt==18 || (cnt>=40 && cnt<=45)) {
01611       macro_grid_data_type qgrid(fvec.bbox());
01612       MacroType *q = (MacroType *)qgrid.databuffer();
01613       for (register int k=mzs; k<=mze; k++)
01614         for (register int j=mys; j<=mye; j++)
01615           for (register int i=mxs; i<=mxe; i++) {
01616             q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
01617             q[base::idx(i,j,k,Nx,Ny)](0)*=R0;
01618             q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
01619             q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
01620             q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
01621           }
01622       
01623       if (skip_ghosts==0) {
01624         switch(cnt) {
01625         case 40:
01626           mxs+=1; mxe-=1; 
01627           break;
01628         case 42:
01629           mys+=1; mye-=1; 
01630           break;
01631         case 45:
01632           mzs+=1; mze-=1; 
01633           break;
01634         case 11:
01635         case 44:
01636           mys+=1; mye-=1; mzs+=1; mze-=1;
01637           break;
01638         case 12:
01639         case 43:
01640           mxs+=1; mxe-=1; mzs+=1; mze-=1;
01641           break;
01642         case 13:
01643         case 41:
01644           mxs+=1; mxe-=1; mys+=1; mye-=1; 
01645           break;
01646         case 14:
01647         case 18:
01648           mxs+=1; mxe-=1; mys+=1; mye-=1; mzs+=1; mze-=1;
01649           break;
01650         }
01651       }
01652 
01653       DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz; 
01654       for (register int k=mzs; k<=mze; k++)
01655         for (register int j=mys; j<=mye; j++)
01656           for (register int i=mxs; i<=mxe; i++) {
01657             if (cnt==18 || cnt==40)
01658               dxux=(q[base::idx(i+1,j,k,Nx,Ny)](1)-q[base::idx(i-1,j,k,Nx,Ny)](1))/(2.*dx(0));
01659             if (cnt==18 || cnt==42)
01660               dyuy=(q[base::idx(i,j+1,k,Nx,Ny)](2)-q[base::idx(i,j-1,k,Nx,Ny)](2))/(2.*dx(1));
01661             if (cnt==18 || cnt==45)
01662               dzuz=(q[base::idx(i,j,k+1,Nx,Ny)](3)-q[base::idx(i,j,k-1,Nx,Ny)](3))/(2.*dx(2));
01663             if (cnt==11 || cnt==14 || cnt==18 || cnt==44) {
01664               dyuz=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(2.*dx(1));
01665               dzuy=(q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(2.*dx(2));
01666             }
01667             if (cnt==12 || cnt==14 || cnt==18 || cnt==43) {
01668               dxuz=(q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(2.*dx(0));
01669               dzux=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(2.*dx(2));
01670             }
01671             if (cnt==13 || cnt==14 || cnt==18 || cnt==41) {         
01672               dxuy=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(2.*dx(0));
01673               dyux=(q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(2.*dx(1));
01674             }
01675             switch(cnt) {
01676             case 11: 
01677               work[base::idx(i,j,k,Nx,Ny)]=dyuz-dzuy;
01678               break;
01679             case 12: 
01680               work[base::idx(i,j,k,Nx,Ny)]=dzux-dxuz;
01681               break;
01682             case 13: 
01683               work[base::idx(i,j,k,Nx,Ny)]=dxuy-dyux;
01684               break;
01685             case 14: {
01686               DataType ox = dyuz-dzuy, oy = dzux-dxuz, oz = dxuy-dyux;
01687               work[base::idx(i,j,k,Nx,Ny)] = std::sqrt(ox*ox+oy*oy+oz*oz);
01688               break;
01689             }
01690             case 18:
01691               work[base::idx(i,j,k,Nx,Ny)] = DataType(-0.5)*(4.*dxux*dxux+4.*dyuy*dyuy+4.*dzuz*dzuz+4.*dxuy*dyux+4.*dxuz*dzux+4.*dyuz*dzuy);
01692               break;
01693             case 40:
01694               work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dxux; 
01695               break;
01696             case 41:
01697               work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuy+dyux); 
01698               break;
01699             case 42:
01700               work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dyuy; 
01701               break;
01702             case 43:
01703               work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuz+dzux); 
01704               break;
01705             case 44:
01706               work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dyuz+dzuy); 
01707               break;
01708             case 45:
01709               work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dzuz; 
01710               break;
01711             }
01712           }
01713     }
01714   }
01715 
01716   virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01717                      const int skip_ghosts=1) const {
01718     int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
01719     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
01720     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
01721     MicroType *f = (MicroType *)fvec.databuffer();
01722     DataType *work = (DataType *)workvec.databuffer();
01723     DCoords dx = base::GH().worldStep(fvec.stepsize());
01724 
01725     for (register int k=mzs; k<=mze; k++)
01726       for (register int j=mys; j<=mye; j++)
01727         for (register int i=mxs; i<=mxe; i++) {
01728           switch(cnt) {
01729           case 1:
01730             f[base::idx(i,j,k,Nx,Ny)](0)=work[base::idx(i,j,k,Nx,Ny)]/R0;
01731             break;
01732           case 2:
01733             f[base::idx(i,j,k,Nx,Ny)](1)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
01734             break;
01735           case 3:
01736             f[base::idx(i,j,k,Nx,Ny)](2)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
01737             break;
01738           case 4:
01739             f[base::idx(i,j,k,Nx,Ny)](3)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
01740             f[base::idx(i,j,k,Nx,Ny)] = Equilibrium((const MacroType &)f[base::idx(i,j,k,Nx,Ny)]);
01741             break;
01742           }
01743         }
01744   }
01745 
01746   virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
01747                     const int verbose) const {
01748     int Nx = fvec.extents(0), Ny = fvec.extents(1);
01749     int b = fvec.bottom();
01750     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01751     BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
01752     int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01753     int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01754     int mzs = std::max(bb.lower(2)/bb.stepsize(2),bbmax.lower(2)/bbmax.stepsize(2));
01755     int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01756     int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01757     int mze = std::min(bb.upper(2)/bb.stepsize(2),bbmax.upper(2)/bbmax.stepsize(2));
01758     MicroType *f = (MicroType *)fvec.databuffer();
01759 
01760     //std::cout << "CHECK Nx=" << Nx << "   Ny=" << Ny << "   Nz=" << Nz << std::endl;
01761     int result = 1;
01762     for (register int k=mzs; k<=mze; k++)
01763       for (register int j=mys; j<=mye; j++)
01764         for (register int i=mxs; i<=mxe; i++)
01765           for (register int l=0; l<base::NMicroVar(); l++)
01766             if (!(f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l)>-1.e37)) {
01767               result = 0;
01768               if (verbose) {
01769                 DCoords lbcorner = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
01770                 DCoords dx = base::GH().worldStep(fvec.stepsize());
01771                 
01772                 std::cerr << "D3Q19-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << "," << k-mdz[l] << ")(" << l << ")="
01773                           << f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l) << " " << fvec.bbox() 
01774                           << " Coords: (" << (i+0.5)*dx(0)+lbcorner(0) << "," << (j+0.5)*dx(1)+lbcorner(1) << "," 
01775                           << (k+0.5)*dx(2)+lbcorner(2) << ")" << std::endl;
01776               }
01777             }
01778     return result;
01779   }
01780 
01781   inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
01782     TensorType P;
01783     if (method[1]==0) {
01784       P(0) = q(0)*q(1)*q(1)-(f(1)+f(2)+f(7)+f(8)+f(9)+f(10)+f(15)+f(16)+f(17)+f(18))+cs2*q(0);
01785       P(1) = q(0)*q(1)*q(2)-(f(7)+f(8)-f(9)-f(10));
01786       P(2) = q(0)*q(2)*q(2)-(f(3)+f(4)+f(7)+f(8)+f(9)+f(10)+f(11)+f(12)+f(13)+f(14))+cs2*q(0);
01787       P(3) = q(0)*q(1)*q(3)-(f(15)+f(16)-f(17)-f(18));
01788       P(4) = q(0)*q(2)*q(3)-(f(11)+f(12)-f(13)-f(14));
01789       P(5) = q(0)*q(3)*q(3)-(f(5)+f(6)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18))+cs2*q(0);
01790       P *= -(DataType(1.0)-DataType(0.5)*om);
01791     }
01792     else {
01793       MicroType feq = Equilibrium(q);     
01794       P = DeviatoricStress(f,feq,om);
01795     }
01796     P(0) -= LatticeBasePressure(q(0));  
01797     P(2) -= LatticeBasePressure(q(0));
01798     P(5) -= LatticeBasePressure(q(0));
01799     return P;
01800   }
01801 
01802   inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
01803     TensorType Sigma;
01804     if (method[2]==0) {
01805       MicroType fneq = feq - f;
01806       Sigma(0) = fneq(1)+fneq(2)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
01807       Sigma(1) = fneq(7)+fneq(8)-fneq(9)-fneq(10);
01808       Sigma(2) = fneq(3)+fneq(4)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(11)+fneq(12)+fneq(13)+fneq(14);
01809       Sigma(3) = fneq(15)+fneq(16)-fneq(17)-fneq(18);
01810       Sigma(4) = fneq(11)+fneq(12)-fneq(13)-fneq(14);
01811       Sigma(5) = fneq(5)+fneq(6)+fneq(11)+fneq(12)+fneq(13)+fneq(14)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
01812       Sigma *= -(DataType(1.0)-DataType(0.5)*om);
01813     }
01814     else {
01815       MacroType q = MacroVariables(f);
01816       Sigma = Stress(f,q,om);
01817       Sigma(0) += LatticeBasePressure(q(0));  
01818       Sigma(2) += LatticeBasePressure(q(0));
01819       Sigma(5) += LatticeBasePressure(q(0));
01820     } 
01821     return Sigma;
01822   }
01823 
01824   virtual int NMethodOrder() const { return 2; }
01825 
01826   inline const DataType& L0() const { return base::LengthScale(); }
01827   inline const DataType& T0() const { return base::TimeScale(); }
01828   inline void SetDensityScale(const DataType r0) { R0 = r0; }
01829   inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
01830   inline void SetSpeedUp(const DataType s0) { S0 = s0; }
01831   inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
01832 
01833   inline const DataType& DensityScale() const { return R0; }
01834   inline const DataType VelocityScale() const { return U0/S0; }
01835   // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
01836   inline const DataType& SpeedUp() const { return S0; }
01837 
01838   // Lattice quantities for the operator
01839   inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
01840   inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
01841   inline virtual DataType LatticeBasePressure(const DataType rho) const {
01842     return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
01843   }
01844 
01845   // Physical quantities for the operator
01846   inline void SetGas(DataType rho, DataType nu, DataType cs) {
01847     rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
01848     SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
01849   }
01850   inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
01851   inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q,
01852                                               const DataType dt) const {
01853     DataType M2 = 0.0;
01854     //x-y
01855     M2 += (f(7) - feq(7))*(f(7) - feq(7));
01856     M2 += (f(8) - feq(8))*(f(8) - feq(8));
01857     M2 += (f(9) - feq(9))*(f(9) - feq(9));
01858     M2 += (f(10) - feq(10))*(f(10) - feq(10));
01859     //y-z
01860     M2 += (f(11) - feq(11))*(f(11) - feq(11));
01861     M2 += (f(12) - feq(12))*(f(12) - feq(12));
01862     M2 += (f(13) - feq(13))*(f(13) - feq(13));
01863     M2 += (f(14) - feq(14))*(f(14) - feq(14));
01864     //x-z
01865     M2 += (f(15) - feq(15))*(f(15) - feq(15));
01866     M2 += (f(16) - feq(16))*(f(16) - feq(16));
01867     M2 += (f(17) - feq(17))*(f(17) - feq(17));
01868     M2 += (f(18) - feq(18))*(f(18) - feq(18));
01869     M2 = std::sqrt(M2);
01870 
01871     DataType tau = DataType(1.0)/Omega(dt);
01872     DataType om_turb =  (DataType(2.0)/( tau
01873                                          + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
01874 
01875     return ( om_turb );
01876   }
01877   inline const DataType Omega_LES_dynamic(const MicroType &f, const MicroType &feq, const MacroType& q,
01878                                           const DataType dt) const {
01879     DataType M2 = 0.0;
01880     //x-y
01881     M2 += (f(7) - feq(7))*(f(7) - feq(7));
01882     M2 += (f(8) - feq(8))*(f(8) - feq(8));
01883     M2 += (f(9) - feq(9))*(f(9) - feq(9));
01884     M2 += (f(10) - feq(10))*(f(10) - feq(10));
01885     //y-z
01886     M2 += (f(11) - feq(11))*(f(11) - feq(11));
01887     M2 += (f(12) - feq(12))*(f(12) - feq(12));
01888     M2 += (f(13) - feq(13))*(f(13) - feq(13));
01889     M2 += (f(14) - feq(14))*(f(14) - feq(14));
01890     //x-z
01891     M2 += (f(15) - feq(15))*(f(15) - feq(15));
01892     M2 += (f(16) - feq(16))*(f(16) - feq(16));
01893     M2 += (f(17) - feq(17))*(f(17) - feq(17));
01894     M2 += (f(18) - feq(18))*(f(18) - feq(18));
01895     M2 = std::sqrt(M2);
01896 
01897     DataType tau = DataType(1.0)/Omega(dt);
01898     DataType vel = sqrt( q(1)*q(1)+q(2)*q(2)+q(3)*q(3) );
01899     DataType Cs_SmagDYN = DataType(0.5)*L0()*L0()*( (vel*vel/csp*(M2-vel)) / ((M2-vel)*(M2-vel))  );
01900     DataType om_turb =  (DataType(2.0)/( tau
01901                                          + sqrt(tau*tau + (DataType(18.)*Cs_SmagDYN*M2)/(R0*q(0)*cs2/S0/S0)) ));
01902 
01903     return ( om_turb );
01904   }
01905   inline const int TurbulenceType() const { return turbulence; }
01906   inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
01907   inline const DataType& GasDensity() const { return rhop; }
01908   inline const DataType& GasViscosity() const { return nup; }
01909   inline const DataType& GasSpeedofSound() const { return csp; }
01910   inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
01911   { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
01912   // Quantities for output only
01913   inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
01914   inline virtual const DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); } 
01915   inline virtual const DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); } 
01916 
01917 protected:
01918   DataType cs2, cs22, cssq;
01919   DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
01920   DataType Cs_Smagorinsky, turbulence;
01921   int method[3], mdx[19], mdy[19], mdz[19];
01922 };
01923 
01924 
01925 #endif