00001
00002
00003
00004
00005
00006 #ifndef LBM_D2Q9_H
00007 #define LBM_D2Q9_H
00008
00016 #include "LBMBase.h"
00017 #include <cfloat>
00018
00019 #define LB_FACTOR 1.0e5
00020
00035 #ifndef NUMPLUS
00036 #define NUMPLUS 0
00037 #endif
00038 #define NUMMICRODIST (9+NUMPLUS)
00039
00040 template <class DataType>
00041 class LBMD2Q9 : public LBMBase<Vector<DataType,NUMMICRODIST>, Vector<DataType,3>, 2> {
00042 typedef LBMBase<Vector<DataType,NUMMICRODIST>, Vector<DataType,3>, 2> base;
00043 public:
00044 typedef typename base::vec_grid_data_type vec_grid_data_type;
00045 typedef typename base::grid_data_type grid_data_type;
00046 typedef typename base::MicroType MicroType;
00047 typedef typename base::MacroType MacroType;
00048 typedef GridData<MacroType,2> macro_grid_data_type;
00049 typedef typename base::SideName SideName;
00050 typedef typename base::point_type point_type;
00051 typedef typename base::MacroType TensorType;
00052
00053 enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
00054 enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall, CharacteristicOutlet, CharacteristicInlet, NoSlipWallEntranceExit, Periodic, NoSlipWallNeq, VanDriest, CSMBC };
00055 enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall, GFMComplex, GFMComplex2 };
00056 enum TurbulenceModel { laminar, LES_Smagorinsky, LES_SmagorinskyMemory, LES_dynamic, LES_dynamicMemory, WALE, CSM };
00057
00058 LBMD2Q9() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.), mfp(1.) {
00059 cs2 = DataType(1.)/DataType(3.);
00060 cs22 = DataType(2.)*cs2;
00061 cssq = DataType(2.)/DataType(9.);
00062 method[0] = 2; method[1] = 0; method[2] = 0;
00063 turbulence = laminar;
00064 Cs_Smagorinsky = DataType(0.2);
00065 mdx[0]=mdx[3]=mdx[6]=0; mdx[1]=mdx[4]=mdx[7]=1; mdx[2]=mdx[5]=mdx[8]=-1;
00066 mdy[0]=mdy[1]=mdy[2]=0; mdy[3]=mdy[4]=mdy[5]=1; mdy[6]=mdy[7]=mdy[8]=-1;
00067 stressPath=0;
00068 }
00069
00070 virtual ~LBMD2Q9() {}
00071
00072 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00073 base::LocCtrl = Ctrl.getSubDevice(prefix+"D2Q9");
00074 RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00075 RegisterAt(base::LocCtrl,"Method(2)",method[1]);
00076 RegisterAt(base::LocCtrl,"Method(3)",method[2]);
00077 RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00078 RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00079 RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00080 RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00081 RegisterAt(base::LocCtrl,"Gas_nu",nup);
00082 RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00083 RegisterAt(base::LocCtrl,"Gas_W",Wp);
00084 RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00085 RegisterAt(base::LocCtrl,"SpeedUp",S0);
00086 RegisterAt(base::LocCtrl,"StressPath",stressPath);
00087 }
00088
00089 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00090 base::SetupData(gh,ghosts);
00091 if (rhop>0.) R0 = rhop;
00092 if (csp>0.) {
00093 cs2p = csp*csp;
00094 SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00095 if (gp==DataType(0.))
00096 mfp = GasViscosity()*std::sqrt(std::asin(1.)*DataType(1.4)/cs2p);
00097 else
00098 mfp = GasViscosity()*std::sqrt(std::asin(1.)*gp/cs2p);
00099 }
00100 std::cout << "D2Q9: Gas_rho=" << rhop << " Gas_Cs=" << csp
00101 << " Gas_nu=" << nup << std::endl;
00102 std::cout << "D2Q9: L0=" << L0() << " T0=" << T0() << " U0=" << U0
00103 << " R0=" << R0 << " Omega=" << Omega(T0())
00104 << std::endl;
00105 if ( TurbulenceType() == laminar ) {
00106 std::cout << "LBM Setup() Laminar" << std::endl;
00107 }
00108 if ( TurbulenceType() == LES_Smagorinsky || TurbulenceType() == LES_SmagorinskyMemory ) {
00109 std::cout << "LBM Setup() LES_Smagorinsky type "<<TurbulenceType()
00110 <<" : Smagorinsky Constant = " << SmagorinskyConstant() << std::endl;
00111 }
00112 if ( TurbulenceType() == LES_dynamic || TurbulenceType() == LES_dynamicMemory ) {
00113 std::cout << "LBM Setup() LES_dynamic type " <<TurbulenceType()<< std::endl;
00114 assert(base::NGhosts()>2);
00115 }
00116 if ( TurbulenceType() == LES_SmagorinskyMemory || TurbulenceType() == LES_dynamicMemory ) {
00117 if (NUMPLUS<3)
00118 std::cerr << "LBM Setup() LES type "<<TurbulenceType()
00119 <<" requires additional data structures. Use LBMD2Q9plus_LES.h or LBMD2Q9plus_DL.h\n";
00120 assert(NUMPLUS>=3);
00121 }
00122 if ( TurbulenceType() == WALE ) {
00123 std::cout << "LBM Setup() WALE" << std::endl;
00124 }
00125 if ( TurbulenceType() == CSM ) {
00126 std::cout << "LBM Setup() CSM mfp="<<mfp << std::endl;
00127 }
00128 WriteInit();
00129 }
00130
00131 virtual void WriteInit() const {
00132 int me = MY_PROC;
00133 if (me == VizServer) {
00134 std::ofstream ofs("chem.dat", std::ios::out);
00135 ofs << "RU " << Rp << std::endl;
00136 ofs << "PA " << BasePressure() << std::endl;
00137 ofs << "Sp(01) Vicosity" << std::endl;
00138 ofs << "W(01) " << Wp << std::endl;
00139 ofs << "Sp(02) |Velocity|" << std::endl;
00140 ofs << "W(02) " << 1. << std::endl;
00141 ofs << "Sp(03) vz" << std::endl;
00142 ofs << "W(03) " << 1. << std::endl;
00143 ofs << "Sp(04) |v|" << std::endl;
00144 ofs << "W(04) " << 1. << std::endl;
00145 ofs << "Sp(05) |Stress|" << std::endl;
00146 ofs << "W(05) " << 1. << std::endl;
00147 ofs << "Sp(06) |DevStress|" << std::endl;
00148 ofs << "W(06) " << 1. << std::endl;
00149 ofs << "Sp(07) StrainRate" << std::endl;
00150 ofs << "W(07) " << 1. << std::endl;
00151 ofs << "Sp(08) qcrit" << std::endl;
00152 ofs << "W(08) " << 1. << std::endl;
00153 ofs << "Sp(09) dxx" << std::endl;
00154 ofs << "W(09) " << 1. << std::endl;
00155 ofs << "Sp(10) dxy" << std::endl;
00156 ofs << "W(10) " << 1. << std::endl;
00157 ofs << "Sp(11) dyy" << std::endl;
00158 ofs << "W(11) " << 1. << std::endl;
00159 ofs.close();
00160 }
00161 }
00162
00163 inline virtual MacroType MacroVariables(const MicroType &f) const {
00164 MacroType q;
00165 q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8);
00166 if (method[0]!=3) {
00167 q(1)=(f(1)-f(2)+f(4)-f(5)+f(7)-f(8))/q(0);
00168 q(2)=(f(3)+f(4)+f(5)-f(6)-f(7)-f(8))/q(0);
00169 }
00170 else {
00171 q(1)=(f(1)-f(2)+f(4)-f(5)+f(7)-f(8))/(rhop/R0);
00172 q(2)=(f(3)+f(4)+f(5)-f(6)-f(7)-f(8))/(rhop/R0);
00173 }
00174 return q;
00175 }
00176
00177 inline virtual MicroType Equilibrium(const MacroType &q) const {
00178 MicroType feq;
00179
00180 DataType rho = q(0);
00181 DataType u = q(1);
00182 DataType v = q(2);
00183
00184 DataType ri=0., rt0=0., rt1=0., rt2=0.;
00185 if (method[0]<3) {
00186 if (method[0]==0) {
00187 ri = DataType(1.);
00188 rt0 = R0 * DataType(4.) / DataType(9.);
00189 rt1 = R0 / DataType(9.);
00190 rt2 = R0 / DataType(36.);
00191 }
00192 if (method[0]==1) {
00193 ri = rho;
00194 rt0 = DataType(4.) / DataType(9.);
00195 rt1 = DataType(1.) / DataType(9.);
00196 rt2 = DataType(1.) / DataType(36.);
00197 }
00198 else if (method[0]==2) {
00199 ri = DataType(1.);
00200 rt0 = rho * DataType(4.) / DataType(9.);
00201 rt1 = rho / DataType(9.);
00202 rt2 = rho / DataType(36.);
00203 }
00204
00205 DataType usq = u * u;
00206 DataType vsq = v * v;
00207 DataType sumsq = (usq + vsq) / cs22;
00208 DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00209 DataType u2 = usq / cssq;
00210 DataType v2 = vsq / cssq;
00211 DataType ui = u / cs2;
00212 DataType vi = v / cs2;
00213 DataType uv = ui * vi;
00214
00215 feq(0) = rt0 * (ri - sumsq);
00216
00217 feq(1) = rt1 * (ri - sumsq + u2 + ui);
00218 feq(2) = rt1 * (ri - sumsq + u2 - ui);
00219 feq(3) = rt1 * (ri - sumsq + v2 + vi);
00220 feq(6) = rt1 * (ri - sumsq + v2 - vi);
00221
00222 feq(4) = rt2 * (ri + sumsq2 + ui + vi + uv);
00223 feq(5) = rt2 * (ri + sumsq2 - ui + vi - uv);
00224 feq(7) = rt2 * (ri + sumsq2 + ui - vi - uv);
00225 feq(8) = rt2 * (ri + sumsq2 - ui - vi + uv);
00226 }
00227 else if (method[0]==3) {
00228 ri = rhop/R0;
00229 rt0 = DataType(4.) / DataType(9.);
00230 rt1 = DataType(1.) / DataType(9.);
00231 rt2 = DataType(1.) / DataType(36.);
00232
00233 DataType usq = u * u;
00234 DataType vsq = v * v;
00235 DataType sumsq = (usq + vsq) / cs22;
00236 DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00237 DataType u2 = usq / cssq;
00238 DataType v2 = vsq / cssq;
00239 DataType ui = u / cs2;
00240 DataType vi = v / cs2;
00241 DataType uv = ui * vi;
00242
00243 feq(0) = rt0 * ( rho + ri*( - sumsq));
00244
00245 feq(1) = rt1 * ( rho + ri*( - sumsq + u2 + ui));
00246 feq(2) = rt1 * ( rho + ri*( - sumsq + u2 - ui));
00247 feq(3) = rt1 * ( rho + ri*( - sumsq + v2 + vi));
00248 feq(6) = rt1 * ( rho + ri*( - sumsq + v2 - vi));
00249
00250 feq(4) = rt2 * ( rho + ri*( + sumsq2 + ui + vi + uv));
00251 feq(5) = rt2 * ( rho + ri*( + sumsq2 - ui + vi - uv));
00252 feq(7) = rt2 * ( rho + ri*( + sumsq2 + ui - vi - uv));
00253 feq(8) = rt2 * ( rho + ri*( + sumsq2 - ui - vi + uv));
00254 }
00255
00256 return feq;
00257 }
00258
00259 inline virtual void Collision(MicroType &f, const DataType dt) const {
00260 MacroType q = MacroVariables(f);
00261 MicroType feq = Equilibrium(q);
00262 DataType omega = 1.;
00263 if (turbulence == laminar ) {
00264 omega = Omega(dt);
00265 }
00266 else
00267 omega = Omega_LES_Smagorinsky(f,feq,q(0),Omega(dt),dt);
00268
00269 for (register int qi=0; qi<9; qi++)
00270 f(qi)=(DataType(1.)-omega)*f(qi) + omega*feq(qi);
00271 }
00272
00273 virtual int IncomingIndices(const int side, int indices[]) const {
00274 switch (side) {
00275 case base::Left:
00276 indices[0] = 1; indices[1] = 4; indices[2] = 7;
00277 break;
00278 case base::Right:
00279 indices[0] = 2; indices[1] = 5; indices[2] = 8;
00280 break;
00281 case base::Bottom:
00282 indices[0] = 3; indices[1] = 4; indices[2] = 5;
00283 break;
00284 case base::Top:
00285 indices[0] = 6; indices[1] = 7; indices[2] = 8;
00286 break;
00287 }
00288 return 3;
00289 }
00290
00291 virtual int OutgoingIndices(const int side, int indices[]) const {
00292 switch (side) {
00293 case base::Left:
00294 indices[0] = 2; indices[1] = 5; indices[2] = 8;
00295 break;
00296 case base::Right:
00297 indices[0] = 1; indices[1] = 4; indices[2] = 7;
00298 break;
00299 case base::Bottom:
00300 indices[0] = 6; indices[1] = 7; indices[2] = 8;
00301 break;
00302 case base::Top:
00303 indices[0] = 3; indices[1] = 4; indices[2] = 5;
00304 break;
00305 }
00306 return 3;
00307 }
00308
00309
00310 virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb,
00311 const int side) const {
00312 int Nx = fvec.extents(0);
00313 int b = fvec.bottom();
00314 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00315 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00316 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00317 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00318 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00319 MicroType *f = (MicroType *)fvec.databuffer();
00320
00321 #ifdef DEBUG_PRINT_FIXUP
00322 ( comm_service::log() << "Reverse streaming at Side: " << side << " of "
00323 << fvec.bbox() << " using " << bb << "\n").flush();
00324 #endif
00325
00326 switch (side) {
00327 case base::Left:
00328 for (register int j=mys; j<=mye; j++) {
00329 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](2);
00330 f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](8);
00331 f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](5);
00332 }
00333 break;
00334 case base::Right:
00335 for (register int j=mys; j<=mye; j++) {
00336 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](1);
00337 f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](7);
00338 f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](4);
00339 }
00340 break;
00341 case base::Bottom:
00342 for (register int i=mxs; i<=mxe; i++) {
00343 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](6);
00344 f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](8);
00345 f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](7);
00346 }
00347 break;
00348 case base::Top:
00349 for (register int i=mxs; i<=mxe; i++) {
00350 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](3);
00351 f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](5);
00352 f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](4);
00353 }
00354 break;
00355 }
00356 }
00357
00358 virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec,
00359 const BBox &bb, const double& dt) const {
00360 int Nx = fvec.extents(0);
00361 int b = fvec.bottom();
00362 MicroType *f = (MicroType *)fvec.databuffer();
00363 MicroType *of = (MicroType *)ovec.databuffer();
00364
00365 #ifdef DEBUG_PRINT_FIXUP
00366 ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox()
00367 << " using " << bb << "\n").flush();
00368 #endif
00369 assert (fvec.extents(0)==ovec.extents(0));
00370 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) &&
00371 fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1));
00372 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00373 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00374 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00375 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00376
00377 if (turbulence != WALE && turbulence != CSM) {
00378 for (register int j=mys; j<=mye; j++)
00379 for (register int i=mxs; i<=mxe; i++) {
00380 for (register int n=0; n<base::NMicroVar(); n++)
00381 f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00382 Collision(f[b+base::idx(i,j,Nx)],dt);
00383 }
00384 }
00385 else if (turbulence == WALE) {
00386 DCoords dx = base::GH().worldStep(fvec.stepsize());
00387 MicroType fxp, fxm, fyp, fym;
00388 for (register int j=mys; j<=mye; j++)
00389 for (register int i=mxs; i<=mxe; i++) {
00390 for (register int n=0; n<base::NMicroVar(); n++) {
00391 f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00392 fxp(n) = of[b+base::idx(i-mdx[n]+1,j-mdy[n],Nx)](n);
00393 fxm(n) = of[b+base::idx(i-mdx[n]-1,j-mdy[n],Nx)](n);
00394 fyp(n) = of[b+base::idx(i-mdx[n],j-mdy[n]+1,Nx)](n);
00395 fym(n) = of[b+base::idx(i-mdx[n],j-mdy[n]-1,Nx)](n);
00396 }
00397 LocalCollisionWALE(f[b+base::idx(i,j,Nx)],
00398 MacroVariables(fxp),MacroVariables(fxm),
00399 MacroVariables(fyp),MacroVariables(fym),dx,dt);
00400 }
00401 }
00402 else if (turbulence == CSM) {
00403 DCoords dx = base::GH().worldStep(fvec.stepsize());
00404 MicroType fxp, fxm, fyp, fym;
00405 for (register int j=mys; j<=mye; j++)
00406 for (register int i=mxs; i<=mxe; i++) {
00407 for (register int n=0; n<base::NMicroVar(); n++) {
00408 f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00409 fxp(n) = of[b+base::idx(i-mdx[n]+1,j-mdy[n],Nx)](n);
00410 fxm(n) = of[b+base::idx(i-mdx[n]-1,j-mdy[n],Nx)](n);
00411 fyp(n) = of[b+base::idx(i-mdx[n],j-mdy[n]+1,Nx)](n);
00412 fym(n) = of[b+base::idx(i-mdx[n],j-mdy[n]-1,Nx)](n);
00413 }
00414 LocalCollisionCSM(f[b+base::idx(i,j,Nx)],
00415 MacroVariables(fxp),MacroVariables(fxm),
00416 MacroVariables(fyp),MacroVariables(fym),dx,dt);
00417 }
00418 }
00419
00420 }
00421
00422 inline virtual macro_grid_data_type Filter(vec_grid_data_type &fvec) const {
00423 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00424 int mxs = base::NGhosts(), mys = base::NGhosts();
00425 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00426 MicroType *fv = (MicroType *)fvec.databuffer();
00427 macro_grid_data_type qgrid(fvec.bbox());
00428 MacroType *qg = (MacroType *)qgrid.databuffer();
00429 for (register int j=mys-2; j<=mye+2; j++)
00430 for (register int i=mxs-2; i<=mxe+2; i++) {
00431 qg[base::idx(i,j,Nx)] = DataType(0.25)*MacroVariables(fv[base::idx(i-1,j,Nx)])
00432 + DataType(0.5)*MacroVariables(fv[base::idx(i,j,Nx)])
00433 + DataType(0.25)*MacroVariables(fv[base::idx(i+1,j,Nx)]);
00434 }
00435 macro_grid_data_type qgrid1(fvec.bbox());
00436 MacroType *sv = (MacroType *)qgrid1.databuffer();
00437 for (register int j=mys-1; j<=mye+1; j++)
00438 for (register int i=mxs-1; i<=mxe+1; i++) {
00439 sv[base::idx(i,j,Nx)] = DataType(0.25)*qg[base::idx(i,j+1,Nx)]
00440 + DataType(0.5)*qg[base::idx(i,j,Nx)]
00441 + DataType(0.25)*qg[base::idx(i,j-1,Nx)];
00442 }
00443 return qgrid1;
00444 }
00445
00446 virtual void CollisionDynamicSmagorinskyLES(vec_grid_data_type &fvec, const double& dt) const {
00456 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00457 int mxs = base::NGhosts(), mys = base::NGhosts();
00458 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00459 MicroType *fv = (MicroType *)fvec.databuffer();
00460 MacroType q;
00461 MicroType feq;
00462 DataType csdyn;
00463 DataType dx2=DataType(1.), dxf2 = DataType(4.), dxf = DataType(2.);
00464 DataType M2, SF2;
00465 DataType sfxx, sfxy, sfyy;
00466 DataType lxx, lxy, lyy;
00467 DataType mxx, mxy, myy;
00468 DataType omega = Omega(dt);
00469 DataType om = omega;
00470 DataType nu_t = DataType(0.);
00471
00472 TensorType Sigma, S_;
00473 DataType S;
00474 DataType clim = DataType(1.0e7);
00475
00476 macro_grid_data_type qgrid1 = Filter(fvec);
00477 MacroType *sv = (MacroType *)qgrid1.databuffer();
00478
00479 for (register int j=mys; j<=mye; j++)
00480 for (register int i=mxs; i<=mxe; i++) {
00481 q = MacroVariables(fv[base::idx(i,j,Nx)]);
00482 feq = Equilibrium(q);
00483 Sigma = DeviatoricStress(fv[base::idx(i,j,Nx)],feq,om)/(DataType(1.0)-DataType(0.5)*om);
00484 S_ = StrainComponents(q(0),Sigma,om,Cs_Smagorinsky);
00485
00486 S = Magnitude(S_);
00487
00488 sfxx = ( sv[base::idx(i+1,j,Nx)](1) - sv[base::idx(i-1,j,Nx)](1) )/dxf;
00489 sfxy = DataType(0.5)*( sv[base::idx(i+1,j,Nx)](2) - sv[base::idx(i-1,j,Nx)](2)
00490 + sv[base::idx(i,j+1,Nx)](1) - sv[base::idx(i,j-1,Nx)](1) )/dxf;
00491 sfyy = ( sv[base::idx(i,j+1,Nx)](2) - sv[base::idx(i,j-1,Nx)](2) )/dxf;
00492 SF2 = std::sqrt( DataType(2.)*(sfxx*sfxx + DataType(2.)*sfxy*sfxy + sfyy*sfyy) );
00493
00494 mxx = (dxf2*SF2*sfxx - dx2*S*S_(0));
00495 mxy = (dxf2*SF2*sfxy - dx2*S*S_(1));
00496 myy = (dxf2*SF2*sfyy - dx2*S*S_(2));
00497 M2 = (mxx*mxx + DataType(2.)*mxy*mxy + myy*myy);
00498
00499 lxx = (q(1)*q(1) - sv[base::idx(i,j,Nx)](1)* sv[base::idx(i,j,Nx)](1) );
00500 lxy = (q(1)*q(2) - sv[base::idx(i,j,Nx)](1)* sv[base::idx(i,j,Nx)](2) );
00501 lyy = (q(2)*q(2) - sv[base::idx(i,j,Nx)](2)* sv[base::idx(i,j,Nx)](2) );
00502
00503 DataType ntmp = (lxx*mxx + DataType(2.)*lxy*mxy + lyy*myy);
00504 if (std::isnan(ntmp)==1) { ntmp = 0.0; }
00505 if (std::isnan(M2)==1) { csdyn = Cs_Smagorinsky*Cs_Smagorinsky; }
00506 else if (M2 < mfp/U0) { csdyn = DataType(-0.5)*ntmp/(mfp/U0); }
00507 else if (M2 > clim) { csdyn = DataType(-0.5)*ntmp/(clim); }
00508 else { csdyn = DataType(-0.5)*ntmp/(M2); }
00509 if (csdyn < DataType(-20.0)*Cs_Smagorinsky) { csdyn = DataType(-20.0)*Cs_Smagorinsky; }
00510 else if (csdyn > DataType(20.0)*Cs_Smagorinsky) { csdyn = DataType(20.0)*Cs_Smagorinsky; }
00511 nu_t = std::fabs(csdyn)*S;
00512 omega = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
00513 if (omega < DataType(1.000001)) { omega = DataType(1.00001); }
00514 else if (omega > DataType(1.999999)) { omega = DataType(1.999999); }
00515 for (register int qi=0; qi<9; qi++)
00516 fv[base::idx(i,j,Nx)](qi) = (DataType(1.)-omega)*fv[base::idx(i,j,Nx)](qi) + omega*feq(qi);
00517 }
00518 }
00519
00520 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00521 const double& t, const double& dt, const int& mpass) const {
00522
00523
00524 DCoords dx = base::GH().worldStep(fvec.stepsize());
00525
00526 if (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR) {
00527 std::cerr << "LBMD2Q9 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00528 << " to be equal." << std::endl << std::flush;
00529 std::exit(-1);
00530 }
00531 if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR*LB_FACTOR*LB_FACTOR) {
00532 std::cerr << "LBMD2Q9 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00533 << " to be equal."
00534 << " dt=" << dt << " TO=" << T0()
00535 << std::endl << std::flush;
00536
00537 }
00538
00539 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00540 MicroType *f = (MicroType *)fvec.databuffer();
00541
00542
00543 for (register int j=Ny-1; j>=1; j--)
00544 for (register int i=0; i<=Nx-2; i++) {
00545 f[base::idx(i,j,Nx)](3) = f[base::idx(i,j-1,Nx)](3);
00546 f[base::idx(i,j,Nx)](5) = f[base::idx(i+1,j-1,Nx)](5);
00547 }
00548
00549 for (register int j=Ny-1; j>=1; j--)
00550 for (register int i=Nx-1; i>=1; i--) {
00551 f[base::idx(i,j,Nx)](1) = f[base::idx(i-1,j,Nx)](1);
00552 f[base::idx(i,j,Nx)](4) = f[base::idx(i-1,j-1,Nx)](4);
00553 }
00554
00555 for (register int j=0; j<=Ny-2; j++)
00556 for (register int i=Nx-1; i>=1; i--) {
00557 f[base::idx(i,j,Nx)](6) = f[base::idx(i,j+1,Nx)](6);
00558 f[base::idx(i,j,Nx)](7) = f[base::idx(i-1,j+1,Nx)](7);
00559 }
00560
00561 for (register int j=0; j<=Ny-2; j++)
00562 for (register int i=0; i<=Nx-2; i++) {
00563 f[base::idx(i,j,Nx)](2) = f[base::idx(i+1,j,Nx)](2);
00564 f[base::idx(i,j,Nx)](8) = f[base::idx(i+1,j+1,Nx)](8);
00565 }
00566
00567
00568 int mxs = base::NGhosts(), mys = base::NGhosts();
00569 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00570
00571 if (turbulence==laminar || turbulence==LES_Smagorinsky || turbulence==LES_SmagorinskyMemory) {
00572 for (register int j=mys; j<=mye; j++)
00573 for (register int i=mxs; i<=mxe; i++)
00574 Collision(f[base::idx(i,j,Nx)],dt);
00575 }
00576 else if ( turbulence == LES_dynamic || turbulence == LES_dynamicMemory ){
00577 CollisionDynamicSmagorinskyLES(fvec,dt);
00578 }
00579 else if ( turbulence == WALE ){
00580 CollisionWALE(fvec,dt);
00581 }
00582 else if ( turbulence == CSM ){
00583 CollisionCSM(fvec,dt);
00584 }
00585
00586 return DataType(1.);
00587 }
00588
00589 virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00590 DataType* aux=0, const int naux=0, const int scaling=0) const {
00591 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00592 int mxs = base::NGhosts(), mys = base::NGhosts();
00593 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00594 MicroType *f = (MicroType *)fvec.databuffer();
00595
00596 if (type==ConstantMicro) {
00597 MicroType g;
00598 for (register int i=0; i<base::NMicroVar(); i++)
00599 g(i) = aux[i];
00600
00601 for (register int j=mys; j<=mye; j++)
00602 for (register int i=mxs; i<=mxe; i++)
00603 f[base::idx(i,j,Nx)] = g;
00604 }
00605
00606 else if (type==ConstantMacro) {
00607 MacroType q(0.);
00608 if (scaling==base::Physical) {
00609 q(0) = aux[0]/R0;
00610 q(1) = aux[1]/U0;
00611 q(2) = aux[2]/U0;
00612 }
00613 else
00614 for (register int i=0; i<base::NMacroVar(); i++)
00615 q(i) = aux[i];
00616 q(1) *= S0; q(2) *= S0;
00617
00618 for (register int j=mys; j<=mye; j++)
00619 for (register int i=mxs; i<=mxe; i++)
00620 f[base::idx(i,j,Nx)] = Equilibrium(q);
00621 }
00622
00623 else if (type==GasAtRest) {
00624 MacroType q;
00625 q(0) = GasDensity()/R0;
00626 q(1) = DataType(0.);
00627 q(2) = DataType(0.);
00628 for (register int j=mys; j<=mye; j++)
00629 for (register int i=mxs; i<=mxe; i++)
00630 f[base::idx(i,j,Nx)] = Equilibrium(q);
00631 }
00632 }
00633
00634
00635 virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00636 DataType* aux=0, const int naux=0, const int scaling=0) const {
00637 int Nx = fvec.extents(0);
00638 int b = fvec.bottom();
00639 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00640 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00641 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00642 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00643 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00644 MicroType *f = (MicroType *)fvec.databuffer();
00645
00646 if (type==Symmetry) {
00647 switch (side) {
00648 case base::Left:
00649 for (register int j=mys; j<=mye; j++) {
00650 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j,Nx)](8);
00651 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00652 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j,Nx)](5);
00653 }
00654 break;
00655 case base::Right:
00656 for (register int j=mys; j<=mye; j++) {
00657 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j,Nx)](7);
00658 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00659 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j,Nx)](4);
00660 }
00661 break;
00662 case base::Bottom:
00663 for (register int i=mxs; i<=mxe; i++) {
00664 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i,mye+1,Nx)](8);
00665 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00666 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i,mye+1,Nx)](7);
00667 }
00668 break;
00669 case base::Top:
00670 for (register int i=mxs; i<=mxe; i++) {
00671 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i,mys-1,Nx)](5);
00672 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00673 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i,mys-1,Nx)](4);
00674 }
00675 break;
00676 }
00677 }
00678
00679 else if (type==SlipWall) {
00680 switch (side) {
00681 case base::Left:
00682 for (register int j=mys; j<=mye; j++) {
00683 if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00684 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00685 if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00686 }
00687 break;
00688 case base::Right:
00689 for (register int j=mys; j<=mye; j++) {
00690 if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00691 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00692 if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00693 }
00694 break;
00695 case base::Bottom:
00696 for (register int i=mxs; i<=mxe; i++) {
00697 if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00698 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00699 if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00700 }
00701 break;
00702 case base::Top:
00703 for (register int i=mxs; i<=mxe; i++) {
00704 if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00705 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00706 if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00707 }
00708 break;
00709 }
00710 }
00711
00712 else if (type==NoSlipWall) {
00713 switch (side) {
00714 case base::Left:
00715 for (register int j=mys; j<=mye; j++) {
00716 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00717 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00718 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00719 }
00720 break;
00721 case base::Right:
00722 for (register int j=mys; j<=mye; j++) {
00723 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00724 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00725 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00726 }
00727 break;
00728 case base::Bottom:
00729 for (register int i=mxs; i<=mxe; i++) {
00730 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00731 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00732 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00733 }
00734 break;
00735 case base::Top:
00736 for (register int i=mxs; i<=mxe; i++) {
00737 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00738 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00739 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00740 }
00741 break;
00742 }
00743 }
00744
00745 else if (type==Inlet) {
00746 if (aux==0 || naux<=0) return;
00747 DataType a0=S0*aux[0], a1=0.;
00748 if (naux>1) a1 = S0*aux[1];
00749 if (scaling==base::Physical) {
00750 a0 /= U0;
00751 a1 /= U0;
00752 }
00753 switch (side) {
00754 case base::Left:
00755 for (register int j=mys; j<=mye; j++) {
00756 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00757 q(1) = a0;
00758 if (naux>1) q(2) = a1;
00759 f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00760 }
00761 break;
00762 case base::Right:
00763 for (register int j=mys; j<=mye; j++) {
00764 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00765 q(1) = a0;
00766 if (naux>1) q(2) = a1;
00767 f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00768 }
00769 break;
00770 case base::Bottom:
00771 for (register int i=mxs; i<=mxe; i++) {
00772 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00773 if (naux==1)
00774 q(2) = a0;
00775 else {
00776 q(1) = a0;
00777 q(2) = a1;
00778 }
00779 f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00780 }
00781 break;
00782 case base::Top:
00783 for (register int i=mxs; i<=mxe; i++) {
00784 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00785 if (naux==1)
00786 q(2) = a0;
00787 else {
00788 q(1) = a0;
00789 q(2) = a1;
00790 }
00791 f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00792 }
00793 break;
00794 }
00795 }
00796
00797 else if (type==Outlet) {
00798 switch (side) {
00799 case base::Left:
00800 for (register int j=mys; j<=mye; j++) {
00801 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00802 if (q(0)>0) f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)]/q(0);
00803 }
00804 break;
00805 case base::Right:
00806 for (register int j=mys; j<=mye; j++) {
00807 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00808 if (q(0)>0) f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)]/q(0);
00809 }
00810 break;
00811 case base::Bottom:
00812 for (register int i=mxs; i<=mxe; i++) {
00813 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00814 if (q(0)>0) f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)]/q(0);
00815 }
00816 break;
00817 case base::Top:
00818 for (register int i=mxs; i<=mxe; i++) {
00819 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00820 if (q(0)>0) f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)]/q(0);
00821 }
00822 break;
00823 }
00824 }
00825
00826
00827 else if (type==Pressure) {
00828 if (aux==0 || naux<=0) return;
00829 DataType a0=aux[0];
00830 if (scaling==base::Physical)
00831 a0 /= R0;
00832 switch (side) {
00833 case base::Left:
00834 for (register int j=mys; j<=mye; j++) {
00835 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00836 q(0) = a0;
00837 f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00838 }
00839 break;
00840 case base::Right:
00841 for (register int j=mys; j<=mye; j++) {
00842 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00843 q(0) = a0;
00844 f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00845 }
00846 break;
00847 case base::Bottom:
00848 for (register int i=mxs; i<=mxe; i++) {
00849 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00850 q(0) = a0;
00851 f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00852 }
00853 break;
00854 case base::Top:
00855 for (register int i=mxs; i<=mxe; i++) {
00856 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00857 q(0) = a0;
00858 f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00859 }
00860 break;
00861 }
00862 }
00863
00864 else if (type==SlidingWall) {
00865 if (aux==0 || naux<=0) return;
00866 DataType a0=S0*aux[0];
00867 if (scaling==base::Physical)
00868 a0 /= U0;
00869 switch (side) {
00870 case base::Left:
00871 for (register int j=mys; j<=mye; j++) {
00872 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00873 DataType rd1 = f[b+base::idx(mxe+1,j,Nx)](5), rd2 = f[b+base::idx(mxe+1,j,Nx)](8);
00874 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00875 DataType pw = DataType(1.)-qw;
00876 if (j<mye) f[b+base::idx(mxe,j+1,Nx)](7) = pw*rd1+qw*rd2;
00877 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00878 if (j>mys) f[b+base::idx(mxe,j-1,Nx)](4) = qw*rd1+pw*rd2;
00879 }
00880 break;
00881 case base::Right:
00882 for (register int j=mys; j<=mye; j++) {
00883 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00884 DataType rd1 = f[b+base::idx(mxs-1,j,Nx)](4), rd2 = f[b+base::idx(mxs-1,j,Nx)](7);
00885 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00886 DataType pw = DataType(1.)-qw;
00887 if (j<mye) f[b+base::idx(mxs,j+1,Nx)](8) = pw*rd1+qw*rd2;
00888 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00889 if (j>mys) f[b+base::idx(mxs,j-1,Nx)](5) = qw*rd1+pw*rd2;
00890 }
00891 break;
00892 case base::Bottom:
00893 for (register int i=mxs; i<=mxe; i++) {
00894 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00895 DataType rd1 = f[b+base::idx(i,mye+1,Nx)](7), rd2 = f[b+base::idx(i,mye+1,Nx)](8);
00896 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00897 DataType pw = DataType(1.)-qw;
00898 if (i<mxe) f[b+base::idx(i+1,mye,Nx)](5) = pw*rd1+qw*rd2;
00899 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00900 if (i>mxs) f[b+base::idx(i-1,mye,Nx)](4) = qw*rd1+pw*rd2;
00901 }
00902 break;
00903 case base::Top:
00904 for (register int i=mxs; i<=mxe; i++) {
00905 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00906 DataType rd1 = f[b+base::idx(i,mys-1,Nx)](4), rd2 = f[b+base::idx(i,mys-1,Nx)](5);
00907 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00908 DataType pw = DataType(1.)-qw;
00909 if (i<mxe) f[b+base::idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00910 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00911 if (i>mxs) f[b+base::idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00912 }
00913 break;
00914 }
00915 }
00916
00920 else if (type==CharacteristicOutlet) {
00921 DataType rhod=DataType(1.0);
00922 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
00923 DCoords dx = base::GH().worldStep(fvec.stepsize());
00924 DataType dt_temp = dx(0)/VelocityScale();
00925 switch (side) {
00926 case base::Left:
00927 for (register int j=mys; j<=mye; j++) {
00928 MicroType ftmp = f[b+base::idx(mxe,j,Nx)];
00929 MacroType q = MacroVariables(ftmp);
00930 if (turbulence!=CSM)
00931 Collision(ftmp,dt_temp);
00932 else {
00933 MacroType qxp=0., qxm=0., qyp=0., qym=0.;
00934 qxp = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00935 qxm = MacroVariables(f[b+base::idx(mxe-1,j,Nx)]);
00936 if (j<mye) qyp = MacroVariables(f[b+base::idx(mxe,j+1,Nx)]);
00937 else qyp=q;
00938 if (j>mys) qym = MacroVariables(f[b+base::idx(mxe,j-1,Nx)]);
00939 else qym=q;
00940 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,dx,dt_temp);
00941 }
00942 MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00943 MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,Nx)]);
00944
00945 MacroType dqdn = NormalDerivative(q,q1,q2);
00946 if (EquilibriumType()!=3) rhod = q(0);
00947 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
00948
00949
00950 MacroType qn;
00951 qn(0) = q(0) - c1oCsSq2*L(0);
00952 qn(1) = q(1) + c1o2cs*L(0)/rhod;
00953 qn(2) = q(2) - L(1);
00954 f[b+base::idx(mxe,j,Nx)] = Equilibrium(qn);
00955 }
00956 break;
00957 case base::Right:
00958 for (register int j=mys; j<=mye; j++) {
00959 MicroType ftmp = f[b+base::idx(mxs,j,Nx)];
00960 MacroType q = MacroVariables(ftmp);
00961 if (turbulence!=CSM)
00962 Collision(ftmp,dt_temp);
00963 else {
00964 MacroType qxp=0., qxm=0., qyp=0., qym=0.;
00965 qxp = MacroVariables(f[b+base::idx(mxs+1,j,Nx)]);
00966 qxm = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00967 if (j<mye) qyp = MacroVariables(f[b+base::idx(mxs,j+1,Nx)]);
00968 else qyp=q;
00969 if (j>mys) qym = MacroVariables(f[b+base::idx(mxs,j-1,Nx)]);
00970 else qym=q;
00971 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,dx,dt_temp);
00972 }
00973 MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00974 MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,Nx)]);
00975
00976 MacroType dqdn = -NormalDerivative(q,q1,q2);
00977 if (EquilibriumType()!=3) rhod = q(0);
00978 Vector<DataType,2> L = WaveAmplitudes(rhod,-q(1),dqdn(0),dqdn(1),dqdn(2));
00979
00980
00981 MacroType qn;
00982 qn(0) = q(0) - c1oCsSq2*L(0);
00983 qn(1) = q(1) - c1o2cs*L(0)/rhod;
00984 qn(2) = q(2) + L(1);
00985 f[b+base::idx(mxs,j,Nx)] = Equilibrium(qn);
00986 }
00987 break;
00988 case base::Bottom:
00989
00990 for (register int i=mxs; i<=mxe; i++) {
00991 MicroType ftmp = f[b+base::idx(i,mye,Nx)];
00992 MacroType q = MacroVariables(ftmp);
00993 if (turbulence!=CSM)
00994 Collision(ftmp,dt_temp);
00995 else {
00996 MacroType qxp=0., qxm=0., qyp=0., qym=0.;
00997 if (i<mxe) qxp = MacroVariables(f[b+base::idx(i+1,mye,Nx)]);
00998 else qxp=q;
00999 if (i>mxs) qxm = MacroVariables(f[b+base::idx(i-1,mye,Nx)]);
01000 else qxm=q;
01001 qyp = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
01002 qym = MacroVariables(f[b+base::idx(i,mye-1,Nx)]);
01003 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,dx,dt_temp);
01004 }
01005 MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
01006 MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,Nx)]);
01007
01008 MacroType dqdn = NormalDerivative(q,q1,q2);
01009 if (EquilibriumType()!=3) rhod = q(0);
01010 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01011
01012
01013 MacroType qn;
01014 qn(0) = q(0) - c1oCsSq2*L(0);
01015 qn(1) = q(1) - L(1);
01016 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01017 f[b+base::idx(i,mye,Nx)] = Equilibrium(qn);
01018 }
01019 break;
01020 case base::Top:
01021 for (register int i=mxs+base::NGhosts(); i<=mxe-base::NGhosts(); i++) {
01022 MicroType ftmp = f[b+base::idx(i,mys,Nx)];
01023 MacroType q = MacroVariables(ftmp);
01024 if (turbulence!=CSM)
01025 Collision(ftmp,dt_temp);
01026 else {
01027 MacroType qxp=0., qxm=0., qyp=0., qym=0.;
01028 if (i<mxe) qxp = MacroVariables(f[b+base::idx(i+1,mys,Nx)]);
01029 else qxp=q;
01030 if (i>mxs) qxm = MacroVariables(f[b+base::idx(i-1,mys,Nx)]);
01031 else qxm=q;
01032 qyp = MacroVariables(f[b+base::idx(i,mys+1,Nx)]);
01033 qym = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
01034 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,dx,dt_temp);
01035 }
01036 MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
01037 MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,Nx)]);
01038
01039 MacroType dqdn = -NormalDerivative(q,q1,q2);
01040 if (EquilibriumType()!=3) rhod = q(0);
01041 Vector<DataType,2> L = WaveAmplitudes(rhod,-q(2),dqdn(0),dqdn(2),dqdn(1));
01042
01043
01044 MacroType qn;
01045 qn(0) = q(0) - c1oCsSq2*L(0);
01046 qn(1) = q(1) - L(1);
01047 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01048 f[b+base::idx(i,mys,Nx)] = Equilibrium(qn);
01049 }
01050 break;
01051 }
01052 }
01053
01057 else if (type==CharacteristicInlet) {
01058 if (aux==0 || naux<=0) return;
01059 DataType rho=aux[0];
01060 DataType a0=S0*aux[1];
01061 DataType a1=S0*aux[2];
01062 if (scaling==base::Physical) {
01063 rho /= R0;
01064 a0 /= U0;
01065 a1 /= U0;
01066 }
01067 MacroType q, q1, q2, qn;
01068 DataType rhod=DataType(1.0);
01069 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
01070 q(0) = rho;
01071 q(1) = a0;
01072 q(2) = a1;
01073 switch (side) {
01074 case base::Left:
01075 for (register int j=mys; j<=mye; j++) {
01076 q1 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
01077 q2 = MacroVariables(f[b+base::idx(mxe+2,j,Nx)]);
01078 if (naux==1) {
01079 q(1) = q1(1); q(2) = q1(2);
01080 MacroType dqdn = NormalDerivative(q,q1,q2);
01081 if (EquilibriumType()!=3)
01082 rhod = q(0);
01083 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01084 qn(0) = q(0) - c1oCsSq2*L(0);
01085 qn(1) = q1(1) + c1o2cs*L(0)/rhod;
01086 qn(2) = q1(2) - L(1);
01087 }
01088 else if (naux==2) {
01089 q(0) = q1(0); q(2) = q1(2);
01090 MacroType dqdn = NormalDerivative(q,q1,q2);
01091 if (EquilibriumType()!=3)
01092 rhod = q(0);
01093 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
01094 qn(0) = q1(0) - c1oCsSq2*L(0);
01095 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01096 qn(2) = q1(2) - L(1);
01097 }
01098 else if (naux==3) {
01099 q(0) = q1(0);
01100 MacroType dqdn = NormalDerivative(q,q1,q2);
01101 if (EquilibriumType()!=3)
01102 rhod = q(0);
01103 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
01104 qn(0) = q1(0) - c1oCsSq2*L(0);
01105 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01106 qn(2) = q(2) - L(1);
01107 }
01108 f[b+base::idx(mxe,j,Nx)] = Equilibrium(qn);
01109 }
01110 break;
01111 case base::Right:
01112 for (register int j=mys; j<=mye; j++) {
01113 q1 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
01114 q2 = MacroVariables(f[b+base::idx(mxs-2,j,Nx)]);
01115 if (naux==1) {
01116 q(1) = q1(1); q(2) = q1(2);
01117 MacroType dqdn = NormalDerivative(q,q1,q2);
01118 if (EquilibriumType()!=3)
01119 rhod = q(0);
01120 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
01121 qn(0) = q(0) - c1oCsSq2*L(0);
01122 qn(1) = q1(1) + c1o2cs*L(0)/rhod;
01123 qn(2) = q1(2) - L(1);
01124 }
01125 else if (naux==2) {
01126 q(0) = q1(0); q(2) = q1(2);
01127 MacroType dqdn = NormalDerivative(q,q1,q2);
01128 if (EquilibriumType()!=3)
01129 rhod = q(0);
01130 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
01131 qn(0) = q1(0) - c1oCsSq2*L(0);
01132 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01133 qn(2) = q1(2) - L(1);
01134 }
01135 else if (naux==3) {
01136 q(0) = q1(0);
01137 MacroType dqdn = NormalDerivative(q,q1,q2);
01138 if (EquilibriumType()!=3)
01139 rhod = q(0);
01140 Vector<DataType,2> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2));
01141 qn(0) = q1(0) - c1oCsSq2*L(0);
01142 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01143 qn(2) = q(2) - L(1);
01144 }
01145 f[b+base::idx(mxs,j,Nx)] = Equilibrium(qn);
01146 }
01147 break;
01148 case base::Bottom:
01149 for (register int i=mxs; i<=mxe; i++) {
01150 q1 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
01151 q2 = MacroVariables(f[b+base::idx(i,mye+2,Nx)]);
01152 if (naux==1) {
01153 q(1) = q1(1); q(2) = q1(2);
01154 MacroType dqdn = NormalDerivative(q,q1,q2);
01155 if (EquilibriumType()!=3)
01156 rhod = q(0);
01157 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01158 qn(0) = q(0) - c1oCsSq2*L(0);
01159 qn(1) = q1(1) - L(1);
01160 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
01161 }
01162 else if (naux==2) {
01163 q(0) = q1(0); q(1) = q1(1);
01164 MacroType dqdn = NormalDerivative(q,q1,q2);
01165 if (EquilibriumType()!=3)
01166 rhod = q(0);
01167 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01168 qn(0) = q1(0) - c1oCsSq2*L(0);
01169 qn(1) = q1(1) - L(1);
01170 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01171 }
01172 else if (naux==3) {
01173 q(0) = q1(0);
01174 MacroType dqdn = NormalDerivative(q,q1,q2);
01175 if (EquilibriumType()!=3)
01176 rhod = q(0);
01177 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01178 qn(0) = q1(0) - c1oCsSq2*L(0);
01179 qn(1) = q(1) - L(1);
01180 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01181 }
01182 f[b+base::idx(i,mye,Nx)] = Equilibrium(qn);
01183 }
01184 break;
01185 case base::Top:
01186 for (register int i=mxs; i<=mxe; i++) {
01187 q1 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
01188 q2 = MacroVariables(f[b+base::idx(i,mys-2,Nx)]);
01189 if (naux==1) {
01190 q(1) = q1(1); q(2) = q1(2);
01191 MacroType dqdn = NormalDerivative(q,q1,q2);
01192 if (EquilibriumType()!=3)
01193 rhod = q(0);
01194 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01195 qn(0) = q(0) - c1oCsSq2*L(0);
01196 qn(1) = q1(1) - L(1);
01197 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
01198 }
01199 else if (naux==2) {
01200 q(0) = q1(0); q(1) = q1(1);
01201 MacroType dqdn = NormalDerivative(q,q1,q2);
01202 if (EquilibriumType()!=3)
01203 rhod = q(0);
01204 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01205 qn(0) = q1(0) - c1oCsSq2*L(0);
01206 qn(1) = q1(1) - L(1);
01207 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01208 }
01209 else if (naux==3) {
01210 q(0) = q1(0);
01211 MacroType dqdn = NormalDerivative(q,q1,q2);
01212 if (EquilibriumType()!=3)
01213 rhod = q(0);
01214 Vector<DataType,2> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1));
01215 qn(0) = q1(0) - c1oCsSq2*L(0);
01216 qn(1) = q(1) - L(1);
01217 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01218 }
01219 f[b+base::idx(i,mys,Nx)] = Equilibrium(qn);
01220 }
01221 break;
01222 }
01223 }
01224
01229 else if (type==NoSlipWallEntranceExit) {
01230 DCoords wc;
01231 int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
01232 DataType WLB[2], WUB[2];
01233 base::GH().wlb(WLB);
01234 base::GH().wub(WUB);
01235 DataType slip = 0.;
01236 DataType lxs, lxe, lys, lye;
01237 DataType min[2], max[2];
01238 if (naux!=4) assert(0);
01239 lxs = aux[0] - WLB[0]; lxe = WUB[0] - aux[1];
01240 lys = aux[2] - WLB[1]; lye = WUB[1] - aux[3];
01241 min[0] = aux[0]; max[0] = aux[1];
01242 min[1] = aux[2]; max[1] = aux[3];
01243
01244 switch (side) {
01245 case base::Left:
01246 for (register int j=mys; j<=mye; j++) {
01247 lc_indx[1] = j*fvec.stepsize(1);
01248 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01249 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01250 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
01251 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
01252 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
01253 }
01254 else {
01255 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01256 slip = DataType(0.);
01257 else if (wc(1)<min[1])
01258 slip = (wc(1) - WLB[1])/lys;
01259 else if (wc(1)>max[1])
01260 slip = (WUB[1] - wc(1))/lye;
01261 if (j>mys && j<mye) {
01262 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
01263 + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
01264 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
01265 + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
01266 }
01267 else if (j==mys) {
01268 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](5)
01269 + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
01270 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
01271 + slip*f[b+base::idx(mxe+1,j,Nx)](5);
01272 }
01273 else if (j==mye) {
01274 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
01275 + slip*f[b+base::idx(mxe+1,j,Nx)](8);
01276 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](8)
01277 + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
01278 }
01279 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
01280 }
01281 }
01282 break;
01283 case base::Right:
01284 for (register int j=mys; j<=mye; j++) {
01285 lc_indx[1] = j*fvec.stepsize(1);
01286 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01287 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01288 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
01289 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
01290 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
01291 }
01292 else {
01293 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01294 slip = DataType(0.);
01295 else if (wc(1)<min[1])
01296 slip = (wc(1) - WLB[1])/lys;
01297 else if (wc(1)>max[1])
01298 slip = (WUB[1] - wc(1))/lye;
01299 if (j>mys && j<mye) {
01300 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
01301 + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
01302 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
01303 + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
01304 }
01305 else if (j==mys) {
01306 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](4)
01307 + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
01308 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
01309 + slip*f[b+base::idx(mxs-1,j,Nx)](4);
01310 }
01311 else if (j==mye) {
01312 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
01313 + slip*f[b+base::idx(mxs-1,j,Nx)](7);
01314 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](7)
01315 + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
01316 }
01317 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
01318 }
01319 }
01320 break;
01321 case base::Bottom:
01322 for (register int i=mxs; i<=mxe; i++) {
01323 lc_indx[0] = i*fvec.stepsize(0);
01324 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01325 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01326 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
01327 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
01328 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
01329 }
01330 else {
01331 if (wc(0)<WLB[0] || wc(0)>WUB[0])
01332 slip = DataType(0.);
01333 else if (wc(0)<min[0])
01334 slip = (wc(0) - WLB[0])/lxs;
01335 else if (wc(0)>max[0])
01336 slip = (WUB[0] - wc(0))/lxe;
01337 if (i>mxs && i<mxe) {
01338 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,Nx)](7)
01339 + slip*f[b+base::idx(i+1,mye+1,Nx)](8);
01340 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,Nx)](8)
01341 + slip*f[b+base::idx(i-1,mye+1,Nx)](7);
01342 }
01343 else if (i==mxs) {
01344 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](7)
01345 + slip*f[b+base::idx(i+1,mye+1,Nx)](8);
01346 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,Nx)](8)
01347 + slip*f[b+base::idx(i,mye+1,Nx)](7);
01348 }
01349 else if (i==mxe) {
01350 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](7)
01351 + slip*f[b+base::idx(i,mye+1,Nx)](8);
01352 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](8)
01353 + slip*f[b+base::idx(i-1,mye+1,Nx)](7);
01354 }
01355 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
01356 }
01357 }
01358 break;
01359 case base::Top:
01360 for (register int i=mxs; i<=mxe; i++) {
01361 lc_indx[0] = i*fvec.stepsize(0);
01362 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01363 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01364 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
01365 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
01366 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
01367 }
01368 else {
01369 if (wc(0)<WLB[0] || wc(0)>WUB[0])
01370 slip = DataType(0.);
01371 else if (wc(0)<min[0])
01372 slip = (wc(0) - WLB[0])/lxs;
01373 else if (wc(0)>max[0])
01374 slip = (WUB[0] - wc(0))/lxe;
01375 if (i>mxs && i<mxe) {
01376 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,Nx)](4)
01377 + slip*f[b+base::idx(i+1,mys-1,Nx)](5);
01378 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,Nx)](5)
01379 + slip*f[b+base::idx(i-1,mys-1,Nx)](4);
01380 }
01381 else if (i==mxs) {
01382 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,Nx)](4)
01383 + slip*f[b+base::idx(i+1,mys-1,Nx)](5);
01384 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,Nx)](5)
01385 + slip*f[b+base::idx(i,mys-1,Nx)](4);
01386 }
01387 else if (i==mxe) {
01388 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,Nx)](4)
01389 + slip*f[b+base::idx(i,mys-1,Nx)](5);
01390 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,Nx)](5)
01391 + slip*f[b+base::idx(i-1,mys-1,Nx)](4);
01392 }
01393 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
01394 }
01395 }
01396 break;
01397 }
01398 }
01399
01400 else if (type==NoSlipWallNeq) {
01401 switch (side) {
01402 case base::Left:
01403 for (register int j=mys; j<=mye; j++) {
01404 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
01405 MacroType qt = q; qt(1) = 0.; qt(2) = 0.;
01406 MicroType feq = Equilibrium(q);
01407 MicroType feqt = Equilibrium(qt);
01408 f[b+base::idx(mxe,j,Nx)] += feqt-feq;
01409
01410 f[b+base::idx(mxe,j,Nx)](7) = DataType(2.)*feqt(7)-feq(7);
01411 f[b+base::idx(mxe,j,Nx)](1) = DataType(2.)*feqt(1)-feq(1);
01412
01413 f[b+base::idx(mxe,j,Nx)](4) = DataType(2.)*feqt(4)-feq(4);
01414 }
01415 break;
01416 case base::Right:
01417 for (register int j=mys; j<=mye; j++) {
01418 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
01419 MacroType qt = q; qt(1) = 0.; qt(2) = 0.;
01420 MicroType feq = Equilibrium(q);
01421 MicroType feqt = Equilibrium(qt);
01422 f[b+base::idx(mxs,j,Nx)] += feqt-feq;
01423
01424 f[b+base::idx(mxs,j,Nx)](8) = DataType(2.)*feqt(8)-feq(8);
01425 f[b+base::idx(mxs,j,Nx)](2) = DataType(2.)*feqt(2)-feq(2);
01426
01427 f[b+base::idx(mxs,j,Nx)](5) = DataType(2.)*feqt(5)-feq(5);
01428 }
01429 break;
01430 case base::Bottom:
01431 for (register int i=mxs; i<=mxe; i++) {
01432 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
01433 MacroType qt = q; qt(1) = 0.; qt(2) = 0.;
01434 MicroType feq = Equilibrium(q);
01435 MicroType feqt = Equilibrium(qt);
01436 f[b+base::idx(i,mye,Nx)] += feqt-feq;
01437
01438 f[b+base::idx(i,mye,Nx)](5) = DataType(2.)*feqt(5)-feq(5);
01439 f[b+base::idx(i,mye,Nx)](3) = DataType(2.)*feqt(3)-feq(3);
01440
01441 f[b+base::idx(i,mye,Nx)](4) = DataType(2.)*feqt(4)-feq(4);
01442 }
01443 break;
01444 case base::Top:
01445 for (register int i=mxs; i<=mxe; i++) {
01446 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
01447 MacroType qt = q; qt(1) = 0.; qt(2) = 0.;
01448 MicroType feq = Equilibrium(q);
01449 MicroType feqt = Equilibrium(qt);
01450 f[b+base::idx(i,mys,Nx)] += feqt-feq;
01451
01452 f[b+base::idx(i,mys,Nx)](8) = DataType(2.)*feqt(8)-feq(8);
01453 f[b+base::idx(i,mys,Nx)](6) = DataType(2.)*feqt(6)-feq(6);
01454
01455 f[b+base::idx(i,mys,Nx)](7) = DataType(2.)*feqt(7)-feq(7);
01456 }
01457 break;
01458 }
01459 }
01460
01461 else if (type==Periodic) {
01462 switch (side) {
01463 case base::Left:
01464 for (register int j=mys; j<=mye; j++) {
01465 f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)];
01466 }
01467 break;
01468 case base::Right:
01469 for (register int j=mys; j<=mye; j++) {
01470 f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)];
01471 }
01472 break;
01473 case base::Bottom:
01474 for (register int i=mxs; i<=mxe; i++) {
01475 f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)];
01476 }
01477 break;
01478 case base::Top:
01479 for (register int i=mxs; i<=mxe; i++) {
01480 f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)];
01481 }
01482 break;
01483 }
01484 }
01485
01486 else if (type==VanDriest) {
01487 DCoords dx = base::GH().worldStep(fvec.stepsize());
01488 DataType dt = dx(0)*S0/U0;
01489 DataType omlam = Omega(dt);
01490 switch (side) {
01491 case base::Left:
01492 for (register int j=mys; j<=mye; j++) {
01493 MacroType q0 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
01494 MicroType feq = Equilibrium(q0);
01495 DataType om = vanDriestRelaxation(f[b+base::idx(mxe+1,j,Nx)],feq,q0(0),q0(2),dx(0),omlam,dt);
01496 if (j>mys)
01497 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.)-om)*f[b+base::idx(mxe+1,j-1,Nx)](5) + om*feq(5);
01498 f[b+base::idx(mxe,j,Nx)](1) = (DataType(1.)-om)*f[b+base::idx(mxe+1,j,Nx)](2) + om*feq(2);
01499 if (j<mye)
01500 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.)-om)*f[b+base::idx(mxe+1,j+1,Nx)](8) + om*feq(8);
01501 }
01502 break;
01503 case base::Right:
01504 for (register int j=mys; j<=mye; j++) {
01505 MacroType q0 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
01506 MicroType feq = Equilibrium(q0);
01507 DataType om = vanDriestRelaxation(f[b+base::idx(mxs-1,j,Nx)],feq,q0(0),q0(2),dx(0),omlam,dt);
01508 if (j>mys)
01509 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.)-om)*f[b+base::idx(mxs-1,j-1,Nx)](4) + om*feq(4);
01510 f[b+base::idx(mxs,j,Nx)](2) = (DataType(1.)-om)*f[b+base::idx(mxs-1,j,Nx)](1) + om*feq(1);
01511 if (j<mye)
01512 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.)-om)*f[b+base::idx(mxs-1,j+1,Nx)](8) + om*feq(8);
01513 }
01514 break;
01515 case base::Bottom:
01516 for (register int i=mxs; i<=mxe; i++) {
01517 MacroType q0 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
01518 MicroType feq = Equilibrium(q0);
01519 DataType om = vanDriestRelaxation(f[b+base::idx(i,mye+1,Nx)],feq,q0(0),q0(1),dx(1),omlam,dt);
01520 if (i>mxs)
01521 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.)-om)*f[b+base::idx(i-1,mye+1,Nx)](7) + om*feq(7);
01522 f[b+base::idx(i,mye,Nx)](3) = (DataType(1.)-om)*f[b+base::idx(i,mye+1,Nx)](6) + om*feq(6);
01523 if (i<mxe)
01524 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.)-om)*f[b+base::idx(i+1,mye+1,Nx)](8) + om*feq(8);
01525 }
01526 break;
01527 case base::Top:
01528 for (register int i=mxs; i<=mxe; i++) {
01529 MacroType q0 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
01530 MicroType feq = Equilibrium(q0);
01531 DataType om = vanDriestRelaxation(f[b+base::idx(i,mys-1,Nx)],feq,q0(0),q0(1),dx(1),omlam,dt);
01532 if (i>mxs)
01533 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.)-om)*f[b+base::idx(i-1,mys-1,Nx)](4) + om*feq(4);
01534 f[b+base::idx(i,mys,Nx)](6) = (DataType(1.)-om)*f[b+base::idx(i,mys-1,Nx)](3) + om*feq(3);
01535 if (i<mxe)
01536 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.)-om)*f[b+base::idx(i+1,mys-1,Nx)](5) + om*feq(5);
01537 }
01538 break;
01539 }
01540 }
01541
01542 else if (type==CSMBC) {
01543 DCoords dx = base::GH().worldStep(fvec.stepsize());
01544 DataType dt_temp = dx(0)/VelocityScale();
01545 switch (side) {
01546 case base::Left:
01547 for (register int j=mys; j<=mye; j++) {
01548
01549 MacroType qxp = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
01550 MacroType qxm = MacroVariables(f[b+base::idx(mxe-1,j,Nx)]);
01551 MacroType qyp = MacroVariables(f[b+base::idx(mxe,j+1,Nx)]);
01552 MacroType qym = MacroVariables(f[b+base::idx(mxe,j-1,Nx)]);
01553 LocalCollisionCSM(f[b+base::idx(mxe,j,Nx)],qxp,qxm,qyp,qym,dx,dt_temp);
01554
01555 }
01556 break;
01557 case base::Right:
01558 for (register int j=mys; j<=mye; j++) {
01559
01560 MacroType qxp = MacroVariables(f[b+base::idx(mxs+1,j,Nx)]);
01561 MacroType qxm = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
01562 MacroType qyp = MacroVariables(f[b+base::idx(mxs,j+1,Nx)]);
01563 MacroType qym = MacroVariables(f[b+base::idx(mxs,j-1,Nx)]);
01564 LocalCollisionCSM(f[b+base::idx(mxs,j,Nx)],qxp,qxm,qyp,qym,dx,dt_temp);
01565
01566 }
01567 break;
01568 case base::Bottom:
01569 for (register int i=mxs+base::NGhosts(); i<=mxe-base::NGhosts(); i++) {
01570
01571 MacroType qxp = MacroVariables(f[b+base::idx(i+1,mye,Nx)]);
01572 MacroType qxm = MacroVariables(f[b+base::idx(i-1,mye,Nx)]);
01573 MacroType qyp = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
01574 MacroType qym = MacroVariables(f[b+base::idx(i,mye-1,Nx)]);
01575 LocalCollisionCSM(f[b+base::idx(i,mye,Nx)],qxp,qxm,qyp,qym,dx,dt_temp);
01576
01577 }
01578 break;
01579 case base::Top:
01580 for (register int i=mxs+base::NGhosts(); i<=mxe-base::NGhosts(); i++) {
01581
01582 MacroType qxp = MacroVariables(f[b+base::idx(i+1,mys,Nx)]);
01583 MacroType qxm = MacroVariables(f[b+base::idx(i-1,mys,Nx)]);
01584 MacroType qyp = MacroVariables(f[b+base::idx(i,mys+1,Nx)]);
01585 MacroType qym = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
01586 LocalCollisionCSM(f[b+base::idx(i,mys,Nx)],qxp,qxm,qyp,qym,dx,dt_temp);
01587
01588 }
01589 break;
01590 }
01591 }
01592
01593 }
01594
01595 virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
01596 const MicroType* f, const point_type* xc, const DataType* distance,
01597 const point_type* normal, DataType* aux=0, const int naux=0,
01598 const int scaling=0) const {
01599
01600 if (type==GFMNoSlipWall || type==GFMSlipWall) {
01601 DataType fac = S0;
01602 if (scaling==base::Physical)
01603 fac /= U0;
01604 for (int n=0; n<nc; n++) {
01605 MacroType q=MacroVariables(f[n]);
01606 DataType u=-q(1);
01607 DataType v=-q(2);
01608
01609
01610 if (naux>=2) {
01611 u += fac*aux[naux*n];
01612 v += fac*aux[naux*n+1];
01613 }
01614 if (type==GFMNoSlipWall) {
01615
01616 q(1) += DataType(2.)*u;
01617 q(2) += DataType(2.)*v;
01618 }
01619 else if (type==GFMSlipWall) {
01620
01621 DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
01622 q(1) += vl*normal[n](0);
01623 q(2) += vl*normal[n](1);
01624 }
01625 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01626 }
01627 }
01628
01629 else if (type==GFMExtrapolation) {
01630 for (int n=0; n<nc; n++) {
01631 MacroType q=MacroVariables(f[n]);
01632 DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1);
01633 q(1) = vl*normal[n](0);
01634 q(2) = vl*normal[n](1);
01635 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01636 }
01637 }
01638
01639 else if (type==GFMComplex) {
01640 DCoords dx = base::GH().worldStep(fvec.stepsize());
01641 DataType h=0, d=dx(0), theta=0.;
01642 DataType fac = S0;
01643 DataType fp1[2], fp2[2], fm[2];
01644 MicroType fw, fwu;
01645 fw(0) = DataType(4.) / DataType(9.);
01646 fw(1) = fw(2) = fw(3) = fw(6) = DataType(1.) / DataType(9.);
01647 fw(4) = fw(5) = fw(7) = fw(8) = DataType(1.) / DataType(36.);
01648 if (scaling==base::Physical)
01649 fac /= U0;
01650 for (int n=0; n<nc; n++) {
01651 DataType nx=std::fabs(normal[n](0));
01652 DataType ny=std::fabs(normal[n](1));
01653 if (nx>ny) { theta = std::atan(ny/nx); }
01654 else { theta = std::atan(nx/ny); }
01655 d = dx(0)/std::cos(theta);
01656
01657
01658 h = std::fabs(distance[n]/d);
01659 fp1[0] = (xc[n](0) - dx(0)*normal[n](0));
01660 fp1[1] = (xc[n](1) - dx(1)*normal[n](1));
01661 fp2[0] = (xc[n](0) - DataType(2.)*dx(0)*normal[n](0));
01662 fp2[1] = (xc[n](1) - DataType(2.)*dx(1)*normal[n](1));
01663 fm[0] = (xc[n](0) + dx(0)*normal[n](0));
01664 fm[1] = (xc[n](1) + dx(1)*normal[n](1));
01665 MacroType q=MacroVariables(f[n]);
01666 MacroType qm=MacroVariables(fvec( base::GH().localCoords((const DataType *) fm) ) );
01667
01668 DataType u=-q(1);
01669 DataType v=-q(2);
01670 DataType um=-qm(1);
01671 DataType vm=-qm(2);
01672
01673
01674 if (naux>=2) {
01675 u += fac*aux[naux*n];
01676 v += fac*aux[naux*n+1];
01677 um += fac*aux[naux*n];
01678 vm += fac*aux[naux*n+1];
01679 }
01680 q(1) += DataType(2.)*u;
01681 q(2) += DataType(2.)*v;
01682 qm(1) += DataType(2.)*vm;
01683 qm(2) += DataType(2.)*vm;
01684
01685
01686
01687
01688
01689
01690
01691 if (h<0.5) {
01692 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = h*(DataType(1.)+2.*h)*f[n]
01693 + (1.-4.*h*h)*fvec( base::GH().localCoords((const DataType *) fp1) )
01694 - h*(1.-2.*h)*fvec( base::GH().localCoords((const DataType *) fp2) )
01695 + 3.*Equilibrium(q);
01696 }
01697 else {
01698 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = 1./(h*(DataType(1.)+2.*h))*f[n]
01699 + (2.*h-1.)/h*fvec( base::GH().localCoords((const DataType *) fp1) )
01700 - (2.*h-1.)/(2.*h+1.)*fvec( base::GH().localCoords((const DataType *) fp2) )
01701 + 3./(h*(2.*h+1))*Equilibrium(q);
01702 }
01703 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01704
01705 }
01706 }
01707
01708 else if (type==GFMComplex2) {
01709 DCoords dx = base::GH().worldStep(fvec.stepsize());
01710 DataType theta=0., h=dx(0), d=0.;
01711 DataType fac = S0;
01712 DataType fp1[2];
01713 if (scaling==base::Physical)
01714 fac /= U0;
01715 for (int n=0; n<nc; n++) {
01716 DataType nx=std::fabs(normal[n](0));
01717 DataType ny=std::fabs(normal[n](1));
01718 if (nx>ny) { theta = std::atan(ny/nx); }
01719 else { theta = std::atan(nx/ny); }
01720 h = DataType(2.)*dx(0)/std::cos(theta);
01721 d = std::fabs(distance[n]);
01722
01723
01724 fp1[0] = (xc[n](0) - dx(0)*normal[n](0));
01725 fp1[1] = (xc[n](1) - dx(1)*normal[n](1));
01726 MacroType q=MacroVariables(f[n]);
01727 MacroType q1=MacroVariables( fvec( base::GH().localCoords((const DataType *) fp1) ) );
01728 DataType u=-q(1);
01729 DataType v=-q(2);
01730
01731
01732
01733 if (naux>=2) {
01734 u += fac*aux[naux*n];
01735 v += fac*aux[naux*n+1];
01736 }
01737 q(1) += DataType(2.)*u + d/h*(fac*aux[naux*n]-q1(1));
01738 q(2) += DataType(2.)*v + d/h*(fac*aux[naux*n+1]-q1(2));
01739 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01740 }
01741 }
01742
01743 }
01744
01745 virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01746 const int skip_ghosts=1) const {
01747 int Nx = fvec.extents(0), Ny = fvec.extents(1);
01748 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01749 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01750 MicroType *f = (MicroType *)fvec.databuffer();
01751 DataType *work = (DataType *)workvec.databuffer();
01752 DCoords dx = base::GH().worldStep(fvec.stepsize());
01753 DataType dt = dx(0)*S0/U0;
01754
01755 if ((cnt>=1 && cnt<=9) || (cnt>=16 && cnt<=18) || (cnt>=23 && cnt<=25)) {
01756 for (register int j=mys; j<=mye; j++)
01757 for (register int i=mxs; i<=mxe; i++) {
01758 MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
01759 switch(cnt) {
01760 case 1:
01761 if (method[0]==0) work[base::idx(i,j,Nx)]=q(0)*R0+rhop;
01762 else
01763 work[base::idx(i,j,Nx)]=q(0)*R0;
01764 break;
01765 case 2:
01766 work[base::idx(i,j,Nx)]=q(1)*U0/S0;
01767 break;
01768 case 3:
01769 work[base::idx(i,j,Nx)]=q(2)*U0/S0;
01770 break;
01771 case 5:
01772 work[base::idx(i,j,Nx)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure());
01773 break;
01774 case 6:
01775 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure();
01776 break;
01777 case 8: {
01778 MicroType feq = Equilibrium(q);
01779 work[base::idx(i,j,Nx)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q(0),Omega(dt),dt),
01780 GasSpeedofSound(),dt);
01781 break;
01782 }
01783 case 9:
01784 work[base::idx(i,j,Nx)]=std::sqrt(q(1)*q(1)+q(2)*q(2))*U0/S0;
01785 break;
01786 case 16:
01787 case 17:
01788 case 18: {
01789 MicroType feq = Equilibrium(q);
01790 DataType omega=Omega(dt);
01791
01792
01793 if (turbulence == LES_Smagorinsky)
01794 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q(0),Omega(dt),dt);
01795 else if (turbulence == WALE) {
01796 DataType dxux(0.), dxuy(0.), dyux(0.), dyuy(0.);
01797 LocalGradVel(f,i,j,Nx,dx,dxux,dxuy,dyux,dyuy);
01798 omega = Omega_WALE(dxux,dxuy,dyux,dyuy,dx,dt);
01799 }
01800 else if (turbulence == CSM) {
01801 DataType dxux(0.), dxuy(0.), dyux(0.), dyuy(0.);
01802 if (skip_ghosts!=1) {
01803 if (i>mxs && i<mxe)
01804 if (j>mys && j<mye)
01805 LocalGradVel(f,i,j,Nx,dx,dxux,dxuy,dyux,dyuy);
01806 }
01807 else
01808 LocalGradVel(f,i,j,Nx,dx,dxux,dxuy,dyux,dyuy);
01809 omega = Omega_CSM(dxux,dxuy,dyux,dyuy,dx,dt);
01810 }
01811 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*DeviatoricStress(f[base::idx(i,j,Nx)],feq,omega)(cnt-16);
01812 break;
01813 }
01814 case 23:
01815 case 24:
01816 case 25: {
01817 DataType omega=Omega(dt);
01818
01819
01820 if (turbulence == LES_Smagorinsky) {
01821 MicroType feq = Equilibrium(q);
01822 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q(0),Omega(dt),dt);
01823 }
01824 else if (turbulence == WALE) {
01825 DataType dxux(0.), dxuy(0.), dyux(0.), dyuy(0.);
01826 LocalGradVel(f,i,j,Nx,dx,dxux,dxuy,dyux,dyuy);
01827 omega = Omega_WALE(dxux,dxuy,dyux,dyuy,dx,dt);
01828 }
01829 else if (turbulence == CSM) {
01830 DataType dxux(0.), dxuy(0.), dyux(0.), dyuy(0.);
01831 if (skip_ghosts!=1) {
01832 if (i>mxs && i<mxe)
01833 if (j>mys && j<mye)
01834 LocalGradVel(f,i,j,Nx,dx,dxux,dxuy,dyux,dyuy);
01835 }
01836 else
01837 LocalGradVel(f,i,j,Nx,dx,dxux,dxuy,dyux,dyuy);
01838 omega = Omega_CSM(dxux,dxuy,dyux,dyuy,dx,dt);
01839 }
01840 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*Stress(f[base::idx(i,j,Nx)],q,omega)(cnt-23);
01841 if (cnt==23 || cnt==25) work[base::idx(i,j,Nx)]-=BasePressure();
01842 break;
01843 }
01844 }
01845 }
01846 }
01847 else if ((cnt>=10 && cnt<=11) || cnt==15 || (cnt>=30 && cnt<=32)) {
01848 macro_grid_data_type qgrid(fvec.bbox());
01849 MacroType *q = (MacroType *)qgrid.databuffer();
01850 for (register int j=mys; j<=mye; j++)
01851 for (register int i=mxs; i<=mxe; i++) {
01852 q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
01853 q[base::idx(i,j,Nx)](0)*=R0;
01854 q[base::idx(i,j,Nx)](1)*=U0/S0;
01855 q[base::idx(i,j,Nx)](2)*=U0/S0;
01856 }
01857
01858 if (skip_ghosts==0) {
01859 switch(cnt) {
01860 case 30:
01861 mxs+=1; mxe-=1;
01862 break;
01863 case 32:
01864 mys+=1; mye-=1;
01865 break;
01866 case 10:
01867 case 11:
01868 case 15:
01869 case 31:
01870 mxs+=1; mxe-=1; mys+=1; mye-=1;
01871 break;
01872 }
01873 }
01874
01875 DataType dxux(0.), dxuy(0.), dyux(0.), dyuy(0.);
01876 for (register int j=mys; j<=mye; j++)
01877 for (register int i=mxs; i<=mxe; i++) {
01878 if (cnt==15 || cnt==30)
01879 dxux=(q[base::idx(i+1,j,Nx)](1)-q[base::idx(i-1,j,Nx)](1))/(2.*dx(0));
01880 if (cnt==15 || cnt==32)
01881 dyuy=(q[base::idx(i,j+1,Nx)](2)-q[base::idx(i,j-1,Nx)](2))/(2.*dx(1));
01882 if (cnt==10 || cnt==11 || cnt==15 || cnt==31) {
01883 dxuy=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(2.*dx(0));
01884 dyux=(q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(2.*dx(1));
01885 }
01886
01887 switch (cnt) {
01888 case 10:
01889 work[base::idx(i,j,Nx)]=dxuy-dyux;
01890 break;
01891 case 11:
01892 work[base::idx(i,j,Nx)]=std::fabs(dxuy-dyux);
01893 break;
01894 case 15:
01895 work[base::idx(i,j,Nx)]=DataType(-0.5)*(4.*dxux*dxux+8.*dyux*dxuy+4.*dyuy*dyuy);
01896 break;
01897 case 30:
01898 work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dxux;
01899 break;
01900 case 31:
01901 work[base::idx(i,j,Nx)]=q[base::idx(i,j,Nx)](0)*GasViscosity()*(dxuy+dyux);
01902 break;
01903 case 32:
01904 work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dyuy;
01905 break;
01906 }
01907 }
01908 }
01909 }
01910
01911 virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01912 const int skip_ghosts=1) const {
01913 int Nx = fvec.extents(0), Ny = fvec.extents(1);
01914 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01915 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01916 MicroType *f = (MicroType *)fvec.databuffer();
01917 DataType *work = (DataType *)workvec.databuffer();
01918
01919 for (register int j=mys; j<=mye; j++)
01920 for (register int i=mxs; i<=mxe; i++) {
01921 switch(cnt) {
01922 case 1:
01923 f[base::idx(i,j,Nx)](0)=work[base::idx(i,j,Nx)]/R0;
01924 break;
01925 case 2:
01926 f[base::idx(i,j,Nx)](1)=work[base::idx(i,j,Nx)]/(U0/S0);
01927 break;
01928 case 3:
01929 f[base::idx(i,j,Nx)](2)=work[base::idx(i,j,Nx)]/(U0/S0);
01930 f[base::idx(i,j,Nx)] = Equilibrium((const MacroType &)f[base::idx(i,j,Nx)]);
01931 break;
01932 }
01933 }
01934 }
01935
01936 virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
01937 const int verbose) const {
01938 int Nx = fvec.extents(0);
01939 int b = fvec.bottom();
01940 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01941 BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
01942 int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01943 int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01944 int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01945 int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01946 MicroType *f = (MicroType *)fvec.databuffer();
01947
01948 int result = 1;
01949 for (register int j=mys; j<=mye; j++)
01950 for (register int i=mxs; i<=mxe; i++)
01951 for (register int l=0; l<base::NMicroVar(); l++)
01952 if (!(f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l)>-1.e37)) {
01953 result = 0;
01954 if (verbose) {
01955 DataType WLB[2];
01956 base::GH().wlb(WLB);
01957 DCoords dx = base::GH().worldStep(fvec.stepsize());
01958 std::cerr << "D2Q9-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << ")(" << l << ")="
01959 << f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l) << " " << fvec.bbox()
01960 << " Coords: (" << (i+0.5)*dx(0)+WLB[0] << "," << (j+0.5)*dx(1)+WLB[1] << ")"<< std::endl;
01961 }
01962 }
01963 return result;
01964 }
01965
01966 inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
01967 TensorType P;
01968 if (stressPath==0) {
01969 if (method[1]==0) {
01970 P(0) = q(0)*q(1)*q(1)-(f(1)+f(2)+f(4)+f(5)+f(7)+f(8))+cs2*q(0);
01971 P(1) = q(0)*q(1)*q(2)-(f(4)-f(5)-f(7)+f(8));
01972 P(2) = q(0)*q(2)*q(2)-(f(3)+f(4)+f(5)+f(6)+f(7)+f(8))+cs2*q(0);
01973 P *= -(DataType(1.0)-DataType(0.5)*om);
01974 }
01975 else {
01976 MicroType feq = Equilibrium(q);
01977 P = DeviatoricStress(f,feq,om);
01978 }
01979 P(0) -= LatticeBasePressure(q(0));
01980 P(2) -= LatticeBasePressure(q(0));
01981 }
01982 else if (stressPath==1)
01983 P = Stress_velocitySpace(f,Equilibrium(q),om);
01984 return P;
01985 }
01986
01987 inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
01988 TensorType Sigma;
01989 if (stressPath==0) {
01990 if (method[2]==0) {
01991 MicroType fneq = feq - f;
01992 Sigma(0) = fneq(1)+fneq(2)+fneq(4)+fneq(5)+fneq(7)+fneq(8);
01993 Sigma(1) = fneq(4)-fneq(5)-fneq(7)+fneq(8);
01994 Sigma(2) = fneq(3)+fneq(4)+fneq(5)+fneq(6)+fneq(7)+fneq(8);
01995 Sigma *= -(DataType(1.0)-DataType(0.5)*om);
01996 }
01997 else {
01998 MacroType q = MacroVariables(f);
01999 Sigma = Stress(f,q,om);
02000 Sigma(0) += LatticeBasePressure(q(0));
02001 Sigma(2) += LatticeBasePressure(q(0));
02002 }
02003 }
02004 else if (stressPath==1)
02005 Sigma = DeviatoricStress_velocitySpace(f,feq,om);
02006 return Sigma;
02007 }
02008
02009
02010 virtual int NMethodOrder() const { return 2; }
02011
02012 inline const DataType& L0() const { return base::LengthScale(); }
02013 inline const DataType& T0() const { return base::TimeScale(); }
02014 inline void SetDensityScale(const DataType r0) { R0 = r0; }
02015 inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
02016 inline void SetSpeedUp(const DataType s0) { S0 = s0; }
02017 inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
02018
02019 inline const DataType& DensityScale() const { return R0; }
02020 inline const DataType VelocityScale() const { return U0/S0; }
02021
02022 inline const DataType& SpeedUp() const { return S0; }
02023
02024
02025 inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
02026 inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
02027 inline virtual DataType LatticeBasePressure(const DataType rho) const {
02028 return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
02029 }
02030
02031
02032 inline void SetGas(DataType rho, DataType nu, DataType cs) {
02033 rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
02034 SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
02035 }
02036 inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
02037 inline const MacroType Stress_velocitySpace(const MicroType &f, const MicroType &feq, const DataType om) const {
02038 if (method[0]==3) {
02039 MacroType P;
02040 MacroType q = MacroVariables(f);
02042 DataType cl = DataType(1.0);
02043 DataType cmu_ = cl - q(1);
02044 DataType cmu2_ = cmu_*cmu_;
02045 DataType cpu_ = cl + q(1);
02046 DataType cpu2_ = cpu_*cpu_;
02047 DataType cmv_ = cl - q(2);
02048 DataType cmv2_ = cmv_*cmv_;
02049 DataType cpv_ = cl + q(2);
02050 DataType cpv2_ = cpv_*cpv_;
02051 DataType u2 = q(1)*q(1);
02052 DataType v2 = q(2)*q(2);
02053
02054 P(0) = (f(0)+f(3)+f(6))*u2 + (f(2)+f(5)+f(8))*cpu2_
02055 + (f(1)+f(4)+f(7))*cmu2_;
02056 P(1) = f(8)*cpu_*cpv_ + f(2)*q(2)*cpu_ + f(6)*q(1)*cpv_ - f(5)*cpu_*cmv_
02057 - f(7)*cpv_*cmu_ + f(0)*q(1)*q(2) - f(1)*q(2)*cmu_ - f(3)*q(1)*cmv_
02058 + f(4)*cmu_*cmv_;
02059 P(2) = (f(6)+f(7)+f(8))*cpv2_ + (f(0)+f(1)+f(2))*v2
02060 + (f(3)+f(4)+f(5))*cmv2_;
02061 return P/std::sqrt(DataType(2.0));
02062 }
02063 else {
02064 MacroType P = DeviatoricStress_velocitySpace(f,feq,om);
02065 MacroType q = MacroVariables(f);
02066 P(0) -= cs2*q(0); P(2) -= cs2*q(0);
02067 return P;
02068 }
02069 }
02070 inline const MacroType DeviatoricStress_velocitySpace(const MicroType &f, const MicroType &feq, const DataType om) const {
02071 if (method[0]==3) {
02072 TensorType Sigma = Stress_velocitySpace(f,feq,om);
02073 MacroType q = MacroVariables(f);
02074 Sigma(0) += cs2*q(0); Sigma(2) += cs2*q(0);
02075 return Sigma;
02076 }
02077 else {
02078 MicroType fneq;
02079 TensorType Sigma;
02080 DataType cl = DataType(1.0);
02081 fneq = (f - feq);
02082 Sigma(0) = (cl*cl-DataType(0.5))*(fneq(1)+fneq(2)+fneq(4)+fneq(5)+fneq(7)+fneq(8))
02083 -DataType(0.5)*(fneq(3)+fneq(6));
02084 Sigma(1) = (cl*cl-DataType(0.5))*(fneq(4)-fneq(5)-fneq(7)+fneq(8));
02085 Sigma(2) = -DataType(0.5)*(fneq(1)+fneq(2))
02086 +(cl*cl-DataType(0.5))*(fneq(3)+fneq(4)+fneq(5)+fneq(6)+fneq(7)+fneq(8));
02087 return (DataType(1.0)-DataType(0.5)*om)*Sigma;
02088 }
02089 }
02090 inline const MacroType StrainComponents(const DataType rho, const MacroType& Sigma, const DataType om_old, const DataType cs_old) const {
02091 DataType omTmp = rho*cs2/om_old;
02092 DataType rhocspcsmTmp = DataType(2.0)*rho*cs2*cs2*cs_old*cs_old*cs_old*cs_old;
02093 DataType sigma2 = Magnitude(Sigma);
02094 DataType stmp = ( (rhocspcsmTmp*sigma2 > (mfp/U0)) ? (-omTmp + std::sqrt( (omTmp)*(omTmp) + rhocspcsmTmp*sigma2 ) )/( rhocspcsmTmp*sigma2 )
02095 : (-omTmp + std::sqrt( (omTmp)*(omTmp) + mfp/U0 ) )/( mfp/U0 ));
02096 MacroType Strain = stmp*Sigma;
02097 return Strain;
02098 }
02099 inline const DataType Strain(const DataType rho, const MacroType& Sigma, const DataType om_old, const DataType cs_old) const {
02100 return Magnitude(StrainComponents(rho,Sigma,om_old,cs_old));
02101 }
02102 inline const DataType StrainLaminar(const DataType rho, const MacroType& Sigma, const DataType om) const {
02103 return om/(DataType(2.)*rho*cs2)*Magnitude(Sigma);
02104 }
02105 inline const DataType Magnitude(const MacroType& A) const {
02106 return std::sqrt(DataType(2.0)*(A(0)*A(0) + DataType(2.0)*A(1)*A(1) + A(2)*A(2)));
02107 }
02108 inline const DataType Omega_LES_Smagorinsky( MicroType &f, const MicroType &feq, const DataType rho,
02109 const DataType om, const DataType dt) const {
02118 TensorType Sigma = DeviatoricStress(f,feq,om)/(DataType(1.0)-DataType(0.5)*om);
02119 DataType S = Strain(rho,Sigma,om,Cs_Smagorinsky);
02120 DataType nu_t = Cs_Smagorinsky*Cs_Smagorinsky*S;
02121 if (nu_t < DataType(0.)) nu_t = DataType(0.);
02122 DataType om_eff = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
02123
02124 return ( om_eff );
02125 }
02126 inline const int EquilibriumType() const { return method[0]; }
02127 inline const int TurbulenceType() const { return turbulence; }
02128 inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
02129 inline const DataType& GasDensity() const { return rhop; }
02130 inline const DataType& GasViscosity() const { return nup; }
02131 inline const DataType& GasSpeedofSound() const { return csp; }
02132 inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
02133 { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
02134 inline virtual const MacroType NormalDerivative(const MacroType qa,const MacroType qb, const MacroType qc) const
02135 { return DataType(-1.5)*qa + DataType(2.)*qb - DataType(0.5)*qc; }
02136 inline virtual Vector<DataType,2> WaveAmplitudes(const DataType rho, const DataType vn, const DataType drhodn, const DataType dvndn, const DataType dvtdn) const {
02137 Vector<DataType,2> L;
02138 L(0) = (vn - LatticeSpeedOfSound())*(cs2*drhodn - rho*LatticeSpeedOfSound()*dvndn);
02139 L(1) = vn*dvtdn;
02140 return L;
02141 }
02142
02143 inline virtual DataType vanDriestRelaxation(const MicroType f, const MicroType feq, const DataType rho, const DataType ut, const DataType dx, const DataType om, const DataType dt) const {
02144 DataType yplus = std::sqrt(dx*(DataType(2.)*std::fabs(ut))/nup );
02145 DataType Csd = DataType(0.17)*dx*(DataType(1.)-std::exp(-yplus/DataType(25.)));
02146 TensorType Sigma = DeviatoricStress(f,feq,om)/(DataType(1.0)-DataType(0.5)*om);
02147 DataType S = Strain(rho,Sigma,om,Cs_Smagorinsky);
02148 DataType nu_t = Csd*Csd*S;
02149 return (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
02150 }
02151
02152 inline virtual DataType Omega_WALE(const DataType dxux, const DataType dxuy,
02153 const DataType dyux, const DataType dyuy,
02154 const DCoords dx, const DataType dt) const {
02159 DataType S2 = DataType(0.5 )*(dxuy + dyux)*(dxuy + dyux) + dxux*dxux + dyuy*dyuy;
02160 if (std::isfinite(S2)==0 ) { S2 = DataType(0.0); }
02161 DataType O2 = DataType(0.5)*(dxuy - dyux)*(dxuy - dyux);
02162 if (std::isfinite(O2)==0 ) { O2 = DataType(0.0); }
02163 DataType IV = -((dxuy - dyux)*(dxuy - dyux) *(DataType(2.)*dxux*dxux + dxuy*dxuy + DataType(2.)*dxuy*dyux + dyux*dyux + DataType(2.)*dyuy*dyuy))/DataType(8.);
02164 if (std::isfinite(IV)==0 ) { IV = DataType(0.0); }
02165 DataType Sdij2 = DataType(1./6.)*(S2*S2+O2*O2)+DataType(2./3.)*S2*O2+DataType(2.)*IV;
02166 DataType eddyVisc = std::pow(Sdij2,DataType(1.5))/( std::pow(S2,DataType(2.5)) + std::pow(Sdij2,DataType(1.25)) );
02167 if (std::isfinite(eddyVisc)==0 ) { eddyVisc = DataType(0.0); }
02168 if (eddyVisc<DataType(0.0)) { (std::cout << "NEG="<<eddyVisc<<std::endl).flush(); eddyVisc=DataType(0.0);}
02169 DataType nu_t = DataType(0.25)*dx(0)*dx(1)*std::sqrt(DataType(2.)*eddyVisc);
02170 DataType omega = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
02171 if (omega < DataType(0.5)) { std::cout <<" nup="<<nup<<" nu_t="<<nu_t<<"om="<<omega << " lower limit\t"; omega = DataType(0.5); }
02172 return omega;
02173 }
02174
02175 inline virtual void LocalCollisionWALE( MicroType &f, const MacroType qxp, const MacroType qxm,
02176 const MacroType qyp, const MacroType qym,
02177 const DCoords dx, const DataType dt) const {
02178 DataType dxux=(qxp(1)-qxm(1))/(2.*dx(0))*U0/S0;
02179 DataType dxuy=(qxp(2)-qxm(2))/(2.*dx(0))*U0/S0;
02180 DataType dyux=(qyp(1)-qym(1))/(2.*dx(1))*U0/S0;
02181 DataType dyuy=(qyp(2)-qym(2))/(2.*dx(1))*U0/S0;
02182 DataType omega = Omega_WALE(dxux,dxuy,dyux,dyuy,dx,dt);
02183 MicroType feq = Equilibrium(MacroVariables(f));
02184 for (register int qi=0; qi<9; qi++)
02185 f(qi) = (DataType(1.)-omega)*f(qi) + omega*feq(qi);
02186 }
02187
02188 virtual void CollisionWALE(vec_grid_data_type &fvec, const double& dt) const {
02189 int Nx = fvec.extents(0), Ny = fvec.extents(1);
02190 int mxs = base::NGhosts(), mys = base::NGhosts();
02191 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
02192 DCoords dx = base::GH().worldStep(fvec.stepsize());
02193 MicroType *f = (MicroType *)fvec.databuffer();
02194 MicroType feq;
02195 DataType omega = Omega(dt);
02196 macro_grid_data_type qgrid(fvec.bbox());
02197 MacroType *q = (MacroType *)qgrid.databuffer();
02198 int x=0, y=0;
02199 for (register int j=mys-1; j<=mye+1; j++) {
02200 if (j==mys-1) y=1;
02201 else if (j==mye+1) y=-1;
02202 else y=0;
02203 for (register int i=mxs-1; i<=mxe+1; i++) {
02204 if (i==mxs-1) x=1;
02205 else if (i==mxe+1) x=-1;
02206 else x=0;
02207 q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i+x,j+y,Nx)]);
02208 q[base::idx(i,j,Nx)](1)*=U0/S0;
02209 q[base::idx(i,j,Nx)](2)*=U0/S0;
02210 }
02211 }
02212 DataType dxux, dxuy, dyux, dyuy;
02213 for (register int j=mys; j<=mye; j++)
02214 for (register int i=mxs; i<=mxe; i++) {
02215 dxux=(q[base::idx(i+1,j,Nx)](1)-q[base::idx(i-1,j,Nx)](1))/(DataType(2.)*dx(0));
02216 dxuy=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(DataType(2.)*dx(0));
02217 dyux=(q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(DataType(2.)*dx(1));
02218 dyuy=(q[base::idx(i,j+1,Nx)](2)-q[base::idx(i,j-1,Nx)](2))/(DataType(2.)*dx(1));
02219 omega = Omega_WALE(dxux,dxuy,dyux,dyuy,dx,dt);
02220 feq = Equilibrium(MacroVariables(f[base::idx(i,j,Nx)]));
02221 for (register int qi=0; qi<9; qi++)
02222 f[base::idx(i,j,Nx)](qi) = (DataType(1.)-omega)*f[base::idx(i,j,Nx)](qi) + omega*feq(qi);
02223 }
02224 }
02225
02226 inline virtual DataType Omega_CSM(const DataType dxux, const DataType dxuy,
02227 const DataType dyux, const DataType dyuy,
02228 const DCoords dx, const DataType dt) const {
02234 DataType S2, Q, E, FCS, C, eddyVisc, nu_t;
02235 S2 = DataType(0.5)*(dxuy + dyux)*(dxuy + dyux) + dxux*dxux + dyuy*dyuy;
02236 Q = DataType(-2.)*(dxux*dxux+dyuy*dyuy)-DataType(4.)*dyux*dxuy;
02237 E = DataType(0.5)*(dxux*dxux+dxuy*dxuy+dyux*dyux+dyuy*dyuy);
02238 if ( std::isfinite(S2)==0 || std::isfinite(Q)==0 || std::isfinite(E)==0 ) {
02239 nu_t = DataType(0.0);
02240 }
02241 else {
02242 FCS = Q/E;
02243 if (std::isfinite(FCS)==0) {
02244 nu_t = DataType(0.0);
02245 }
02246 else {
02247 if (FCS<DataType(-1.0)) { FCS = DataType(-1.0); }
02248 else if (FCS>DataType(1.0)) { FCS = DataType(1.0); }
02249 C = DataType(1./22.)*std::pow(std::fabs(FCS),DataType(1.5))*(DataType(1.)-FCS);
02250 eddyVisc = dx(0)*dx(1)*std::sqrt(DataType(2.)*S2);
02251 if (std::isfinite(eddyVisc)==0 ) { std::cout <<" eddyVisc="<< eddyVisc << " reset local\t"; eddyVisc = DataType(0.0); }
02252 if (std::isfinite(C)==0 ) { std::cout <<" C="<< C << " reset local\t"; C = DataType(0.0); }
02253 nu_t = C*eddyVisc;
02254 }
02255 }
02256 DataType omega = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
02257 return omega;
02258 }
02259
02260 inline virtual void LocalCollisionCSM( MicroType &f, const MacroType qxp, const MacroType qxm,
02261 const MacroType qyp, const MacroType qym,
02262 const DCoords dx, const DataType dt) const {
02263 DataType dxux=(qxp(1)-qxm(1))/(DataType(2.)*dx(0))*U0/S0;
02264 DataType dxuy=(qxp(2)-qxm(2))/(DataType(2.)*dx(0))*U0/S0;
02265 DataType dyux=(qyp(1)-qym(1))/(DataType(2.)*dx(1))*U0/S0;
02266 DataType dyuy=(qyp(2)-qym(2))/(DataType(2.)*dx(1))*U0/S0;
02267 DataType omega = Omega_CSM(dxux,dxuy,dyux,dyuy,dx,dt);
02268 MicroType feq = Equilibrium(MacroVariables(f));
02269 for (register int qi=0; qi<9; qi++)
02270 f(qi) = (DataType(1.)-omega)*f(qi) + omega*feq(qi);
02271 }
02272
02273 virtual void CollisionCSM(vec_grid_data_type &fvec, const double& dt) const {
02274 int Nx = fvec.extents(0), Ny = fvec.extents(1);
02275 int mxs = base::NGhosts(), mys = base::NGhosts();
02276 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
02277 DCoords dx = base::GH().worldStep(fvec.stepsize());
02278 MicroType *f = (MicroType *)fvec.databuffer();
02279 MicroType feq;
02280 DataType omega = Omega(dt);
02281 macro_grid_data_type qgrid(fvec.bbox());
02282 MacroType *q = (MacroType *)qgrid.databuffer();
02283 for (register int j=mys-1; j<=mye+1; j++) {
02284 for (register int i=mxs-1; i<=mxe+1; i++) {
02285 q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
02286 q[base::idx(i,j,Nx)](1)*=U0/S0;
02287 q[base::idx(i,j,Nx)](2)*=U0/S0;
02288 }
02289 }
02290 DataType dxux, dxuy, dyux, dyuy;
02291 for (register int j=mys; j<=mye; j++)
02292 for (register int i=mxs; i<=mxe; i++) {
02293 dxux=(q[base::idx(i+1,j,Nx)](1)-q[base::idx(i-1,j,Nx)](1))/(DataType(2.)*dx(0));
02294 dxuy=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(DataType(2.)*dx(0));
02295 dyux=(q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(DataType(2.)*dx(1));
02296 dyuy=(q[base::idx(i,j+1,Nx)](2)-q[base::idx(i,j-1,Nx)](2))/(DataType(2.)*dx(1));
02297 omega = Omega_CSM(dxux,dxuy,dyux,dyuy,dx,dt);
02298 feq = Equilibrium(MacroVariables(f[base::idx(i,j,Nx)]));
02299 for (register int qi=0; qi<9; qi++)
02300 f[base::idx(i,j,Nx)](qi) = (DataType(1.)-omega)*f[base::idx(i,j,Nx)](qi) + omega*feq(qi);
02301 }
02302 }
02303
02304 inline void LocalGradVel(const MicroType *f, const int i, const int j, const int Nx, const DCoords dx,
02305 DataType &dxux, DataType &dxuy, DataType &dyux, DataType &dyuy) const {
02306 MacroType fxp=MacroVariables(f[base::idx(i+1,j,Nx)]);
02307 MacroType fxm=MacroVariables(f[base::idx(i-1,j,Nx)]);
02308 MacroType fyp=MacroVariables(f[base::idx(i,j+1,Nx)]);
02309 MacroType fym=MacroVariables(f[base::idx(i,j-1,Nx)]);
02310 dxux = (fxp(1)-fxm(1))/(DataType(2.)*dx(0));
02311 dxuy = (fxp(2)-fxm(2))/(DataType(2.)*dx(0));
02312 dyux = (fyp(1)-fym(1))/(DataType(2.)*dx(1));
02313 dyuy = (fyp(2)-fym(2))/(DataType(2.)*dx(1));
02314 }
02315
02316
02317 inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
02318 inline virtual DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); }
02319 inline virtual DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); }
02320
02321 protected:
02322 DataType cs2, cs22, cssq;
02323 DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
02324 DataType Cs_Smagorinsky, turbulence, mfp;
02325 int method[3], mdx[NUMMICRODIST], mdy[NUMMICRODIST];
02326 int stressPath;
02327 };
02328
02329
02330 #endif