amroc/lbm/LBMD2Q9.h

Go to the documentation of this file.
00001 // -*- C++ -*-
00002 
00003 // Copyright (C) 2012 Oak Ridge National Laboratory
00004 // Ralf Deiterding, deiterdingr@ornl.gov
00005 
00006 #ifndef LBM_D2Q9_H
00007 #define LBM_D2Q9_H
00008 
00016 #include "LBMBase.h"
00017 #include <cfloat>
00018 
00019 #define LB_FACTOR 1.0e5
00020 
00035 template <class DataType>
00036 class LBMD2Q9 : public LBMBase<Vector<DataType,9>, Vector<DataType,3>, 2> {
00037   typedef LBMBase<Vector<DataType,9>, Vector<DataType,3>, 2> base;
00038 public:
00039   typedef typename base::vec_grid_data_type vec_grid_data_type;
00040   typedef typename base::grid_data_type grid_data_type;
00041   typedef typename base::MicroType MicroType;
00042   typedef typename base::MacroType MacroType;
00043   typedef GridData<MacroType,2> macro_grid_data_type;
00044   typedef typename base::SideName SideName;
00045   typedef typename base::point_type point_type;
00046   typedef typename base::MacroType TensorType;
00047 
00048   enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro }; 
00049   enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall,
00050                       CharacteristicOutlet, CharacteristicInlet, NoSlipWallEntranceExit }; 
00051   enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall };
00052   enum TurbulenceModel { laminar, LES_Smagorinsky };
00053 
00054   LBMD2Q9() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.) {
00055     cs2  = DataType(1.)/DataType(3.);
00056     cs22 = DataType(2.)*cs2;
00057     cssq = DataType(2.)/DataType(9.);
00058     method[0] = 2; method[1] = 0; method[2] = 0; method[3] = 0;
00059     turbulence = laminar;
00060     mdx[0]=mdx[3]=mdx[6]=0; mdx[1]=mdx[4]=mdx[7]=1; mdx[2]=mdx[5]=mdx[8]=-1;
00061     mdy[0]=mdy[1]=mdy[2]=0; mdy[3]=mdy[4]=mdy[5]=1; mdy[6]=mdy[7]=mdy[8]=-1;
00062   }
00063  
00064   virtual ~LBMD2Q9() {}
00065 
00066   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00067     base::LocCtrl = Ctrl.getSubDevice(prefix+"D2Q9");
00068     RegisterAt(base::LocCtrl,"Method(1)",method[0]); 
00069     RegisterAt(base::LocCtrl,"Method(2)",method[1]); 
00070     RegisterAt(base::LocCtrl,"Method(3)",method[2]); 
00071     RegisterAt(base::LocCtrl,"Method(4)",method[3]); 
00072     RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00073     RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00074     RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00075     RegisterAt(base::LocCtrl,"Gas_Cs",csp); 
00076     RegisterAt(base::LocCtrl,"Gas_nu",nup); 
00077     RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00078     RegisterAt(base::LocCtrl,"Gas_W",Wp);
00079     RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00080     RegisterAt(base::LocCtrl,"SpeedUp",S0);     
00081    }
00082 
00083   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
00084     base::SetupData(gh,ghosts);
00085     if (rhop>0.) R0 = rhop;
00086     if (csp>0.) {
00087       cs2p = csp*csp;
00088       SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00089     }
00090     std::cout << "D2Q9: Gas_rho=" << rhop << "   Gas_Cs=" << csp 
00091               << "   Gas_nu=" << nup << std::endl; 
00092     std::cout << "D2Q9: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0 
00093               << "   R0=" << R0 << std::endl; 
00094     WriteInit();
00095   }
00096 
00097   virtual void WriteInit() const {
00098     int me = MY_PROC; 
00099     if (me == VizServer) {
00100       std::ofstream ofs("chem.dat", std::ios::out);
00101       ofs << "RU         " << Rp << std::endl;
00102       ofs << "PA         " << BasePressure() << std::endl;
00103       ofs << "Sp(01)     Vicosity" << std::endl;
00104       ofs << "W(01)      " << Wp << std::endl;
00105       ofs << "Sp(02)     |Velocity|" << std::endl;
00106       ofs << "W(02)      " << 1.0 << std::endl;
00107       ofs << "Sp(03)     Vorticity" << std::endl;
00108       ofs << "W(03)      " << 1.0 << std::endl;
00109       ofs << "Sp(04)     |Vorticity|" << std::endl;
00110       ofs << "W(04)      " << 1.0 << std::endl;
00111      ofs.close();
00112     }    
00113   }
00114 
00115   inline virtual MacroType MacroVariables(const MicroType &f) const {
00116     MacroType q;
00117     q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8);
00118     q(1)=(f(1)-f(2)+f(4)-f(5)+f(7)-f(8))/q(0);
00119     q(2)=(f(3)+f(4)+f(5)-f(6)-f(7)-f(8))/q(0);
00120     return q;
00121   }
00122 
00123   inline virtual MicroType Equilibrium(const MacroType &q) const {
00124     MicroType feq;
00125 
00126     DataType rho = q(0);
00127     DataType u   = q(1);
00128     DataType v   = q(2);
00129 
00130     DataType ri, rt0, rt1, rt2;
00131     if (method[0]==0) {
00132       ri = DataType(1.);
00133       rt0 = R0 * DataType(4.) / DataType(9.);
00134       rt1 = R0 / DataType(9.);
00135       rt2 = R0 / DataType(36.);
00136     }
00137     if (method[0]==1) {
00138       ri = rho;
00139       rt0 = DataType(4.) / DataType(9.);
00140       rt1 = DataType(1.) / DataType(9.);
00141       rt2 = DataType(1.) / DataType(36.);
00142     }
00143     else if (method[0]==2) {
00144       ri = DataType(1.);
00145       rt0 = rho * DataType(4.) / DataType(9.);
00146       rt1 = rho / DataType(9.);
00147       rt2 = rho / DataType(36.);
00148     }
00149                
00150     DataType usq = u * u; 
00151     DataType vsq = v * v;
00152     DataType sumsq = (usq + vsq) / cs22;
00153     DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00154     DataType u2 = usq / cssq; 
00155     DataType v2 = vsq / cssq;
00156     DataType ui = u / cs2;
00157     DataType vi = v / cs2;
00158     DataType uv = ui * vi;
00159 
00160     feq(0) = rt0 * (ri - sumsq);
00161                
00162     feq(1) = rt1 * (ri - sumsq + u2 + ui);
00163     feq(2) = rt1 * (ri - sumsq + u2 - ui);
00164     feq(3) = rt1 * (ri - sumsq + v2 + vi);
00165     feq(6) = rt1 * (ri - sumsq + v2 - vi);
00166                
00167     feq(4) = rt2 * (ri + sumsq2 + ui + vi + uv);
00168     feq(5) = rt2 * (ri + sumsq2 - ui + vi - uv);
00169     feq(7) = rt2 * (ri + sumsq2 + ui - vi - uv);
00170     feq(8) = rt2 * (ri + sumsq2 - ui - vi + uv);    
00171 
00172     // Third order 
00173     if (method[3]>0) {
00174       DataType cs4 = DataType(1.)/(DataType(4.)*cs2*cs2);
00175       DataType upv = (u+v)*cs4*(usq-(u+v)*(u+v)+vsq);
00176       DataType umv = (u-v)*cs4*(usq-(u-v)*(u-v)+vsq);
00177 
00178       feq(1) -= rt1 * u*vsq*cs4;
00179       feq(2) += rt1 * u*vsq*cs4;
00180       feq(3) -= rt1 * usq*v*cs4;
00181       feq(6) += rt1 * usq*v*cs4;
00182 
00183       feq(4) -= rt2 * upv;
00184       feq(5) += rt2 * umv;
00185       feq(7) -= rt2 * umv;
00186       feq(8) += rt2 * upv;
00187     }
00188     return feq;
00189   }
00190 
00191   inline virtual void Collision(MicroType &f, const DataType dt) const {
00192     MacroType q = MacroVariables(f);
00193     MicroType feq = Equilibrium(q);
00194 
00195     DataType omega;
00196     if (turbulence == laminar) 
00197       omega = Omega(dt);
00198     else if (turbulence == LES_Smagorinsky) 
00199       omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00200 
00201     f=(DataType(1.)-omega)*f + omega*feq;
00202   }     
00203 
00204   virtual int IncomingIndices(const int side, int indices[]) const {
00205     switch (side) {
00206     case base::Left:
00207       indices[0] = 1; indices[1] = 4; indices[2] = 7;  
00208       break;
00209     case base::Right:
00210       indices[0] = 2; indices[1] = 5; indices[2] = 8;  
00211       break;
00212     case base::Bottom:
00213       indices[0] = 3; indices[1] = 4; indices[2] = 5;  
00214       break;
00215     case base::Top:
00216       indices[0] = 6; indices[1] = 7; indices[2] = 8;  
00217       break;
00218     }
00219     return 3;
00220   }
00221 
00222   virtual int OutgoingIndices(const int side, int indices[]) const {
00223     switch (side) {
00224     case base::Left:
00225       indices[0] = 2; indices[1] = 5; indices[2] = 8;  
00226       break;
00227     case base::Right:
00228       indices[0] = 1; indices[1] = 4; indices[2] = 7;  
00229       break;
00230     case base::Bottom:
00231       indices[0] = 6; indices[1] = 7; indices[2] = 8;  
00232       break;
00233     case base::Top:
00234       indices[0] = 3; indices[1] = 4; indices[2] = 5;  
00235       break;
00236     }
00237     return 3;
00238   }
00239 
00240   // Revert just one layer of cells at the boundary 
00241   virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
00242                              const int side) const {
00243     int Nx = fvec.extents(0);
00244     int b = fvec.bottom();
00245     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00246     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00247     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00248     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00249     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00250     MicroType *f = (MicroType *)fvec.databuffer();
00251 
00252 #ifdef DEBUG_PRINT_FIXUP
00253     ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
00254       << fvec.bbox() << " using " << bb << "\n").flush();
00255 #endif
00256     
00257     switch (side) {
00258     case base::Left:
00259       for (register int j=mys; j<=mye; j++) {
00260         f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](2);
00261         f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](8);
00262         f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](5);
00263       }
00264       break;
00265     case base::Right:
00266       for (register int j=mys; j<=mye; j++) {
00267         f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](1);
00268         f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](7);
00269         f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](4);
00270       }
00271       break;
00272     case base::Bottom:
00273       for (register int i=mxs; i<=mxe; i++) {
00274         f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](6);
00275         f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](8);
00276         f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](7);
00277       }
00278       break;
00279     case base::Top:
00280       for (register int i=mxs; i<=mxe; i++) {
00281         f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](3);
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     if (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR) {
00324       std::cerr << "LBMD2Q9 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00325                 << " to be equal." << std::endl << std::flush;
00326       std::exit(-1);
00327     }
00328     if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
00329       std::cerr << "LBMD2Q9 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0() 
00330                 << " to be equal." 
00331                 << "   dt=" << dt << "   TO=" << T0()
00332                 << std::endl << std::flush;
00333       std::exit(-1);
00334     }
00335 
00336     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00337     MicroType *f = (MicroType *)fvec.databuffer();
00338 
00339     // Streaming - update all cells
00340     for (register int j=Ny-1; j>=1; j--)
00341       for (register int i=0; i<=Nx-2; i++) {
00342         f[base::idx(i,j,Nx)](3) = f[base::idx(i,j-1,Nx)](3);
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)](4) = f[base::idx(i-1,j-1,Nx)](4);
00350       }
00351 
00352     for (register int j=0; j<=Ny-2; j++)
00353       for (register int i=Nx-1; i>=1; i--) {
00354         f[base::idx(i,j,Nx)](6) = f[base::idx(i,j+1,Nx)](6);
00355         f[base::idx(i,j,Nx)](7) = f[base::idx(i-1,j+1,Nx)](7);
00356       }
00357     for (register int j=0; j<=Ny-2; j++)
00358       for (register int i=0; i<=Nx-2; i++) {
00359         f[base::idx(i,j,Nx)](2) = f[base::idx(i+1,j,Nx)](2);
00360         f[base::idx(i,j,Nx)](8) = f[base::idx(i+1,j+1,Nx)](8);
00361       }    
00362 
00363     // Collision - only in the interior region
00364     int mxs = base::NGhosts(), mys = base::NGhosts();
00365     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00366     for (register int j=mys; j<=mye; j++)
00367       for (register int i=mxs; i<=mxe; i++)
00368         Collision(f[base::idx(i,j,Nx)],dt);
00369 
00370     return DataType(1.);
00371   }
00372 
00373    virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00374                            DataType* aux=0, const int naux=0, const int scaling=0) const {
00375     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00376     int mxs = base::NGhosts(), mys = base::NGhosts();
00377     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00378     MicroType *f = (MicroType *)fvec.databuffer();
00379 
00380     if (type==ConstantMicro) {
00381       MicroType g;
00382       for (register int i=0; i<base::NMicroVar(); i++)
00383         g(i) = aux[i];
00384     
00385       for (register int j=mys; j<=mye; j++)
00386         for (register int i=mxs; i<=mxe; i++)
00387           f[base::idx(i,j,Nx)] = g;
00388     }
00389     
00390     else if (type==ConstantMacro) {
00391       MacroType q;
00392       if (scaling==base::Physical) {
00393         q(0) = aux[0]/R0; 
00394         q(1) = aux[1]/U0; 
00395         q(2) = aux[2]/U0;
00396       }
00397       else 
00398         for (register int i=0; i<base::NMacroVar(); i++)
00399           q(i) = aux[i];
00400       q(1) *= S0; q(2) *= S0;
00401     
00402       for (register int j=mys; j<=mye; j++)
00403         for (register int i=mxs; i<=mxe; i++)
00404           f[base::idx(i,j,Nx)] = Equilibrium(q);
00405     }
00406     
00407     else if (type==GasAtRest) {
00408       MacroType q;
00409       q(0) = GasDensity()/R0; 
00410       q(1) = DataType(0.); 
00411       q(2) = DataType(0.);
00412       for (register int j=mys; j<=mye; j++)
00413         for (register int i=mxs; i<=mxe; i++)
00414           f[base::idx(i,j,Nx)] = Equilibrium(q);
00415     }       
00416   }
00417 
00418   // Set just the one layer of ghost cells at the boundary 
00419   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00420                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00421     int Nx = fvec.extents(0);
00422     int b = fvec.bottom();
00423     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00424     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00425     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00426     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00427     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00428     MicroType *f = (MicroType *)fvec.databuffer();
00429     
00430     if (type==Symmetry) {
00431       switch (side) {
00432       case base::Left:
00433         for (register int j=mys; j<=mye; j++) {
00434           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j,Nx)](8);
00435           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00436           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j,Nx)](5);
00437         }
00438         break;
00439       case base::Right:
00440         for (register int j=mys; j<=mye; j++) {
00441           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j,Nx)](7);
00442           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00443           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j,Nx)](4);
00444         }
00445         break;
00446       case base::Bottom:
00447         for (register int i=mxs; i<=mxe; i++) {
00448           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i,mye+1,Nx)](8);
00449           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00450           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i,mye+1,Nx)](7);
00451         }
00452         break;
00453       case base::Top:
00454         for (register int i=mxs; i<=mxe; i++) {
00455           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i,mys-1,Nx)](5);
00456           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00457           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i,mys-1,Nx)](4);
00458         }
00459         break;
00460       }
00461     }
00462     
00463     else if (type==SlipWall) {
00464       switch (side) {
00465       case base::Left:
00466         for (register int j=mys; j<=mye; j++) {
00467           if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5); 
00468           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00469           if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00470         }
00471         break;
00472       case base::Right:
00473         for (register int j=mys; j<=mye; j++) {
00474           if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4); 
00475           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00476           if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00477         }
00478         break;
00479       case base::Bottom:
00480         for (register int i=mxs; i<=mxe; i++) {
00481           if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00482           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00483           if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00484         }
00485         break;
00486       case base::Top:
00487         for (register int i=mxs; i<=mxe; i++) {
00488           if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00489           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00490           if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00491         }
00492         break;
00493       }
00494     }
00495 
00496     else if (type==NoSlipWall) {
00497       switch (side) {
00498       case base::Left:
00499         for (register int j=mys; j<=mye; j++) {
00500           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00501           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00502           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00503         }
00504         break;
00505       case base::Right:
00506         for (register int j=mys; j<=mye; j++) {
00507           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00508           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00509           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00510         }
00511         break;  
00512       case base::Bottom:
00513         for (register int i=mxs; i<=mxe; i++) {
00514           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00515           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00516           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00517         }
00518         break;  
00519       case base::Top:
00520         for (register int i=mxs; i<=mxe; i++) {
00521           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00522           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00523           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00524         }
00525         break;
00526       }
00527     }
00528 
00529     else if (type==Inlet) {
00530       if (aux==0 || naux<=0) return;
00531       DataType a0=S0*aux[0], a1=0.;
00532       if (naux>1) a1 = S0*aux[1];
00533       if (scaling==base::Physical) {
00534         a0 /= U0;
00535         a1 /= U0;
00536       }
00537       switch (side) {
00538       case base::Left:
00539         for (register int j=mys; j<=mye; j++) {
00540           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00541           q(1) = a0;
00542           if (naux>1) q(2) = a1;
00543           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00544         }
00545         break;
00546       case base::Right:
00547         for (register int j=mys; j<=mye; j++) {
00548           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00549           q(1) = a0;
00550           if (naux>1) q(2) = a1;
00551           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00552         }
00553         break;
00554       case base::Bottom:
00555         for (register int i=mxs; i<=mxe; i++) {
00556           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00557           if (naux==1) 
00558             q(2) = a0;
00559           else {
00560             q(1) = a0;
00561             q(2) = a1;
00562           }
00563           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00564         }
00565         break;
00566       case base::Top:
00567         for (register int i=mxs; i<=mxe; i++) {
00568           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00569           if (naux==1) 
00570             q(2) = a0;
00571           else {
00572             q(1) = a0;
00573             q(2) = a1;
00574           }
00575           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00576         }
00577         break;
00578       }
00579     }
00580 
00581     else if (type==Outlet) {
00582       switch (side) {
00583       case base::Left:
00584         for (register int j=mys; j<=mye; j++) {
00585           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);       
00586           if (q(0)>0) f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)]/q(0);
00587         }
00588         break;
00589       case base::Right:
00590         for (register int j=mys; j<=mye; j++) {
00591           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00592           if (q(0)>0) f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)]/q(0);
00593         }
00594         break;
00595       case base::Bottom:
00596         for (register int i=mxs; i<=mxe; i++) {
00597           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00598           if (q(0)>0) f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)]/q(0);
00599         }
00600         break;
00601       case base::Top:
00602         for (register int i=mxs; i<=mxe; i++) {
00603           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00604           if (q(0)>0) f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)]/q(0);
00605         }
00606         break;
00607       }
00608     }
00609     
00610     // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)    
00611     else if (type==Pressure) {
00612       if (aux==0 || naux<=0) return;
00613       DataType a0=aux[0];
00614       if (scaling==base::Physical)
00615         a0 /= R0;
00616       switch (side) {
00617       case base::Left:
00618         for (register int j=mys; j<=mye; j++) {
00619           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00620           q(0) = a0;
00621           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00622         }
00623         break;
00624       case base::Right:
00625         for (register int j=mys; j<=mye; j++) {
00626           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00627           q(0) = a0;
00628           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00629         }
00630         break;
00631       case base::Bottom:
00632         for (register int i=mxs; i<=mxe; i++) {
00633           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00634           q(0) = a0;
00635           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00636         }
00637         break;
00638       case base::Top:
00639         for (register int i=mxs; i<=mxe; i++) {
00640           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00641           q(0) = a0;
00642           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00643         }
00644         break;
00645       }
00646     }
00647     
00648     else if (type==SlidingWall) {
00649       if (aux==0 || naux<=0) return;
00650       if (aux==0 || naux<=0) return;
00651       DataType a0=S0*aux[0];
00652       if (scaling==base::Physical)
00653         a0 /= U0;
00654       switch (side) {
00655       case base::Left:
00656         for (register int j=mys; j<=mye; j++) {
00657           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00658           DataType rd1 = f[b+base::idx(mxe+1,j,Nx)](5), rd2 = f[b+base::idx(mxe+1,j,Nx)](8);
00659           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00660           DataType pw = DataType(1.)-qw;
00661           if (j<mye) f[b+base::idx(mxe,j+1,Nx)](7) = pw*rd1+qw*rd2;
00662           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00663           if (j>mys) f[b+base::idx(mxe,j-1,Nx)](4) = qw*rd1+pw*rd2;
00664         }
00665         break;
00666       case base::Right:
00667         for (register int j=mys; j<=mye; j++) {
00668           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00669           DataType rd1 = f[b+base::idx(mxs-1,j,Nx)](4), rd2 = f[b+base::idx(mxs-1,j,Nx)](7);
00670           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00671           DataType pw = DataType(1.)-qw;
00672           if (j<mye) f[b+base::idx(mxs,j+1,Nx)](8) = pw*rd1+qw*rd2;
00673           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00674           if (j>mys) f[b+base::idx(mxs,j-1,Nx)](5) = qw*rd1+pw*rd2;
00675         }
00676         break;
00677       case base::Bottom:
00678         for (register int i=mxs; i<=mxe; i++) {
00679           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00680           DataType rd1 = f[b+base::idx(i,mye+1,Nx)](7), rd2 = f[b+base::idx(i,mye+1,Nx)](8);
00681           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00682           DataType pw = DataType(1.)-qw;
00683           if (i<mxe) f[b+base::idx(i+1,mye,Nx)](5) = pw*rd1+qw*rd2;
00684           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00685           if (i>mxs) f[b+base::idx(i-1,mye,Nx)](4) = qw*rd1+pw*rd2;
00686         }
00687         break;
00688       case base::Top:
00689         for (register int i=mxs; i<=mxe; i++) {
00690           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00691           DataType rd1 = f[b+base::idx(i,mys-1,Nx)](4), rd2 = f[b+base::idx(i,mys-1,Nx)](5);
00692           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00693           DataType pw = DataType(1.)-qw;
00694           if (i<mxe) f[b+base::idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00695           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00696           if (i>mxs) f[b+base::idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00697         }
00698         break;
00699       }
00700     }
00701     
00705     else if (type==CharacteristicOutlet) {
00706       DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
00707       DCoords dx = base::GH().worldStep(fvec.stepsize());
00708       DataType dt_temp = dx(0)/VelocityScale();
00709       switch (side) {
00710       case base::Left:
00711         for (register int j=mys; j<=mye; j++) {
00712           MicroType ftmp = f[b+base::idx(mxe,j,Nx)];
00713           MacroType q = MacroVariables(ftmp);
00714           Collision(ftmp,dt_temp);
00715           MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00716           MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,Nx)]);
00717 
00718           MacroType dqdn = NormalDerivative(q,q1,q2);
00719           Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00720 
00721           // Solve LODI equations :
00722           MacroType qn;
00723           qn(0) = q(0) - c1oCsSq2*L(0);
00724           qn(1) = q(1) + c1o2cs*L(0)/q(0);
00725           qn(2) = q(2) - L(1);
00726           f[b+base::idx(mxe,j,Nx)] = Equilibrium(qn);
00727         }
00728         break;
00729       case base::Right:
00730         for (register int j=mys; j<=mye; j++) {
00731           MicroType ftmp = f[b+base::idx(mxs,j,Nx)];
00732           MacroType q = MacroVariables(ftmp);
00733           Collision(ftmp,dt_temp);
00734           MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00735           MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,Nx)]);
00736 
00737           MacroType dqdn = -NormalDerivative(q,q1,q2);
00738           if (method[0]!=3) q(0) = q(0);
00739           Vector<DataType,2> L = WaveAmplitudes(q(0),-q(1),dqdn(0),dqdn(1),dqdn(2));
00740 
00741           // Solve LODI equations :
00742           MacroType qn;
00743           qn(0) = q(0) - c1oCsSq2*L(0);
00744           qn(1) = q(1) - c1o2cs*L(0)/q(0);
00745           qn(2) = q(2) + L(1);
00746           f[b+base::idx(mxs,j,Nx)] = Equilibrium(qn);
00747         }
00748         break;
00749       case base::Bottom:
00750         //for (register int i=mxs+base::NGhosts(); i<=mxe-base::NGhosts(); i++) {
00751         for (register int i=mxs; i<=mxe; i++) {
00752           MicroType ftmp = f[b+base::idx(i,mye,Nx)];
00753           MacroType q  = MacroVariables(ftmp);
00754           Collision(ftmp,dt_temp);
00755           MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00756           MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,Nx)]);
00757 
00758           MacroType dqdn = NormalDerivative(q,q1,q2);
00759           if (method[0]!=3) q(0) = q(0);
00760           Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00761 
00762           // Solve LODI equations :
00763           MacroType qn;
00764           qn(0) = q(0) - c1oCsSq2*L(0);
00765           qn(1) = q(1) - L(1);
00766           qn(2) = q(2) + c1o2cs*L(0)/q(0);
00767           f[b+base::idx(i,mye,Nx)] = Equilibrium(qn);
00768         }
00769         break;
00770       case base::Top:
00771         for (register int i=mxs+base::NGhosts(); i<=mxe-base::NGhosts(); i++) {
00772           MicroType ftmp = f[b+base::idx(i,mys,Nx)];
00773           MacroType q  = MacroVariables(ftmp);
00774           Collision(ftmp,dt_temp);
00775           MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00776           MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,Nx)]);
00777 
00778           MacroType dqdn = -NormalDerivative(q,q1,q2);
00779           if (method[0]!=3) q(0) = q(0);
00780           Vector<DataType,2> L = WaveAmplitudes(q(0),-q(2),dqdn(0),dqdn(2),dqdn(1));
00781 
00782           // Solve LODI equations :
00783           MacroType qn;
00784           qn(0) = q(0) - c1oCsSq2*L(0);
00785           qn(1) = q(1) - L(1);
00786           qn(2) = q(2) + c1o2cs*L(0)/q(0);
00787           f[b+base::idx(i,mys,Nx)] = Equilibrium(qn);
00788         }
00789         break;
00790       }
00791     }
00792 
00796     else if (type==CharacteristicInlet) {
00797       if (aux==0 || naux<=0) return;
00798       DataType rho=aux[0];
00799       DataType a0=S0*aux[1];
00800       DataType a1=S0*aux[2];
00801       if (scaling==base::Physical) {
00802         rho /= R0;
00803         a0 /= U0;
00804         a1 /= U0;
00805       }
00806       MacroType q, q1, q2, qn;
00807       DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
00808       q(0) = rho;
00809       q(1) = a0;
00810       q(2) = a1;
00811       switch (side) {
00812       case base::Left:
00813         for (register int j=mys; j<=mye; j++) {
00814           q1 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00815           q2 = MacroVariables(f[b+base::idx(mxe+2,j,Nx)]);
00816           if (naux==1) { // pressure specified
00817             q(1) = q1(1);  q(2) = q1(2);
00818             MacroType dqdn = NormalDerivative(q,q1,q2);
00819             Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00820             qn(0) = q(0) - c1oCsSq2*L(0);
00821             qn(1) = q1(1) + c1o2cs*L(0)/q(0);
00822             qn(2) = q1(2) - L(1);
00823           }
00824           else if (naux==2) { // normal velocity component specified
00825             q(0) = q1(0);  q(2) = q1(2);
00826             MacroType dqdn = NormalDerivative(q,q1,q2);
00827             Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00828             qn(0) = q1(0) - c1oCsSq2*L(0);
00829             qn(1) = q(1) + c1o2cs*L(0)/q(0);
00830             qn(2) = q1(2) - L(1);
00831           }
00832           else if (naux==3) { // normal and transverse velocity components specified
00833             q(0) = q1(0);
00834             MacroType dqdn = NormalDerivative(q,q1,q2);
00835             Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00836             qn(0) = q1(0) - c1oCsSq2*L(0);
00837             qn(1) = q(1) + c1o2cs*L(0)/q(0);
00838             qn(2) = q(2) - L(1);
00839           }
00840           f[b+base::idx(mxe,j,Nx)] = Equilibrium(qn);
00841         }
00842         break;
00843       case base::Right:
00844         for (register int j=mys; j<=mye; j++) {
00845           q1 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00846           q2 = MacroVariables(f[b+base::idx(mxs-2,j,Nx)]);
00847           if (naux==1) {
00848             q(1) = q1(1);  q(2) = q1(2);
00849             MacroType dqdn = NormalDerivative(q,q1,q2);
00850             Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00851             qn(0) = q(0) - c1oCsSq2*L(0);
00852             qn(1) = q1(1) + c1o2cs*L(0)/q(0);
00853             qn(2) = q1(2) - L(1);
00854           }
00855           else if (naux==2) {
00856             q(0) = q1(0);  q(2) = q1(2);
00857             MacroType dqdn = NormalDerivative(q,q1,q2);
00858             Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00859             qn(0) = q1(0) - c1oCsSq2*L(0);
00860             qn(1) = q(1) + c1o2cs*L(0)/q(0);
00861             qn(2) = q1(2) - L(1);
00862           }
00863           else if (naux==3) {
00864             q(0) = q1(0);
00865             MacroType dqdn = NormalDerivative(q,q1,q2);
00866             Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00867             qn(0) = q1(0) - c1oCsSq2*L(0);
00868             qn(1) = q(1) + c1o2cs*L(0)/q(0);
00869             qn(2) = q(2) - L(1);
00870           }
00871           f[b+base::idx(mxs,j,Nx)] = Equilibrium(qn);
00872         }
00873         break;
00874       case base::Bottom:
00875         for (register int i=mxs; i<=mxe; i++) {
00876           q1 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00877           q2 = MacroVariables(f[b+base::idx(i,mye+2,Nx)]);
00878           if (naux==1) {
00879             q(1) = q1(1);  q(2) = q1(2);
00880             MacroType dqdn = NormalDerivative(q,q1,q2);
00881             Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00882             qn(0) = q(0) - c1oCsSq2*L(0);
00883             qn(1) = q1(1) - L(1);
00884             qn(2) = q1(2) + c1o2cs*L(0)/q(0);
00885           }
00886           else if (naux==2) {
00887             q(0) = q1(0);  q(1) = q1(1);
00888             MacroType dqdn = NormalDerivative(q,q1,q2);
00889             Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00890             qn(0) = q1(0) - c1oCsSq2*L(0);
00891             qn(1) = q1(1) - L(1);
00892             qn(2) = q(2) + c1o2cs*L(0)/q(0);
00893           }
00894           else if (naux==3) {
00895             q(0) = q1(0);
00896             MacroType dqdn = NormalDerivative(q,q1,q2);
00897             Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00898             qn(0) = q1(0) - c1oCsSq2*L(0);
00899             qn(1) = q(1) - L(1);
00900             qn(2) = q(2) + c1o2cs*L(0)/q(0);
00901           }
00902           f[b+base::idx(i,mye,Nx)] = Equilibrium(qn);
00903         }
00904         break;
00905       case base::Top:
00906         for (register int i=mxs; i<=mxe; i++) {
00907           q1 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00908           q2 = MacroVariables(f[b+base::idx(i,mys-2,Nx)]);
00909           if (naux==1) {
00910             q(1) = q1(1);  q(2) = q1(2);
00911             MacroType dqdn = NormalDerivative(q,q1,q2);
00912             Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00913             qn(0) = q(0) - c1oCsSq2*L(0);
00914             qn(1) = q1(1) - L(1);
00915             qn(2) = q1(2) + c1o2cs*L(0)/q(0);
00916           }
00917           else if (naux==2) {
00918             q(0) = q1(0);  q(1) = q1(1);
00919             MacroType dqdn = NormalDerivative(q,q1,q2);
00920             Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00921             qn(0) = q1(0) - c1oCsSq2*L(0);
00922             qn(1) = q1(1) - L(1);
00923             qn(2) = q(2) + c1o2cs*L(0)/q(0);
00924           }
00925           else if (naux==3) {
00926             q(0) = q1(0);
00927             MacroType dqdn = NormalDerivative(q,q1,q2);
00928             Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00929             qn(0) = q1(0) - c1oCsSq2*L(0);
00930             qn(1) = q(1) - L(1);
00931             qn(2) = q(2) + c1o2cs*L(0)/q(0);
00932           }
00933           f[b+base::idx(i,mys,Nx)] = Equilibrium(qn);
00934         }
00935         break;
00936       }
00937     }
00938 
00943     else if (type==NoSlipWallEntranceExit) {
00944       DCoords wc;
00945       int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
00946       DataType WLB[2],  WUB[2];
00947       base::GH().wlb(WLB);
00948       base::GH().wub(WUB);
00949       DataType slip = 0.;
00950       DataType lxs, lxe, lys, lye;
00951       DataType min[2], max[2];
00952       if (naux!=4) assert(0);
00953       lxs = aux[0] - WLB[0];  lxe = WUB[0] - aux[1];
00954       lys = aux[2] - WLB[1];  lye = WUB[1] - aux[3];
00955       min[0] = aux[0]; max[0] = aux[1];
00956       min[1] = aux[2]; max[1] = aux[3];
00957 
00958       switch (side) {
00959       case base::Left:
00960         for (register int j=mys; j<=mye; j++) {
00961           lc_indx[1] = j*fvec.stepsize(1);
00962           wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00963           if (wc(1)>=min[1] && wc(1)<=max[1]) {
00964             if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00965             f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00966             if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00967           }
00968           else {
00969             if (wc(1)<WLB[1] || wc(1) >WUB[1])
00970               slip = DataType(0.);
00971             else if (wc(1)<min[1])
00972               slip = (wc(1) - WLB[1])/lys;
00973             else if (wc(1)>max[1])
00974               slip = (WUB[1] - wc(1))/lye;
00975             if (j>mys && j<mye) {
00976               f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
00977                   + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
00978               f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
00979                   + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
00980             }
00981             else if (j==mys) {
00982               f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](5)
00983                   + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
00984               f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
00985                   + slip*f[b+base::idx(mxe+1,j,Nx)](5);
00986             }
00987             else if (j==mye) {
00988               f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
00989                   + slip*f[b+base::idx(mxe+1,j,Nx)](8);
00990               f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](8)
00991                   + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
00992             }
00993             f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00994           }
00995         }
00996         break;
00997       case base::Right:
00998         for (register int j=mys; j<=mye; j++) {
00999           lc_indx[1] = j*fvec.stepsize(1);
01000           wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01001           if (wc(1)>=min[1] && wc(1)<=max[1]) {
01002             if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
01003             f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
01004             if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
01005           }
01006           else {
01007             if (wc(1)<WLB[1] || wc(1) >WUB[1])
01008               slip = DataType(0.);
01009             else if (wc(1)<min[1])
01010               slip = (wc(1) - WLB[1])/lys;
01011             else if (wc(1)>max[1])
01012               slip = (WUB[1] - wc(1))/lye;
01013             if (j>mys && j<mye) {
01014               f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
01015                   + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
01016               f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
01017                   + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
01018             }
01019             else if (j==mys) {
01020               f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](4)
01021                   + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
01022               f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
01023                   + slip*f[b+base::idx(mxs-1,j,Nx)](4);
01024             }
01025             else if (j==mye) {
01026               f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
01027                   + slip*f[b+base::idx(mxs-1,j,Nx)](7);
01028               f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](7)
01029                   + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
01030             }
01031             f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
01032           }
01033         }
01034         break;
01035       case base::Bottom:
01036         for (register int i=mxs; i<=mxe; i++) {
01037           lc_indx[0] = i*fvec.stepsize(0);
01038           wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01039           if (wc(0)>=min[0] && wc(0)<=max[0]) {
01040             if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
01041             f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
01042             if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
01043           }
01044           else {
01045             if (wc(0)<WLB[0] || wc(0)>WUB[0])
01046               slip = DataType(0.);
01047             else if (wc(0)<min[0])
01048               slip = (wc(0) - WLB[0])/lxs;
01049             else if (wc(0)>max[0])
01050               slip = (WUB[0] - wc(0))/lxe;
01051             if (i>mxs && i<mxe) {
01052               f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,Nx)](7)
01053                   + slip*f[b+base::idx(i+1,mye+1,Nx)](8);
01054               f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,Nx)](8)
01055                   + slip*f[b+base::idx(i-1,mye+1,Nx)](7);
01056             }
01057             else if (i==mxs) {
01058               f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](7)
01059                   + slip*f[b+base::idx(i+1,mye+1,Nx)](8);
01060               f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,Nx)](8)
01061                   + slip*f[b+base::idx(i,mye+1,Nx)](7);
01062             }
01063             else if (i==mxe) {
01064               f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](7)
01065                   + slip*f[b+base::idx(i,mye+1,Nx)](8);
01066               f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](8)
01067                   + slip*f[b+base::idx(i-1,mye+1,Nx)](7);
01068             }
01069             f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
01070           }
01071         }
01072         break;
01073       case base::Top:
01074         for (register int i=mxs; i<=mxe; i++) {
01075           lc_indx[0] = i*fvec.stepsize(0);
01076           wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01077           if (wc(0)>=min[0] && wc(0)<=max[0]) {
01078             if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
01079             f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
01080             if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
01081           }
01082           else {
01083             if (wc(0)<WLB[0] || wc(0)>WUB[0])
01084               slip = DataType(0.);
01085             else if (wc(0)<min[0])
01086               slip = (wc(0) - WLB[0])/lxs;
01087             else if (wc(0)>max[0])
01088               slip = (WUB[0] - wc(0))/lxe;
01089             if (i>mxs && i<mxe) {
01090               f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,Nx)](4)
01091                   + slip*f[b+base::idx(i+1,mys-1,Nx)](5);
01092               f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,Nx)](5)
01093                   + slip*f[b+base::idx(i-1,mys-1,Nx)](4);
01094             }
01095             else if (i==mxs) {
01096               f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,Nx)](4)
01097                   + slip*f[b+base::idx(i+1,mys-1,Nx)](5);
01098               f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,Nx)](5)
01099                   + slip*f[b+base::idx(i,mys-1,Nx)](4);
01100             }
01101             else if (i==mxe) {
01102               f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,Nx)](4)
01103                   + slip*f[b+base::idx(i,mys-1,Nx)](5);
01104               f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,Nx)](5)
01105                   + slip*f[b+base::idx(i-1,mys-1,Nx)](4);
01106             }
01107             f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
01108           }
01109         }
01110         break;
01111       }
01112     }
01113   }
01114     
01115   virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
01116                              const MicroType* f, const point_type* xc, const DataType* distance, 
01117                              const point_type* normal, DataType* aux=0, const int naux=0, 
01118                              const int scaling=0) const {    
01119     
01120     if (type==GFMNoSlipWall || type==GFMSlipWall) { 
01121       DataType fac = S0;
01122       if (scaling==base::Physical)
01123         fac /= U0;      
01124       for (int n=0; n<nc; n++) {
01125         MacroType q=MacroVariables(f[n]);
01126         DataType u=-q(1); 
01127         DataType v=-q(2);
01128         
01129         // Add boundary velocities if available
01130         if (naux>=2) {
01131           u += fac*aux[naux*n]; 
01132           v += fac*aux[naux*n+1];
01133         } 
01134         if (type==GFMNoSlipWall) {
01135           // Prescribe entire velocity vector in ghost cells
01136           q(1) += DataType(2.)*u;
01137           q(2) += DataType(2.)*v;
01138         }
01139         else if (type==GFMSlipWall) {
01140           // Prescribe only normal velocity vector in ghost cells
01141           DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
01142           q(1) += vl*normal[n](0);
01143           q(2) += vl*normal[n](1);
01144         }
01145         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01146       } 
01147     }
01148 
01149     else if (type==GFMExtrapolation) {
01150       for (int n=0; n<nc; n++) {
01151         MacroType q=MacroVariables(f[n]);
01152         DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1);
01153         q(1) = vl*normal[n](0);
01154         q(2) = vl*normal[n](1);
01155         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01156       }
01157     }
01158   }
01159     
01160   inline virtual Vector<DataType,2> WaveAmplitudes(const DataType rho, const DataType vn, const DataType drhodn,
01161                                                    const DataType dvndn, const DataType dvtdn) const { 
01162     Vector<DataType,2> L;
01163     L(0) = (vn - LatticeSpeedOfSound())*(cs2*drhodn - rho*LatticeSpeedOfSound()*dvndn);
01164     L(1) = vn*dvtdn;
01165     return L;
01166   }
01167   inline virtual const MacroType NormalDerivative(const MacroType qa, const MacroType qb, const MacroType qc) const 
01168   { return DataType(-1.5)*qa + DataType(2.)*qb - DataType(0.5)*qc; }
01169 
01170   virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01171                       const int skip_ghosts=1) const {
01172     int Nx = fvec.extents(0), Ny = fvec.extents(1);
01173     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01174     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01175     MicroType *f = (MicroType *)fvec.databuffer();
01176     DataType *work = (DataType *)workvec.databuffer();
01177     DCoords dx = base::GH().worldStep(fvec.stepsize());
01178     DataType dt = dx(0)*S0/U0;
01179 
01180     if ((cnt>=1 && cnt<=9) || (cnt>=16 && cnt<=18) || (cnt>=23 && cnt<=25)) { 
01181       for (register int j=mys; j<=mye; j++)
01182         for (register int i=mxs; i<=mxe; i++) {
01183           MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
01184           switch(cnt) {
01185           case 1:
01186             if (method[0]==0) work[base::idx(i,j,Nx)]=q(0)*R0+rhop;
01187             else 
01188               work[base::idx(i,j,Nx)]=q(0)*R0;
01189             break;
01190           case 2:
01191             work[base::idx(i,j,Nx)]=q(1)*U0/S0;
01192             break;
01193           case 3:
01194             work[base::idx(i,j,Nx)]=q(2)*U0/S0;
01195             break;
01196           case 5:
01197             work[base::idx(i,j,Nx)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure()); 
01198             break;
01199           case 6:
01200             work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure(); 
01201             break;
01202           case 8: {
01203             MicroType feq = Equilibrium(q);       
01204             work[base::idx(i,j,Nx)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt),
01205                                                  GasSpeedofSound(),dt);
01206             break;
01207           }
01208           case 9:
01209             work[base::idx(i,j,Nx)]=std::sqrt(q(1)*q(1)+q(2)*q(2))*U0/S0;
01210             break;
01211           case 16:
01212           case 17:
01213           case 18: {
01214             MicroType feq = Equilibrium(q);
01215             DataType omega;
01216             if (turbulence == laminar) 
01217               omega = Omega(dt);
01218             else if (turbulence == LES_Smagorinsky) 
01219               omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt);
01220             work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*DeviatoricStress(f[base::idx(i,j,Nx)],feq,omega)(cnt-16);
01221             break;
01222           }
01223           case 23: 
01224           case 24:
01225           case 25: {
01226             DataType omega;
01227             if (turbulence == laminar) 
01228               omega = Omega(dt);
01229             else if (turbulence == LES_Smagorinsky) {
01230               MicroType feq = Equilibrium(q);
01231               omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt);
01232             }
01233             work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*Stress(f[base::idx(i,j,Nx)],q,omega)(cnt-23);
01234             if (cnt==23 || cnt==25) work[base::idx(i,j,Nx)]-=BasePressure();
01235             break;
01236           }
01237           }
01238         }
01239     }
01240     else if ((cnt>=10 && cnt<=11) || cnt==15 || (cnt>=30 && cnt<=32)) {
01241       macro_grid_data_type qgrid(fvec.bbox());
01242       MacroType *q = (MacroType *)qgrid.databuffer();
01243       for (register int j=mys; j<=mye; j++)
01244         for (register int i=mxs; i<=mxe; i++) {
01245           q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
01246           q[base::idx(i,j,Nx)](0)*=R0;
01247           q[base::idx(i,j,Nx)](1)*=U0/S0;
01248           q[base::idx(i,j,Nx)](2)*=U0/S0;
01249         }
01250       
01251       if (skip_ghosts==0) {
01252         switch(cnt) {
01253         case 30:
01254           mxs+=1; mxe-=1; 
01255           break;
01256         case 32:
01257           mys+=1; mye-=1; 
01258           break;
01259         case 10:
01260         case 11:
01261         case 15:
01262         case 31:
01263           mxs+=1; mxe-=1; mys+=1; mye-=1; 
01264           break;
01265         }
01266       }
01267 
01268       DataType dxux, dxuy, dyux, dyuy;
01269       for (register int j=mys; j<=mye; j++)
01270         for (register int i=mxs; i<=mxe; i++) {
01271           if (cnt==15 || cnt==30)
01272             dxux=(q[base::idx(i+1,j,Nx)](1)-q[base::idx(i-1,j,Nx)](1))/(2.*dx(0));
01273           if (cnt==15 || cnt==32)
01274             dyuy=(q[base::idx(i,j+1,Nx)](2)-q[base::idx(i,j-1,Nx)](2))/(2.*dx(1));
01275           if (cnt==10 || cnt==11 || cnt==15 || cnt==31) {
01276             dxuy=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(2.*dx(0));
01277             dyux=(q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(2.*dx(1));
01278           }
01279 
01280           switch (cnt) {
01281           case 10:
01282             work[base::idx(i,j,Nx)]=dxuy-dyux;
01283             break;
01284           case 11: 
01285             work[base::idx(i,j,Nx)]=std::abs(dxuy-dyux);
01286             break;
01287           case 15: 
01288             work[base::idx(i,j,Nx)]=DataType(-0.5)*(4.*dxux*dxux+8.*dyux*dxuy+4.*dyuy*dyuy); 
01289             break;
01290           case 30:
01291             work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dxux; 
01292             break;
01293           case 31:
01294             work[base::idx(i,j,Nx)]=q[base::idx(i,j,Nx)](0)*GasViscosity()*(dxuy+dyux); 
01295             break;
01296           case 32:
01297             work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dyuy; 
01298             break;
01299           }
01300         }
01301     }
01302   }
01303 
01304   virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01305                      const int skip_ghosts=1) const {
01306     int Nx = fvec.extents(0), Ny = fvec.extents(1);
01307     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01308     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01309     MicroType *f = (MicroType *)fvec.databuffer();
01310     DataType *work = (DataType *)workvec.databuffer();
01311 
01312     for (register int j=mys; j<=mye; j++)
01313       for (register int i=mxs; i<=mxe; i++) {
01314         switch(cnt) {
01315         case 1:
01316           f[base::idx(i,j,Nx)](0)=work[base::idx(i,j,Nx)]/R0;
01317           break;
01318         case 2:
01319           f[base::idx(i,j,Nx)](1)=work[base::idx(i,j,Nx)]/(U0/S0);
01320           break;
01321         case 3:
01322           f[base::idx(i,j,Nx)](2)=work[base::idx(i,j,Nx)]/(U0/S0);
01323           f[base::idx(i,j,Nx)] = Equilibrium((const MacroType &)f[base::idx(i,j,Nx)]);
01324           break;
01325         }
01326       }
01327   }
01328 
01329   virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time, 
01330                     const int verbose) const {
01331     int Nx = fvec.extents(0);
01332     int b = fvec.bottom();
01333     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01334     BBox bbmax = grow(fvec.bbox(),-base::NGhosts());    
01335     int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01336     int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01337     int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01338     int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01339     MicroType *f = (MicroType *)fvec.databuffer();
01340 
01341     int result = 1;
01342     for (register int j=mys; j<=mye; j++)
01343       for (register int i=mxs; i<=mxe; i++)
01344         for (register int l=0; l<base::NMicroVar(); l++) 
01345           if (!(f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l)>-1.e37)) {
01346             result = 0;
01347             if (verbose) 
01348               std::cerr << "D2Q9-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << ")(" << l << ")=" 
01349                         << f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l) << " " << fvec.bbox() << std::endl;
01350           }
01351     return result;
01352   }  
01353         
01354   inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
01355     TensorType P;
01356     if (method[1]==0) {
01357       P(0) = q(0)*q(1)*q(1)-(f(1)+f(2)+f(4)+f(5)+f(7)+f(8))+cs2*q(0);
01358       P(1) = q(0)*q(1)*q(2)-(f(4)-f(5)-f(7)+f(8));
01359       P(2) = q(0)*q(2)*q(2)-(f(3)+f(4)+f(5)+f(6)+f(7)+f(8))+cs2*q(0);
01360       P *= -(DataType(1.0)-DataType(0.5)*om);
01361     }
01362     else {
01363       MicroType feq = Equilibrium(q);     
01364       P = DeviatoricStress(f,feq,om);
01365     }
01366     P(0) -= LatticeBasePressure(q(0));  
01367     P(2) -= LatticeBasePressure(q(0));
01368     return P;
01369   }
01370 
01371   inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
01372     TensorType Sigma;
01373     if (method[2]==0) {
01374       MicroType fneq = feq - f;
01375       Sigma(0) = fneq(1)+fneq(2)+fneq(4)+fneq(5)+fneq(7)+fneq(8);
01376       Sigma(1) = fneq(4)-fneq(5)-fneq(7)+fneq(8);
01377       Sigma(2) = fneq(3)+fneq(4)+fneq(5)+fneq(6)+fneq(7)+fneq(8);
01378       Sigma *= -(DataType(1.0)-DataType(0.5)*om);
01379     }
01380     else {
01381       MacroType q = MacroVariables(f);
01382       Sigma = Stress(f,q,om);
01383       Sigma(0) += LatticeBasePressure(q(0));  
01384       Sigma(2) += LatticeBasePressure(q(0));
01385     } 
01386     return Sigma;
01387   }
01388 
01389   virtual int NMethodOrder() const { return 2; } 
01390 
01391   inline const DataType& L0() const { return base::LengthScale(); }
01392   inline const DataType& T0() const { return base::TimeScale(); }
01393   inline void SetDensityScale(const DataType r0) { R0 = r0; }
01394   inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
01395   inline void SetSpeedUp(const DataType s0) { S0 = s0; }
01396   inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
01397 
01398   inline const DataType& DensityScale() const { return R0; }
01399   inline const DataType VelocityScale() const { return U0/S0; }
01400   // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
01401   inline const DataType& SpeedUp() const { return S0; }
01402 
01403   // Lattice quantities for the operator
01404   inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
01405   inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
01406   inline virtual DataType LatticeBasePressure(const DataType rho) const {
01407     return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
01408   }
01409 
01410   // Physical quantities for the operator
01411   inline void SetGas(DataType rho, DataType nu, DataType cs) { 
01412     rhop = rho; nup = nu; csp = cs; cs2p = cs*cs; 
01413     SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
01414   }
01415   inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
01416   inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q, 
01417                                               const DataType dt) const {
01418     DataType M2 = 0.0;
01419 
01420     M2 += (f(4) - feq(4))*(f(4) - feq(4));
01421     M2 += (f(5) - feq(5))*(f(5) - feq(5));
01422     M2 += (f(7) - feq(7))*(f(7) - feq(7));
01423     M2 += (f(8) - feq(8))*(f(8) - feq(8));
01424     M2 = std::sqrt(M2);
01425 
01426     DataType tau = DataType(1.0)/Omega(dt);
01427     DataType om_turb =  (DataType(2.0)/( tau
01428              + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
01429 
01430     return ( om_turb );
01431   }
01432   inline const int TurbulenceType() const { return turbulence; }
01433   inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
01434   inline const DataType& GasDensity() const { return rhop; }
01435   inline const DataType& GasViscosity() const { return nup; }
01436   inline const DataType& GasSpeedofSound() const { return csp; }
01437   inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
01438   { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
01439   // Quantities for output only
01440   inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
01441   inline virtual DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); } 
01442   inline virtual DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); } 
01443 
01444 protected:
01445   DataType cs2, cs22, cssq;
01446   DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
01447   DataType Cs_Smagorinsky, turbulence;
01448   int method[4], mdx[9], mdy[9];
01449 };
01450 
01451 
01452 #endif