amroc/lbm/LBMD2Q9Thermal.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 
00044 template <class DataType>
00045 class LBMD2Q9Thermal : public LBMBase<Vector<DataType,14>, Vector<DataType,4>, 2> {
00046   typedef LBMBase<Vector<DataType,14>, Vector<DataType,4>, 2> base;
00047 public:
00048   typedef typename base::vec_grid_data_type vec_grid_data_type;
00049   typedef typename base::grid_data_type grid_data_type;
00050   typedef typename base::MicroType MicroType;
00051   typedef typename base::MacroType MacroType;
00052   typedef GridData<MacroType,2> macro_grid_data_type;
00053   typedef typename base::SideName SideName;
00054   typedef typename base::point_type point_type;
00055 
00056   enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro }; 
00057   enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall, ExtrapolationEquilibrium, SlipWallTemperature, NoSlipWallTemperature,InletTemp }; 
00058   enum GFMPredefined { GFMExtrapolation, GFMSlipWallTemp, GFMNoSlipWallTemp };
00059   enum TurbulenceModel { laminar, LES_Smagorinsky };
00060 
00061   LBMD2Q9Thermal() : base(), P0(1.), R0(1.), U0(1.), S0(1.), Temp0(1.), rhop(0.), Tpmin(0.), Tpmax(0.), Tref(0.5), csp(0.), 
00062                      cs2p(0.), nup(0.), nutp(0.) , prp(0.) {
00063     cs2  = DataType(1.)/DataType(3.);
00064     cs22 = DataType(2.)*cs2;
00065     cssq = DataType(2.)/DataType(9.);
00066     ds2  = DataType(1.)/DataType(2.);
00067     dcs2 = ds2/cs2;
00068     method[0] = 2;
00069     turbulence = laminar;
00070 
00071     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;
00072     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;
00073     
00074   }
00075  
00076   virtual ~LBMD2Q9Thermal() {}
00077 
00078   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00079     base::LocCtrl = Ctrl.getSubDevice(prefix+"D2Q9");
00080     RegisterAt(base::LocCtrl,"Method(1)",method[0]); 
00081     RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00082     RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00083     RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00084     RegisterAt(base::LocCtrl,"Gas_Cs",csp); 
00085     RegisterAt(base::LocCtrl,"Gas_nu",nup); 
00086     RegisterAt(base::LocCtrl,"Gas_nuT",nutp);
00087     RegisterAt(base::LocCtrl,"Gas_pr",prp);  
00088     RegisterAt(base::LocCtrl,"Gas_Tmin",Tpmin); 
00089     RegisterAt(base::LocCtrl,"Gas_Tmax",Tpmax); 
00090     RegisterAt(base::LocCtrl,"SpeedUp",S0);  
00091   }
00092 
00093   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
00094     base::SetupData(gh,ghosts);
00095     if (rhop>0.) R0 = rhop;
00096     if (prp>0.) P0 = prp;
00097     if (Tpmin!=0. || Tpmax!=1.) SetTemperatureScale(Tpmin,Tpmax);
00098     if (csp>0.) {
00099       cs2p = csp*csp;     
00100       SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00101     }   
00102     std::cout << "D2Q9Thermal: Gas_rho=" << rhop << "   Gas_Cs=" << csp 
00103               << "   Gas_nu=" << nup << "   Gas_nut=" << nutp << "   Gas_pr=" << prp << std::endl; 
00104     std::cout << "D2Q9Thermal: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0 
00105               << "   R0=" << R0 << "   Temp0=" << Temp0 << "   S0=" << S0 << "   P0=" << P0 << std::endl;
00106 
00107   }
00108 
00109   inline virtual MacroType MacroVariables(const MicroType &f) const {
00110     MacroType q;
00111     DataType part1, part2;
00112 
00113     q(1)=f(1)-f(2)+f(4)-f(5)+f(7)-f(8);
00114     q(2)=f(3)+f(4)+f(5)-f(6)-f(7)-f(8);
00115 
00116     part1 = f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8);
00117     part2 = cs22*(std::pow(q(1),2)+std::pow(q(2),2));
00118     q(0)= (DataType(3.)/DataType(5.))*(part1 - part2);
00119 
00120     q(3)=f(10)+f(11)+f(12)+f(13);
00121 
00122     return q;
00123   }
00124 
00125   inline virtual MicroType Equilibrium(const MacroType &q) const {
00126     MicroType feq;
00127 
00128     DataType p = q(0);
00129     DataType u = q(1);
00130     DataType v = q(2);
00131     DataType T = q(3);
00132 
00133     DataType ri, rt0, rt1, rt2, pre0, pre1, pre2, Temp;
00134 
00135     ri = DataType(0.);
00136     rt0 = DataType(4.) / DataType(9.);
00137     rt1 = DataType(1.) / DataType(9.);
00138     rt2 = DataType(1.) / DataType(36.);
00139 
00140     pre0 = -(DataType(5.) / DataType(3.)) * p;
00141     pre1 = (DataType(1.) / DataType(3.)) * p;
00142     pre2 = (DataType(1.) / DataType(12.)) * p;
00143 
00144     Temp = T/DataType(4.);
00145                
00146     DataType usq = u * u; 
00147     DataType vsq = v * v;
00148     DataType sumsq = (usq + vsq) / cs22;
00149     DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00150     DataType u2 = usq / cssq; 
00151     DataType v2 = vsq / cssq;
00152     DataType ui = u / cs2;
00153     DataType vi = v / cs2;
00154     DataType uv = ui * vi;
00155     DataType uT = DataType(2.)*u;
00156     DataType vT = DataType(2.)*v;
00157 
00158     feq(0) = pre0 + rt0 * (ri - sumsq);
00159                
00160     feq(1) = pre1 + rt1 * (ri - sumsq + u2 + ui);
00161     feq(2) = pre1 + rt1 * (ri - sumsq + u2 - ui);
00162     feq(3) = pre1 + rt1 * (ri - sumsq + v2 + vi);
00163     feq(6) = pre1 + rt1 * (ri - sumsq + v2 - vi);
00164                
00165     feq(4) = pre2 + rt2 * (ri + sumsq2 + ui + vi + uv);
00166     feq(5) = pre2 + rt2 * (ri + sumsq2 - ui + vi - uv);
00167     feq(7) = pre2 + rt2 * (ri + sumsq2 + ui - vi - uv);
00168     feq(8) = pre2 + rt2 * (ri + sumsq2 - ui - vi + uv);   
00169 
00170     feq(9) = DataType(0.);
00171     feq(10) = Temp*(DataType(1.)+uT);
00172     feq(11) = Temp*(DataType(1.)-uT);
00173     feq(12) = Temp*(DataType(1.)+vT);
00174     feq(13) = Temp*(DataType(1.)-vT);
00175 
00176     return feq;
00177   }
00178 
00179   inline virtual void Collision(MicroType &f, const DataType dt) const {
00180     MacroType q = MacroVariables(f);
00181     MicroType feq = Equilibrium(q);
00182 
00183     DataType force = DataType(0.5)*dt/S0*gbeta*(q(3)-Tref);
00184 
00185     DataType omega = Omega(dt), omegaT = OmegaT(dt);
00186     for (int i=0; i<=8;i++){
00187         f(i)=(DataType(1.)-omega)*f(i) + omega*feq(i);
00188         if(i==3) f(i)+= force;
00189         if(i==6) f(i)-= force;
00190     }
00191     for (int i=9; i<=13;i++){
00192         f(i)=(DataType(1.)-omegaT)*f(i) + omegaT*feq(i);
00193     }
00194    
00195   }     
00196 
00197   virtual int IncomingIndices(const int side, int indices[]) const {
00198     switch (side) {
00199     case base::Left:
00200       indices[0] = 1; indices[1] = 4; indices[2] = 7; indices[3] = 10;
00201       break;
00202     case base::Right:
00203       indices[0] = 2; indices[1] = 5; indices[2] = 8; indices[3] = 11;
00204       break;
00205     case base::Bottom:
00206       indices[0] = 3; indices[1] = 4; indices[2] = 5; indices[3] = 12; 
00207       break;
00208     case base::Top:
00209       indices[0] = 6; indices[1] = 7; indices[2] = 8; indices[3] = 13; 
00210       break;
00211     }
00212     return 4;
00213   }
00214 
00215   virtual int OutgoingIndices(const int side, int indices[]) const {
00216     switch (side) {
00217     case base::Left:
00218       indices[0] = 2; indices[1] = 5; indices[2] = 8; indices[3] = 11; 
00219       break;
00220     case base::Right:
00221       indices[0] = 1; indices[1] = 4; indices[2] = 7; indices[3] = 10; 
00222       break;
00223     case base::Bottom:
00224       indices[0] = 6; indices[1] = 7; indices[2] = 8; indices[3] = 13; 
00225       break;
00226     case base::Top:
00227       indices[0] = 3; indices[1] = 4; indices[2] = 5; indices[3] = 12; 
00228       break;
00229     }
00230     return 4;
00231   }
00232 
00233   // Revert just one layer of cells at the boundary 
00234   virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
00235                              const int side) const {
00236     int Nx = fvec.extents(0);
00237     int b = fvec.bottom();
00238     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00239     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00240     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00241     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00242     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00243     MicroType *f = (MicroType *)fvec.databuffer();
00244 
00245 #ifdef DEBUG_PRINT_FIXUP
00246     ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
00247       << fvec.bbox() << " using " << bb << "\n").flush();
00248 #endif
00249     
00250     switch (side) {
00251     case base::Left:
00252       for (register int j=mys; j<=mye; j++) {
00253         f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](2);
00254         f[b+base::idx(mxs,j,Nx)](11)= f[b+base::idx(mxs-1,j,Nx)](11);
00255         f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](8);
00256         f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](5);
00257       }
00258       break;
00259     case base::Right:
00260       for (register int j=mys; j<=mye; j++) {
00261         f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](1);
00262         f[b+base::idx(mxe,j,Nx)](10)= f[b+base::idx(mxe+1,j,Nx)](10);
00263         f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](7);
00264         f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](4);
00265       }
00266       break;
00267     case base::Bottom:
00268       for (register int i=mxs; i<=mxe; i++) {
00269         f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](6);
00270         f[b+base::idx(i,mys,Nx)](13)= f[b+base::idx(i,mys-1,Nx)](13);
00271         f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](8);
00272         f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](7);
00273       }
00274       break;
00275     case base::Top:
00276       for (register int i=mxs; i<=mxe; i++) {
00277         f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](3);
00278         f[b+base::idx(i,mye,Nx)](12)= f[b+base::idx(i,mye+1,Nx)](12);
00279         f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](5);
00280         f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](4);
00281       }
00282       break;
00283     }
00284   }   
00285 
00286   virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec, 
00287                          const BBox &bb, const double& dt) const {
00288     int Nx = fvec.extents(0);
00289     int b = fvec.bottom();
00290     MicroType *f = (MicroType *)fvec.databuffer();
00291     MicroType *of = (MicroType *)ovec.databuffer();
00292 
00293 #ifdef DEBUG_PRINT_FIXUP
00294     ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox() 
00295       << " using " << bb << "\n").flush();
00296 #endif
00297     
00298     assert (fvec.extents(0)==ovec.extents(0));
00299     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) &&
00300             fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1));
00301     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00302     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00303     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00304     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00305 
00306     for (register int j=mys; j<=mye; j++) 
00307       for (register int i=mxs; i<=mxe; i++) {
00308         for (register int n=0; n<base::NMicroVar(); n++) 
00309           f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00310         Collision(f[b+base::idx(i,j,Nx)],dt);
00311       }
00312   }  
00313 
00314   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00315                       const double& t, const double& dt, const int& mpass) const {
00316     // Make sure that the physical length and times are scaling consistently
00317     // into lattice quantities
00318     DCoords dx = base::GH().worldStep(fvec.stepsize());
00319 
00320 
00321     if (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR) {
00322       std::cerr << "LBMD2Q9 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00323                 << " to be equal." << std::endl << std::flush;
00324       std::exit(-1);
00325     }
00326     if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
00327       std::cerr << "LBMD2Q9 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0() 
00328                 << " to be equal." 
00329                 << "   dt=" << dt << "   TO=" << T0()
00330                 << std::endl << std::flush;
00331       std::exit(-1);
00332     }
00333 
00334     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00335     MicroType *f = (MicroType *)fvec.databuffer();
00336     
00337 
00338     // Streaming - update all cells
00339     for (register int j=Ny-1; j>=1; j--)
00340       for (register int i=0; i<=Nx-2; i++) {
00341         f[base::idx(i,j,Nx)](3) = f[base::idx(i,j-1,Nx)](3);
00342         f[base::idx(i,j,Nx)](12) = f[base::idx(i,j-1,Nx)](12);
00343         f[base::idx(i,j,Nx)](5) = f[base::idx(i+1,j-1,Nx)](5);  
00344       }
00345 
00346     for (register int j=Ny-1; j>=1; j--)
00347       for (register int i=Nx-1; i>=1; i--) {
00348         f[base::idx(i,j,Nx)](1) = f[base::idx(i-1,j,Nx)](1);
00349         f[base::idx(i,j,Nx)](10) = f[base::idx(i-1,j,Nx)](10);
00350         f[base::idx(i,j,Nx)](4) = f[base::idx(i-1,j-1,Nx)](4);
00351       }
00352 
00353     for (register int j=0; j<=Ny-2; j++)
00354       for (register int i=Nx-1; i>=1; i--) {
00355         f[base::idx(i,j,Nx)](6) = f[base::idx(i,j+1,Nx)](6);
00356         f[base::idx(i,j,Nx)](13) = f[base::idx(i,j+1,Nx)](13);
00357         f[base::idx(i,j,Nx)](7) = f[base::idx(i-1,j+1,Nx)](7);
00358       }
00359 
00360     for (register int j=0; j<=Ny-2; j++)
00361       for (register int i=0; i<=Nx-2; i++) {
00362         f[base::idx(i,j,Nx)](2) = f[base::idx(i+1,j,Nx)](2);
00363         f[base::idx(i,j,Nx)](11) = f[base::idx(i+1,j,Nx)](11);
00364         f[base::idx(i,j,Nx)](8) = f[base::idx(i+1,j+1,Nx)](8);
00365       }    
00366 
00367     // Collision - only in the interior region
00368     int mxs = base::NGhosts(), mys = base::NGhosts();
00369     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00370     for (register int j=mys; j<=mye; j++)
00371       for (register int i=mxs; i<=mxe; i++)
00372         Collision(f[base::idx(i,j,Nx)],dt);
00373     
00374     return DataType(1.);
00375   }
00376 
00377    virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00378                            DataType* aux=0, const int naux=0, const int scaling=0) const {
00379     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00380     int mxs = base::NGhosts(), mys = base::NGhosts();
00381     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00382     MicroType *f = (MicroType *)fvec.databuffer();
00383 
00384     if (type==ConstantMicro) {
00385       MicroType g;
00386       for (register int i=0; i<base::NMicroVar(); i++)
00387         g(i) = aux[i];
00388     
00389       for (register int j=mys; j<=mye; j++)
00390         for (register int i=mxs; i<=mxe; i++)
00391           f[base::idx(i,j,Nx)] = g;
00392     }
00393     
00394     else if (type==ConstantMacro) {
00395       MacroType q;
00396       if (scaling==base::Physical) {
00397         q(0) = aux[0]/P0; 
00398         q(1) = aux[1]/U0; 
00399         q(2) = aux[2]/U0;
00400         q(3) = (aux[3]-Tpmin)/Temp0;
00401       }
00402       else 
00403         for (register int i=0; i<base::NMacroVar(); i++)
00404           q(i) = aux[i];
00405       q(1) *= S0; q(2) *= S0;
00406     
00407       for (register int j=mys; j<=mye; j++)
00408         for (register int i=mxs; i<=mxe; i++)
00409           f[base::idx(i,j,Nx)] = Equilibrium(q);
00410     }
00411     
00412     else if (type==GasAtRest) {
00413       MacroType q;
00414       q(0) = GasDensity()/R0; 
00415       q(1) = DataType(0.); 
00416       q(2) = DataType(0.);
00417       for (register int j=mys; j<=mye; j++)
00418         for (register int i=mxs; i<=mxe; i++)
00419           f[base::idx(i,j,Nx)] = Equilibrium(q);
00420     }       
00421   }
00422 
00423   // Set just the one layer of ghost cells at the boundary 
00424   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00425                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00426     int Nx = fvec.extents(0);
00427     int b = fvec.bottom();
00428 
00429     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00430     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00431     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00432     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00433     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00434     MicroType *f = (MicroType *)fvec.databuffer();
00435  
00436     if (type==Symmetry) {
00437       switch (side) {
00438       case base::Left:
00439         for (register int j=mys; j<=mye; j++) {
00440           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j,Nx)](8);
00441           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00442           f[b+base::idx(mxe,j,Nx)](10) = f[b+base::idx(mxe+1,j,Nx)](11);
00443           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j,Nx)](5);
00444         }
00445         break;
00446       case base::Right:
00447         for (register int j=mys; j<=mye; j++) {
00448           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j,Nx)](7);
00449           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00450           f[b+base::idx(mxs,j,Nx)](11) = f[b+base::idx(mxs-1,j,Nx)](10);
00451           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j,Nx)](4);
00452         }
00453         break;
00454       case base::Bottom:
00455         for (register int i=mxs; i<=mxe; i++) {
00456           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i,mye+1,Nx)](8);
00457           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00458           f[b+base::idx(i,mye,Nx)](12) = f[b+base::idx(i,mye+1,Nx)](13);
00459           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i,mye+1,Nx)](7);
00460         }
00461         break;
00462       case base::Top:
00463         for (register int i=mxs; i<=mxe; i++) {
00464           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i,mys-1,Nx)](5);
00465           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00466           f[b+base::idx(i,mys,Nx)](13) = f[b+base::idx(i,mys-1,Nx)](12);
00467           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i,mys-1,Nx)](4);
00468         }
00469         break;
00470       }
00471     }
00472     
00473     else if (type==SlipWall) {
00474       switch (side) {
00475       case base::Left:
00476         for (register int j=mys; j<=mye; j++) {
00477           if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5); 
00478           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00479           f[b+base::idx(mxe,j,Nx)](10) = f[b+base::idx(mxe+1,j,Nx)](11);
00480           if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00481         }
00482         break;
00483       case base::Right:
00484         for (register int j=mys; j<=mye; j++) {
00485           if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4); 
00486           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00487           f[b+base::idx(mxs,j,Nx)](11) = f[b+base::idx(mxs-1,j,Nx)](10);
00488           if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00489         }
00490         break;
00491       case base::Bottom:
00492         for (register int i=mxs; i<=mxe; i++) {
00493           if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00494           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00495           f[b+base::idx(i,mye,Nx)](12) = f[b+base::idx(i,mye+1,Nx)](13);
00496           if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00497         }
00498         break;
00499       case base::Top:
00500         for (register int i=mxs; i<=mxe; i++) {
00501           if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00502           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00503           f[b+base::idx(i,mys,Nx)](13) = f[b+base::idx(i,mys-1,Nx)](12);
00504           if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00505         }
00506         break;
00507       }
00508     }
00509    
00510     else if (type==NoSlipWall) {
00511       switch (side) {
00512       case base::Left:
00513         for (register int j=mys; j<=mye; j++) {
00514           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00515           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00516           f[b+base::idx(mxe,j,Nx)](10) = f[b+base::idx(mxe+1,j,Nx)](11);
00517           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00518         }
00519         break;
00520       case base::Right:
00521         for (register int j=mys; j<=mye; j++) {
00522           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00523           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00524           f[b+base::idx(mxs,j,Nx)](11) = f[b+base::idx(mxs-1,j,Nx)](10);
00525           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00526         }
00527         break;  
00528       case base::Bottom:
00529         for (register int i=mxs; i<=mxe; i++) {
00530           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00531           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00532           f[b+base::idx(i,mye,Nx)](12) = f[b+base::idx(i,mye+1,Nx)](13);
00533           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00534         }
00535         break;  
00536       case base::Top:
00537         for (register int i=mxs; i<=mxe; i++) {
00538           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00539           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00540           f[b+base::idx(i,mys,Nx)](13) = f[b+base::idx(i,mys-1,Nx)](12);
00541           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00542         }
00543         break;
00544       }
00545     }
00546     // 
00547     else if (type==Inlet) {
00548       if (aux==0 || naux<=0) return;
00549       DataType a0=S0*aux[0], a1=0.;
00550       if (naux>1) a1 = S0*aux[1];
00551       if (scaling==base::Physical) {
00552         a0 /= U0;
00553         a1 /= U0;
00554       }
00555       switch (side) {
00556       case base::Left:
00557         for (register int j=mys; j<=mye; j++) {
00558           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00559           q(1) = a0;
00560           if (naux>1) q(2) = a1;
00561           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00562         }
00563         break;
00564       case base::Right:
00565         for (register int j=mys; j<=mye; j++) {
00566           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00567           q(1) = a0;
00568           if (naux>1) q(2) = a1;
00569           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00570         }
00571         break;
00572       case base::Bottom:
00573         for (register int i=mxs; i<=mxe; i++) {
00574           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00575           if (naux==1) 
00576             q(2) = a0;
00577           else {
00578             q(1) = a0;
00579             q(2) = a1;
00580           }
00581           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00582         }
00583         break;
00584       case base::Top:
00585         for (register int i=mxs; i<=mxe; i++) {
00586           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00587           if (naux==1)
00588             q(2) = a0;
00589           else {
00590             q(1) = a0;
00591             q(2) = a1;
00592           }
00593           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00594         }
00595         break;
00596       }
00597     }
00598     // Without temperature model
00599     else if (type==Outlet) {
00600       switch (side) {
00601       case base::Left:
00602         for (register int j=mys; j<=mye; j++) { 
00603                 f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)];
00604         }
00605         break;
00606       case base::Right:
00607         for (register int j=mys; j<=mye; j++) {
00608                 f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)];
00609 
00610                 //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);
00611                 //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);
00612                 //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);
00613                 //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);                
00614           }
00615         break;
00616       case base::Bottom:
00617         for (register int i=mxs; i<=mxe; i++) {
00618                 f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)];
00619         }
00620         break;
00621       case base::Top:
00622         for (register int i=mxs; i<=mxe; i++) {
00623                 f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)];
00624         }
00625         break;
00626       }
00627     }
00628     
00629     // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)    
00630     // Without temperature model
00631     else if (type==Pressure) {
00632       if (aux==0 || naux<=0) return;
00633       DataType a0=aux[0];
00634       if (scaling==base::Physical)
00635         a0 /= R0;
00636       switch (side) {
00637       case base::Left:
00638         for (register int j=mys; j<=mye; j++) {
00639           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00640           q(0) = a0;
00641           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00642         }
00643         break;
00644       case base::Right:
00645         for (register int j=mys; j<=mye; j++) {
00646           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00647           q(0) = a0;
00648           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00649         }
00650         break;
00651       case base::Bottom:
00652         for (register int i=mxs; i<=mxe; i++) {
00653           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00654           q(0) = a0;
00655           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00656         }
00657         break;
00658       case base::Top:
00659         for (register int i=mxs; i<=mxe; i++) {
00660           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00661           q(0) = a0;
00662           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00663         }
00664         break;
00665       }
00666     }
00667     //Without temperature model
00668     else if (type==SlidingWall) {
00669       if (aux==0 || naux<=0) return;
00670       if (aux==0 || naux<=0) return;
00671       DataType a0=S0*aux[0];
00672       if (scaling==base::Physical)
00673         a0 /= U0;
00674       switch (side) {
00675       case base::Left:
00676         for (register int j=mys; j<=mye; j++) {
00677           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00678           DataType rd1 = f[b+base::idx(mxe+1,j,Nx)](5), rd2 = f[b+base::idx(mxe+1,j,Nx)](8);
00679           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00680           DataType pw = DataType(1.)-qw;
00681           if (j<mye) f[b+base::idx(mxe,j+1,Nx)](7) = pw*rd1+qw*rd2;
00682           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00683           if (j>mys) f[b+base::idx(mxe,j-1,Nx)](4) = qw*rd1+pw*rd2;
00684         }
00685         break;
00686       case base::Right:
00687         for (register int j=mys; j<=mye; j++) {
00688           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00689           DataType rd1 = f[b+base::idx(mxs-1,j,Nx)](4), rd2 = f[b+base::idx(mxs-1,j,Nx)](7);
00690           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00691           DataType pw = DataType(1.)-qw;
00692           if (j<mye) f[b+base::idx(mxs,j+1,Nx)](8) = pw*rd1+qw*rd2;
00693           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00694           if (j>mys) f[b+base::idx(mxs,j-1,Nx)](5) = qw*rd1+pw*rd2;
00695         }
00696         break;
00697       case base::Bottom:
00698         for (register int i=mxs; i<=mxe; i++) {
00699           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00700           DataType rd1 = f[b+base::idx(i,mye+1,Nx)](7), rd2 = f[b+base::idx(i,mye+1,Nx)](8);
00701           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00702           DataType pw = DataType(1.)-qw;
00703           if (i<mxe) f[b+base::idx(i+1,mye,Nx)](5) = pw*rd1+qw*rd2;
00704           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00705           if (i>mxs) f[b+base::idx(i-1,mye,Nx)](4) = qw*rd1+pw*rd2;
00706         }
00707         break;
00708       case base::Top:
00709         for (register int i=mxs; i<=mxe; i++) {
00710           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00711           DataType rd1 = f[b+base::idx(i,mys-1,Nx)](4), rd2 = f[b+base::idx(i,mys-1,Nx)](5);
00712           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00713           DataType pw = DataType(1.)-qw;
00714           if (i<mxe) f[b+base::idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00715           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00716           if (i>mxs) f[b+base::idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00717         }
00718         break;
00719       }
00720     }
00721 
00722     // (Z.Guo, B. Shi and C.Zheng - A coupled lattice BGK model for the boussinesq equations)    
00723     else if (type==ExtrapolationEquilibrium) {
00724       switch (side) {
00725       case base::Left:
00726         for (register int j=mys; j<=mye; j++) {
00727           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00728           MacroType qBC = q;
00729           if (aux==0 || naux<=0) return;
00730           for(int k=0; k<=3; k++){
00731                 if(aux[k] != 0){
00732                         if(k==0) qBC(k) = 2.*aux[k+4]/P0-q(k);
00733                         //else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);
00734                         //else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00735                         else if(k>0 && k<3) qBC(k) = 2.*aux[k+4]/(U0/S0)-q(k);
00736                         else qBC(k)=2.*((aux[k+4]-Tpmin)/Temp0)-q(k);
00737                 }
00738           }
00739           f[b+base::idx(mxe,j,Nx)] = Equilibrium(qBC)+f[b+base::idx(mxe+1,j,Nx)]-Equilibrium(q);
00740         }
00741         break;
00742 
00743      case base::Right:
00744         for (register int j=mys; j<=mye; j++) {
00745           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00746           MacroType qBC = q;
00747           if (aux==0 || naux<=0) return;
00748           for(int k=0; k<=3; k++){
00749                 if(aux[k] != 0){
00750                         if(k==0) qBC(k) = 2.*aux[k+4]/P0-q(k);
00751                         //else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);        
00752                         //else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00753                         else if(k>0 && k<3) qBC(k) = 2.*aux[k+4]/(U0/S0)-q(k);
00754                         else qBC(k)=2.*((aux[k+4]-Tpmin)/Temp0)-q(k);
00755                 }
00756           }
00757           f[b+base::idx(mxs,j,Nx)] = Equilibrium(qBC)+f[b+base::idx(mxs-1,j,Nx)]-Equilibrium(q);
00758         }
00759         break;
00760 
00761      case base::Bottom:
00762         for (register int i=mxs; i<=mxe; i++) {
00763           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00764           MacroType qBC = q;
00765           if (aux==0 || naux<=0) return;
00766           for(int k=0; k<=3; k++){
00767                 if(aux[k] != 0){
00768                         if(k==0) qBC(k) = 2.*aux[k+4]/P0-q(k);
00769                         //else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);        
00770                         //else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00771                         else if(k>0 && k<3) qBC(k) = 2.*aux[k+4]/(U0/S0)-q(k);
00772                         else qBC(k)=2.*((aux[k+4]-Tpmin)/Temp0)-q(k);
00773                 }
00774           }
00775           f[b+base::idx(i,mye,Nx)] = Equilibrium(qBC)+f[b+base::idx(i,mye+1,Nx)]-Equilibrium(q);
00776         }
00777         break;
00778 
00779       case base::Top:
00780         for (register int i=mxs; i<=mxe; i++) {
00781           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00782           MacroType qBC = q;
00783           if (aux==0 || naux<=0) return;
00784           for(int k=0; k<=3; k++){
00785                 if(aux[k] != 0){
00786                         if(k==0) qBC(k) = 2.*aux[k+4]/P0-q(k);
00787                         //else if(k>0 && k<3) qBC(k) = aux[k+4]/(U0/S0);        
00788                         //else qBC(k)=(aux[k+4]-Tpmin)/Temp0;
00789                         else if(k>0 && k<3) qBC(k) = 2.*aux[k+4]/(U0/S0)-q(k);
00790                         else qBC(k)=2.*((aux[k+4]-Tpmin)/Temp0)-q(k);
00791                 }
00792           }
00793          f[b+base::idx(i,mys,Nx)] = Equilibrium(qBC)+f[b+base::idx(i,mys-1,Nx)]-Equilibrium(q);
00794         }
00795         break;
00796 
00797       }
00798     }
00799 
00800 else if (type==SlipWallTemperature) {
00801       switch (side) {
00802       case base::Left:
00803         for (register int j=mys; j<=mye; j++) {
00804           MacroType qBC;
00805           qBC(3) = (aux[0]-Tpmin)/Temp0;
00806           MicroType gBC = Equilibrium(qBC);
00807           if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5); 
00808           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00809           f[b+base::idx(mxe,j,Nx)](10) = gBC(10);
00810           if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00811         }
00812         break;
00813       case base::Right:
00814         for (register int j=mys; j<=mye; j++) {
00815           MacroType qBC;
00816           qBC(3) = (aux[0]-Tpmin)/Temp0;
00817           MicroType gBC = Equilibrium(qBC);
00818           if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4); 
00819           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00820           f[b+base::idx(mxs,j,Nx)](11) = gBC(11);
00821           if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00822         }
00823         break;
00824       case base::Bottom:
00825         for (register int i=mxs; i<=mxe; i++) {
00826           MacroType qBC;
00827           qBC(3) = (aux[0]-Tpmin)/Temp0;
00828           MicroType gBC = Equilibrium(qBC);
00829           if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00830           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00831           f[b+base::idx(i,mye,Nx)](12) = gBC(12);
00832           if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00833         }
00834         break;
00835       case base::Top:
00836         for (register int i=mxs; i<=mxe; i++) {
00837           MacroType qBC;
00838           qBC(3) = (aux[0]-Tpmin)/Temp0;
00839           MicroType gBC = Equilibrium(qBC);
00840           if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00841           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00842           f[b+base::idx(i,mys,Nx)](13) = gBC(13);
00843           if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00844         }
00845         break;
00846       }
00847     }
00848    
00849     else if (type==NoSlipWallTemperature) {
00850       switch (side) {
00851       case base::Left:
00852         for (register int j=mys; j<=mye; j++) {
00853           MacroType qBC;
00854           qBC(3) = (aux[0]-Tpmin)/Temp0;
00855           MicroType gBC = Equilibrium(qBC);
00856           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00857           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00858           f[b+base::idx(mxe,j,Nx)](10) = gBC(10);
00859           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00860         }
00861         break;
00862       case base::Right:
00863         for (register int j=mys; j<=mye; j++) {
00864           MacroType qBC;
00865           qBC(3) = (aux[0]-Tpmin)/Temp0;
00866           MicroType gBC = Equilibrium(qBC);
00867           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00868           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00869           f[b+base::idx(mxs,j,Nx)](11) = gBC(11);
00870           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00871         }
00872         break;  
00873       case base::Bottom:
00874         for (register int i=mxs; i<=mxe; i++) {
00875           MacroType qBC;
00876           qBC(3) = (aux[0]-Tpmin)/Temp0;
00877           MicroType gBC = Equilibrium(qBC);
00878           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00879           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00880           f[b+base::idx(i,mye,Nx)](12) = gBC(12);
00881           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00882         }
00883         break;  
00884       case base::Top:
00885         for (register int i=mxs; i<=mxe; i++) {
00886           MacroType qBC;
00887           qBC(3) = (aux[0]-Tpmin)/Temp0;
00888           MicroType gBC = Equilibrium(qBC);
00889           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00890           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00891           f[b+base::idx(i,mys,Nx)](13) = gBC(13);
00892           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00893         }
00894         break;
00895       }
00896     }
00897 
00898    else if (type==InletTemp) {
00899       if (aux==0 || naux<=0) return;
00900       //a0: Vel, a1: Temp, a2: Vel
00901       DataType a0=S0*aux[0], a1=aux[1], a2=0.;
00902       if (naux>2) a2 = S0*aux[2];
00903       if (scaling==base::Physical) {
00904         a0 /= U0;
00905         a1 = (a1-Tpmin)/Temp0;
00906         a2 /= U0;
00907       }
00908       switch (side) {
00909       case base::Left:
00910         for (register int j=mys; j<=mye; j++) {
00911           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00912           q(1) = a0;
00913           q(3) = a1;
00914           if (naux>2) q(2) = a2;
00915           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00916         }
00917         break;
00918       case base::Right:
00919         for (register int j=mys; j<=mye; j++) {
00920           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00921           q(1) = a0;
00922           q(3) = a1;
00923           if (naux>2) q(2) = a2;
00924           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00925         }
00926         break;
00927       case base::Bottom:
00928         for (register int i=mxs; i<=mxe; i++) {
00929           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00930           if (naux==2){ 
00931             q(2) = a0;
00932             q(3) = a1;
00933           }
00934           else {
00935             q(1) = a0;
00936             q(3) = a1;
00937             q(2) = a2;
00938           }
00939           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00940         }
00941         break;
00942       case base::Top:
00943         for (register int i=mxs; i<=mxe; i++) {
00944           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00945           if (naux==2){
00946             q(2) = a0;
00947             q(3) = a1;
00948           }
00949           else {
00950             q(1) = a0;
00951             q(3) = a1;
00952             q(2) = a2;
00953           }
00954           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00955         }
00956         break;
00957       }
00958     }
00959 
00960   }
00961 
00962   virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
00963                              const MicroType* f, const point_type* xc, const DataType* distance, 
00964                              const point_type* normal, DataType* aux=0, const int naux=0, 
00965                              const int scaling=0) const {    
00966     
00967     if (type==GFMNoSlipWallTemp || type==GFMSlipWallTemp) { 
00968       DataType fac = S0;
00969       if (scaling==base::Physical)
00970         fac /= U0;      
00971       for (int n=0; n<nc; n++) {
00972         MacroType q=MacroVariables(f[n]);
00973         DataType u=-q(1); 
00974         DataType v=-q(2);
00975         DataType T= q(3);
00976         
00977         // Add boundary velocities if available
00978         if (naux>=3) {
00979           u += fac*aux[naux*n]; 
00980           v += fac*aux[naux*n+1];
00981           T = (aux[naux*n+2]-Tpmin)/(Tpmax-Tpmin);
00982         } 
00983         if (type==GFMNoSlipWallTemp) {
00984           // Prescribe entire velocity vector in ghost cells
00985           q(1) += DataType(2.)*u;
00986           q(2) += DataType(2.)*v;
00987           q(3) = T;
00988         }
00989         else if (type==GFMSlipWallTemp) {
00990           // Prescribe only normal velocity vector in ghost cells
00991           DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
00992           q(1) += vl*normal[n](0);
00993           q(2) += vl*normal[n](1);
00994           q(3) = T;
00995         }
00996         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
00997       } 
00998     }
00999 
01000     else if (type==GFMExtrapolation) {
01001       for (int n=0; n<nc; n++) {
01002         MacroType q=MacroVariables(f[n]);
01003         DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1);
01004         q(1) = vl*normal[n](0);
01005         q(2) = vl*normal[n](1);
01006         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01007       }
01008     }
01009   }
01010     
01011   virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01012                       const int skip_ghosts=1) const {
01013     int Nx = fvec.extents(0), Ny = fvec.extents(1);
01014     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01015     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01016     MicroType *f = (MicroType *)fvec.databuffer();
01017     DataType *work = (DataType *)workvec.databuffer();
01018     DCoords dx = base::GH().worldStep(fvec.stepsize());
01019     DataType dt_temp = dx(0)*S0/U0;
01020 
01021     if (cnt>=1 && cnt<=9) {
01022       for (register int j=mys; j<=mye; j++)
01023         for (register int i=mxs; i<=mxe; i++) {
01024           MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
01025           // 1:r, 2:u, 3:v, 4:E, 5:T, 6:p, 7:gamma, >=8, variable: \nu
01026           switch(cnt) {   
01027           case 1:
01028             work[base::idx(i,j,Nx)]=0;
01029             break;
01030           case 2:
01031             work[base::idx(i,j,Nx)]=q(1)*U0/S0;
01032             break;
01033           case 3:
01034             work[base::idx(i,j,Nx)]=q(2)*U0/S0;
01035             break;
01036           case 5:
01037             work[base::idx(i,j,Nx)]=q(3)*(Tpmax-Tpmin)+Tpmin;
01038             break;
01039           case 6:
01040             work[base::idx(i,j,Nx)]=q(0)*P0;
01041             break;
01042           case 8: {
01043             MicroType feq = Equilibrium(q);       
01044             work[base::idx(i,j,Nx)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt_temp),
01045                                                GasSpeedofSound(),dt_temp );
01046             break;
01047           }
01048           case 9:
01049             work[base::idx(i,j,Nx)]=std::sqrt(q(1)*q(1)+q(2)*q(2))*U0/S0;
01050             break;
01051           }       
01052         }
01053     }
01054     else if (cnt>=10 && cnt<=11) {
01055       macro_grid_data_type qgrid(fvec.bbox());
01056       MacroType *q = (MacroType *)qgrid.databuffer();
01057       for (register int j=mys; j<=mye; j++)
01058         for (register int i=mxs; i<=mxe; i++) {
01059           q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
01060           q[base::idx(i,j,Nx)](1)*=U0/S0;
01061           q[base::idx(i,j,Nx)](2)*=U0/S0;
01062         }
01063 
01064       if (skip_ghosts==0) {
01065         mxs+=1; mxe-=1; mys+=1; mye-=1;
01066       }
01067 
01068       for (register int j=mys; j<=mye; j++)
01069         for (register int i=mxs; i<=mxe; i++) {
01070           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))-
01071             (q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(2.*dx(1));
01072           if (cnt==11) work[base::idx(i,j,Nx)]=std::abs(work[base::idx(i,j,Nx)]);
01073         }
01074     }
01075   }
01076 
01077   virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01078                      const int skip_ghosts=1) const {
01079     int Nx = fvec.extents(0), Ny = fvec.extents(1);
01080     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01081     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01082     MicroType *f = (MicroType *)fvec.databuffer();
01083     DataType *work = (DataType *)workvec.databuffer();
01084 
01085     for (register int j=mys; j<=mye; j++)
01086       for (register int i=mxs; i<=mxe; i++) {
01087         switch(cnt) {
01088         case 2:
01089           f[base::idx(i,j,Nx)](1)=work[base::idx(i,j,Nx)]/(U0/S0);
01090           break;
01091         case 3:
01092           f[base::idx(i,j,Nx)](2)=work[base::idx(i,j,Nx)]/(U0/S0);
01093           f[base::idx(i,j,Nx)] = Equilibrium((const MacroType &)f[base::idx(i,j,Nx)]);
01094           break;
01095         case 5:
01096           f[base::idx(i,j,Nx)](3)=(work[base::idx(i,j,Nx)]-Tpmin)/(Tpmax-Tpmin);
01097           break;
01098         case 6:
01099           f[base::idx(i,j,Nx)](0)=work[base::idx(i,j,Nx)]/P0;
01100           break;
01101         }
01102       }
01103   }
01104 
01105   virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time, 
01106                     const int verbose) const {
01107     int Nx = fvec.extents(0);
01108     int b = fvec.bottom();
01109     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01110     BBox bbmax = grow(fvec.bbox(),-base::NGhosts());    
01111     int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01112     int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01113     int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01114     int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01115     MicroType *f = (MicroType *)fvec.databuffer();
01116     
01117     int result = 1;
01118     for (register int j=mys; j<=mye; j++)
01119       for (register int i=mxs; i<=mxe; i++)
01120         for (register int l=0; l<base::NMicroVar(); l++)          
01121           if (!(f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l)>-1.e37)) {
01122             result = 0;
01123             if (verbose) 
01124               std::cerr << "D2Q9Thermal-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << ")(" << l << ")=" 
01125                         << f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l) << " " << fvec.bbox() << std::endl;
01126           }
01127     return result;
01128   }  
01129         
01130   virtual int NMethodOrder() const { return 2; } 
01131 
01132   inline const DataType& L0() const { return base::LengthScale(); }
01133   inline const DataType& T0() const { return base::TimeScale(); }
01134   inline void SetDensityScale(const DataType r0) { R0 = r0; }
01135   inline void SetPressureScale(const DataType p0) { P0 = p0; }
01136   inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
01137   inline void SetSpeedUp(const DataType s0) { S0 = s0; }
01138   inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
01139   inline void SetTemperatureScale(const DataType tpmin, const DataType tpmax) { 
01140     Tpmin = tpmin; Tpmax = tpmax; 
01141     Temp0 = tpmax-tpmin;
01142     Tref = DataType(0.5)*(Tpmax+Tpmin)/Temp0;
01143   }
01144 
01145   inline const DataType& DensityScale() const { return R0; }
01146   inline const DataType& PressureScale() const { return P0; }
01147   inline const DataType VelocityScale() const { return U0/S0; }
01148   inline const DataType& TemperatureScale() const { return Temp0; }
01149 
01150   // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
01151   inline const DataType& SpeedUp() const { return S0; }
01152 
01153   // Lattice quantities for the operator
01154   inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
01155   inline DataType LatticeViscosityT(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*ds2); }
01156   inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
01157   inline DataType LatticeSpeedOfSoundT() const { return (std::sqrt(ds2)); }
01158 
01159   // Physical quantities for the operator
01160   inline void SetGas(const DataType pr, const DataType rho, const DataType nu, const DataType cs) { 
01161     rhop = rho; nup = nu; csp = cs; cs2p = cs*cs; prp = pr; 
01162     SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
01163     SetPressureScale(prp);
01164   }
01165   //inline void SetThermalGas(const DataType tmin, const DataType tmax, const DataType diff, const DataType force, const DataType width) {
01166   inline void SetThermalGas(const DataType tmin, const DataType tmax, const DataType diff, const DataType gp, const DataType betap) {        
01167     nutp = diff; 
01168     SetTemperatureScale(tmin,tmax);
01169 
01170     //DataType H = width;
01171     DataType g = gp/(L0()/T0()/T0()); //g in Lattice Coordinates
01172     DataType beta = betap*Temp0; //beta in Lattice Coordinates
01173     gbeta = g*beta;
01174 
01175     //Force input is the Rayleigh Number -> Calculating in Lattice directly
01176     //gbeta = force*LatticeViscosity(Omega(T0()))*LatticeViscosityT(OmegaT(T0()))/(H/L0())/(H/L0())/(H/L0());
01177 
01178   }
01179   
01180   inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
01181   inline virtual const DataType OmegaT(const DataType dt) const { return (dcs2*cs2p*dt/S0/(nutp*S0+dcs2*cs2p*dt/S0*DataType(0.5))); }
01182 
01183   inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q, 
01184                                               const DataType dt) const {
01185     DataType M2 = 0.0;
01186 
01187     M2 += (f(4) - feq(4))*(f(4) - feq(4));
01188     M2 += (f(5) - feq(5))*(f(5) - feq(5));
01189     M2 += (f(7) - feq(7))*(f(7) - feq(7));
01190     M2 += (f(8) - feq(8))*(f(8) - feq(8));
01191     M2 = std::sqrt(M2);
01192 
01193     DataType tau = DataType(1.0)/Omega(dt);
01194     DataType om_turb =  (DataType(2.0)/( tau
01195              + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
01196 
01197     return ( om_turb );
01198   }
01199   inline const int TurbulenceType() const { return turbulence; }
01200   inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
01201 
01202   inline const DataType& GasDensity() const { return rhop; }
01203   inline const DataType& GasPressure() const { return prp; }
01204   inline const DataType& GasViscosity() const { return nup; }
01205   inline const DataType& GasViscosityT() const { return nutp; }
01206   inline const DataType& GasSpeedofSound() const { return csp; }
01207   inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
01208   { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
01209   inline const DataType GasViscosityT(const DataType omegat, const DataType cs, const DataType dt) const
01210   { return ((DataType(1.)/omegat-DataType(0.5))*dcs2*cs*cs*dt/S0/S0); }
01211 
01212 protected:
01213   DataType cs2, cs22, cssq, ds2, dcs2;
01214   DataType P0, R0, U0, S0, Temp0, rhop, Tpmin, Tpmax, Tref, csp, cs2p, nup, nutp, prp, gbeta;
01215   DataType Cs_Smagorinsky, turbulence;
01216   int method[1], mdx[14], mdy[14];
01217 };
01218 
01219 
01220 #endif