amroc/lbm/LBMD2Q9Thermal_BC.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Working still in progress: Incompressible LBM operator
00004 
00005 // Copyright (C) 2012 Oak Ridge National Laboratory
00006 // Ralf Deiterding, deiterdingr@ornl.gov
00007 
00008 #ifndef LBM_D2Q9_H
00009 #define LBM_D2Q9_H
00010 
00018 #include "LBMBase.h"
00019 #include <cfloat>
00020 
00021 #define LB_FACTOR 1.0e5
00022 
00047 template <class DataType>
00048 class LBMD2Q9Thermal : public LBMBase<Vector<DataType,14>, Vector<DataType,4>, 2> {
00049   typedef LBMBase<Vector<DataType,14>, Vector<DataType,4>, 2> base;
00050 public:
00051   typedef typename base::vec_grid_data_type vec_grid_data_type;
00052   typedef typename base::grid_data_type grid_data_type;
00053   typedef typename base::MicroType MicroType;
00054   typedef typename base::MacroType MacroType;
00055   typedef GridData<MacroType,2> macro_grid_data_type;
00056   typedef typename base::SideName SideName;
00057   typedef typename base::point_type point_type;
00058 
00059   enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro }; 
00060   enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall, ExtrapolationEquilibrium, SlipWallTemperature, NoSlipWallTemperature,InletTemp }; 
00061   enum GFMPredefined { GFMExtrapolation, GFMSlipWallTemp, GFMNoSlipWallTemp, GFMBounceBack };
00062   enum TurbulenceModel { laminar, LES_Smagorinsky };
00063 
00064   LBMD2Q9Thermal() : base(), P0(1.), R0(1.), U0(1.), S0(1.), Temp0(1.), rhop(0.), Tpmin(0.), Tpmax(0.), Tref(0.5), csp(0.), 
00065                      cs2p(0.), nup(0.), nutp(0.) , prp(0.) {
00066     cs2  = DataType(1.)/DataType(3.);
00067     cs22 = DataType(2.)*cs2;
00068     cssq = DataType(2.)/DataType(9.);
00069     ds2  = DataType(1.)/DataType(2.);
00070     dcs2 = ds2/cs2;
00071     method[0] = 2;
00072     turbulence = laminar;
00073 
00074     mdx[0]=mdx[3]=mdx[6]=mdx[9]=mdx[12]=mdx[13]=0; mdx[1]=mdx[4]=mdx[7]=mdx[10]=1; mdx[2]=mdx[5]=mdx[8]=mdx[11]=-1;
00075     mdy[0]=mdy[1]=mdy[2]=mdy[9]=mdy[10]=mdy[11]=0; mdy[3]=mdy[4]=mdy[5]=mdy[12]=1; mdy[6]=mdy[7]=mdy[8]=mdy[13]=-1;
00076     
00077   }
00078  
00079   virtual ~LBMD2Q9Thermal() {}
00080 
00081   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00082     base::LocCtrl = Ctrl.getSubDevice(prefix+"D2Q9");
00083     RegisterAt(base::LocCtrl,"Method(1)",method[0]); 
00084     RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00085     RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00086     RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00087     RegisterAt(base::LocCtrl,"Gas_Cs",csp); 
00088     RegisterAt(base::LocCtrl,"Gas_nu",nup); 
00089     RegisterAt(base::LocCtrl,"Gas_nuT",nutp);
00090     RegisterAt(base::LocCtrl,"Gas_pr",prp);  
00091     RegisterAt(base::LocCtrl,"Gas_Tmin",Tpmin); 
00092     RegisterAt(base::LocCtrl,"Gas_Tmax",Tpmax); 
00093     RegisterAt(base::LocCtrl,"SpeedUp",S0);  
00094   }
00095 
00096   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
00097     base::SetupData(gh,ghosts);
00098     if (rhop>0.) R0 = rhop;
00099     if (prp>0.) P0 = prp;
00100     if (Tpmin!=0. || Tpmax!=1.) SetTemperatureScale(Tpmin,Tpmax);
00101     if (csp>0.) {
00102       cs2p = csp*csp;     
00103       SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00104     }   
00105     std::cout << "D2Q9Thermal: Gas_rho=" << rhop << "   Gas_Cs=" << csp 
00106               << "   Gas_nu=" << nup << "   Gas_nut=" << nutp << "   Gas_pr=" << prp << std::endl; 
00107     std::cout << "D2Q9Thermal: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0 
00108               << "   R0=" << R0 << "   Temp0=" << Temp0 << "   S0=" << S0 << "   P0=" << P0 << std::endl;
00109 
00110   }
00111 
00112   inline virtual MacroType MacroVariables(const MicroType &f) const {
00113     MacroType q;
00114     DataType part1, part2;
00115 
00116     q(1)=f(1)-f(2)+f(4)-f(5)+f(7)-f(8);
00117     q(2)=f(3)+f(4)+f(5)-f(6)-f(7)-f(8);
00118 
00119     part1 = f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8);
00120     part2 = cs22*(std::pow(q(1),2)+std::pow(q(2),2));
00121     q(0)= (DataType(3.)/DataType(5.))*(part1 - part2);
00122 
00123     q(3)=f(10)+f(11)+f(12)+f(13);
00124 
00125     return q;
00126   }
00127 
00128   inline virtual MicroType Equilibrium(const MacroType &q) const {
00129     MicroType feq;
00130 
00131     DataType p = q(0);
00132     DataType u = q(1);
00133     DataType v = q(2);
00134     DataType T = q(3);
00135 
00136     DataType ri, rt0, rt1, rt2, pre0, pre1, pre2, Temp;
00137 
00138     ri = DataType(0.);
00139     rt0 = DataType(4.) / DataType(9.);
00140     rt1 = DataType(1.) / DataType(9.);
00141     rt2 = DataType(1.) / DataType(36.);
00142 
00143     pre0 = -(DataType(5.) / DataType(3.)) * p;
00144     pre1 = (DataType(1.) / DataType(3.)) * p;
00145     pre2 = (DataType(1.) / DataType(12.)) * p;
00146 
00147     Temp = T/DataType(4.);
00148                
00149     DataType usq = u * u; 
00150     DataType vsq = v * v;
00151     DataType sumsq = (usq + vsq) / cs22;
00152     DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00153     DataType u2 = usq / cssq; 
00154     DataType v2 = vsq / cssq;
00155     DataType ui = u / cs2;
00156     DataType vi = v / cs2;
00157     DataType uv = ui * vi;
00158     DataType uT = DataType(2.)*u;
00159     DataType vT = DataType(2.)*v;
00160 
00161     feq(0) = pre0 + rt0 * (ri - sumsq);
00162                
00163     feq(1) = pre1 + rt1 * (ri - sumsq + u2 + ui);
00164     feq(2) = pre1 + rt1 * (ri - sumsq + u2 - ui);
00165     feq(3) = pre1 + rt1 * (ri - sumsq + v2 + vi);
00166     feq(6) = pre1 + rt1 * (ri - sumsq + v2 - vi);
00167                
00168     feq(4) = pre2 + rt2 * (ri + sumsq2 + ui + vi + uv);
00169     feq(5) = pre2 + rt2 * (ri + sumsq2 - ui + vi - uv);
00170     feq(7) = pre2 + rt2 * (ri + sumsq2 + ui - vi - uv);
00171     feq(8) = pre2 + rt2 * (ri + sumsq2 - ui - vi + uv);   
00172 
00173     feq(9) = DataType(0.);
00174     feq(10) = Temp*(DataType(1.)+uT);
00175     feq(11) = Temp*(DataType(1.)-uT);
00176     feq(12) = Temp*(DataType(1.)+vT);
00177     feq(13) = Temp*(DataType(1.)-vT);
00178 
00179     return feq;
00180   }
00181 
00182   inline virtual void Collision(MicroType &f, const DataType dt) const {
00183     MacroType q = MacroVariables(f);
00184     MicroType feq = Equilibrium(q);
00185 
00186     DataType force = DataType(0.5)*dt/S0*gbeta*(q(3)-Tref);
00187 
00188     DataType omega = Omega(dt), omegaT = OmegaT(dt);
00189     for (int i=0; i<=8;i++){
00190         f(i)=(DataType(1.)-omega)*f(i) + omega*feq(i);
00191         if(i==3) f(i)+= force;
00192         if(i==6) f(i)-= force;
00193     }
00194     for (int i=9; i<=13;i++){
00195         f(i)=(DataType(1.)-omegaT)*f(i) + omegaT*feq(i);
00196     }
00197    
00198   }     
00199 
00200   virtual int IncomingIndices(const int side, int indices[]) const {
00201     switch (side) {
00202     case base::Left:
00203       indices[0] = 1; indices[1] = 4; indices[2] = 7; indices[3] = 10;
00204       break;
00205     case base::Right:
00206       indices[0] = 2; indices[1] = 5; indices[2] = 8; indices[3] = 11;
00207       break;
00208     case base::Bottom:
00209       indices[0] = 3; indices[1] = 4; indices[2] = 5; indices[3] = 12; 
00210       break;
00211     case base::Top:
00212       indices[0] = 6; indices[1] = 7; indices[2] = 8; indices[3] = 13; 
00213       break;
00214     }
00215     return 4;
00216   }
00217 
00218   virtual int OutgoingIndices(const int side, int indices[]) const {
00219     switch (side) {
00220     case base::Left:
00221       indices[0] = 2; indices[1] = 5; indices[2] = 8; indices[3] = 11; 
00222       break;
00223     case base::Right:
00224       indices[0] = 1; indices[1] = 4; indices[2] = 7; indices[3] = 10; 
00225       break;
00226     case base::Bottom:
00227       indices[0] = 6; indices[1] = 7; indices[2] = 8; indices[3] = 13; 
00228       break;
00229     case base::Top:
00230       indices[0] = 3; indices[1] = 4; indices[2] = 5; indices[3] = 12; 
00231       break;
00232     }
00233     return 4;
00234   }
00235 
00236   // Revert just one layer of cells at the boundary 
00237   virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
00238                              const int side) const {
00239     int Nx = fvec.extents(0);
00240     int b = fvec.bottom();
00241     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00242     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00243     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00244     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00245     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00246     MicroType *f = (MicroType *)fvec.databuffer();
00247 
00248 #ifdef DEBUG_PRINT_FIXUP
00249     ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
00250       << fvec.bbox() << " using " << bb << "\n").flush();
00251 #endif
00252     
00253     switch (side) {
00254     case base::Left:
00255       for (register int j=mys; j<=mye; j++) {
00256         f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](2);
00257         f[b+base::idx(mxs,j,Nx)](11)= f[b+base::idx(mxs-1,j,Nx)](11);
00258         f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](8);
00259         f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](5);
00260       }
00261       break;
00262     case base::Right:
00263       for (register int j=mys; j<=mye; j++) {
00264         f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](1);
00265         f[b+base::idx(mxe,j,Nx)](10)= f[b+base::idx(mxe+1,j,Nx)](10);
00266         f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](7);
00267         f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](4);
00268       }
00269       break;
00270     case base::Bottom:
00271       for (register int i=mxs; i<=mxe; i++) {
00272         f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](6);
00273         f[b+base::idx(i,mys,Nx)](13)= f[b+base::idx(i,mys-1,Nx)](13);
00274         f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](8);
00275         f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](7);
00276       }
00277       break;
00278     case base::Top:
00279       for (register int i=mxs; i<=mxe; i++) {
00280         f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](3);
00281         f[b+base::idx(i,mye,Nx)](12)= f[b+base::idx(i,mye+1,Nx)](12);
00282         f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](5);
00283         f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](4);
00284       }
00285       break;
00286     }
00287   }   
00288 
00289   virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec, 
00290                          const BBox &bb, const double& dt) const {
00291     int Nx = fvec.extents(0);
00292     int b = fvec.bottom();
00293     MicroType *f = (MicroType *)fvec.databuffer();
00294     MicroType *of = (MicroType *)ovec.databuffer();
00295 
00296 #ifdef DEBUG_PRINT_FIXUP
00297     ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox() 
00298       << " using " << bb << "\n").flush();
00299 #endif
00300     
00301     assert (fvec.extents(0)==ovec.extents(0));
00302     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) &&
00303             fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1));
00304     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00305     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00306     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00307     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00308 
00309     for (register int j=mys; j<=mye; j++) 
00310       for (register int i=mxs; i<=mxe; i++) {
00311         for (register int n=0; n<base::NMicroVar(); n++) 
00312           f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00313         Collision(f[b+base::idx(i,j,Nx)],dt);
00314       }
00315   }  
00316 
00317   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00318                       const double& t, const double& dt, const int& mpass) const {
00319     // Make sure that the physical length and times are scaling consistently
00320     // into lattice quantities
00321     DCoords dx = base::GH().worldStep(fvec.stepsize());
00322 
00323 
00324     if (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR) {
00325       std::cerr << "LBMD2Q9 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00326                 << " to be equal." << std::endl << std::flush;
00327       std::exit(-1);
00328     }
00329     if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
00330       std::cerr << "LBMD2Q9 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0() 
00331                 << " to be equal." 
00332                 << "   dt=" << dt << "   TO=" << T0()
00333                 << std::endl << std::flush;
00334       std::exit(-1);
00335     }
00336 
00337     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00338     MicroType *f = (MicroType *)fvec.databuffer();
00339     
00340 
00341     // Streaming - update all cells
00342     for (register int j=Ny-1; j>=1; j--)
00343       for (register int i=0; i<=Nx-2; i++) {
00344         f[base::idx(i,j,Nx)](3) = f[base::idx(i,j-1,Nx)](3);
00345         f[base::idx(i,j,Nx)](12) = f[base::idx(i,j-1,Nx)](12);
00346         f[base::idx(i,j,Nx)](5) = f[base::idx(i+1,j-1,Nx)](5);  
00347       }
00348 
00349     for (register int j=Ny-1; j>=1; j--)
00350       for (register int i=Nx-1; i>=1; i--) {
00351         f[base::idx(i,j,Nx)](1) = f[base::idx(i-1,j,Nx)](1);
00352         f[base::idx(i,j,Nx)](10) = f[base::idx(i-1,j,Nx)](10);
00353         f[base::idx(i,j,Nx)](4) = f[base::idx(i-1,j-1,Nx)](4);
00354       }
00355 
00356     for (register int j=0; j<=Ny-2; j++)
00357       for (register int i=Nx-1; i>=1; i--) {
00358         f[base::idx(i,j,Nx)](6) = f[base::idx(i,j+1,Nx)](6);
00359         f[base::idx(i,j,Nx)](13) = f[base::idx(i,j+1,Nx)](13);
00360         f[base::idx(i,j,Nx)](7) = f[base::idx(i-1,j+1,Nx)](7);
00361       }
00362 
00363     for (register int j=0; j<=Ny-2; j++)
00364       for (register int i=0; i<=Nx-2; i++) {
00365         f[base::idx(i,j,Nx)](2) = f[base::idx(i+1,j,Nx)](2);
00366         f[base::idx(i,j,Nx)](11) = f[base::idx(i+1,j,Nx)](11);
00367         f[base::idx(i,j,Nx)](8) = f[base::idx(i+1,j+1,Nx)](8);
00368       }    
00369 
00370     // Collision - only in the interior region
00371     int mxs = base::NGhosts(), mys = base::NGhosts();
00372     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00373     for (register int j=mys; j<=mye; j++)
00374       for (register int i=mxs; i<=mxe; i++)
00375         Collision(f[base::idx(i,j,Nx)],dt);
00376     
00377     return DataType(1.);
00378   }
00379 
00380    virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00381                            DataType* aux=0, const int naux=0, const int scaling=0) const {
00382     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00383     int mxs = base::NGhosts(), mys = base::NGhosts();
00384     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00385     MicroType *f = (MicroType *)fvec.databuffer();
00386 
00387     if (type==ConstantMicro) {
00388       MicroType g;
00389       for (register int i=0; i<base::NMicroVar(); i++)
00390         g(i) = aux[i];
00391     
00392       for (register int j=mys; j<=mye; j++)
00393         for (register int i=mxs; i<=mxe; i++)
00394           f[base::idx(i,j,Nx)] = g;
00395     }
00396     
00397     else if (type==ConstantMacro) {
00398       MacroType q;
00399       if (scaling==base::Physical) {
00400         q(0) = aux[0]/P0; 
00401         q(1) = aux[1]/U0; 
00402         q(2) = aux[2]/U0;
00403         q(3) = (aux[3]-Tpmin)/Temp0;
00404       }
00405       else 
00406         for (register int i=0; i<base::NMacroVar(); i++)
00407           q(i) = aux[i];
00408       q(1) *= S0; q(2) *= S0;
00409     
00410       for (register int j=mys; j<=mye; j++)
00411         for (register int i=mxs; i<=mxe; i++)
00412           f[base::idx(i,j,Nx)] = Equilibrium(q);
00413     }
00414     
00415     else if (type==GasAtRest) {
00416       MacroType q;
00417       q(0) = GasDensity()/R0; 
00418       q(1) = DataType(0.); 
00419       q(2) = DataType(0.);
00420       for (register int j=mys; j<=mye; j++)
00421         for (register int i=mxs; i<=mxe; i++)
00422           f[base::idx(i,j,Nx)] = Equilibrium(q);
00423     }       
00424   }
00425 
00426   // Set just the one layer of ghost cells at the boundary 
00427   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00428                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00429     int Nx = fvec.extents(0);
00430     int b = fvec.bottom();
00431 
00432     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00433     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00434     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00435     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00436     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00437     MicroType *f = (MicroType *)fvec.databuffer();
00438  
00439     if (type==Symmetry) {
00440       switch (side) {
00441       case base::Left:
00442         for (register int j=mys; j<=mye; j++) {
00443           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j,Nx)](8);
00444           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00445           f[b+base::idx(mxe,j,Nx)](10) = f[b+base::idx(mxe+1,j,Nx)](11);
00446           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j,Nx)](5);
00447         }
00448         break;
00449       case base::Right:
00450         for (register int j=mys; j<=mye; j++) {
00451           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j,Nx)](7);
00452           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00453           f[b+base::idx(mxs,j,Nx)](11) = f[b+base::idx(mxs-1,j,Nx)](10);
00454           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j,Nx)](4);
00455         }
00456         break;
00457       case base::Bottom:
00458         for (register int i=mxs; i<=mxe; i++) {
00459           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i,mye+1,Nx)](8);
00460           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00461           f[b+base::idx(i,mye,Nx)](12) = f[b+base::idx(i,mye+1,Nx)](13);
00462           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i,mye+1,Nx)](7);
00463         }
00464         break;
00465       case base::Top:
00466         for (register int i=mxs; i<=mxe; i++) {
00467           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i,mys-1,Nx)](5);
00468           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00469           f[b+base::idx(i,mys,Nx)](13) = f[b+base::idx(i,mys-1,Nx)](12);
00470           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i,mys-1,Nx)](4);
00471         }
00472         break;
00473       }
00474     }
00475     
00476     else if (type==SlipWall) {
00477       switch (side) {
00478       case base::Left:
00479         for (register int j=mys; j<=mye; j++) {
00480           if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5); 
00481           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00482           f[b+base::idx(mxe,j,Nx)](10) = f[b+base::idx(mxe+1,j,Nx)](11);
00483           if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00484         }
00485         break;
00486       case base::Right:
00487         for (register int j=mys; j<=mye; j++) {
00488           if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4); 
00489           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00490           f[b+base::idx(mxs,j,Nx)](11) = f[b+base::idx(mxs-1,j,Nx)](10);
00491           if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00492         }
00493         break;
00494       case base::Bottom:
00495         for (register int i=mxs; i<=mxe; i++) {
00496           if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00497           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00498           f[b+base::idx(i,mye,Nx)](12) = f[b+base::idx(i,mye+1,Nx)](13);
00499           if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00500         }
00501         break;
00502       case base::Top:
00503         for (register int i=mxs; i<=mxe; i++) {
00504           if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00505           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00506           f[b+base::idx(i,mys,Nx)](13) = f[b+base::idx(i,mys-1,Nx)](12);
00507           if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00508         }
00509         break;
00510       }
00511     }
00512    
00513     else if (type==NoSlipWall) {
00514       switch (side) {
00515       case base::Left:
00516         for (register int j=mys; j<=mye; j++) {
00517           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00518           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00519           f[b+base::idx(mxe,j,Nx)](10) = f[b+base::idx(mxe+1,j,Nx)](11);
00520           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00521         }
00522         break;
00523       case base::Right:
00524         for (register int j=mys; j<=mye; j++) {
00525           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00526           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00527           f[b+base::idx(mxs,j,Nx)](11) = f[b+base::idx(mxs-1,j,Nx)](10);
00528           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00529         }
00530         break;  
00531       case base::Bottom:
00532         for (register int i=mxs; i<=mxe; i++) {
00533           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00534           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00535           f[b+base::idx(i,mye,Nx)](12) = f[b+base::idx(i,mye+1,Nx)](13);
00536           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00537         }
00538         break;  
00539       case base::Top:
00540         for (register int i=mxs; i<=mxe; i++) {
00541           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00542           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00543           f[b+base::idx(i,mys,Nx)](13) = f[b+base::idx(i,mys-1,Nx)](12);
00544           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00545         }
00546         break;
00547       }
00548     }
00549     // 
00550     else if (type==Inlet) {
00551       if (aux==0 || naux<=0) return;
00552       DataType a0=S0*aux[0], a1=0.;
00553       if (naux>1) a1 = S0*aux[1];
00554       if (scaling==base::Physical) {
00555         a0 /= U0;
00556         a1 /= U0;
00557       }
00558       switch (side) {
00559       case base::Left:
00560         for (register int j=mys; j<=mye; j++) {
00561           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00562           q(1) = a0;
00563           if (naux>1) q(2) = a1;
00564           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00565         }
00566         break;
00567       case base::Right:
00568         for (register int j=mys; j<=mye; j++) {
00569           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00570           q(1) = a0;
00571           if (naux>1) q(2) = a1;
00572           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00573         }
00574         break;
00575       case base::Bottom:
00576         for (register int i=mxs; i<=mxe; i++) {
00577           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00578           if (naux==1) 
00579             q(2) = a0;
00580           else {
00581             q(1) = a0;
00582             q(2) = a1;
00583           }
00584           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00585         }
00586         break;
00587       case base::Top:
00588         for (register int i=mxs; i<=mxe; i++) {
00589           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00590           if (naux==1)
00591             q(2) = a0;
00592           else {
00593             q(1) = a0;
00594             q(2) = a1;
00595           }
00596           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00597         }
00598         break;
00599       }
00600     }
00601     // Without temperature model
00602     else if (type==Outlet) {
00603       switch (side) {
00604       case base::Left:
00605         for (register int j=mys; j<=mye; j++) { 
00606                 f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)];
00607         }
00608         break;
00609       case base::Right:
00610         for (register int j=mys; j<=mye; j++) {
00611                 f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)];
00612 
00613                 //f[b+base::idx(mxs,j,Nx)](1) = 2.*f[b+base::idx(mxs-1,j,Nx)](1)-f[b+base::idx(mxs-2,j,Nx)](1);
00614                 //f[b+base::idx(mxs,j,Nx)](4) = 2.*f[b+base::idx(mxs-1,j,Nx)](4)-f[b+base::idx(mxs-2,j,Nx)](4);
00615                 //f[b+base::idx(mxs,j,Nx)](7) = 2.*f[b+base::idx(mxs-1,j,Nx)](7)-f[b+base::idx(mxs-2,j,Nx)](7);
00616                 //f[b+base::idx(mxs,j,Nx)](10) = 2.*f[b+base::idx(mxs-1,j,Nx)](10)-f[b+base::idx(mxs-2,j,Nx)](10);                
00617           }
00618         break;
00619       case base::Bottom:
00620         for (register int i=mxs; i<=mxe; i++) {
00621                 f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)];
00622         }
00623         break;
00624       case base::Top:
00625         for (register int i=mxs; i<=mxe; i++) {
00626                 f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)];
00627         }
00628         break;
00629       }
00630     }
00631     
00632     // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)    
00633     // Without temperature model
00634     else if (type==Pressure) {
00635       if (aux==0 || naux<=0) return;
00636       DataType a0=aux[0];
00637       if (scaling==base::Physical)
00638         a0 /= R0;
00639       switch (side) {
00640       case base::Left:
00641         for (register int j=mys; j<=mye; j++) {
00642           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00643           q(0) = a0;
00644           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00645         }
00646         break;
00647       case base::Right:
00648         for (register int j=mys; j<=mye; j++) {
00649           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00650           q(0) = a0;
00651           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00652         }
00653         break;
00654       case base::Bottom:
00655         for (register int i=mxs; i<=mxe; i++) {
00656           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00657           q(0) = a0;
00658           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00659         }
00660         break;
00661       case base::Top:
00662         for (register int i=mxs; i<=mxe; i++) {
00663           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00664           q(0) = a0;
00665           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00666         }
00667         break;
00668       }
00669     }
00670     //Without temperature model
00671     else if (type==SlidingWall) {
00672       if (aux==0 || naux<=0) return;
00673       if (aux==0 || naux<=0) return;
00674       DataType a0=S0*aux[0];
00675       if (scaling==base::Physical)
00676         a0 /= U0;
00677       switch (side) {
00678       case base::Left:
00679         for (register int j=mys; j<=mye; j++) {
00680           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00681           DataType rd1 = f[b+base::idx(mxe+1,j,Nx)](5), rd2 = f[b+base::idx(mxe+1,j,Nx)](8);
00682           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00683           DataType pw = DataType(1.)-qw;
00684           if (j<mye) f[b+base::idx(mxe,j+1,Nx)](7) = pw*rd1+qw*rd2;
00685           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00686           if (j>mys) f[b+base::idx(mxe,j-1,Nx)](4) = qw*rd1+pw*rd2;
00687         }
00688         break;
00689       case base::Right:
00690         for (register int j=mys; j<=mye; j++) {
00691           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00692           DataType rd1 = f[b+base::idx(mxs-1,j,Nx)](4), rd2 = f[b+base::idx(mxs-1,j,Nx)](7);
00693           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00694           DataType pw = DataType(1.)-qw;
00695           if (j<mye) f[b+base::idx(mxs,j+1,Nx)](8) = pw*rd1+qw*rd2;
00696           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00697           if (j>mys) f[b+base::idx(mxs,j-1,Nx)](5) = qw*rd1+pw*rd2;
00698         }
00699         break;
00700       case base::Bottom:
00701         for (register int i=mxs; i<=mxe; i++) {
00702           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00703           DataType rd1 = f[b+base::idx(i,mye+1,Nx)](7), rd2 = f[b+base::idx(i,mye+1,Nx)](8);
00704           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00705           DataType pw = DataType(1.)-qw;
00706           if (i<mxe) f[b+base::idx(i+1,mye,Nx)](5) = pw*rd1+qw*rd2;
00707           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00708           if (i>mxs) f[b+base::idx(i-1,mye,Nx)](4) = qw*rd1+pw*rd2;
00709         }
00710         break;
00711       case base::Top:
00712         for (register int i=mxs; i<=mxe; i++) {
00713           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00714           DataType rd1 = f[b+base::idx(i,mys-1,Nx)](4), rd2 = f[b+base::idx(i,mys-1,Nx)](5);
00715           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00716           DataType pw = DataType(1.)-qw;
00717           if (i<mxe) f[b+base::idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00718           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00719           if (i>mxs) f[b+base::idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00720         }
00721         break;
00722       }
00723     }
00724 
00725     // (Z.Guo, B. Shi and C.Zheng - A coupled lattice BGK model for the boussinesq equations)    
00726     else if (type==ExtrapolationEquilibrium) {
00727       switch (side) {
00728       case base::Left:
00729         for (register int j=mys; j<=mye; j++) {
00730           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00731           MacroType qBC = q;
00732           if (aux==0 || naux<=0) return;
00733           for(int k=0; k<=3; k++){
00734                 if(aux[k] != 0){
00735                   if(k==0) qBC(k) = aux[k+4]/P0;
00736                         else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);
00737                         else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00738                 }
00739           }
00740           //f[b+base::idx(mxe,j,Nx)] = Equilibrium(qBC)+f[b+base::idx(mxe+1,j,Nx)]-Equilibrium(q);
00741 
00742           DataType Force = DataType(2.)/DataType(3.);
00743           DataType ForceT = (DataType(1.)/DataType(2.))*qBC(3);
00744           DataType Vel = qBC(1)+qBC(2);
00745 
00746 
00747           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2)+Force*qBC(1);
00748           if(j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8)+Force/DataType(4.)*Vel;
00749           if(j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5)+Force/DataType(4.)*Vel;
00750           //MINUS HINZUFÜGEN!!!!
00751           f[b+base::idx(mxe,j,Nx)](10) = -f[b+base::idx(mxe+1,j,Nx)](11)+ForceT;
00752 
00753         }
00754         break;
00755 
00756      case base::Right:
00757         for (register int j=mys; j<=mye; j++) {
00758           
00759           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00760           MacroType qBC = q;
00761           if (aux==0 || naux<=0) return;
00762           for(int k=0; k<=3; k++){
00763                 if(aux[k] != 0){
00764                         if(k==0) qBC(k) = aux[k+4]/P0;
00765                         else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);  
00766                         else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00767                 }
00768           }
00769           //f[b+base::idx(mxs,j,Nx)] = Equilibrium(qBC)+f[b+base::idx(mxs-1,j,Nx)]-Equilibrium(q);
00770 
00771           DataType Force = DataType(2.)/DataType(3.);
00772           DataType ForceT = (DataType(1.)/DataType(2.))*qBC(3);
00773           DataType Vel = qBC(1)+qBC(2);
00774           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1)+Force*qBC(1);
00775           if(j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7)+Force/DataType(4.)*Vel;
00776           if(j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4)+Force/DataType(4.)*Vel;
00777           f[b+base::idx(mxs,j,Nx)](11) = -f[b+base::idx(mxs-1,j,Nx)](10)+ForceT;
00778           
00779         }
00780         break;
00781 
00782      case base::Bottom:
00783         for (register int i=mxs; i<=mxe; i++) {
00784           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00785           MacroType qBC = q;
00786           if (aux==0 || naux<=0) return;
00787           for(int k=0; k<=3; k++){
00788                 if(aux[k] != 0){
00789                         if(k==0) qBC(k) = aux[k+4]/P0;
00790                         else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);  
00791                         else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00792                 }
00793           }
00794           //f[b+base::idx(i,mye,Nx)] = Equilibrium(qBC)+f[b+base::idx(i,mye+1,Nx)]-Equilibrium(q);
00795 
00796        
00797           DataType Force = DataType(2.)/DataType(3.);
00798           DataType ForceT = (DataType(1.)/DataType(2.))*qBC(3);
00799           DataType Vel = qBC(1)+qBC(2);
00800       
00801           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6)+Force*qBC(1);
00802           if(i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7)+Force/DataType(4.)*Vel;
00803           if(i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8)+Force/DataType(4.)*Vel;
00804           f[b+base::idx(i,mye,Nx)](12) = -f[b+base::idx(i,mye+1,Nx)](13)+ForceT;
00805     
00806 
00807         }
00808         break;
00809 
00810       case base::Top:
00811         for (register int i=mxs; i<=mxe; i++) { 
00812 
00813           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00814           MacroType qBC = q;
00815           if (aux==0 || naux<=0) return;
00816           for(int k=0; k<=3; k++){
00817                 if(aux[k] != 0){
00818                         if(k==0) qBC(k) = aux[k+4]/P0;
00819                         else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);  
00820                         else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00821                 }
00822           }
00823           //f[b+base::idx(i,mys,Nx)] = Equilibrium(qBC)+f[b+base::idx(i,mys-1,Nx)]-Equilibrium(q);
00824 
00825           DataType Force = DataType(2.)/DataType(3.);
00826           DataType ForceT = (DataType(1.)/DataType(2.))*qBC(3);
00827           DataType Vel = qBC(1)+qBC(2);
00828       
00829           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3)+Force*qBC(1);
00830           if(i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4)+Force/DataType(4.)*Vel;
00831           if(i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5)+Force/DataType(4.)*Vel;
00832           f[b+base::idx(i,mys,Nx)](13) = -f[b+base::idx(i,mys-1,Nx)](12)+ForceT;
00833         }
00834         break;
00835 
00836       }
00837     }
00838 
00839 else if (type==SlipWallTemperature) {
00840       switch (side) {
00841       case base::Left:
00842         for (register int j=mys; j<=mye; j++) {
00843           MacroType qBC;
00844           qBC(3) = (aux[0]-Tpmin)/Temp0;
00845           MicroType gBC = Equilibrium(qBC);
00846           if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5); 
00847           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00848           f[b+base::idx(mxe,j,Nx)](10) = gBC(10);
00849           if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00850         }
00851         break;
00852       case base::Right:
00853         for (register int j=mys; j<=mye; j++) {
00854           MacroType qBC;
00855           qBC(3) = (aux[0]-Tpmin)/Temp0;
00856           MicroType gBC = Equilibrium(qBC);
00857           if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4); 
00858           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00859           f[b+base::idx(mxs,j,Nx)](11) = gBC(11);
00860           if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00861         }
00862         break;
00863       case base::Bottom:
00864         for (register int i=mxs; i<=mxe; i++) {
00865           MacroType qBC;
00866           qBC(3) = (aux[0]-Tpmin)/Temp0;
00867           MicroType gBC = Equilibrium(qBC);
00868           if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00869           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00870           f[b+base::idx(i,mye,Nx)](12) = gBC(12);
00871           if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00872         }
00873         break;
00874       case base::Top:
00875         for (register int i=mxs; i<=mxe; i++) {
00876           MacroType qBC;
00877           qBC(3) = (aux[0]-Tpmin)/Temp0;
00878           MicroType gBC = Equilibrium(qBC);
00879           if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00880           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00881           f[b+base::idx(i,mys,Nx)](13) = gBC(13);
00882           if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00883         }
00884         break;
00885       }
00886     }
00887    
00888     else if (type==NoSlipWallTemperature) {
00889       switch (side) {
00890       case base::Left:
00891         for (register int j=mys; j<=mye; j++) {
00892           //MacroType qBC;
00893           MacroType qBC = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00894           qBC(3) = (aux[0]-Tpmin)/Temp0;
00895           MicroType gBC = Equilibrium(qBC);
00896           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00897           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00898           f[b+base::idx(mxe,j,Nx)](10) = gBC(10);
00899           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00900 
00901           //Method1:
00902           //f[b+base::idx(mxe,j,Nx)](10)=DataType(2.)/DataType(4.)*qBC(3)-f[b+base::idx(mxe+1,j,Nx)](11);
00903           //Method2: 
00904           //f[b+base::idx(mxe,j,Nx)](10)=qBC(3)*(DataType(1.)/DataType(2.)-f[b+base::idx(mxe+1,j,Nx)](11));
00905 
00906 
00907         }
00908         break;
00909       case base::Right:
00910         for (register int j=mys; j<=mye; j++) {
00911           //MacroType qBC;
00912           MacroType qBC = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00913           qBC(3) = (aux[0]-Tpmin)/Temp0;
00914           MicroType gBC = Equilibrium(qBC);
00915           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00916           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00917           f[b+base::idx(mxs,j,Nx)](11) = gBC(11);
00918           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);         
00919            //Method1:
00920           //f[b+base::idx(mxs,j,Nx)](11)=DataType(2.)/DataType(4.)*qBC(3)-f[b+base::idx(mxs-1,j,Nx)](10);
00921           //Method2: 
00922           //f[b+base::idx(mxs,j,Nx)](11)=qBC(3)*(DataType(1.)/DataType(2.)-f[b+base::idx(mxs-1,j,Nx)](10));
00923 
00924         }
00925         break;  
00926       case base::Bottom:
00927         for (register int i=mxs; i<=mxe; i++) {
00928           MacroType qBC;
00929           qBC(3) = (aux[0]-Tpmin)/Temp0;
00930           MicroType gBC = Equilibrium(qBC);
00931           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00932           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00933           f[b+base::idx(i,mye,Nx)](12) = gBC(12);
00934           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00935         }
00936         break;  
00937       case base::Top:
00938         for (register int i=mxs; i<=mxe; i++) {
00939           MacroType qBC;
00940           qBC(3) = (aux[0]-Tpmin)/Temp0;
00941           MicroType gBC = Equilibrium(qBC);
00942           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00943           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00944           f[b+base::idx(i,mys,Nx)](13) = gBC(13);
00945           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00946         }
00947         break;
00948       }
00949     }
00950 
00951    else if (type==InletTemp) {
00952       if (aux==0 || naux<=0) return;
00953       //a0: Vel, a1: Temp, a2: Vel
00954       DataType a0=S0*aux[0], a1=aux[1], a2=0.;
00955       if (naux>2) a2 = S0*aux[2];
00956       if (scaling==base::Physical) {
00957         a0 /= U0;
00958         a1 = (a1-Tpmin)/Temp0;
00959         a2 /= U0;
00960       }
00961       switch (side) {
00962       case base::Left:
00963         for (register int j=mys; j<=mye; j++) {
00964           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00965           q(1) = a0;
00966           q(3) = a1;
00967           if (naux>2) q(2) = a2;
00968           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00969         }
00970         break;
00971       case base::Right:
00972         for (register int j=mys; j<=mye; j++) {
00973           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00974           q(1) = a0;
00975           q(3) = a1;
00976           if (naux>2) q(2) = a2;
00977           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00978         }
00979         break;
00980       case base::Bottom:
00981         for (register int i=mxs; i<=mxe; i++) {
00982           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00983           if (naux==2){ 
00984             q(2) = a0;
00985             q(3) = a1;
00986           }
00987           else {
00988             q(1) = a0;
00989             q(3) = a1;
00990             q(2) = a2;
00991           }
00992           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00993         }
00994         break;
00995       case base::Top:
00996         for (register int i=mxs; i<=mxe; i++) {
00997           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00998           if (naux==2){
00999             q(2) = a0;
01000             q(3) = a1;
01001           }
01002           else {
01003             q(1) = a0;
01004             q(3) = a1;
01005             q(2) = a2;
01006           }
01007           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
01008         }
01009         break;
01010       }
01011     }
01012 
01013   }
01014 
01015   virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
01016                              const MicroType* f, const point_type* xc, const DataType* distance, 
01017                              const point_type* normal, DataType* aux=0, const int naux=0, 
01018                              const int scaling=0) const {    
01019     
01020     if (type==GFMNoSlipWallTemp || type==GFMSlipWallTemp) { 
01021       DataType fac = S0;
01022       if (scaling==base::Physical)
01023         fac /= U0;      
01024       for (int n=0; n<nc; n++) {
01025         MacroType q=MacroVariables(f[n]);
01026         DataType u=-q(1); 
01027         DataType v=-q(2);
01028         DataType T= q(3);
01029         
01030         // Add boundary velocities if available
01031         if (naux>=3) {
01032           u += fac*aux[naux*n]; 
01033           v += fac*aux[naux*n+1];
01034           T = (aux[naux*n+2]-Tpmin)/(Tpmax-Tpmin);
01035         } 
01036         if (type==GFMNoSlipWallTemp) {
01037           // Prescribe entire velocity vector in ghost cells
01038           q(1) += DataType(2.)*u;
01039           q(2) += DataType(2.)*v;
01040           q(3) = T;
01041         }
01042         else if (type==GFMSlipWallTemp) {
01043           // Prescribe only normal velocity vector in ghost cells
01044           DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
01045           q(1) += vl*normal[n](0);
01046           q(2) += vl*normal[n](1);
01047           q(3) = T;
01048         }
01049         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01050       } 
01051     }
01052 
01053     else if (type==GFMExtrapolation) {
01054       for (int n=0; n<nc; n++) {
01055         MacroType q=MacroVariables(f[n]);
01056         DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1);
01057         q(1) = vl*normal[n](0);
01058         q(2) = vl*normal[n](1);
01059         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01060       }
01061     }
01062 
01063     else if (type==GFMBounceBack) {
01064       DCoords dx = base::GH().worldStep(fvec.stepsize());
01065       DataType fac = S0;
01066       if(scaling==base::Physical) fac/= U0;
01067 
01068       DataType dt_temp = dx(0)*S0;
01069    
01070       //xc : Cell midpoints of the ghost cells
01071       //nc : Number of ghost cells
01072       //normal : normalized normal vectors, pointed inside
01073       //distance: Level-Set distance between ghost cell midpoint and boundary
01074        
01075       DataType distB[2], distBScalar, cosa, leng, d_dist, qNorm, mdx2, mdy2, zaehler, nenner, lengXi;
01076         
01077       for (int n=0; n<nc; n++){
01078         distB[0] = (double &)distance[n]*normal[n](0);//*dx(0);
01079         distB[1] = (double &)distance[n]*normal[n](1);//*dx(1);
01080 
01081         //std::ofstream deb("debugOut.txt", std::ios::app);
01082 
01083 
01084         //alpha: Angle between discrete velocity direction (Ghostcell) and inner normal direction (Boundary)
01085         //leng: Distance between ghost cell midpoint and curved boundary along discrete velocity direction
01086         //d_dist: Distance between fluid-cellmidpoint and boundary along discrete velocity direction
01087         //qNorm: Normalized length from fluid-cellmidpoint to boundary along discrete velocity direction
01088  
01089         for(int ind=0; ind<base::NMicroVar(); ind++){
01090                 
01091                 mdx2 = mdx[ind]*mdx[ind];
01092                 mdy2 = mdy[ind]*mdy[ind];       
01093                 distBScalar = std::sqrt(std::pow(distB[0],2)+std::pow(distB[1],2));
01094                 leng = distBScalar*distBScalar*std::sqrt(mdx2+mdy2)/(distB[0]*mdx[ind]+distB[1]*mdy[ind]);
01095                 lengXi = std::sqrt(mdx2*dx(0)*dx(0)+mdy2*dx(1)*dx(1));
01096                 d_dist = lengXi-leng;
01097                 qNorm = d_dist/lengXi;
01098                 
01099                 //deb << ind << " " << dx(0) << " " << dx(1) << " " << lengXi << " " << d_dist << " " << qNorm << " " << xc[n](0) << " " << xc[n](1) << std::endl;      
01100 
01101                 if(qNorm>0 && qNorm < 1.){ //std::sqrt(2.)){
01102                 //The positive lengths are important only       
01103                         
01104                         //Two cases: qNorm greater equal or lower than 0.5
01105                         int idxA, idxB;
01106                         if(ind==0)      idxA=0;
01107                         else if(ind==1) idxA=2;
01108                         else if(ind==2) idxA=1;
01109                         else if(ind==3) idxA=6;
01110                         else if(ind==4) idxA=8;
01111                         else if(ind==5) idxA=7;
01112                         else if(ind==6) idxA=3;
01113                         else if(ind==7) idxA=5;
01114                         else if(ind==8) idxA=4;
01115 
01116                         //Temperature distributions starts here
01117                         else if(ind==9) idxA=9;
01118                         else if(ind==10)idxA=11;
01119                         else if(ind==11)idxA=10;
01120                         else if(ind==12)idxA=13;
01121                         else            idxA=12;
01122 
01123                         idxB=idxA;
01124                         if(qNorm>=DataType(0.5)) idxB=ind;
01125 
01126                         //Calculate Coordinates of needed fluidcell
01127                         //Coords cn,cnA,cnB;
01128                         //cn = Coords(base::Dim(),idx+base::Dim()*n); //Coords of ghostcell
01129                         
01130                         DataType cn[2], cnA[2], cnB[2];
01131                         cn[0] = xc[n](0);
01132                         cn[1] = xc[n](1);
01133 
01134                         //cnA = cn; //Coords of fluidcell
01135                         cnA[0] = cn[0]+mdx[ind]*dx(0);
01136                         cnA[1] = cn[1]+mdy[ind]*dx(1);
01137 
01138 
01139                         //cnB = cnA; //Coords of 2nd fluidcell
01140                         cnB[0] = cnA[0];
01141                         cnB[1] = cnA[1];
01142 
01143                         if(qNorm<DataType(0.5) && qNorm > DataType(0.)){
01144                                 cnB[0] += mdx[ind]*dx(0);
01145                                 cnB[1] += mdy[ind]*dx(1);
01146                         }
01147 
01148                         //Getting the corresponding partial probability distribution functions
01149                         DataType distFuncA, distFuncB;
01150 
01151                         //distFuncA = fvec(cnA)[idxA];
01152                         //distFuncB = fvec(cnB)[idxB];
01153 
01154                         //std::cout << "cn = " << cn[0] << " " << cn[1] << ", cnA = " << cnA[0] << " " << cnB[1] << ", cnB = " << cnB[0] << " " << cnB[1] << ", dx= " << dx(0) << " " << dx(1) << std::endl;
01155         
01156         
01157                         distFuncA = fvec(base::GH().localCoords((const DataType*) cnA))(idxA);
01158                         distFuncB = fvec(base::GH().localCoords((const DataType*) cnB))(idxB);
01159 
01160                         DataType intp1,intp2,intp,q2;
01161                         q2 = DataType(2.)*qNorm;
01162 
01163                         //Using formular from Eitel-Amor for density distribution functions
01164                         if(ind<9){      
01165                                 if(qNorm<DataType(0.5) && qNorm >= DataType(0.)){
01166                                         intp1=q2*distFuncA;
01167                                         intp2=(DataType(1.)-q2)*distFuncB;
01168                                 }
01169                                 else{
01170                                         if(q2!=0){
01171                                                 intp1 = (DataType(1.)/q2)*distFuncA;
01172                                                 intp2 = ((q2-DataType(1.))/q2)*distFuncB;
01173                                         }
01174                                         else{
01175                                         intp1 = DataType(0.);
01176                                         intp2 = intp1;
01177                                         }
01178                                 }
01179                         }
01180 
01181                         //Using formular from L. Li et al. for temperature distribution functions
01182                         else{           
01183                                 if(qNorm<=DataType(0.5) && qNorm >= DataType(0.)){
01184                                         intp1 = -q2*distFuncA;
01185                                         intp2 = (q2-DataType(1.))*distFuncB;
01186                                 }
01187                                 else{
01188                                         if(q2!=0){
01189                                                 intp1 = -(DataType(1.)/q2)*distFuncA;
01190                                                 intp2 = (DataType(1.)-(DataType(1.)/q2))*distFuncB;
01191                                         }
01192                                         else{
01193                                                 intp1 = DataType(0.);
01194                                                 intp2 = intp1;
01195                                         }
01196                                 }
01197                         }
01198 
01199                         intp = intp1+intp2;
01200                         
01201                         //fvec(Coords(base::Dim(),idx+base::Dim()*n))[ind] = intp;
01202 
01203                         fvec(base::GH().localCoords((const DataType*) cn))(ind) = intp;
01204 
01205                         //ADDITION PART FOR EXTERNAL FORCE LIKE VELOCITY OR TEMPERATURE
01206                         //From: Momentum transfer of a Boltzmann-lattice fluid with boundaries - M. Bouzidi et al. (2001) Physics of Fluids
01207 
01208                         if(naux>=2 && ind<9){
01209                                 DataType u = fac*aux[naux*n];
01210                                 DataType v = fac*aux[naux*n+1];
01211         
01212                                 DataType weight;
01213                                 weight = DataType(1.)/DataType(3);
01214                                 if(ind==0) weight *= 0;
01215                                 else if(ind==4 || ind ==5 || ind==7 || ind==8) weight *= DataType(1.)/DataType(4.);
01216         
01217                                 DataType velForce;
01218                                         
01219                                 if(qNorm<DataType(0.5)) velForce = DataType(2.)*weight;
01220                                 else velForce = (DataType(1.)/qNorm)*weight;
01221                                         
01222                                 velForce *= ((mdx[ind]*u)+(mdy[ind]*v));
01223                                 //fvec(Coords(base::Dim(),idx+base::Dim()*n))[ind] +=velForce;
01224 
01225                                 fvec(base::GH().localCoords((const DataType*) cn))(ind) +=velForce;
01226                         }
01227 
01228                         //From: Boundary conditions for thermal lattice Boltzmann equation method - L. Li et al. (2013) Journal of Computational Physics 237
01229                         if(naux>=3 && ind>=9){
01230                                 DataType T = (aux[naux*n+2]-Tpmin)/(Tpmax-Tpmin);
01231                                 DataType TForce = DataType(1.)/DataType(2.);
01232 
01233                                 if(qNorm>DataType(0.5)) TForce /= q2;
01234 
01235                                 //fvec(Coords(base::Dim(),idx+base::Dim()*n))[ind] += TForce*T;         
01236                                 fvec(base::GH().localCoords((const DataType*) cn))(ind) +=TForce*T;
01237                         
01238                         }       
01239                 }       
01240         }
01241       }
01242     }
01243   }
01244     
01245   virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01246                       const int skip_ghosts=1) const {
01247     int Nx = fvec.extents(0), Ny = fvec.extents(1);
01248     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01249     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01250     MicroType *f = (MicroType *)fvec.databuffer();
01251     DataType *work = (DataType *)workvec.databuffer();
01252     DCoords dx = base::GH().worldStep(fvec.stepsize());
01253     DataType dt_temp = dx(0)*S0/U0;
01254 
01255 //    std::cout << "mxs = " << mxs << ", mys = " << mys << ", mxe = " << mxe << ", mye = " << mye << std::endl;
01256 
01257     if (cnt>=1 && cnt<=9) {
01258       for (register int j=mys; j<=mye; j++)
01259         for (register int i=mxs; i<=mxe; i++) {
01260           MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
01261           // 1:r, 2:u, 3:v, 4:E, 5:T, 6:p, 7:gamma, >=8, variable: \nu
01262           switch(cnt) {   
01263           case 1:
01264             work[base::idx(i,j,Nx)]=0;
01265             break;
01266           case 2:
01267             work[base::idx(i,j,Nx)]=q(1)*U0/S0;
01268             break;
01269           case 3:
01270             work[base::idx(i,j,Nx)]=q(2)*U0/S0;
01271             break;
01272           case 5:
01273             work[base::idx(i,j,Nx)]=q(3)*(Tpmax-Tpmin)+Tpmin;
01274             break;
01275           case 6:
01276             work[base::idx(i,j,Nx)]=q(0)*P0;
01277             break;
01278           case 8: {
01279             MicroType feq = Equilibrium(q);       
01280             work[base::idx(i,j,Nx)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt_temp),
01281                                                GasSpeedofSound(),dt_temp );
01282             break;
01283           }
01284           case 9:
01285             work[base::idx(i,j,Nx)]=std::sqrt(q(1)*q(1)+q(2)*q(2))*U0/S0;
01286             break;
01287           }       
01288         }
01289     }
01290     else if (cnt>=10 && cnt<=11) {
01291       macro_grid_data_type qgrid(fvec.bbox());
01292       MacroType *q = (MacroType *)qgrid.databuffer();
01293       for (register int j=mys; j<=mye; j++)
01294         for (register int i=mxs; i<=mxe; i++) {
01295           q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
01296           q[base::idx(i,j,Nx)](1)*=U0/S0;
01297           q[base::idx(i,j,Nx)](2)*=U0/S0;
01298         }
01299 
01300       if (skip_ghosts==0) {
01301         mxs+=1; mxe-=1; mys+=1; mye-=1;
01302       }
01303 
01304       for (register int j=mys; j<=mye; j++)
01305         for (register int i=mxs; i<=mxe; i++) {
01306           work[base::idx(i,j,Nx)]=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(2.*dx(0))-
01307             (q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(2.*dx(1));
01308           if (cnt==11) work[base::idx(i,j,Nx)]=std::abs(work[base::idx(i,j,Nx)]);
01309         }
01310     }
01311 
01312     else if (cnt>=90 && cnt <=99){
01313         switch(cnt){
01314         case 90:
01315         //std::cout << "mxs = " << mxs << ", mys = " << mys << ", mxe = " << mxe << ", mye = " << mye << std::endl;
01316         for (register int j=mys; j<=mye; j++)
01317           for (register int i=mxs; i<=mxe; i++) {
01318                 MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
01319                 MacroType q2=MacroVariables(f[base::idx(i+1,j,Nx)]);
01320                 work[base::idx(i,j,Nx)] = q2(3)-q(3);           
01321           }
01322         
01323         break;
01324         case 91:
01325         case 92:
01326         case 93:
01327         case 94:
01328         case 95:
01329         case 96:
01330         case 97:
01331         case 98:
01332         case 99:
01333         break;
01334         }
01335     }
01336   }
01337 
01338   virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01339                      const int skip_ghosts=1) const {
01340     int Nx = fvec.extents(0), Ny = fvec.extents(1);
01341     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01342     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01343     MicroType *f = (MicroType *)fvec.databuffer();
01344     DataType *work = (DataType *)workvec.databuffer();
01345 
01346     for (register int j=mys; j<=mye; j++)
01347       for (register int i=mxs; i<=mxe; i++) {
01348         switch(cnt) {
01349         case 2:
01350           f[base::idx(i,j,Nx)](1)=work[base::idx(i,j,Nx)]/(U0/S0);
01351           break;
01352         case 3:
01353           f[base::idx(i,j,Nx)](2)=work[base::idx(i,j,Nx)]/(U0/S0);
01354           f[base::idx(i,j,Nx)] = Equilibrium((const MacroType &)f[base::idx(i,j,Nx)]);
01355           break;
01356         case 5:
01357           f[base::idx(i,j,Nx)](3)=(work[base::idx(i,j,Nx)]-Tpmin)/(Tpmax-Tpmin);
01358           break;
01359         case 6:
01360           f[base::idx(i,j,Nx)](0)=work[base::idx(i,j,Nx)]/P0;
01361           break;
01362         }
01363       }
01364   }
01365 
01366   virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time, 
01367                     const int verbose) const {
01368     int Nx = fvec.extents(0);
01369     int b = fvec.bottom();
01370     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01371     BBox bbmax = grow(fvec.bbox(),-base::NGhosts());    
01372     int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01373     int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01374     int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01375     int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01376     MicroType *f = (MicroType *)fvec.databuffer();
01377     
01378     int result = 1;
01379     for (register int j=mys; j<=mye; j++)
01380       for (register int i=mxs; i<=mxe; i++)
01381         for (register int l=0; l<base::NMicroVar(); l++)          
01382           if (!(f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l)>-1.e37)) {
01383             result = 0;
01384             if (verbose) 
01385               std::cerr << "D2Q9Thermal-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << ")(" << l << ")=" 
01386                         << f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l) << " " << fvec.bbox() << std::endl;
01387           }
01388     return result;
01389   }  
01390         
01391   virtual int NMethodOrder() const { return 2; } 
01392 
01393   inline const DataType& L0() const { return base::LengthScale(); }
01394   inline const DataType& T0() const { return base::TimeScale(); }
01395   inline void SetDensityScale(const DataType r0) { R0 = r0; }
01396   inline void SetPressureScale(const DataType p0) { P0 = p0; }
01397   inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
01398   inline void SetSpeedUp(const DataType s0) { S0 = s0; }
01399   inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
01400   inline void SetTemperatureScale(const DataType tpmin, const DataType tpmax) { 
01401     Tpmin = tpmin; Tpmax = tpmax; 
01402     Temp0 = tpmax-tpmin;
01403     Tref = DataType(0.5)*(Tpmax+Tpmin)/Temp0;
01404   }
01405 
01406   inline const DataType& DensityScale() const { return R0; }
01407   inline const DataType& PressureScale() const { return P0; }
01408   inline const DataType VelocityScale() const { return U0/S0; }
01409   inline const DataType& TemperatureScale() const { return Temp0; }
01410 
01411   // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
01412   inline const DataType& SpeedUp() const { return S0; }
01413 
01414   // Lattice quantities for the operator
01415   inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
01416   inline DataType LatticeViscosityT(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*ds2); }
01417   inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
01418   inline DataType LatticeSpeedOfSoundT() const { return (std::sqrt(ds2)); }
01419 
01420   // Physical quantities for the operator
01421   inline void SetGas(const DataType pr, const DataType rho, const DataType nu, const DataType cs) { 
01422     rhop = rho; nup = nu; csp = cs; cs2p = cs*cs; prp = pr; 
01423     SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
01424     SetPressureScale(prp);
01425   }
01426   //inline void SetThermalGas(const DataType tmin, const DataType tmax, const DataType diff, const DataType force, const DataType width) {
01427   inline void SetThermalGas(const DataType tmin, const DataType tmax, const DataType diff, const DataType gp, const DataType betap) {        
01428     nutp = diff; 
01429     SetTemperatureScale(tmin,tmax);
01430 
01431     DataType g = gp/(L0()/T0()/T0()); //g in Lattice Coordinates
01432     DataType beta = betap*Temp0; //beta in Lattice Coordinates
01433     gbeta = g*beta;
01434 
01435     //Force input is the Rayleigh Number -> Calculating in Lattice directly
01436     //gbeta = force*LatticeViscosity(Omega(T0()))*LatticeViscosityT(OmegaT(T0()))/(H/L0())/(H/L0())/(H/L0());
01437 
01438   }
01439   
01440   inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
01441   inline virtual const DataType OmegaT(const DataType dt) const { return (dcs2*cs2p*dt/S0/(nutp*S0+dcs2*cs2p*dt/S0*DataType(0.5))); }
01442 
01443   inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q, 
01444                                               const DataType dt) const {
01445     DataType M2 = 0.0;
01446 
01447     M2 += (f(4) - feq(4))*(f(4) - feq(4));
01448     M2 += (f(5) - feq(5))*(f(5) - feq(5));
01449     M2 += (f(7) - feq(7))*(f(7) - feq(7));
01450     M2 += (f(8) - feq(8))*(f(8) - feq(8));
01451     M2 = std::sqrt(M2);
01452 
01453     DataType tau = DataType(1.0)/Omega(dt);
01454     DataType om_turb =  (DataType(2.0)/( tau
01455              + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
01456 
01457     return ( om_turb );
01458   }
01459   inline const int TurbulenceType() const { return turbulence; }
01460   inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
01461 
01462   inline const DataType& GasDensity() const { return rhop; }
01463   inline const DataType& GasPressure() const { return prp; }
01464   inline const DataType& GasViscosity() const { return nup; }
01465   inline const DataType& GasViscosityT() const { return nutp; }
01466   inline const DataType& GasSpeedofSound() const { return csp; }
01467   inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
01468   { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
01469   inline const DataType GasViscosityT(const DataType omegat, const DataType cs, const DataType dt) const
01470   { return ((DataType(1.)/omegat-DataType(0.5))*dcs2*cs*cs*dt/S0/S0); }
01471 
01472 protected:
01473   DataType cs2, cs22, cssq, ds2, dcs2;
01474   DataType P0, R0, U0, S0, Temp0, rhop, Tpmin, Tpmax, Tref, csp, cs2p, nup, nutp, prp, gbeta;
01475   DataType Cs_Smagorinsky, turbulence;
01476   int method[1], mdx[14], mdy[14];
01477 };
01478 
01479 
01480 #endif
01481