amroc/lbm/LBMD2Q9_DR.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_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]) { // left
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]) { // right
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]) { // bottom
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]) { // top
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) { // left no-slip wall
00300             if (wc(0)-WLB[0]<DRminThk[1] && wc(1)-WLB[1]<DRminThk[1]) { // bottom
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]) { // top
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) { // right no-slip wall
00309             if (WUB[0]-wc(0)<DRminThk[1] && wc(1)-WLB[1]<=DRminThk[1]) { // bottom
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]) { // top
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) { // bottom no-slip wall
00318             if (wc(0)-WLB[0]<DRminThk[0] && wc(1)-WLB[1]<DRminThk[0]) { // left
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]) { // right
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) { // top no-slip wall
00334             if (wc(0)-WLB[0]<DRminThk[0] && WUB[1]-wc(1)<DRminThk[0]) { // left
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]) { // right
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     // Damping Region
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) { // pressure specified
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) { // normal velocity component specified
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) { // normal and transverse velocity components specified
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