00001
00002
00003
00004
00005
00006 #ifndef LBM_D2Q9_DR_H
00007 #define LBM_D2Q9_DR_H
00008
00017 #include "LBMD2Q9.h"
00018 #include <cfloat>
00019
00020 #define LB_FACTOR 1.0e5
00021
00037 typedef LBMD2Q9<DataType> LBMBaseType;
00038
00039 template <class DataType>
00040 class LBMD2Q9DR : public LBMBaseType {
00041 typedef LBMBaseType base;
00042 public:
00043 typedef typename base::vec_grid_data_type vec_grid_data_type;
00044 typedef typename base::grid_data_type grid_data_type;
00045 typedef typename base::MicroType MicroType;
00046 typedef typename base::MacroType MacroType;
00047 typedef GridData<MacroType,2> macro_grid_data_type;
00048 typedef typename base::SideName SideName;
00049 typedef typename base::point_type point_type;
00050
00051 LBMD2Q9DR() : base() {
00052 DRCoeff=0.005;
00053 DRFlow(0)=1.;
00054 DRFlow(1)=0.;
00055 DRFlow(2)=0.;
00056 fmeq = Equilibrium(DRFlow);
00057 Cs_Smagorinsky = DataType(0.2);
00058 }
00059
00060 virtual ~LBMD2Q9DR() {}
00061
00062 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00063 base::register_at(Ctrl,prefix);
00064 RegisterAt(base::LocCtrl,"DR_damping",DRCoeff);
00065 RegisterAt(base::LocCtrl,"DR_rho",DRFlow(0));
00066 RegisterAt(base::LocCtrl,"DR_u",DRFlow(1));
00067 RegisterAt(base::LocCtrl,"DR_v",DRFlow(2));
00068
00069 RegisterAt(base::LocCtrl,"DR_xmin",DRmin[0]);
00070 RegisterAt(base::LocCtrl,"DR_xmax",DRmax[0]);
00071 RegisterAt(base::LocCtrl,"DR_ymin",DRmin[1]);
00072 RegisterAt(base::LocCtrl,"DR_ymax",DRmax[1]);
00073
00074 RegisterAt(base::LocCtrl,"DR_Left",DRbc[0]);
00075 RegisterAt(base::LocCtrl,"DR_Right",DRbc[1]);
00076 RegisterAt(base::LocCtrl,"DR_Bottom",DRbc[2]);
00077 RegisterAt(base::LocCtrl,"DR_Top",DRbc[3]);
00078 }
00079
00080 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00081 base::SetupData(gh,ghosts);
00082 base::GH().wlb(WLB);
00083 base::GH().wub(WUB);
00084 DRminThk[0] = DRmin[0]-WLB[0];
00085 DRminThk[1] = DRmin[1]-WLB[1];
00086 DRmaxThk[0] = WUB[0]-DRmax[0];
00087 DRmaxThk[1] = WUB[1]-DRmax[1];
00088 for (register int i=0; i<2; i++) {
00089 if (DRminThk[i]<DataType(0.)) DRminThk[i] = DataType(0.);
00090 if (DRmaxThk[i]<DataType(0.)) DRmaxThk[i] = DataType(0.);
00091 }
00092 std::cout << "wholebbox " << base::GH().wholebbox() << WLB[0] <<","<<WLB[1] <<" " <<WUB[0]<<","<<WUB[1] <<")"<< std::endl;
00093 std::cout << "DR surrounds rectangle : min=(" << DRmin[0] <<","<< DRmin[1] <<") max=(" <<DRmax[0] <<","<<DRmax[1]<< ") [m]"<<std::endl;
00094 std::cout << "DR thickness : min=(" << DRminThk[0] <<","<< DRminThk[1] <<") max=(" <<DRmaxThk[0] <<","<<DRmaxThk[1]<<") [m]"<<std::endl;
00095 std::cout << "DR thickness : min=(" << DRminThk[0]/L0() <<","<< DRminThk[1]/L0() <<") max=(" <<DRmaxThk[0]/L0() <<","<<DRmaxThk[1]/L0()<<") lattice lengths"<<std::endl;
00096 std::cout << "DR type : min=(" << DRbc[0] <<","<< DRbc[2] <<") max=(" <<DRbc[1] <<","<<DRbc[3]<< ")"<<std::endl;
00097 (std::cout << "DR Damping Coefficient " << DRCoeff <<std::endl
00098 << "DR mean flow " << DRFlow << std::endl).flush();
00099 DRFlow(0) /= R0;
00100 DRFlow(1) /= U0/S0;
00101 DRFlow(2) /= U0/S0;
00102 fmeq = Equilibrium(DRFlow);
00103 }
00104
00105 virtual int DRInside(const DCoords lower, const DCoords upper) const {
00110 int position = 100;
00111 if ( ( lower(0) >= DRmin[0] ) && ( lower(0) <= DRmax[0] ) ) {
00112 if ( ( lower(1) >= DRmin[1] ) && ( lower(1) <= DRmax[1] ) ) {
00113 if ( ( upper(0) >= DRmin[0] ) && ( upper(0) <= DRmax[0] ) ) {
00114 if ( ( upper(1) >= DRmin[1] ) && ( upper(1) <= DRmax[1] ) ) {
00115 position = 0;
00116 }
00117 }
00118 }
00119 }
00120 return position;
00121 }
00122
00123 inline virtual void DRCollision(MicroType &f, const MicroType fmeql, const DataType damp, const DataType om) const {
00125 MacroType q = MacroVariables(f);
00126 MicroType feq = Equilibrium(q);
00127 MacroType Sigma = DeviatoricStress(f,feq,om);
00128 DataType S = Strain(q(0),Sigma,om,Cs_Smagorinsky);
00129 #if NUMPLUS >= 3
00130 if ( (turbulence==LES_SmagorinskyMemory || turbulence==LES_dynamicMemory) && f(10)>DataType(1.0) ) {
00131 Sigma = DeviatoricStress(f,feq,f(10));
00132 if (turbulence==LES_SmagorinskyMemory)
00133 S = Strain(q(0),Sigma,f(10),Cs_Smagorinsky);
00134 else
00135 S = Strain(q(0),Sigma,f(10),f(11));
00136 }
00137 #endif
00138 DataType strainDamp = damp*S/S0;
00139 for (register int qi=0; qi<9; qi++) {
00140 f(qi) -= strainDamp*(f(qi)-feq(qi)) + damp*(feq(qi)-fmeql(qi));
00141 }
00142 }
00143
00144 virtual void DRDamp(vec_grid_data_type &fvec) const {
00151 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00152 int mxs = base::NGhosts(), mys = base::NGhosts();
00153 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00154 MicroType *f = (MicroType *)fvec.databuffer();
00155 DCoords dx = base::GH().worldStep(fvec.stepsize());
00156 DataType dt_temp = dx(0)/VelocityScale();
00157 DataType om = Omega(dt_temp);
00158 DCoords wc;
00159 int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
00160 DataType WLB[2], WUB[2];
00161 base::GH().wlb(WLB);
00162 base::GH().wub(WUB);
00163 MacroType qtmp;
00164
00165 for (register int j=mys; j<=mye; j++) {
00166 lc_indx[1] = j*fvec.stepsize(1)+fvec.lower(1);
00167 for (register int i=mxs; i<=mxe; i++) {
00168 lc_indx[0] = i*fvec.stepsize(0)+fvec.lower(0);
00169 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00170 DataType d1 = (DRmaxThk[0]>DataType(0.) ? (wc(0)-DRmax[0])/DRmaxThk[0] : DataType(0.));
00171 DataType d2 = (DRminThk[0]>DataType(0.) ? (DRmin[0]-wc(0))/DRminThk[0] : DataType(0.));
00172 DataType d3 = (DRmaxThk[1]>DataType(0.) ? (wc(1)-DRmax[1])/DRmaxThk[1] : DataType(0.));
00173 DataType d4 = (DRminThk[1]>DataType(0.) ? (DRmin[1]-wc(1))/DRminThk[1] : DataType(0.));
00174 DataType dist = std::max(d1,std::max(d2,std::max(d3,d4)));
00175 if (dist>DataType(0.)) {
00176 DataType dampVal = DRCoeff*dist*dist;
00177 MicroType feqtmp = fmeq;
00178 if (DRbc[0]>1 || DRbc[1]>1 || DRbc[2]>1 || DRbc[3]>1) {
00179 qtmp = MacroVariables(f[base::idx(i,j,Nx)]);
00180 if (wc(0)-WLB[0]<DRminThk[0]) {
00181 switch(DRbc[0]) {
00182 case 0:
00183 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00184 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00185 qtmp(0) = DRFlow(0);
00186 break;
00187 case 2:
00188 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00189 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00190 break;
00191 case 3:
00192 qtmp(0) = DRFlow(0);
00193 break;
00194 case 4:
00195 case 5:
00196 if (DRbc[0]==4) {
00197 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00198 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00199 }
00200 else if (DRbc[0]==5)
00201 qtmp(0) = DRFlow(0);
00202 DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,Nx)]), MacroVariables(f[base::idx(i-1,j,Nx)]),
00203 MacroVariables(f[base::idx(i,j+1,Nx)]), MacroVariables(f[base::idx(i,j-1,Nx)]), dx);
00204 if (Qcrit>DataType(0.)) {
00205
00206 }
00207 else {
00208 dampVal *= DataType(0.5);
00209 }
00210 break;
00211 }
00212 }
00213 else if (WUB[0]-wc(0)<=DRmaxThk[0]) {
00214 switch(DRbc[1]) {
00215 case 2:
00216 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00217 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00218 break;
00219 case 3:
00220 qtmp(0) = DRFlow(0);
00221 break;
00222 case 4:
00223 case 5:
00224 if (DRbc[1]==4) {
00225 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00226 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00227 }
00228 else if (DRbc[1]==5)
00229 qtmp(0) = DRFlow(0);
00230 DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,Nx)]), MacroVariables(f[base::idx(i-1,j,Nx)]),
00231 MacroVariables(f[base::idx(i,j+1,Nx)]), MacroVariables(f[base::idx(i,j-1,Nx)]), dx);
00232 if (Qcrit>DataType(0.)) {
00233
00234 }
00235 else {
00236 dampVal *= DataType(0.5);
00237 }
00238 break;
00239 }
00240 }
00241 else if (wc(1)-WLB[1]<DRminThk[1]) {
00242 switch(DRbc[2]) {
00243 case 2:
00244 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00245 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00246 break;
00247 case 3:
00248 qtmp(0) = DRFlow(0);
00249 break;
00250 case 4:
00251 case 5:
00252 if (DRbc[2]==4) {
00253 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00254 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00255 }
00256 else if (DRbc[2]==5)
00257 qtmp(0) = DRFlow(0);
00258 DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,Nx)]), MacroVariables(f[base::idx(i-1,j,Nx)]),
00259 MacroVariables(f[base::idx(i,j+1,Nx)]), MacroVariables(f[base::idx(i,j-1,Nx)]), dx);
00260 if (Qcrit>DataType(0.)) {
00261
00262 }
00263 else {
00264 dampVal *= DataType(0.5);
00265 }
00266 break;
00267 }
00268 }
00269 else if (WUB[1]-wc(1)<=DRmaxThk[1]) {
00270 switch(DRbc[3]) {
00271 case 2:
00272 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00273 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00274 break;
00275 case 3:
00276 qtmp(0) = DRFlow(0);
00277 break;
00278 case 4:
00279 case 5:
00280 if (DRbc[3]==4) {
00281 qtmp(1) = DRFlow(1)*qtmp(0)/DRFlow(0);
00282 qtmp(2) = DRFlow(2)*qtmp(0)/DRFlow(0);
00283 }
00284 else if (DRbc[3]==5)
00285 qtmp(0) = DRFlow(0);
00286 DataType Qcrit = QcritLocal(MacroVariables(f[base::idx(i+1,j,Nx)]), MacroVariables(f[base::idx(i-1,j,Nx)]),
00287 MacroVariables(f[base::idx(i,j+1,Nx)]), MacroVariables(f[base::idx(i,j-1,Nx)]), dx);
00288 if (Qcrit>DataType(0.)) {
00289
00290 }
00291 else {
00292 dampVal *= DataType(0.5);
00293 }
00294 break;
00295 }
00296 }
00297 feqtmp = Equilibrium(qtmp);
00298 }
00299 if (DRbc[0]==1) {
00300 if (wc(0)-WLB[0]<DRminThk[1] && wc(1)-WLB[1]<DRminThk[1]) {
00301 dampVal *= (wc(0)-WLB[0])*(wc(0)-WLB[0])/DRminThk[1]/DRminThk[1];
00302 }
00303 else if (wc(0)-WLB[0]<DRminThk[1] && WUB[1]-wc(1)<DRmaxThk[1]) {
00304 dampVal *= (wc(0)-WLB[0])*(wc(0)-WLB[0])/DRmaxThk[1]/DRmaxThk[1];
00305 }
00306 feqtmp = Equilibrium(qtmp);
00307 }
00308 if (DRbc[1]==1) {
00309 if (WUB[0]-wc(0)<DRminThk[1] && wc(1)-WLB[1]<=DRminThk[1]) {
00310 dampVal *= (WUB[0]-wc(0))*(WUB[0]-wc(0))/DRminThk[1]/DRminThk[1];
00311 }
00312 else if (WUB[0]-wc(0)<DRmaxThk[1] && WUB[1]-wc(1)<=DRmaxThk[1]) {
00313 dampVal *= (WUB[0]-wc(0))*(WUB[0]-wc(0))/DRmaxThk[1]/DRmaxThk[1];
00314 }
00315 feqtmp = Equilibrium(qtmp);
00316 }
00317 if (DRbc[2]==1) {
00318 if (wc(0)-WLB[0]<DRminThk[0] && wc(1)-WLB[1]<DRminThk[0]) {
00319 DataType miny = WLB[1]+DRminThk[0];
00320 DataType px = DataType(1.)-(wc(0)-WLB[0])/DRminThk[0];
00321 DataType py = (miny-wc(1))/DRminThk[0];
00322 dampVal = std::max(DRCoeff*dist*dist, DRCoeff*py*py );
00323 qtmp = MacroVariables(f[base::idx(i,j,Nx)]);
00324 qtmp(0) = DRFlow(0);
00325 qtmp(1) = px*DRFlow(1);
00326 qtmp(2) = px*DRFlow(2);
00327 }
00328 else if (WUB[0]-wc(0)<DRmaxThk[0] && wc(1)-WLB[1]<DRmaxThk[0]) {
00329 dampVal *= (wc(1)-WLB[1])*(wc(1)-WLB[1])/DRmaxThk[0]/DRmaxThk[0];
00330 }
00331 feqtmp = Equilibrium(qtmp);
00332 }
00333 if (DRbc[3]==1) {
00334 if (wc(0)-WLB[0]<DRminThk[0] && WUB[1]-wc(1)<DRminThk[0]) {
00335 DataType miny = WUB[1]-DRminThk[0];
00336 DataType px = DataType(1.)-(wc(0)-WLB[0])/DRminThk[0];
00337 DataType py = (wc(1)-miny)/DRminThk[0];
00338 dampVal = std::max(DRCoeff*dist*dist, DRCoeff*py*py );
00339 qtmp = MacroVariables(f[base::idx(i,j,Nx)]);
00340 qtmp(0) = DRFlow(0);
00341 qtmp(1) = px*DRFlow(1);
00342 qtmp(2) = px*DRFlow(2);
00343 }
00344 else if (WUB[0]-wc(0)<=DRmaxThk[0] && WUB[1]-wc(1)<=DRmaxThk[0]) {
00345 dampVal *= (WUB[1]+dx(1)-wc(1))*(WUB[1]+dx(1)-wc(1))/DRmaxThk[0]/DRmaxThk[0];
00346 }
00347 feqtmp = Equilibrium(qtmp);
00348 }
00349
00350 if (dampVal>DataType(0.)) {
00351 DRCollision(f[base::idx(i,j,Nx)],feqtmp,dampVal,om);
00352 }
00353 }
00354 }
00355 }
00356 }
00357
00358 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00359 const double& t, const double& dt, const int& mpass) const {
00360 base::Step(fvec,ovec,Flux,t,dt,mpass);
00361
00362
00363 DCoords lower = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
00364 DCoords upper = base::GH().worldCoords(fvec.upper(), fvec.stepsize());
00365 if (DRInside(lower,upper)!=0) {
00366 DRDamp(fvec);
00367 }
00368
00369 return DataType(1.);
00370 }
00371
00372 virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00373 DataType* aux=0, const int naux=0, const int scaling=0) const {
00374 int Nx = fvec.extents(0);
00375 int b = fvec.bottom();
00376 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00377 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00378 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00379 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00380 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00381 MicroType *f = (MicroType *)fvec.databuffer();
00382
00386 if (type==80) {
00387 DCoords wc;
00388 int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
00389 DataType WLB[2], WUB[2];
00390 base::GH().wlb(WLB);
00391 base::GH().wub(WUB);
00392 DataType slip = 0.;
00393 DataType lxs, lxe, lys, lye;
00394 DataType min[2], max[2];
00395 lxs = DRmin[0] - WLB[0]; lxe = WUB[0] - DRmax[0];
00396 lys = DRmin[1] - WLB[1]; lye = WUB[1] - DRmax[1];
00397 min[0] = DRmin[0]; max[0] = DRmax[0];
00398 min[1] = DRmin[1]; max[1] = DRmax[1];
00399
00400 if (aux==0 || naux<=0) return;
00401 DataType rho=aux[0];
00402 DataType a0=S0*aux[1];
00403 DataType a1=S0*aux[2];
00404 if (scaling==base::Physical) {
00405 rho /= R0;
00406 a0 /= U0;
00407 a1 /= U0;
00408 }
00409 MacroType q, q1, q2, qn;
00410 DataType rhod=DataType(1.0);
00411 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
00412 q(0) = rho;
00413 q(1) = a0;
00414 q(2) = a1;
00415 switch (side) {
00416 case base::Left:
00417 for (register int j=mys; j<=mye; j++) {
00418 lc_indx[1] = j*fvec.stepsize(1);
00419 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00420 q1 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00421 q2 = MacroVariables(f[b+base::idx(mxe+2,j,Nx)]);
00422 if (naux==1) {
00423 q(1) = q1(1); q(2) = q1(2);
00424 MacroType dqdn = NormalDerivative(q,q1,q2);
00425 if (EquilibriumType()!=3)
00426 rhod = q(0);
00427 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
00428 qn(0) = q(0) - c1oCsSq2*L(0);
00429 qn(1) = q1(1) + c1o2cs*L(0)/rhod;
00430 qn(2) = q1(2) - L(1);
00431 }
00432 else if (naux==2) {
00433 q(0) = q1(0); q(2) = q1(2);
00434 MacroType dqdn = NormalDerivative(q,q1,q2);
00435 if (EquilibriumType()!=3)
00436 rhod = q(0);
00437 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
00438 qn(0) = q1(0) - c1oCsSq2*L(0);
00439 qn(1) = q(1) + c1o2cs*L(0)/rhod;
00440 qn(2) = q1(2) - L(1);
00441 }
00442 else if (naux==3) {
00443 q(0) = q1(0);
00444 MacroType dqdn = NormalDerivative(q,q1,q2);
00445 if (EquilibriumType()!=3)
00446 rhod = q(0);
00447 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
00448 qn(0) = q1(0) - c1oCsSq2*L(0);
00449 qn(1) = q(1) + c1o2cs*L(0)/rhod;
00450 qn(2) = q(2) - L(1);
00451 }
00452 DataType d=(WUB[1]-wc(1));
00453 if (d<=DRminThk[0]) {
00454 qn(0) = rho;
00455 }
00456 d=(wc(1)-WLB[1]);
00457 if (d<=DRminThk[0]) {
00458 qn(0) = rho;
00459 }
00460
00461 f[b+base::idx(mxe,j,Nx)] = Equilibrium(qn);
00462 }
00463 break;
00464 case base::Right:
00465 for (register int j=mys; j<=mye; j++) {
00466 q1 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00467 q2 = MacroVariables(f[b+base::idx(mxs-2,j,Nx)]);
00468 if (naux==1) {
00469 q(1) = q1(1); q(2) = q1(2);
00470 MacroType dqdn = NormalDerivative(q,q1,q2);
00471 if (EquilibriumType()!=3)
00472 rhod = q(0);
00473 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
00474 qn(0) = q(0) - c1oCsSq2*L(0);
00475 qn(1) = q1(1) + c1o2cs*L(0)/rhod;
00476 qn(2) = q1(2) - L(1);
00477 }
00478 else if (naux==2) {
00479 q(0) = q1(0); q(2) = q1(2);
00480 MacroType dqdn = NormalDerivative(q,q1,q2);
00481 if (EquilibriumType()!=3)
00482 rhod = q(0);
00483 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
00484 qn(0) = q1(0) - c1oCsSq2*L(0);
00485 qn(1) = q(1) + c1o2cs*L(0)/rhod;
00486 qn(2) = q1(2) - L(1);
00487 }
00488 else if (naux==3) {
00489 q(0) = q1(0);
00490 MacroType dqdn = NormalDerivative(q,q1,q2);
00491 if (EquilibriumType()!=3)
00492 rhod = q(0);
00493 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
00494 qn(0) = q1(0) - c1oCsSq2*L(0);
00495 qn(1) = q(1) + c1o2cs*L(0)/rhod;
00496 qn(2) = q(2) - L(1);
00497 }
00498 f[b+base::idx(mxs,j,Nx)] = Equilibrium(qn);
00499 }
00500 break;
00501 case base::Bottom:
00502 for (register int i=mxs; i<=mxe; i++) {
00503 q1 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00504 q2 = MacroVariables(f[b+base::idx(i,mye+2,Nx)]);
00505 if (naux==1) {
00506 q(1) = q1(1); q(2) = q1(2);
00507 MacroType dqdn = NormalDerivative(q,q1,q2);
00508 if (EquilibriumType()!=3)
00509 rhod = q(0);
00510 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
00511 qn(0) = q(0) - c1oCsSq2*L(0);
00512 qn(1) = q1(1) - L(1);
00513 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
00514 }
00515 else if (naux==2) {
00516 q(0) = q1(0); q(1) = q1(1);
00517 MacroType dqdn = NormalDerivative(q,q1,q2);
00518 if (EquilibriumType()!=3)
00519 rhod = q(0);
00520 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
00521 qn(0) = q1(0) - c1oCsSq2*L(0);
00522 qn(1) = q1(1) - L(1);
00523 qn(2) = q(2) + c1o2cs*L(0)/rhod;
00524 }
00525 else if (naux==3) {
00526 q(0) = q1(0);
00527 MacroType dqdn = NormalDerivative(q,q1,q2);
00528 if (EquilibriumType()!=3)
00529 rhod = q(0);
00530 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
00531 qn(0) = q1(0) - c1oCsSq2*L(0);
00532 qn(1) = q(1) - L(1);
00533 qn(2) = q(2) + c1o2cs*L(0)/rhod;
00534 }
00535 f[b+base::idx(i,mye,Nx)] = Equilibrium(qn);
00536 }
00537 break;
00538 case base::Top:
00539 for (register int i=mxs; i<=mxe; i++) {
00540 q1 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00541 q2 = MacroVariables(f[b+base::idx(i,mys-2,Nx)]);
00542 if (naux==1) {
00543 q(1) = q1(1); q(2) = q1(2);
00544 MacroType dqdn = NormalDerivative(q,q1,q2);
00545 if (EquilibriumType()!=3)
00546 rhod = q(0);
00547 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
00548 qn(0) = q(0) - c1oCsSq2*L(0);
00549 qn(1) = q1(1) - L(1);
00550 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
00551 }
00552 else if (naux==2) {
00553 q(0) = q1(0); q(1) = q1(1);
00554 MacroType dqdn = NormalDerivative(q,q1,q2);
00555 if (EquilibriumType()!=3)
00556 rhod = q(0);
00557 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
00558 qn(0) = q1(0) - c1oCsSq2*L(0);
00559 qn(1) = q1(1) - L(1);
00560 qn(2) = q(2) + c1o2cs*L(0)/rhod;
00561 }
00562 else if (naux==3) {
00563 q(0) = q1(0);
00564 MacroType dqdn = NormalDerivative(q,q1,q2);
00565 if (EquilibriumType()!=3)
00566 rhod = q(0);
00567 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
00568 qn(0) = q1(0) - c1oCsSq2*L(0);
00569 qn(1) = q(1) - L(1);
00570 qn(2) = q(2) + c1o2cs*L(0)/rhod;
00571 }
00572 f[b+base::idx(i,mys,Nx)] = Equilibrium(qn);
00573 }
00574 break;
00575 }
00576 }
00577
00582 else if (type==NoSlipWallEntranceExit) {
00583 const DataType pi = DataType(2.)*std::asin(DataType(1.));
00584 DCoords wc;
00585 int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
00586 DataType WLB[2], WUB[2];
00587 base::GH().wlb(WLB);
00588 base::GH().wub(WUB);
00589 DataType slip = 0.;
00590 DataType lxs, lxe, lys, lye;
00591 DataType min[2], max[2];
00592 if (naux==4) {
00593 lxs = aux[0] - WLB[0]; lxe = WUB[0] - aux[1];
00594 lys = aux[2] - WLB[1]; lye = WUB[1] - aux[3];
00595 min[0] = aux[0]; max[0] = aux[1];
00596 min[1] = aux[2]; max[1] = aux[3];
00597 }
00598 else {
00599 lxs = DRmin[0] - WLB[0]; lxe = WUB[0] - DRmax[0];
00600 lys = DRmin[1] - WLB[1]; lye = WUB[1] - DRmax[1];
00601 min[0] = DRmin[0]; max[0] = DRmax[0];
00602 min[1] = DRmin[1]; max[1] = DRmax[1];
00603 }
00604
00605 switch (side) {
00606 case base::Left:
00607 for (register int j=mys; j<=mye; j++) {
00608 lc_indx[1] = j*fvec.stepsize(1);
00609 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00610 if (wc(1)>=min[1] && wc(1)<=max[1]) {
00611 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00612 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00613 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00614 }
00615 else {
00616 if (wc(1)<WLB[1] || wc(1) >WUB[1])
00617 slip = DataType(0.);
00618 else if (wc(1)<min[1])
00619 slip = (wc(1) - WLB[1])/lys;
00620 else if (wc(1)>max[1])
00621 slip = (WUB[1] - wc(1))/lye;
00622 if (j>mys && j<mye) {
00623 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
00624 + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
00625 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
00626 + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
00627 }
00628 else if (j==mys) {
00629 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](5)
00630 + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
00631 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
00632 + slip*f[b+base::idx(mxe+1,j,Nx)](5);
00633 }
00634 else if (j==mye) {
00635 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
00636 + slip*f[b+base::idx(mxe+1,j,Nx)](8);
00637 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](8)
00638 + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
00639 }
00640 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00641 }
00642 }
00643 break;
00644 case base::Right:
00645 for (register int j=mys; j<=mye; j++) {
00646 lc_indx[1] = j*fvec.stepsize(1);
00647 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00648 if (wc(1)>=min[1] && wc(1)<=max[1]) {
00649 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00650 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00651 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00652 }
00653 else {
00654 if (wc(1)<WLB[1] || wc(1) >WUB[1])
00655 slip = DataType(0.);
00656 else if (wc(1)<min[1])
00657 slip = (wc(1) - WLB[1])/lys;
00658 else if (wc(1)>max[1])
00659 slip = (WUB[1] - wc(1))/lye;
00660 if (j>mys && j<mye) {
00661 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
00662 + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
00663 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
00664 + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
00665 }
00666 else if (j==mys) {
00667 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](4)
00668 + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
00669 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
00670 + slip*f[b+base::idx(mxs-1,j,Nx)](4);
00671 }
00672 else if (j==mye) {
00673 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
00674 + slip*f[b+base::idx(mxs-1,j,Nx)](7);
00675 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](7)
00676 + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
00677 }
00678 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00679 }
00680 }
00681 break;
00682 case base::Bottom:
00683 for (register int i=mxs; i<=mxe; i++) {
00684 lc_indx[0] = i*fvec.stepsize(0);
00685 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00686 if (wc(0)>=min[0] && wc(0)<=max[0]) {
00687 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00688 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00689 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00690 }
00691 else {
00692 if (wc(0)<=WLB[0] || wc(0)>=WUB[0])
00693 slip = DataType(0.);
00694 else if (wc(0)<min[0]) {
00695 DataType x = (min[0]-wc(0))/lxs;
00696 slip = DataType(0.5)*(std::sin(x*pi-pi*DataType(0.5))+DataType(1.));
00697 }
00698 else if (wc(0)>max[0])
00699 slip = (WUB[0] - wc(0))/lxe;
00700 if (i>=mxs && i<=mxe) {
00701 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,Nx)](7)
00702 + slip*f[b+base::idx(i+1,mye+1,Nx)](8);
00703 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,Nx)](8)
00704 + slip*f[b+base::idx(i-1,mye+1,Nx)](7);
00705 }
00706 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00707 MacroType q=MacroVariables(f[b+base::idx(i,mye,Nx)]);
00708 f[b+base::idx(i,mye,Nx)] *= rhop/(q(0)*R0);
00709 }
00710 }
00711 break;
00712 case base::Top:
00713 for (register int i=mxs; i<=mxe; i++) {
00714 lc_indx[0] = i*fvec.stepsize(0);
00715 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00716 if (wc(0)>=min[0] && wc(0)<=max[0]) {
00717 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00718 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00719 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00720 }
00721 else {
00722 if (wc(0)<=WLB[0] || wc(0)>=WUB[0])
00723 slip = DataType(0.);
00724 else if (wc(0)<min[0]) {
00725 DataType x = (min[0]-wc(0))/lxs;
00726 slip = DataType(0.5)*(std::sin(x*pi-pi*DataType(0.5))+DataType(1.));
00727 }
00728 else if (wc(0)>max[0])
00729 slip = (wc(0) - max[0])/lxe;
00730 if (slip<DataType(0.))
00731 slip = DataType(0.);
00732 if (i>=mxs && i<=mxe) {
00733 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,Nx)](5)
00734 + slip*f[b+base::idx(i-1,mys-1,Nx)](4);
00735 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,Nx)](4)
00736 + slip*f[b+base::idx(i+1,mys-1,Nx)](5);
00737 }
00738 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00739 MacroType q=MacroVariables(f[b+base::idx(i,mys,Nx)]);
00740 f[b+base::idx(i,mys,Nx)] *= rhop/(q(0)*R0);
00741 }
00742 }
00743 break;
00744 }
00745 }
00746 else
00747 base::BCStandard(fvec,bb,type,side,aux,naux,scaling);
00748
00749 }
00750
00751 inline DataType QcritLocal(const MacroType qxp, const MacroType qxm,
00752 const MacroType qyp, const MacroType qym, const DCoords dx) const {
00753 DataType dxux=(qxp(1)-qxm(1))/(2.*dx(0));
00754 DataType dxuy=(qxp(2)-qxm(2))/(2.*dx(0));
00755 DataType dyux=(qyp(1)-qym(1))/(2.*dx(1));
00756 DataType dyuy=(qyp(2)-qym(2))/(2.*dx(1));
00757 return DataType(-0.5)*(4.*dxux*dxux+8.*dyux*dxuy+4.*dyuy*dyuy);
00758 }
00759
00760 inline virtual const MacroType DRflow() const { return DRFlow; }
00761
00762 protected:
00763
00764 int WLBI[2], WUBI[2], DRbc[4];
00765 DataType WLB[2], WUB[2];
00766 DataType DRmin[2], DRmax[2];
00767 DataType DRminThk[2], DRmaxThk[2];
00768 DataType DRCoeff;
00769 MacroType DRFlow;
00770 MicroType fmeq;
00771 };
00772
00773
00774 #endif