amroc/lbm/LBMD2Q9Smag.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_SMAGORINSKY_H
00007 #define LBM_D2Q9_SMAGORINSKY_H
00008 
00016 #include "LBMD2Q9.h"
00017 #include <cfloat>
00018 
00019 #define LB_FACTOR 1.0e5
00020 
00035 template <class DataType>
00036 class LBMD2Q9Smag : public LBMD2Q9<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 
00047   enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro }; 
00048   enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall }; 
00049   enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall };
00050   enum TurbulenceModel { laminar, LES_Smagorinsky };
00051 
00052   LBMD2Q9() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.) {
00053     cs2  = DataType(1.)/DataType(3.);
00054     cs22 = DataType(2.)*cs2;
00055     cssq = DataType(2.)/DataType(9.);
00056     method[0] = 2;
00057     turbulence = laminar;
00058     mdx[0]=mdx[3]=mdx[6]=0; mdx[1]=mdx[4]=mdx[7]=1; mdx[2]=mdx[5]=mdx[8]=-1;
00059     mdy[0]=mdy[1]=mdy[2]=0; mdy[3]=mdy[4]=mdy[5]=1; mdy[6]=mdy[7]=mdy[8]=-1;
00060   }
00061  
00062   virtual ~LBMD2Q9() {}
00063 
00064   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00065     base::LocCtrl = Ctrl.getSubDevice(prefix+"D2Q9");
00066     RegisterAt(base::LocCtrl,"Method(1)",method[0]); 
00067     RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00068     RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00069     RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00070     RegisterAt(base::LocCtrl,"Gas_Cs",csp); 
00071     RegisterAt(base::LocCtrl,"Gas_nu",nup); 
00072     RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00073     RegisterAt(base::LocCtrl,"Gas_W",Wp);
00074     RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00075     RegisterAt(base::LocCtrl,"SpeedUp",S0);     
00076    }
00077 
00078   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {  
00079     base::SetupData(gh,ghosts);
00080     if (rhop>0.) R0 = rhop;
00081     if (csp>0.) {
00082       cs2p = csp*csp;
00083       SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00084     }
00085     std::cout << "D2Q9: Gas_rho=" << rhop << "   Gas_Cs=" << csp 
00086               << "   Gas_nu=" << nup << std::endl; 
00087     std::cout << "D2Q9: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0 
00088               << "   R0=" << R0 << std::endl; 
00089     WriteInit();
00090   }
00091 
00092   virtual void WriteInit() const {
00093     int me = MY_PROC; 
00094     if (me == VizServer) {
00095       std::ofstream ofs("chem.dat", std::ios::out);
00096       ofs << "RU         " << Rp << std::endl;
00097       ofs << "PA         " << BasePressure() << std::endl;
00098       ofs << "Sp(01)     Vicosity" << std::endl;
00099       ofs << "W(01)      " << Wp << std::endl;
00100       ofs << "Sp(02)     |Velocity|" << std::endl;
00101       ofs << "W(02)      " << 1.0 << std::endl;
00102       ofs << "Sp(03)     Vorticity" << std::endl;
00103       ofs << "W(03)      " << 1.0 << std::endl;
00104       ofs << "Sp(04)     |Vorticity|" << std::endl;
00105       ofs << "W(04)      " << 1.0 << std::endl;
00106      ofs.close();
00107     }    
00108   }
00109 
00110   inline virtual MacroType MacroVariables(const MicroType &f) const {
00111     MacroType q;
00112     q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8);
00113     q(1)=(f(1)-f(2)+f(4)-f(5)+f(7)-f(8))/q(0);
00114     q(2)=(f(3)+f(4)+f(5)-f(6)-f(7)-f(8))/q(0);
00115     return q;
00116   }
00117 
00118   inline virtual MicroType Equilibrium(const MacroType &q) const {
00119     MicroType feq;
00120 
00121     DataType rho = q(0);
00122     DataType u   = q(1);
00123     DataType v   = q(2);
00124 
00125     DataType ri, rt0, rt1, rt2;
00126     if (method[0]==0) {
00127       ri = DataType(1.);
00128       rt0 = R0 * DataType(4.) / DataType(9.);
00129       rt1 = R0 / DataType(9.);
00130       rt2 = R0 / DataType(36.);
00131     }
00132     if (method[0]==1) {
00133       ri = rho;
00134       rt0 = DataType(4.) / DataType(9.);
00135       rt1 = DataType(1.) / DataType(9.);
00136       rt2 = DataType(1.) / DataType(36.);
00137     }
00138     else if (method[0]==2) {
00139       ri = DataType(1.);
00140       rt0 = rho * DataType(4.) / DataType(9.);
00141       rt1 = rho / DataType(9.);
00142       rt2 = rho / DataType(36.);
00143     }
00144                
00145     DataType usq = u * u; 
00146     DataType vsq = v * v;
00147     DataType sumsq = (usq + vsq) / cs22;
00148     DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00149     DataType u2 = usq / cssq; 
00150     DataType v2 = vsq / cssq;
00151     DataType ui = u / cs2;
00152     DataType vi = v / cs2;
00153     DataType uv = ui * vi;
00154 
00155     feq(0) = rt0 * (ri - sumsq);
00156                
00157     feq(1) = rt1 * (ri - sumsq + u2 + ui);
00158     feq(2) = rt1 * (ri - sumsq + u2 - ui);
00159     feq(3) = rt1 * (ri - sumsq + v2 + vi);
00160     feq(6) = rt1 * (ri - sumsq + v2 - vi);
00161                
00162     feq(4) = rt2 * (ri + sumsq2 + ui + vi + uv);
00163     feq(5) = rt2 * (ri + sumsq2 - ui + vi - uv);
00164     feq(7) = rt2 * (ri + sumsq2 + ui - vi - uv);
00165     feq(8) = rt2 * (ri + sumsq2 - ui - vi + uv);    
00166 
00167     return feq;
00168   }
00169 
00170   inline virtual void Collision(MicroType &f, const DataType dt) const {
00171     MacroType q = MacroVariables(f);
00172     MicroType feq = Equilibrium(q);
00173 
00174     DataType omega;
00175     if (turbulence == laminar) 
00176       omega = Omega(dt);
00177     else if (turbulence == LES_Smagorinsky) 
00178       omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00179 
00180     f=(DataType(1.)-omega)*f + omega*feq;
00181   }     
00182 
00183   virtual int IncomingIndices(const int side, int indices[]) const {
00184     switch (side) {
00185     case base::Left:
00186       indices[0] = 1; indices[1] = 4; indices[2] = 7;  
00187       break;
00188     case base::Right:
00189       indices[0] = 2; indices[1] = 5; indices[2] = 8;  
00190       break;
00191     case base::Bottom:
00192       indices[0] = 3; indices[1] = 4; indices[2] = 5;  
00193       break;
00194     case base::Top:
00195       indices[0] = 6; indices[1] = 7; indices[2] = 8;  
00196       break;
00197     }
00198     return 3;
00199   }
00200 
00201   virtual int OutgoingIndices(const int side, int indices[]) const {
00202     switch (side) {
00203     case base::Left:
00204       indices[0] = 2; indices[1] = 5; indices[2] = 8;  
00205       break;
00206     case base::Right:
00207       indices[0] = 1; indices[1] = 4; indices[2] = 7;  
00208       break;
00209     case base::Bottom:
00210       indices[0] = 6; indices[1] = 7; indices[2] = 8;  
00211       break;
00212     case base::Top:
00213       indices[0] = 3; indices[1] = 4; indices[2] = 5;  
00214       break;
00215     }
00216     return 3;
00217   }
00218 
00219   // Revert just one layer of cells at the boundary 
00220   virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
00221                              const int side) const {
00222     int Nx = fvec.extents(0);
00223     int b = fvec.bottom();
00224     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00225     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00226     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00227     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00228     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00229     MicroType *f = (MicroType *)fvec.databuffer();
00230 
00231 #ifdef DEBUG_PRINT_FIXUP
00232     ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
00233       << fvec.bbox() << " using " << bb << "\n").flush();
00234 #endif
00235     
00236     switch (side) {
00237     case base::Left:
00238       for (register int j=mys; j<=mye; j++) {
00239         f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](2);
00240         f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](8);
00241         f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](5);
00242       }
00243       break;
00244     case base::Right:
00245       for (register int j=mys; j<=mye; j++) {
00246         f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](1);
00247         f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](7);
00248         f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](4);
00249       }
00250       break;
00251     case base::Bottom:
00252       for (register int i=mxs; i<=mxe; i++) {
00253         f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](6);
00254         f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](8);
00255         f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](7);
00256       }
00257       break;
00258     case base::Top:
00259       for (register int i=mxs; i<=mxe; i++) {
00260         f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](3);
00261         f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](5);
00262         f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](4);
00263       }
00264       break;
00265     }
00266   }   
00267 
00268   virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec, 
00269                          const BBox &bb, const double& dt) const {
00270     int Nx = fvec.extents(0);
00271     int b = fvec.bottom();
00272     MicroType *f = (MicroType *)fvec.databuffer();
00273     MicroType *of = (MicroType *)ovec.databuffer();
00274 
00275 #ifdef DEBUG_PRINT_FIXUP
00276     ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox() 
00277       << " using " << bb << "\n").flush();
00278 #endif
00279     
00280     assert (fvec.extents(0)==ovec.extents(0));
00281     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) &&
00282             fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1));
00283     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00284     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00285     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00286     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00287 
00288     for (register int j=mys; j<=mye; j++) 
00289       for (register int i=mxs; i<=mxe; i++) {
00290         for (register int n=0; n<base::NMicroVar(); n++) 
00291           f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00292         Collision(f[b+base::idx(i,j,Nx)],dt);
00293       }
00294   }  
00295 
00296   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00297                       const double& t, const double& dt, const int& mpass) const {
00298     // Make sure that the physical length and times are scaling consistently
00299     // into lattice quantities
00300     DCoords dx = base::GH().worldStep(fvec.stepsize());
00301     
00302     if (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR) {
00303       std::cerr << "LBMD2Q9 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00304                 << " to be equal." << std::endl << std::flush;
00305       std::exit(-1);
00306     }
00307     if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
00308       std::cerr << "LBMD2Q9 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0() 
00309                 << " to be equal." 
00310                 << "   dt=" << dt << "   TO=" << T0()
00311                 << std::endl << std::flush;
00312       std::exit(-1);
00313     }
00314 
00315     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00316     MicroType *f = (MicroType *)fvec.databuffer();
00317 
00318     // Streaming - update all cells
00319     for (register int j=Ny-1; j>=1; j--)
00320       for (register int i=0; i<=Nx-2; i++) {
00321         f[base::idx(i,j,Nx)](3) = f[base::idx(i,j-1,Nx)](3);
00322         f[base::idx(i,j,Nx)](5) = f[base::idx(i+1,j-1,Nx)](5);
00323       }
00324 
00325     for (register int j=Ny-1; j>=1; j--)
00326       for (register int i=Nx-1; i>=1; i--) {
00327         f[base::idx(i,j,Nx)](1) = f[base::idx(i-1,j,Nx)](1);
00328         f[base::idx(i,j,Nx)](4) = f[base::idx(i-1,j-1,Nx)](4);
00329       }
00330 
00331     for (register int j=0; j<=Ny-2; j++)
00332       for (register int i=Nx-1; i>=1; i--) {
00333         f[base::idx(i,j,Nx)](6) = f[base::idx(i,j+1,Nx)](6);
00334         f[base::idx(i,j,Nx)](7) = f[base::idx(i-1,j+1,Nx)](7);
00335       }
00336 
00337     for (register int j=0; j<=Ny-2; j++)
00338       for (register int i=0; i<=Nx-2; i++) {
00339         f[base::idx(i,j,Nx)](2) = f[base::idx(i+1,j,Nx)](2);
00340         f[base::idx(i,j,Nx)](8) = f[base::idx(i+1,j+1,Nx)](8);
00341       }    
00342 
00343     // Collision - only in the interior region
00344     int mxs = base::NGhosts(), mys = base::NGhosts();
00345     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00346     for (register int j=mys; j<=mye; j++)
00347       for (register int i=mxs; i<=mxe; i++)
00348         Collision(f[base::idx(i,j,Nx)],dt);
00349 
00350     return DataType(1.);
00351   }
00352 
00353    virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00354                            DataType* aux=0, const int naux=0, const int scaling=0) const {
00355     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00356     int mxs = base::NGhosts(), mys = base::NGhosts();
00357     int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00358     MicroType *f = (MicroType *)fvec.databuffer();
00359 
00360     if (type==ConstantMicro) {
00361       MicroType g;
00362       for (register int i=0; i<base::NMicroVar(); i++)
00363         g(i) = aux[i];
00364     
00365       for (register int j=mys; j<=mye; j++)
00366         for (register int i=mxs; i<=mxe; i++)
00367           f[base::idx(i,j,Nx)] = g;
00368     }
00369     
00370     else if (type==ConstantMacro) {
00371       MacroType q;
00372       if (scaling==base::Physical) {
00373         q(0) = aux[0]/R0; 
00374         q(1) = aux[1]/U0; 
00375         q(2) = aux[2]/U0;
00376       }
00377       else 
00378         for (register int i=0; i<base::NMacroVar(); i++)
00379           q(i) = aux[i];
00380       q(1) *= S0; q(2) *= S0;
00381     
00382       for (register int j=mys; j<=mye; j++)
00383         for (register int i=mxs; i<=mxe; i++)
00384           f[base::idx(i,j,Nx)] = Equilibrium(q);
00385     }
00386     
00387     else if (type==GasAtRest) {
00388       MacroType q;
00389       q(0) = GasDensity()/R0; 
00390       q(1) = DataType(0.); 
00391       q(2) = DataType(0.);
00392       for (register int j=mys; j<=mye; j++)
00393         for (register int i=mxs; i<=mxe; i++)
00394           f[base::idx(i,j,Nx)] = Equilibrium(q);
00395     }       
00396   }
00397 
00398   // Set just the one layer of ghost cells at the boundary 
00399   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00400                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00401     int Nx = fvec.extents(0);
00402     int b = fvec.bottom();
00403     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00404     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00405     int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00406     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00407     int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00408     MicroType *f = (MicroType *)fvec.databuffer();
00409     
00410     if (type==Symmetry) {
00411       switch (side) {
00412       case base::Left:
00413         for (register int j=mys; j<=mye; j++) {
00414           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j,Nx)](8);
00415           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00416           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j,Nx)](5);
00417         }
00418         break;
00419       case base::Right:
00420         for (register int j=mys; j<=mye; j++) {
00421           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j,Nx)](7);
00422           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00423           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j,Nx)](4);
00424         }
00425         break;
00426       case base::Bottom:
00427         for (register int i=mxs; i<=mxe; i++) {
00428           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i,mye+1,Nx)](8);
00429           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00430           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i,mye+1,Nx)](7);
00431         }
00432         break;
00433       case base::Top:
00434         for (register int i=mxs; i<=mxe; i++) {
00435           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i,mys-1,Nx)](5);
00436           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00437           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i,mys-1,Nx)](4);
00438         }
00439         break;
00440       }
00441     }
00442     
00443     else if (type==SlipWall) {
00444       switch (side) {
00445       case base::Left:
00446         for (register int j=mys; j<=mye; j++) {
00447           if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5); 
00448           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00449           if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00450         }
00451         break;
00452       case base::Right:
00453         for (register int j=mys; j<=mye; j++) {
00454           if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4); 
00455           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00456           if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00457         }
00458         break;
00459       case base::Bottom:
00460         for (register int i=mxs; i<=mxe; i++) {
00461           if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00462           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00463           if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00464         }
00465         break;
00466       case base::Top:
00467         for (register int i=mxs; i<=mxe; i++) {
00468           if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00469           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00470           if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00471         }
00472         break;
00473       }
00474     }
00475 
00476     else if (type==NoSlipWall) {
00477       switch (side) {
00478       case base::Left:
00479         for (register int j=mys; j<=mye; j++) {
00480           if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00481           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00482           if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00483         }
00484         break;
00485       case base::Right:
00486         for (register int j=mys; j<=mye; j++) {
00487           if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00488           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00489           if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00490         }
00491         break;  
00492       case base::Bottom:
00493         for (register int i=mxs; i<=mxe; i++) {
00494           if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00495           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00496           if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00497         }
00498         break;  
00499       case base::Top:
00500         for (register int i=mxs; i<=mxe; i++) {
00501           if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00502           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00503           if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00504         }
00505         break;
00506       }
00507     }
00508 
00509     else if (type==Inlet) {
00510       if (aux==0 || naux<=0) return;
00511       DataType a0=S0*aux[0], a1=0.;
00512       if (naux>1) a1 = S0*aux[1];
00513       if (scaling==base::Physical) {
00514         a0 /= U0;
00515         a1 /= U0;
00516       }
00517       switch (side) {
00518       case base::Left:
00519         for (register int j=mys; j<=mye; j++) {
00520           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00521           q(1) = a0;
00522           if (naux>1) q(2) = a1;
00523           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00524         }
00525         break;
00526       case base::Right:
00527         for (register int j=mys; j<=mye; j++) {
00528           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00529           q(1) = a0;
00530           if (naux>1) q(2) = a1;
00531           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00532         }
00533         break;
00534       case base::Bottom:
00535         for (register int i=mxs; i<=mxe; i++) {
00536           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00537           if (naux==1) 
00538             q(2) = a0;
00539           else {
00540             q(1) = a0;
00541             q(2) = a1;
00542           }
00543           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00544         }
00545         break;
00546       case base::Top:
00547         for (register int i=mxs; i<=mxe; i++) {
00548           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00549           if (naux==1) 
00550             q(2) = a0;
00551           else {
00552             q(1) = a0;
00553             q(2) = a1;
00554           }
00555           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00556         }
00557         break;
00558       }
00559     }
00560 
00561     else if (type==Outlet) {
00562       switch (side) {
00563       case base::Left:
00564         for (register int j=mys; j<=mye; j++) {
00565           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);       
00566           if (q(0)>0) f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)]/q(0);
00567         }
00568         break;
00569       case base::Right:
00570         for (register int j=mys; j<=mye; j++) {
00571           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00572           if (q(0)>0) f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)]/q(0);
00573         }
00574         break;
00575       case base::Bottom:
00576         for (register int i=mxs; i<=mxe; i++) {
00577           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00578           if (q(0)>0) f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)]/q(0);
00579         }
00580         break;
00581       case base::Top:
00582         for (register int i=mxs; i<=mxe; i++) {
00583           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00584           if (q(0)>0) f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)]/q(0);
00585         }
00586         break;
00587       }
00588     }
00589     
00590     // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)    
00591     else if (type==Pressure) {
00592       if (aux==0 || naux<=0) return;
00593       DataType a0=aux[0];
00594       if (scaling==base::Physical)
00595         a0 /= R0;
00596       switch (side) {
00597       case base::Left:
00598         for (register int j=mys; j<=mye; j++) {
00599           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00600           q(0) = a0;
00601           f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00602         }
00603         break;
00604       case base::Right:
00605         for (register int j=mys; j<=mye; j++) {
00606           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00607           q(0) = a0;
00608           f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00609         }
00610         break;
00611       case base::Bottom:
00612         for (register int i=mxs; i<=mxe; i++) {
00613           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00614           q(0) = a0;
00615           f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00616         }
00617         break;
00618       case base::Top:
00619         for (register int i=mxs; i<=mxe; i++) {
00620           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00621           q(0) = a0;
00622           f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00623         }
00624         break;
00625       }
00626     }
00627     
00628     else if (type==SlidingWall) {
00629       if (aux==0 || naux<=0) return;
00630       if (aux==0 || naux<=0) return;
00631       DataType a0=S0*aux[0];
00632       if (scaling==base::Physical)
00633         a0 /= U0;
00634       switch (side) {
00635       case base::Left:
00636         for (register int j=mys; j<=mye; j++) {
00637           MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00638           DataType rd1 = f[b+base::idx(mxe+1,j,Nx)](5), rd2 = f[b+base::idx(mxe+1,j,Nx)](8);
00639           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00640           DataType pw = DataType(1.)-qw;
00641           if (j<mye) f[b+base::idx(mxe,j+1,Nx)](7) = pw*rd1+qw*rd2;
00642           f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00643           if (j>mys) f[b+base::idx(mxe,j-1,Nx)](4) = qw*rd1+pw*rd2;
00644         }
00645         break;
00646       case base::Right:
00647         for (register int j=mys; j<=mye; j++) {
00648           MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00649           DataType rd1 = f[b+base::idx(mxs-1,j,Nx)](4), rd2 = f[b+base::idx(mxs-1,j,Nx)](7);
00650           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00651           DataType pw = DataType(1.)-qw;
00652           if (j<mye) f[b+base::idx(mxs,j+1,Nx)](8) = pw*rd1+qw*rd2;
00653           f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00654           if (j>mys) f[b+base::idx(mxs,j-1,Nx)](5) = qw*rd1+pw*rd2;
00655         }
00656         break;
00657       case base::Bottom:
00658         for (register int i=mxs; i<=mxe; i++) {
00659           MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00660           DataType rd1 = f[b+base::idx(i,mye+1,Nx)](7), rd2 = f[b+base::idx(i,mye+1,Nx)](8);
00661           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00662           DataType pw = DataType(1.)-qw;
00663           if (i<mxe) f[b+base::idx(i+1,mye,Nx)](5) = pw*rd1+qw*rd2;
00664           f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00665           if (i>mxs) f[b+base::idx(i-1,mye,Nx)](4) = qw*rd1+pw*rd2;
00666         }
00667         break;
00668       case base::Top:
00669         for (register int i=mxs; i<=mxe; i++) {
00670           MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00671           DataType rd1 = f[b+base::idx(i,mys-1,Nx)](4), rd2 = f[b+base::idx(i,mys-1,Nx)](5);
00672           DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00673           DataType pw = DataType(1.)-qw;
00674           if (i<mxe) f[b+base::idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00675           f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00676           if (i>mxs) f[b+base::idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00677         }
00678         break;
00679       }
00680     }
00681   }
00682 
00683   virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
00684                              const MicroType* f, const point_type* xc, const DataType* distance, 
00685                              const point_type* normal, DataType* aux=0, const int naux=0, 
00686                              const int scaling=0) const {    
00687     
00688     if (type==GFMNoSlipWall || type==GFMSlipWall) { 
00689       DataType fac = S0;
00690       if (scaling==base::Physical)
00691         fac /= U0;      
00692       for (int n=0; n<nc; n++) {
00693         MacroType q=MacroVariables(f[n]);
00694         DataType u=-q(1); 
00695         DataType v=-q(2);
00696         
00697         // Add boundary velocities if available
00698         if (naux>=2) {
00699           u += fac*aux[naux*n]; 
00700           v += fac*aux[naux*n+1];
00701         } 
00702         if (type==GFMNoSlipWall) {
00703           // Prescribe entire velocity vector in ghost cells
00704           q(1) += DataType(2.)*u;
00705           q(2) += DataType(2.)*v;
00706         }
00707         else if (type==GFMSlipWall) {
00708           // Prescribe only normal velocity vector in ghost cells
00709           DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
00710           q(1) += vl*normal[n](0);
00711           q(2) += vl*normal[n](1);
00712         }
00713         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
00714       } 
00715     }
00716 
00717     else if (type==GFMExtrapolation) {
00718       for (int n=0; n<nc; n++) {
00719         MacroType q=MacroVariables(f[n]);
00720         DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1);
00721         q(1) = vl*normal[n](0);
00722         q(2) = vl*normal[n](1);
00723         fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
00724       }
00725     }
00726   }
00727     
00728   virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
00729                       const int skip_ghosts=1) const {
00730     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00731     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
00732     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
00733     MicroType *f = (MicroType *)fvec.databuffer();
00734     DataType *work = (DataType *)workvec.databuffer();
00735     DCoords dx = base::GH().worldStep(fvec.stepsize());
00736     DataType dt_temp = dx(0)*S0/U0;
00737 
00738     if (cnt>=1 && cnt<=9) { 
00739       for (register int j=mys; j<=mye; j++)
00740         for (register int i=mxs; i<=mxe; i++) {
00741           MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
00742           switch(cnt) {
00743           case 1:
00744             work[base::idx(i,j,Nx)]=q(0)*R0;
00745             break;
00746           case 2:
00747             work[base::idx(i,j,Nx)]=q(1)*U0/S0;
00748             break;
00749           case 3:
00750             work[base::idx(i,j,Nx)]=q(2)*U0/S0;
00751             break;
00752           case 5:
00753             work[base::idx(i,j,Nx)]=TempEquation(U0/S0*U0/S0*R0*cs2*(q(0)-rhop/R0)+BasePressure()); 
00754             break;
00755           case 6:
00756             work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*cs2*(q(0)-rhop/R0)+BasePressure(); 
00757             break;
00758           case 8: {
00759             MicroType feq = Equilibrium(q);       
00760             work[base::idx(i,j,Nx)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt_temp),
00761                                                  GasSpeedofSound(),dt_temp );
00762             break;
00763           }
00764           case 9:
00765             work[base::idx(i,j,Nx)]=std::sqrt(q(1)*q(1)+q(2)*q(2))*U0/S0;
00766             break;
00767           }       
00768         }
00769     }
00770     else if (cnt>=10 && cnt<=11) {
00771       macro_grid_data_type qgrid(fvec.bbox());
00772       MacroType *q = (MacroType *)qgrid.databuffer();
00773       for (register int j=mys; j<=mye; j++)
00774         for (register int i=mxs; i<=mxe; i++) {
00775           q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
00776           q[base::idx(i,j,Nx)](1)*=U0/S0;
00777           q[base::idx(i,j,Nx)](2)*=U0/S0;
00778         }
00779       
00780       if (skip_ghosts==0) {
00781         mxs+=1; mxe-=1; mys+=1; mye-=1; 
00782       }
00783 
00784       for (register int j=mys; j<=mye; j++)
00785         for (register int i=mxs; i<=mxe; i++) {
00786           work[base::idx(i,j,Nx)]=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(2.*dx(0))-
00787             (q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(2.*dx(1));
00788           if (cnt==11) work[base::idx(i,j,Nx)]=std::abs(work[base::idx(i,j,Nx)]);
00789         }
00790     }
00791   }
00792 
00793   virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
00794                      const int skip_ghosts=1) const {
00795     int Nx = fvec.extents(0), Ny = fvec.extents(1);
00796     int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
00797     int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
00798     MicroType *f = (MicroType *)fvec.databuffer();
00799     DataType *work = (DataType *)workvec.databuffer();
00800 
00801     for (register int j=mys; j<=mye; j++)
00802       for (register int i=mxs; i<=mxe; i++) {
00803         switch(cnt) {
00804         case 1:
00805           f[base::idx(i,j,Nx)](0)=work[base::idx(i,j,Nx)]/R0;
00806           break;
00807         case 2:
00808           f[base::idx(i,j,Nx)](1)=work[base::idx(i,j,Nx)]/(U0/S0);
00809           break;
00810         case 3:
00811           f[base::idx(i,j,Nx)](2)=work[base::idx(i,j,Nx)]/(U0/S0);
00812           f[base::idx(i,j,Nx)] = Equilibrium((const MacroType &)f[base::idx(i,j,Nx)]);
00813           break;
00814         }
00815       }
00816   }
00817 
00818   virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time, 
00819                     const int verbose) const {
00820     int Nx = fvec.extents(0);
00821     int b = fvec.bottom();
00822     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00823     BBox bbmax = grow(fvec.bbox(),-base::NGhosts());    
00824     int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
00825     int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
00826     int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
00827     int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
00828     MicroType *f = (MicroType *)fvec.databuffer();
00829 
00830     int result = 1;
00831     for (register int j=mys; j<=mye; j++)
00832       for (register int i=mxs; i<=mxe; i++)
00833         for (register int l=0; l<base::NMicroVar(); l++) 
00834           if (!(f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l)>-1.e37)) {
00835             result = 0;
00836             if (verbose) 
00837               std::cerr << "D2Q9-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << ")(" << l << ")=" 
00838                         << f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l) << " " << fvec.bbox() << std::endl;
00839           }
00840     return result;
00841   }  
00842         
00843   virtual int NMethodOrder() const { return 2; } 
00844 
00845   inline const DataType& L0() const { return base::LengthScale(); }
00846   inline const DataType& T0() const { return base::TimeScale(); }
00847   inline void SetDensityScale(const DataType r0) { R0 = r0; }
00848   inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
00849   inline void SetSpeedUp(const DataType s0) { S0 = s0; }
00850   inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
00851 
00852   inline const DataType& DensityScale() const { return R0; }
00853   inline const DataType VelocityScale() const { return U0/S0; }
00854   // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
00855   inline const DataType& SpeedUp() const { return S0; }
00856 
00857   // Lattice quantities for the operator
00858   inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
00859   inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
00860 
00861   // Physical quantities for the operator
00862   inline void SetGas(DataType rho, DataType nu, DataType cs) { 
00863     rhop = rho; nup = nu; csp = cs; cs2p = cs*cs; 
00864     SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
00865   }
00866   inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
00867   inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q, 
00868                                               const DataType dt) const {
00869     DataType M2 = 0.0;
00870 
00871     M2 += (f(4) - feq(4))*(f(4) - feq(4));
00872     M2 += (f(5) - feq(5))*(f(5) - feq(5));
00873     M2 += (f(7) - feq(7))*(f(7) - feq(7));
00874     M2 += (f(8) - feq(8))*(f(8) - feq(8));
00875     M2 = std::sqrt(M2);
00876 
00877     DataType tau = DataType(1.0)/Omega(dt);
00878     DataType om_turb =  (DataType(2.0)/( tau
00879              + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
00880 
00881     return ( om_turb );
00882   }
00883   inline const int TurbulenceType() const { return turbulence; }
00884   inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
00885   inline const DataType& GasDensity() const { return rhop; }
00886   inline const DataType& GasViscosity() const { return nup; }
00887   inline const DataType& GasSpeedofSound() const { return csp; }
00888   inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
00889   { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
00890   // Quantities for output only
00891   inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
00892   inline virtual DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); } 
00893   inline virtual DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); } 
00894 
00895 protected:
00896   DataType cs2, cs22, cssq;
00897   DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
00898   DataType Cs_Smagorinsky, turbulence;
00899   int method[1], mdx[9], mdy[9];
00900 };
00901 
00902 
00903 #endif