• VTF
  • FSI
  • AMROC
  • SFC
  • Motion
  • STLIB
  • Main Page
  • Related Pages
  • Classes
  • Files
  • File List
  • File Members

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