amroc/lbm/LBMD2Q9SupplementalBC.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_SUPPLEMENTALBOUNDARY_H
00007 #define LBM_SUPPLEMENTALBOUNDARY_H
00008 
00022 
00023 
00024 
00025 
00026 
00027 template <typename T> int sgn(T val) {
00028   return (T(0) < val) - (val < T(0));
00029 }
00030 
00031 inline DataType sgnl2(DataType a)  {
00032   if (a < 0) return -1.;
00033   else return 1.;
00034 }
00035 
00036 class BoundaryConditionsSupplemental : public LBMBoundaryConditions<LBMType,DIM> {
00037   typedef LBMBoundaryConditions<LBMType,DIM> base;
00038 public:
00039   typedef base::MicroType MicroType;
00040   typedef base::MacroType MacroType;
00041   BoundaryConditionsSupplemental(LBMType &lbm) : base(lbm) {}
00042   enum BCSupp { EqInletRamp=int(LBMType::NoSlipWallEntranceExit+1), oneSeventhInletRamp, parabolicInletRamp, slidingWallRamp };
00043 
00044   virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00045     base::SetupData(gh,ghosts);
00046     const double* bndry = base::GH().wholebndry();
00047     const int nbndry = base::GH().nbndry();
00048     DataType l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00049     if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00050     DataType omega = LBM().Omega(LBM().TimeScale());
00051     DataType u_p_avg=Side(0).aux[0];
00052     double Re_p=u_p_avg*l0_p/LBM().GasViscosity();
00053     double Re_lbm=(u_p_avg/LBM().VelocityScale())*(l0_p/LBM().LengthScale())/LBM().LatticeViscosity(omega);
00054     ( comm_service::log() << "BoundaryConditionsSupplemental\nl0_p=" << l0_p << "   u_p_avg=" << u_p_avg
00055                           << "   Re_p=" << Re_p << "   Re_lbm=" << Re_lbm
00056                           << "\n   omega=" << omega
00057                           << "   SpeedUp=" << LBM().SpeedUp()
00058                           << "   u_avg_LBM=" << u_p_avg/LBM().VelocityScale()
00059                           << "   SpeedRatio=" << u_p_avg/LBM().VelocityScale()*std::sqrt(DataType(3.))
00060                           << "   LBM().TO()=" << LBM().T0()
00061                           << "\n   Re_p-Re_lbm=" << Re_p-Re_lbm
00062                           << std::endl ).flush();
00063     assert(std::fabs(Re_p-Re_lbm)<DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR);
00064   }
00065 
00066   virtual void SetBndry(vec_grid_data_type &fvec, const int& level, const BBox &bb,
00067                         const int &dir, const double& time) {
00068 
00069     BBox wholebbox = base::GH().wholebbox();
00070     BBox inside = bb*coarsen(wholebbox,wholebbox.stepsize()/wholebbox.stepsize());
00071     if (!inside.empty()) {
00072       LBM().BCStandard(fvec,bb,LBMType::NoSlipWall,dir);
00073     }
00074     else if ( (dir==LBMType::Left && Side(0).Type<EqInletRamp) ||
00075               (dir==LBMType::Right && Side(1).Type<EqInletRamp) ||
00076               (dir==LBMType::Bottom && Side(2).Type<EqInletRamp) ||
00077               (dir==LBMType::Top && Side(3).Type<EqInletRamp) ) {
00078       base::SetBndry(fvec,level,bb,dir,time);
00079     }
00080     else {
00081       int Nx = fvec.extents(0);
00082       int b = fvec.bottom();
00083       assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00084       int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00085       int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00086       int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00087       int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00088       MicroType *f = (MicroType *)fvec.databuffer();
00089       MacroType q;
00090 
00092       if (dir==LBMType::Left && bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0) &&
00093           Side(0).Type==EqInletRamp) {
00094 
00095         DataType rho=Side(0).aux[0];
00096         DataType a0=Side(0).aux[1];
00097         DataType a1=Side(0).aux[2];
00098         DataType T_ramp=Side(0).aux[3];
00099 
00100         MacroType q, q_;
00101         DataType ramp;
00102 
00103         if (Side(0).Scaling==LBMType::Physical) {
00104           rho /= LBM().DensityScale();
00105           a0 /= LBM().VelocityScale();
00106           a1 /= LBM().VelocityScale();
00107         }
00108 
00109         q_(0) = rho;
00110         q_(1) = a0;
00111         q_(2) = a1;
00112 
00113         if ( time < T_ramp ) {
00114           ramp = time/T_ramp;
00115           q(0) = q_(0);
00116           q(1) = q_(1)*ramp;
00117           q(2) = q_(2)*ramp;
00118         }
00119         else
00120           q = q_;
00121 
00122         for (register int j=mys; j<=mye; j++) {
00123           q_ = LBM().MacroVariables(f[b+LBM().idx(mxe+1,j,Nx)]);
00124           q(0) = q_(0);
00125           f[b+LBM().idx(mxe,j,Nx)] = LBM().Equilibrium(q);
00126         }
00127       }
00128 
00130       else if (dir==LBMType::Left && bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0) &&
00131                Side(0).Type==oneSeventhInletRamp) {
00132         DataType rho=Side(0).aux[0];
00133         DataType a0=Side(0).aux[1];
00134         DataType a1=Side(0).aux[2];
00135         DataType elevation=Side(0).aux[3];
00136 
00137         MacroType q, q_;
00138         DataType profile, height;
00139         DCoords dx = base::GH().worldStep(fvec.stepsize()), wc;
00140         int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
00141 
00142         if (Side(0).Scaling==LBMType::Physical) {
00143           rho /= LBM().DensityScale();
00144           a0 /= LBM().VelocityScale();
00145           a1 /= LBM().VelocityScale();
00146         }
00147 
00148         q_(0) = rho;
00149         q_(1) = a0;
00150         q_(2) = a1;
00151 
00152         for (register int j=mys; j<=mye; j++) {
00153           lc_indx[1] = j*fvec.stepsize(1);
00154           wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00155           height = wc(1)+dx(1)/2.;
00156           if ( height < 0.0) height = 0.0;
00157           if ( height < elevation ) {
00158             profile = std::pow((height/elevation),(1./7.));
00159             q(0) = q_(0);
00160             q(1) = q_(1)*profile;
00161             q(2) = q_(2)*profile;
00162           }
00163           else
00164             q = q_;
00165           f[b+LBM().idx(mxe,j,Nx)] = LBM().Equilibrium(q);
00166         }
00167 
00168       }
00169 
00171       else if (dir==LBMType::Left && bb.upper(0)+bb.stepsize(0)==wholebbox.lower(0) &&
00172                Side(0).Type==parabolicInletRamp) {
00173 
00174         DataType rho=Side(0).aux[0];
00175         DataType a0=Side(0).aux[1];
00176         DataType a1=Side(0).aux[2];
00177         DataType H=Side(0).aux[3];
00178         DataType T_ramp=Side(0).aux[4];
00179 
00180         MacroType qp, q_;
00181         DataType ramp;
00182 
00183         if (Side(0).Scaling==LBMType::Physical) {
00184           rho /= LBM().DensityScale();
00185           a0 /= LBM().VelocityScale();
00186           a1 /= LBM().VelocityScale();
00187         }
00188 
00189         q_(0) = rho;
00190         q_(1) = a0;
00191         q_(2) = a1;
00192 
00193         if ( time < T_ramp ) {
00194           ramp = time/T_ramp;
00195           qp(0) = q_(0);
00196           qp(1) = q_(1)*ramp;
00197           qp(2) = q_(2)*ramp;
00198         }
00199         else
00200           qp = q_;
00201 
00202         const double* bndry = base::GH().wholebndry();
00203         const int nbndry = base::GH().nbndry();
00204         DataType l0_p=bndry[nbndry*(1+2*1)]-bndry[nbndry*(0+2*1)];
00205         if (nbndry>1) l0_p=bndry[nbndry*(1+2*1)]-bndry[1+nbndry*(1+2*1)];
00206         H = l0_p;
00207 
00208         DataType profile, height;
00209         DCoords dx = base::GH().worldStep(fvec.stepsize()), wc;
00210         int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
00211 
00212         for (register int j=mys; j<=mye; j++) {
00213           q = LBM().MacroVariables(f[b+LBM().idx(mxe+1,j,Nx)]);
00214           lc_indx[1] = j*fvec.stepsize(1);
00215           wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00216           height = wc(1)+dx(1)/2.+H/2.;
00217           profile = 4.*height*(H-height)/(H*H);
00218           //q(0) = rho;
00219           q(1) = qp(1)*profile;
00220           q(2) = qp(2)*profile;
00221           f[b+LBM().idx(mxe,j,Nx)] = LBM().Equilibrium(q);
00222         }
00223 
00224       }
00225 
00227       else if (dir==LBMType::Top && bb.lower(1)-wholebbox.stepsize(1)==wholebbox.upper(1) &&
00228                Side(3).Type==slidingWallRamp) {
00229 
00230         DataType rho=Side(3).aux[0];
00231         DataType a0=Side(3).aux[1];
00232         DataType a1=Side(3).aux[2];
00233         DataType T_ramp=Side(3).aux[3];
00234         MacroType q, q_;
00235         DataType ramp;
00236 
00237         if (Side(0).Scaling==LBMType::Physical) {
00238           rho /= LBM().DensityScale();
00239           a0 /= LBM().VelocityScale();
00240           a1 /= LBM().VelocityScale();
00241         }
00242 
00243         q_(0) = rho;
00244         q_(1) = a0;
00245         q_(2) = a1;
00246 
00247         if ( time < T_ramp ) {
00248           ramp = time/T_ramp;
00249           q(0) = q_(0);
00250           q(1) = q_(1)*ramp;
00251           q(2) = q_(2)*ramp;
00252         }
00253         else
00254           q = q_;
00255 
00256         for (register int i=mxs; i<=mxe; i++) {
00257           MacroType qs = LBM().MacroVariables(f[b+LBM().idx(i,mys-1,Nx)]);
00258           DataType rd1 = f[b+LBM().idx(i,mys-1,Nx)](4), rd2 = f[b+LBM().idx(i,mys-1,Nx)](5);
00259           DataType qw = (qs(0)/DataType(6.)*q(1))/(rd1-rd2);
00260           DataType pw = DataType(1.)-qw;
00261           if (i<mxe) f[b+LBM().idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00262           f[b+LBM().idx(i,mys,Nx)](6) = f[b+LBM().idx(i,mys-1,Nx)](3);
00263           if (i>mxs) f[b+LBM().idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00264         }
00265       }
00266     }
00267 
00268   }
00269 
00270   inline DataType sgn(DataType a)  {
00271     if (a >= 0.) return 1.0;
00272     else return -1.0;
00273     assert(0);
00274     return -2.0;
00275   }
00276 
00277 };
00278 
00279 #endif