amroc/lbm/LBMD2Q9_Test.h

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