00001
00002
00003
00004
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
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