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