amroc/lbm/LBMD1Q3.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_D1Q3_H
00007 #define LBM_D1Q3_H
00008 
00016 #include "LBMBase.h"
00017 #include <cfloat>
00018 
00019 #define LB_FACTOR 1.0e5
00020 
00031 template <class DataType>
00032 class LBMD1Q3 : public LBMBase<Vector<DataType,3>, Vector<DataType,2>, 1> {
00033   typedef LBMBase<Vector<DataType,3>, Vector<DataType,2>, 1> base;
00034 public:
00035   typedef typename base::vec_grid_data_type vec_grid_data_type;
00036   typedef typename base::grid_data_type grid_data_type;
00037   typedef typename base::MicroType MicroType;
00038   typedef typename base::MacroType MacroType;
00039   typedef typename base::SideName SideName;
00040   typedef typename base::point_type point_type;
00041 
00042   enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
00043   enum BCPredefined { NoSlipWall, Inlet, Outlet, Outlet2, Pressure };
00044 
00045   LBMD1Q3() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.) {
00046     cs2  = DataType(1.)/DataType(3.);
00047     cs22 = DataType(2.)*cs2;
00048     cssq = DataType(2.)/DataType(9.);
00049     method[0] = 2;
00050     mdx[0]=0; mdx[1]=1; mdx[2]=-1;
00051   }
00052 
00053   virtual ~LBMD1Q3() {}
00054 
00055   virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00056     base::LocCtrl = Ctrl.getSubDevice(prefix+"D1Q3");
00057     RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00058     RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00059     RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00060     RegisterAt(base::LocCtrl,"Gas_nu",nup);
00061     RegisterAt(base::LocCtrl,"SpeedUp",S0);
00062   }
00063 
00064   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00065     base::SetupData(gh,ghosts);
00066     if (rhop>0.) R0 = rhop;
00067     if (csp>0.) {
00068       cs2p = csp*csp;
00069       SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00070     }
00071     std::cout << "D1Q3: Gas_rho=" << rhop << "   Gas_Cs=" << csp
00072               << "   Gas_nu=" << nup << std::endl;
00073     std::cout << "D1Q3: L0=" << L0() << "   T0=" << T0() << "   U0=" << U0
00074               << "   R0=" << R0 << std::endl;
00075   }
00076 
00077   inline virtual MacroType MacroVariables(const MicroType &f) const {
00078     MacroType q;
00079     q(0)=f(0)+f(1)+f(2);
00080     q(1)=(f(1)-f(2))/q(0);
00081     return q;
00082   }
00083 
00084   inline virtual MicroType Equilibrium(const MacroType &q) const {
00085     MicroType feq;
00086 
00087     DataType rho = q(0);
00088     DataType u   = q(1);
00089 
00090     DataType ri, rt0, rt1;
00091     if (method[0]==0) {
00092       ri = DataType(1.);
00093       rt0 = R0 * DataType(4.) / DataType(6.);
00094       rt1 = R0 / DataType(6.);
00095     }
00096     if (method[0]==1) {
00097       ri = rho;
00098       rt0 = DataType(4.) / DataType(6.);
00099       rt1 = DataType(1.) / DataType(6.);
00100     }
00101     else if (method[0]==2) {
00102       ri = DataType(1.);
00103       rt0 = rho * DataType(4.) / DataType(6.);
00104       rt1 = rho / DataType(6.);
00105     }
00106 
00107     DataType usq = u * u;
00108     DataType sumsq = usq / cs22;
00109     DataType u2 = usq / cs2;
00110     DataType ui = u / cs2;
00111 
00112     feq(0) = rt0 * (ri - sumsq);
00113 
00114     feq(1) = rt1 * (ri + u2 + ui);
00115     feq(2) = rt1 * (ri + u2 - ui);
00116 
00117     if (method[0]==3) {
00118       DataType Cr = u*T0()/L0();
00119       DataType Cr2 = Cr*Cr;
00120       DataType Cs2_ = cs2;//p;//*VelocityScale();//p;
00121       DataType Dif = Cs2_*T0()*(1.0/Omega(T0())-DataType(0.5));//*VelocityScale();
00122       DataType Pes = std::fabs(u)*L0()/Dif*(1.0/Omega(T0())-DataType(0.5));
00123 
00124       ri = Cr/Pes;
00125       rt0 = rho * DataType(4.) / DataType(6.);
00126       rt1 = rho / DataType(6.);
00127       usq = u * u;
00128       sumsq = usq / cs22;
00129       u2 = usq / cs2;
00130       ui = u / cs2;
00131 
00132       //rho = R0;
00133       feq(0) = rho *(DataType(1.) - ( ri + Cr2 ));
00134       feq(1) = rho *(DataType(0.5) * (ri + Cr2 + Cr));
00135       feq(2) = rho *(DataType(0.5) * (ri + Cr2 - Cr));
00136     }
00137 
00138     return feq;
00139   }
00140 
00141   inline virtual void Collision(MicroType &f, const DataType dt) const {
00142     MacroType q = MacroVariables(f);
00143     MicroType feq = Equilibrium(q);
00144 
00145     DataType omega = Omega(dt);
00146     f=(DataType(1.)-omega)*f + omega*feq;
00147   }     
00148 
00149   virtual int IncomingIndices(const int side, int indices[]) const {
00150     switch (side) {
00151     case base::Left:
00152       indices[0] = 1;
00153       break;
00154     case base::Right:
00155       indices[0] = 2;
00156       break;
00157     }
00158     return 1;
00159   }
00160 
00161   virtual int OutgoingIndices(const int side, int indices[]) const {
00162     switch (side) {
00163     case base::Left:
00164       indices[0] = 2;  
00165       break;
00166     case base::Right:
00167       indices[0] = 1; 
00168       break;
00169     }
00170     return 1;
00171   }
00172 
00173   // Revert just one layer of cells at the boundary 
00174   virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb, 
00175                              const int side) const {
00176     int b = fvec.bottom();
00177     assert (fvec.stepsize(0)==bb.stepsize(0));
00178     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00179     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00180     MicroType *f = (MicroType *)fvec.databuffer();
00181 
00182 #ifdef DEBUG_PRINT_FIXUP
00183     ( comm_service::log() << "Reverse streaming at Side: " << side << " of " 
00184       << fvec.bbox() << " using " << bb << "\n").flush();
00185 #endif
00186     
00187     switch (side) {
00188     case base::Left:
00189       f[b+mxs](2) = f[b+mxs-1](2);
00190       break;
00191     case base::Right:
00192       f[b+mxe](1) = f[b+mxe+1](1);
00193       break;
00194     }
00195   }   
00196 
00197   virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec, 
00198                          const BBox &bb, const double& dt) const {
00199     int b = fvec.bottom();
00200     MicroType *f = (MicroType *)fvec.databuffer();
00201     MicroType *of = (MicroType *)ovec.databuffer();
00202 
00203 #ifdef DEBUG_PRINT_FIXUP
00204     ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox() 
00205       << " using " << bb << "\n").flush();
00206 #endif
00207     
00208     assert (fvec.extents(0)==ovec.extents(0));
00209     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(0)==ovec.stepsize(0));
00210     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00211     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00212 
00213     for (register int i=mxs; i<=mxe; i++) {
00214       for (register int n=0; n<base::NMicroVar(); n++) 
00215         f[b+i](n) = of[b+i-mdx[n]](n);
00216       Collision(f[b+i],dt);
00217     }
00218   }  
00219 
00220   virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00221                       const double& t, const double& dt, const int& mpass) const {
00222     // Make sure that the physical length and times are scaling consistently
00223     // into lattice quantities
00224     DCoords dx = base::GH().worldStep(fvec.stepsize());
00225 
00226     if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR)   {
00227       std::cerr << "LBMD1Q3 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00228                 << " to be equal." << std::endl << std::flush;
00229       std::exit(-1);
00230     }
00231 
00232     int Nx = fvec.extents(0);
00233     MicroType *f = (MicroType *)fvec.databuffer();
00234 
00235     // Streaming - update all cells
00236     for (register int i=Nx-1; i>=1; i--)
00237       f[i](1) = f[i-1](1);
00238     for (register int i=0; i<=Nx-2; i++)
00239       f[i](2) = f[i+1](2);
00240     
00241     // Collision - only in the interior region
00242     int mxs = base::NGhosts();
00243     int mxe = Nx-base::NGhosts()-1;
00244     for (register int i=mxs; i<=mxe; i++)
00245       Collision(f[i],dt);
00246       
00247     return DataType(1.);
00248   }
00249 
00250   virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00251                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00252     int Nx = fvec.extents(0);
00253     int mxs = base::NGhosts();
00254     int mxe = Nx-base::NGhosts()-1;
00255     MicroType *f = (MicroType *)fvec.databuffer();
00256 
00257     if (type==ConstantMicro) {
00258       MicroType g;
00259       for (register int i=0; i<base::NMicroVar(); i++)
00260         g(i) = aux[i];
00261 
00262       for (register int i=mxs; i<=mxe; i++)
00263         f[i] = g;
00264     }
00265     
00266     else if (type==ConstantMacro) {
00267       MacroType q;
00268       if (scaling==base::Physical) {
00269         q(0) = aux[0]/R0;
00270         q(1) = aux[1]/U0;
00271       }
00272       else
00273         for (register int i=0; i<base::NMacroVar(); i++)
00274           q(i) = aux[i];
00275       q(1) *= S0;
00276 
00277       for (register int i=mxs; i<=mxe; i++)
00278         f[i] = Equilibrium(q);
00279     }
00280     
00281     else if (type==GasAtRest) {
00282       MacroType q;
00283       q(0) = GasDensity()/R0;
00284       q(1) = DataType(0.);
00285       for (register int i=mxs; i<=mxe; i++)
00286         f[i] = Equilibrium(q);
00287     }
00288   }
00289 
00290   // Set just the one layer of ghost cells at the boundary
00291   virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00292                           DataType* aux=0, const int naux=0, const int scaling=0) const {
00293     int b = fvec.bottom();
00294     assert (fvec.stepsize(0)==bb.stepsize(0));
00295     int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00296     int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00297     MicroType *f = (MicroType *)fvec.databuffer();
00298     
00299     if (type==NoSlipWall) {
00300       switch (side) {
00301       case base::Left:
00302         f[b+mxe](1) = f[b+mxe+1](2);
00303         break;
00304       case base::Right:
00305         f[b+mxs](2) = f[b+mxs-1](1);
00306         break;
00307       }
00308     }
00309 
00310     else if (type==Inlet) {
00311       if (aux==0 || naux<=0) return;
00312       DataType a0=S0*aux[0], a1=0.;
00313       MacroType q;
00314       if (naux>1) a1 = S0*aux[1];
00315       if (scaling==base::Physical) {
00316         a0 /= U0;
00317         a1 /= U0;
00318       }
00319       switch (side) {
00320       case base::Left:
00321         q = MacroVariables(f[b+mxe+1]);
00322         q(1) = a0;
00323         f[b+mxe] = Equilibrium(q);
00324         break;
00325       case base::Right:
00326         q = MacroVariables(f[b+mxs-1]);
00327         q(1) = a0;
00328         f[b+mxs] = Equilibrium(q);
00329         break;
00330       }
00331     }
00332 
00333     else if (type==Outlet) {
00334       MacroType q;
00335       switch (side) {
00336       case base::Left:
00337         q = MacroVariables(f[b+mxe+1]);
00338         if (q(0)>0) f[b+mxe] = f[b+mxe+1]/q(0);
00339         break;
00340       case base::Right:
00341         q = MacroVariables(f[b+mxs-1]);
00342         if (q(0)>0) f[b+mxs] = f[b+mxs-1]/q(0);
00343         break;
00344       }
00345     }
00346     
00347     else if (type==Outlet2) {
00348       switch (side) {
00349       case base::Left:
00350         f[b+mxe] = -DataType(1.)/DataType(3.)*f[b+mxe+2]
00351             + DataType(4.)/DataType(3.)*f[b+mxe+1];
00352         break;
00353       case base::Right:
00354         f[b+mxs] = DataType(4.)/DataType(3.)*f[b+mxs-1]
00355             - DataType(1.)/DataType(3.)*f[b+mxs-2];
00356         break;
00357       }
00358     }
00359 
00360     // (Succi pg. 90 and Chen_Martinez_Mei_PhysFluids_8_2527)
00361     else if (type==Pressure) {
00362       if (aux==0 || naux<=0) return;
00363       DataType a0=aux[0];
00364       MacroType q;
00365       if (scaling==base::Physical)
00366         a0 /= R0;
00367       switch (side) {
00368       case base::Left:
00369         q = MacroVariables(f[b+mxe+1]);
00370         q(0) = a0;
00371         f[b+mxe] = Equilibrium(q);
00372         break;
00373       case base::Right:
00374         q = MacroVariables(f[b+mxs-1]);
00375         q(0) = a0;
00376         f[b+mxs] = Equilibrium(q);
00377         break;
00378       }
00379     }
00380   }
00381 
00382   virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
00383                              const MicroType* f, const point_type* xc, const DataType* distance,
00384                              const point_type* normal, DataType* aux=0, const int naux=0,
00385                              const int scaling=0) const {
00386     DataType fac = S0;
00387     if (scaling==base::Physical)
00388       fac /= U0;
00389     
00390     for (int n=0; n<nc; n++) {
00391       MacroType q=MacroVariables(f[n]);
00392       DataType u=-q(1);
00393       //DataType v=-q(2);
00394       
00395       // Add boundary velocities if available
00396       if (naux>=1) {
00397         u += fac*aux[naux*n];
00398       }
00399 
00400       // Construct normal velocity vector
00401       // Tangential velocity remains unchanged
00402       DataType vl = DataType(2.)*(u*normal[n](0) );//+v*normal[n](1));
00403       q(1) += vl*normal[n](0);
00404 
00405       fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
00406     }
00407   }
00408 
00409   virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
00410                       const int skip_ghosts=1) const {
00411     int Nx = fvec.extents(0);
00412     int mxs = base::NGhosts()*skip_ghosts;
00413     int mxe = Nx-base::NGhosts()*skip_ghosts-1;
00414     MicroType *f = (MicroType *)fvec.databuffer();
00415     DataType *work = (DataType *)workvec.databuffer();
00416 
00417     for (register int i=mxs; i<=mxe; i++) {
00418       MacroType q=MacroVariables(f[i]);
00419       switch(cnt) {
00420       case 1:
00421         work[i]=q(0)*R0;
00422         break;
00423       case 2:
00424         work[i]=q(1)*U0/S0;
00425         break;
00426       case 5:
00427         work[i]=U0/S0*U0/S0*R0*cs2*(q(0)-rhop/R0)+rhop*cs2p;
00428         break;
00429       }
00430     }
00431   }
00432 
00433   virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
00434                      const int skip_ghosts=1) const {
00435     int Nx = fvec.extents(0);
00436     int mxs = base::NGhosts()*skip_ghosts;
00437     int mxe = Nx-base::NGhosts()*skip_ghosts-1;
00438     MicroType *f = (MicroType *)fvec.databuffer();
00439     DataType *work = (DataType *)workvec.databuffer();
00440 
00441     for (register int i=mxs; i<=mxe; i++) {
00442       MacroType q=MacroVariables(f[i]);
00443       switch(cnt) {
00444       case 1:
00445         f[i](0)=work[i]/R0;
00446         break;
00447       case 2:
00448         f[i](1)=work[i]/(U0/S0);
00449         f[i] = Equilibrium((const MacroType &)f[i]);
00450         break;
00451       }
00452     }
00453   }
00454   
00455   virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
00456                     const int verbose) const {
00457     int b = fvec.bottom();
00458     assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00459     BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
00460     int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
00461     int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
00462     MicroType *f = (MicroType *)fvec.databuffer();
00463 
00464     int result = 1;
00465     for (register int i=mxs; i<=mxe; i++)
00466       for (register int l=0; l<base::NMicroVar(); l++)
00467         if (!(f[b+i-mdx[l]](l)>-1.e37)) {
00468           result = 0;
00469           if (verbose)
00470             std::cerr << "D1Q3-Check(): f(" << i-mdx[l] << ")(" << l << ")="
00471                       << f[b+i-mdx[l]](l) << " " << fvec.bbox() << std::endl;
00472         }
00473     return result;
00474   }
00475 
00476   virtual int NMethodOrder() const { return 2; }
00477 
00478   inline const DataType& L0() const { return base::LengthScale(); }
00479   inline const DataType& T0() const { return base::TimeScale(); }
00480   inline void SetDensityScale(const DataType r0) { R0 = r0; }
00481   inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
00482   inline void SetSpeedUp(const DataType s0) { S0 = s0; }
00483   inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
00484 
00485   inline const DataType& DensityScale() const { return R0; }
00486   inline const DataType VelocityScale() const { return U0/S0; }
00487   // Multiply all velocities and viscosity by S0 to speed up flow artificially, remove this scaling in output
00488   inline const DataType& SpeedUp() const { return S0; }
00489 
00490   // Lattice quantities for the operator
00491   inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
00492   inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
00493 
00494   // Physical quantities for the operator
00495   inline void SetGas(DataType rho, DataType nu, DataType cs) {
00496     rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
00497     SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
00498   }
00499   inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
00500   inline const DataType& GasDensity() const { return rhop; }
00501   inline const DataType& GasViscosity() const { return nup; }
00502   inline const DataType& GasSpeedofSound() const { return csp; }
00503   inline const DataType& GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
00504   { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
00505 
00506     protected:
00507   DataType cs2, cs22, cssq;
00508   DataType R0, U0, S0, rhop, csp, cs2p, nup;
00509   int method[1], mdx[3], mdy[3];
00510 };
00511 
00512 
00513 #endif