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