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 template <class DataType>
00036 class LBMD2Q9 : public LBMBase<Vector<DataType,9>, Vector<DataType,3>, 2> {
00037 typedef LBMBase<Vector<DataType,9>, Vector<DataType,3>, 2> base;
00038 public:
00039 typedef typename base::vec_grid_data_type vec_grid_data_type;
00040 typedef typename base::grid_data_type grid_data_type;
00041 typedef typename base::MicroType MicroType;
00042 typedef typename base::MacroType MacroType;
00043 typedef GridData<MacroType,2> macro_grid_data_type;
00044 typedef typename base::SideName SideName;
00045 typedef typename base::point_type point_type;
00046 typedef typename base::MacroType TensorType;
00047
00048 enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
00049 enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall,
00050 CharacteristicOutlet, CharacteristicInlet, NoSlipWallEntranceExit };
00051 enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall };
00052 enum TurbulenceModel { laminar, LES_Smagorinsky };
00053
00054 LBMD2Q9() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.) {
00055 cs2 = DataType(1.)/DataType(3.);
00056 cs22 = DataType(2.)*cs2;
00057 cssq = DataType(2.)/DataType(9.);
00058 method[0] = 2; method[1] = 0; method[2] = 0; method[3] = 0;
00059 turbulence = laminar;
00060 mdx[0]=mdx[3]=mdx[6]=0; mdx[1]=mdx[4]=mdx[7]=1; mdx[2]=mdx[5]=mdx[8]=-1;
00061 mdy[0]=mdy[1]=mdy[2]=0; mdy[3]=mdy[4]=mdy[5]=1; mdy[6]=mdy[7]=mdy[8]=-1;
00062 }
00063
00064 virtual ~LBMD2Q9() {}
00065
00066 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00067 base::LocCtrl = Ctrl.getSubDevice(prefix+"D2Q9");
00068 RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00069 RegisterAt(base::LocCtrl,"Method(2)",method[1]);
00070 RegisterAt(base::LocCtrl,"Method(3)",method[2]);
00071 RegisterAt(base::LocCtrl,"Method(4)",method[3]);
00072 RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00073 RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00074 RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00075 RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00076 RegisterAt(base::LocCtrl,"Gas_nu",nup);
00077 RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00078 RegisterAt(base::LocCtrl,"Gas_W",Wp);
00079 RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00080 RegisterAt(base::LocCtrl,"SpeedUp",S0);
00081 }
00082
00083 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00084 base::SetupData(gh,ghosts);
00085 if (rhop>0.) R0 = rhop;
00086 if (csp>0.) {
00087 cs2p = csp*csp;
00088 SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00089 }
00090 std::cout << "D2Q9: Gas_rho=" << rhop << " Gas_Cs=" << csp
00091 << " Gas_nu=" << nup << std::endl;
00092 std::cout << "D2Q9: L0=" << L0() << " T0=" << T0() << " U0=" << U0
00093 << " R0=" << R0 << std::endl;
00094 WriteInit();
00095 }
00096
00097 virtual void WriteInit() const {
00098 int me = MY_PROC;
00099 if (me == VizServer) {
00100 std::ofstream ofs("chem.dat", std::ios::out);
00101 ofs << "RU " << Rp << std::endl;
00102 ofs << "PA " << BasePressure() << std::endl;
00103 ofs << "Sp(01) Vicosity" << std::endl;
00104 ofs << "W(01) " << Wp << std::endl;
00105 ofs << "Sp(02) |Velocity|" << std::endl;
00106 ofs << "W(02) " << 1.0 << std::endl;
00107 ofs << "Sp(03) Vorticity" << std::endl;
00108 ofs << "W(03) " << 1.0 << std::endl;
00109 ofs << "Sp(04) |Vorticity|" << std::endl;
00110 ofs << "W(04) " << 1.0 << std::endl;
00111 ofs.close();
00112 }
00113 }
00114
00115 inline virtual MacroType MacroVariables(const MicroType &f) const {
00116 MacroType q;
00117 q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8);
00118 q(1)=(f(1)-f(2)+f(4)-f(5)+f(7)-f(8))/q(0);
00119 q(2)=(f(3)+f(4)+f(5)-f(6)-f(7)-f(8))/q(0);
00120 return q;
00121 }
00122
00123 inline virtual MicroType Equilibrium(const MacroType &q) const {
00124 MicroType feq;
00125
00126 DataType rho = q(0);
00127 DataType u = q(1);
00128 DataType v = q(2);
00129
00130 DataType ri, rt0, rt1, rt2;
00131 if (method[0]==0) {
00132 ri = DataType(1.);
00133 rt0 = R0 * DataType(4.) / DataType(9.);
00134 rt1 = R0 / DataType(9.);
00135 rt2 = R0 / DataType(36.);
00136 }
00137 if (method[0]==1) {
00138 ri = rho;
00139 rt0 = DataType(4.) / DataType(9.);
00140 rt1 = DataType(1.) / DataType(9.);
00141 rt2 = DataType(1.) / DataType(36.);
00142 }
00143 else if (method[0]==2) {
00144 ri = DataType(1.);
00145 rt0 = rho * DataType(4.) / DataType(9.);
00146 rt1 = rho / DataType(9.);
00147 rt2 = rho / DataType(36.);
00148 }
00149
00150 DataType usq = u * u;
00151 DataType vsq = v * v;
00152 DataType sumsq = (usq + vsq) / cs22;
00153 DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00154 DataType u2 = usq / cssq;
00155 DataType v2 = vsq / cssq;
00156 DataType ui = u / cs2;
00157 DataType vi = v / cs2;
00158 DataType uv = ui * vi;
00159
00160 feq(0) = rt0 * (ri - sumsq);
00161
00162 feq(1) = rt1 * (ri - sumsq + u2 + ui);
00163 feq(2) = rt1 * (ri - sumsq + u2 - ui);
00164 feq(3) = rt1 * (ri - sumsq + v2 + vi);
00165 feq(6) = rt1 * (ri - sumsq + v2 - vi);
00166
00167 feq(4) = rt2 * (ri + sumsq2 + ui + vi + uv);
00168 feq(5) = rt2 * (ri + sumsq2 - ui + vi - uv);
00169 feq(7) = rt2 * (ri + sumsq2 + ui - vi - uv);
00170 feq(8) = rt2 * (ri + sumsq2 - ui - vi + uv);
00171
00172
00173 if (method[3]>0) {
00174 DataType cs4 = DataType(1.)/(DataType(4.)*cs2*cs2);
00175 DataType upv = (u+v)*cs4*(usq-(u+v)*(u+v)+vsq);
00176 DataType umv = (u-v)*cs4*(usq-(u-v)*(u-v)+vsq);
00177
00178 feq(1) -= rt1 * u*vsq*cs4;
00179 feq(2) += rt1 * u*vsq*cs4;
00180 feq(3) -= rt1 * usq*v*cs4;
00181 feq(6) += rt1 * usq*v*cs4;
00182
00183 feq(4) -= rt2 * upv;
00184 feq(5) += rt2 * umv;
00185 feq(7) -= rt2 * umv;
00186 feq(8) += rt2 * upv;
00187 }
00188 return feq;
00189 }
00190
00191 inline virtual void Collision(MicroType &f, const DataType dt) const {
00192 MacroType q = MacroVariables(f);
00193 MicroType feq = Equilibrium(q);
00194
00195 DataType omega;
00196 if (turbulence == laminar)
00197 omega = Omega(dt);
00198 else if (turbulence == LES_Smagorinsky)
00199 omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00200
00201 f=(DataType(1.)-omega)*f + omega*feq;
00202 }
00203
00204 virtual int IncomingIndices(const int side, int indices[]) const {
00205 switch (side) {
00206 case base::Left:
00207 indices[0] = 1; indices[1] = 4; indices[2] = 7;
00208 break;
00209 case base::Right:
00210 indices[0] = 2; indices[1] = 5; indices[2] = 8;
00211 break;
00212 case base::Bottom:
00213 indices[0] = 3; indices[1] = 4; indices[2] = 5;
00214 break;
00215 case base::Top:
00216 indices[0] = 6; indices[1] = 7; indices[2] = 8;
00217 break;
00218 }
00219 return 3;
00220 }
00221
00222 virtual int OutgoingIndices(const int side, int indices[]) const {
00223 switch (side) {
00224 case base::Left:
00225 indices[0] = 2; indices[1] = 5; indices[2] = 8;
00226 break;
00227 case base::Right:
00228 indices[0] = 1; indices[1] = 4; indices[2] = 7;
00229 break;
00230 case base::Bottom:
00231 indices[0] = 6; indices[1] = 7; indices[2] = 8;
00232 break;
00233 case base::Top:
00234 indices[0] = 3; indices[1] = 4; indices[2] = 5;
00235 break;
00236 }
00237 return 3;
00238 }
00239
00240
00241 virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb,
00242 const int side) const {
00243 int Nx = fvec.extents(0);
00244 int b = fvec.bottom();
00245 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00246 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00247 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00248 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00249 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00250 MicroType *f = (MicroType *)fvec.databuffer();
00251
00252 #ifdef DEBUG_PRINT_FIXUP
00253 ( comm_service::log() << "Reverse streaming at Side: " << side << " of "
00254 << fvec.bbox() << " using " << bb << "\n").flush();
00255 #endif
00256
00257 switch (side) {
00258 case base::Left:
00259 for (register int j=mys; j<=mye; j++) {
00260 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](2);
00261 f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](8);
00262 f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](5);
00263 }
00264 break;
00265 case base::Right:
00266 for (register int j=mys; j<=mye; j++) {
00267 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](1);
00268 f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](7);
00269 f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](4);
00270 }
00271 break;
00272 case base::Bottom:
00273 for (register int i=mxs; i<=mxe; i++) {
00274 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](6);
00275 f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](8);
00276 f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](7);
00277 }
00278 break;
00279 case base::Top:
00280 for (register int i=mxs; i<=mxe; i++) {
00281 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](3);
00282 f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](5);
00283 f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](4);
00284 }
00285 break;
00286 }
00287 }
00288
00289 virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec,
00290 const BBox &bb, const double& dt) const {
00291 int Nx = fvec.extents(0);
00292 int b = fvec.bottom();
00293 MicroType *f = (MicroType *)fvec.databuffer();
00294 MicroType *of = (MicroType *)ovec.databuffer();
00295
00296 #ifdef DEBUG_PRINT_FIXUP
00297 ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox()
00298 << " using " << bb << "\n").flush();
00299 #endif
00300
00301 assert (fvec.extents(0)==ovec.extents(0));
00302 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) &&
00303 fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1));
00304 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00305 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00306 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00307 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00308
00309 for (register int j=mys; j<=mye; j++)
00310 for (register int i=mxs; i<=mxe; i++) {
00311 for (register int n=0; n<base::NMicroVar(); n++)
00312 f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00313 Collision(f[b+base::idx(i,j,Nx)],dt);
00314 }
00315 }
00316
00317 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00318 const double& t, const double& dt, const int& mpass) const {
00319
00320
00321 DCoords dx = base::GH().worldStep(fvec.stepsize());
00322
00323 if (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR) {
00324 std::cerr << "LBMD2Q9 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00325 << " to be equal." << std::endl << std::flush;
00326 std::exit(-1);
00327 }
00328 if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR) {
00329 std::cerr << "LBMD2Q9 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00330 << " to be equal."
00331 << " dt=" << dt << " TO=" << T0()
00332 << std::endl << std::flush;
00333 std::exit(-1);
00334 }
00335
00336 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00337 MicroType *f = (MicroType *)fvec.databuffer();
00338
00339
00340 for (register int j=Ny-1; j>=1; j--)
00341 for (register int i=0; i<=Nx-2; i++) {
00342 f[base::idx(i,j,Nx)](3) = f[base::idx(i,j-1,Nx)](3);
00343 f[base::idx(i,j,Nx)](5) = f[base::idx(i+1,j-1,Nx)](5);
00344 }
00345
00346 for (register int j=Ny-1; j>=1; j--)
00347 for (register int i=Nx-1; i>=1; i--) {
00348 f[base::idx(i,j,Nx)](1) = f[base::idx(i-1,j,Nx)](1);
00349 f[base::idx(i,j,Nx)](4) = f[base::idx(i-1,j-1,Nx)](4);
00350 }
00351
00352 for (register int j=0; j<=Ny-2; j++)
00353 for (register int i=Nx-1; i>=1; i--) {
00354 f[base::idx(i,j,Nx)](6) = f[base::idx(i,j+1,Nx)](6);
00355 f[base::idx(i,j,Nx)](7) = f[base::idx(i-1,j+1,Nx)](7);
00356 }
00357 for (register int j=0; j<=Ny-2; j++)
00358 for (register int i=0; i<=Nx-2; i++) {
00359 f[base::idx(i,j,Nx)](2) = f[base::idx(i+1,j,Nx)](2);
00360 f[base::idx(i,j,Nx)](8) = f[base::idx(i+1,j+1,Nx)](8);
00361 }
00362
00363
00364 int mxs = base::NGhosts(), mys = base::NGhosts();
00365 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00366 for (register int j=mys; j<=mye; j++)
00367 for (register int i=mxs; i<=mxe; i++)
00368 Collision(f[base::idx(i,j,Nx)],dt);
00369
00370 return DataType(1.);
00371 }
00372
00373 virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00374 DataType* aux=0, const int naux=0, const int scaling=0) const {
00375 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00376 int mxs = base::NGhosts(), mys = base::NGhosts();
00377 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00378 MicroType *f = (MicroType *)fvec.databuffer();
00379
00380 if (type==ConstantMicro) {
00381 MicroType g;
00382 for (register int i=0; i<base::NMicroVar(); i++)
00383 g(i) = aux[i];
00384
00385 for (register int j=mys; j<=mye; j++)
00386 for (register int i=mxs; i<=mxe; i++)
00387 f[base::idx(i,j,Nx)] = g;
00388 }
00389
00390 else if (type==ConstantMacro) {
00391 MacroType q;
00392 if (scaling==base::Physical) {
00393 q(0) = aux[0]/R0;
00394 q(1) = aux[1]/U0;
00395 q(2) = aux[2]/U0;
00396 }
00397 else
00398 for (register int i=0; i<base::NMacroVar(); i++)
00399 q(i) = aux[i];
00400 q(1) *= S0; q(2) *= S0;
00401
00402 for (register int j=mys; j<=mye; j++)
00403 for (register int i=mxs; i<=mxe; i++)
00404 f[base::idx(i,j,Nx)] = Equilibrium(q);
00405 }
00406
00407 else if (type==GasAtRest) {
00408 MacroType q;
00409 q(0) = GasDensity()/R0;
00410 q(1) = DataType(0.);
00411 q(2) = DataType(0.);
00412 for (register int j=mys; j<=mye; j++)
00413 for (register int i=mxs; i<=mxe; i++)
00414 f[base::idx(i,j,Nx)] = Equilibrium(q);
00415 }
00416 }
00417
00418
00419 virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00420 DataType* aux=0, const int naux=0, const int scaling=0) 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 if (type==Symmetry) {
00431 switch (side) {
00432 case base::Left:
00433 for (register int j=mys; j<=mye; j++) {
00434 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j,Nx)](8);
00435 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00436 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j,Nx)](5);
00437 }
00438 break;
00439 case base::Right:
00440 for (register int j=mys; j<=mye; j++) {
00441 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j,Nx)](7);
00442 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00443 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j,Nx)](4);
00444 }
00445 break;
00446 case base::Bottom:
00447 for (register int i=mxs; i<=mxe; i++) {
00448 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i,mye+1,Nx)](8);
00449 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00450 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i,mye+1,Nx)](7);
00451 }
00452 break;
00453 case base::Top:
00454 for (register int i=mxs; i<=mxe; i++) {
00455 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i,mys-1,Nx)](5);
00456 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00457 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i,mys-1,Nx)](4);
00458 }
00459 break;
00460 }
00461 }
00462
00463 else if (type==SlipWall) {
00464 switch (side) {
00465 case base::Left:
00466 for (register int j=mys; j<=mye; j++) {
00467 if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00468 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00469 if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00470 }
00471 break;
00472 case base::Right:
00473 for (register int j=mys; j<=mye; j++) {
00474 if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00475 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00476 if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00477 }
00478 break;
00479 case base::Bottom:
00480 for (register int i=mxs; i<=mxe; i++) {
00481 if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00482 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00483 if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00484 }
00485 break;
00486 case base::Top:
00487 for (register int i=mxs; i<=mxe; i++) {
00488 if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00489 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00490 if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00491 }
00492 break;
00493 }
00494 }
00495
00496 else if (type==NoSlipWall) {
00497 switch (side) {
00498 case base::Left:
00499 for (register int j=mys; j<=mye; j++) {
00500 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00501 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00502 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00503 }
00504 break;
00505 case base::Right:
00506 for (register int j=mys; j<=mye; j++) {
00507 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00508 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00509 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00510 }
00511 break;
00512 case base::Bottom:
00513 for (register int i=mxs; i<=mxe; i++) {
00514 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00515 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00516 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00517 }
00518 break;
00519 case base::Top:
00520 for (register int i=mxs; i<=mxe; i++) {
00521 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00522 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00523 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00524 }
00525 break;
00526 }
00527 }
00528
00529 else if (type==Inlet) {
00530 if (aux==0 || naux<=0) return;
00531 DataType a0=S0*aux[0], a1=0.;
00532 if (naux>1) a1 = S0*aux[1];
00533 if (scaling==base::Physical) {
00534 a0 /= U0;
00535 a1 /= U0;
00536 }
00537 switch (side) {
00538 case base::Left:
00539 for (register int j=mys; j<=mye; j++) {
00540 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00541 q(1) = a0;
00542 if (naux>1) q(2) = a1;
00543 f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00544 }
00545 break;
00546 case base::Right:
00547 for (register int j=mys; j<=mye; j++) {
00548 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00549 q(1) = a0;
00550 if (naux>1) q(2) = a1;
00551 f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00552 }
00553 break;
00554 case base::Bottom:
00555 for (register int i=mxs; i<=mxe; i++) {
00556 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00557 if (naux==1)
00558 q(2) = a0;
00559 else {
00560 q(1) = a0;
00561 q(2) = a1;
00562 }
00563 f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00564 }
00565 break;
00566 case base::Top:
00567 for (register int i=mxs; i<=mxe; i++) {
00568 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00569 if (naux==1)
00570 q(2) = a0;
00571 else {
00572 q(1) = a0;
00573 q(2) = a1;
00574 }
00575 f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00576 }
00577 break;
00578 }
00579 }
00580
00581 else if (type==Outlet) {
00582 switch (side) {
00583 case base::Left:
00584 for (register int j=mys; j<=mye; j++) {
00585 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00586 if (q(0)>0) f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)]/q(0);
00587 }
00588 break;
00589 case base::Right:
00590 for (register int j=mys; j<=mye; j++) {
00591 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00592 if (q(0)>0) f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)]/q(0);
00593 }
00594 break;
00595 case base::Bottom:
00596 for (register int i=mxs; i<=mxe; i++) {
00597 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00598 if (q(0)>0) f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)]/q(0);
00599 }
00600 break;
00601 case base::Top:
00602 for (register int i=mxs; i<=mxe; i++) {
00603 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00604 if (q(0)>0) f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)]/q(0);
00605 }
00606 break;
00607 }
00608 }
00609
00610
00611 else if (type==Pressure) {
00612 if (aux==0 || naux<=0) return;
00613 DataType a0=aux[0];
00614 if (scaling==base::Physical)
00615 a0 /= R0;
00616 switch (side) {
00617 case base::Left:
00618 for (register int j=mys; j<=mye; j++) {
00619 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00620 q(0) = a0;
00621 f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00622 }
00623 break;
00624 case base::Right:
00625 for (register int j=mys; j<=mye; j++) {
00626 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00627 q(0) = a0;
00628 f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00629 }
00630 break;
00631 case base::Bottom:
00632 for (register int i=mxs; i<=mxe; i++) {
00633 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00634 q(0) = a0;
00635 f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00636 }
00637 break;
00638 case base::Top:
00639 for (register int i=mxs; i<=mxe; i++) {
00640 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00641 q(0) = a0;
00642 f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00643 }
00644 break;
00645 }
00646 }
00647
00648 else if (type==SlidingWall) {
00649 if (aux==0 || naux<=0) return;
00650 if (aux==0 || naux<=0) return;
00651 DataType a0=S0*aux[0];
00652 if (scaling==base::Physical)
00653 a0 /= U0;
00654 switch (side) {
00655 case base::Left:
00656 for (register int j=mys; j<=mye; j++) {
00657 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00658 DataType rd1 = f[b+base::idx(mxe+1,j,Nx)](5), rd2 = f[b+base::idx(mxe+1,j,Nx)](8);
00659 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00660 DataType pw = DataType(1.)-qw;
00661 if (j<mye) f[b+base::idx(mxe,j+1,Nx)](7) = pw*rd1+qw*rd2;
00662 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00663 if (j>mys) f[b+base::idx(mxe,j-1,Nx)](4) = qw*rd1+pw*rd2;
00664 }
00665 break;
00666 case base::Right:
00667 for (register int j=mys; j<=mye; j++) {
00668 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00669 DataType rd1 = f[b+base::idx(mxs-1,j,Nx)](4), rd2 = f[b+base::idx(mxs-1,j,Nx)](7);
00670 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00671 DataType pw = DataType(1.)-qw;
00672 if (j<mye) f[b+base::idx(mxs,j+1,Nx)](8) = pw*rd1+qw*rd2;
00673 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00674 if (j>mys) f[b+base::idx(mxs,j-1,Nx)](5) = qw*rd1+pw*rd2;
00675 }
00676 break;
00677 case base::Bottom:
00678 for (register int i=mxs; i<=mxe; i++) {
00679 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00680 DataType rd1 = f[b+base::idx(i,mye+1,Nx)](7), rd2 = f[b+base::idx(i,mye+1,Nx)](8);
00681 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00682 DataType pw = DataType(1.)-qw;
00683 if (i<mxe) f[b+base::idx(i+1,mye,Nx)](5) = pw*rd1+qw*rd2;
00684 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00685 if (i>mxs) f[b+base::idx(i-1,mye,Nx)](4) = qw*rd1+pw*rd2;
00686 }
00687 break;
00688 case base::Top:
00689 for (register int i=mxs; i<=mxe; i++) {
00690 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00691 DataType rd1 = f[b+base::idx(i,mys-1,Nx)](4), rd2 = f[b+base::idx(i,mys-1,Nx)](5);
00692 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00693 DataType pw = DataType(1.)-qw;
00694 if (i<mxe) f[b+base::idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00695 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00696 if (i>mxs) f[b+base::idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00697 }
00698 break;
00699 }
00700 }
00701
00705 else if (type==CharacteristicOutlet) {
00706 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
00707 DCoords dx = base::GH().worldStep(fvec.stepsize());
00708 DataType dt_temp = dx(0)/VelocityScale();
00709 switch (side) {
00710 case base::Left:
00711 for (register int j=mys; j<=mye; j++) {
00712 MicroType ftmp = f[b+base::idx(mxe,j,Nx)];
00713 MacroType q = MacroVariables(ftmp);
00714 Collision(ftmp,dt_temp);
00715 MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00716 MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,Nx)]);
00717
00718 MacroType dqdn = NormalDerivative(q,q1,q2);
00719 Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00720
00721
00722 MacroType qn;
00723 qn(0) = q(0) - c1oCsSq2*L(0);
00724 qn(1) = q(1) + c1o2cs*L(0)/q(0);
00725 qn(2) = q(2) - L(1);
00726 f[b+base::idx(mxe,j,Nx)] = Equilibrium(qn);
00727 }
00728 break;
00729 case base::Right:
00730 for (register int j=mys; j<=mye; j++) {
00731 MicroType ftmp = f[b+base::idx(mxs,j,Nx)];
00732 MacroType q = MacroVariables(ftmp);
00733 Collision(ftmp,dt_temp);
00734 MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00735 MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,Nx)]);
00736
00737 MacroType dqdn = -NormalDerivative(q,q1,q2);
00738 if (method[0]!=3) q(0) = q(0);
00739 Vector<DataType,2> L = WaveAmplitudes(q(0),-q(1),dqdn(0),dqdn(1),dqdn(2));
00740
00741
00742 MacroType qn;
00743 qn(0) = q(0) - c1oCsSq2*L(0);
00744 qn(1) = q(1) - c1o2cs*L(0)/q(0);
00745 qn(2) = q(2) + L(1);
00746 f[b+base::idx(mxs,j,Nx)] = Equilibrium(qn);
00747 }
00748 break;
00749 case base::Bottom:
00750
00751 for (register int i=mxs; i<=mxe; i++) {
00752 MicroType ftmp = f[b+base::idx(i,mye,Nx)];
00753 MacroType q = MacroVariables(ftmp);
00754 Collision(ftmp,dt_temp);
00755 MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00756 MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,Nx)]);
00757
00758 MacroType dqdn = NormalDerivative(q,q1,q2);
00759 if (method[0]!=3) q(0) = q(0);
00760 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00761
00762
00763 MacroType qn;
00764 qn(0) = q(0) - c1oCsSq2*L(0);
00765 qn(1) = q(1) - L(1);
00766 qn(2) = q(2) + c1o2cs*L(0)/q(0);
00767 f[b+base::idx(i,mye,Nx)] = Equilibrium(qn);
00768 }
00769 break;
00770 case base::Top:
00771 for (register int i=mxs+base::NGhosts(); i<=mxe-base::NGhosts(); i++) {
00772 MicroType ftmp = f[b+base::idx(i,mys,Nx)];
00773 MacroType q = MacroVariables(ftmp);
00774 Collision(ftmp,dt_temp);
00775 MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00776 MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,Nx)]);
00777
00778 MacroType dqdn = -NormalDerivative(q,q1,q2);
00779 if (method[0]!=3) q(0) = q(0);
00780 Vector<DataType,2> L = WaveAmplitudes(q(0),-q(2),dqdn(0),dqdn(2),dqdn(1));
00781
00782
00783 MacroType qn;
00784 qn(0) = q(0) - c1oCsSq2*L(0);
00785 qn(1) = q(1) - L(1);
00786 qn(2) = q(2) + c1o2cs*L(0)/q(0);
00787 f[b+base::idx(i,mys,Nx)] = Equilibrium(qn);
00788 }
00789 break;
00790 }
00791 }
00792
00796 else if (type==CharacteristicInlet) {
00797 if (aux==0 || naux<=0) return;
00798 DataType rho=aux[0];
00799 DataType a0=S0*aux[1];
00800 DataType a1=S0*aux[2];
00801 if (scaling==base::Physical) {
00802 rho /= R0;
00803 a0 /= U0;
00804 a1 /= U0;
00805 }
00806 MacroType q, q1, q2, qn;
00807 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
00808 q(0) = rho;
00809 q(1) = a0;
00810 q(2) = a1;
00811 switch (side) {
00812 case base::Left:
00813 for (register int j=mys; j<=mye; j++) {
00814 q1 = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00815 q2 = MacroVariables(f[b+base::idx(mxe+2,j,Nx)]);
00816 if (naux==1) {
00817 q(1) = q1(1); q(2) = q1(2);
00818 MacroType dqdn = NormalDerivative(q,q1,q2);
00819 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00820 qn(0) = q(0) - c1oCsSq2*L(0);
00821 qn(1) = q1(1) + c1o2cs*L(0)/q(0);
00822 qn(2) = q1(2) - L(1);
00823 }
00824 else if (naux==2) {
00825 q(0) = q1(0); q(2) = q1(2);
00826 MacroType dqdn = NormalDerivative(q,q1,q2);
00827 Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00828 qn(0) = q1(0) - c1oCsSq2*L(0);
00829 qn(1) = q(1) + c1o2cs*L(0)/q(0);
00830 qn(2) = q1(2) - L(1);
00831 }
00832 else if (naux==3) {
00833 q(0) = q1(0);
00834 MacroType dqdn = NormalDerivative(q,q1,q2);
00835 Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00836 qn(0) = q1(0) - c1oCsSq2*L(0);
00837 qn(1) = q(1) + c1o2cs*L(0)/q(0);
00838 qn(2) = q(2) - L(1);
00839 }
00840 f[b+base::idx(mxe,j,Nx)] = Equilibrium(qn);
00841 }
00842 break;
00843 case base::Right:
00844 for (register int j=mys; j<=mye; j++) {
00845 q1 = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00846 q2 = MacroVariables(f[b+base::idx(mxs-2,j,Nx)]);
00847 if (naux==1) {
00848 q(1) = q1(1); q(2) = q1(2);
00849 MacroType dqdn = NormalDerivative(q,q1,q2);
00850 Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00851 qn(0) = q(0) - c1oCsSq2*L(0);
00852 qn(1) = q1(1) + c1o2cs*L(0)/q(0);
00853 qn(2) = q1(2) - L(1);
00854 }
00855 else if (naux==2) {
00856 q(0) = q1(0); q(2) = q1(2);
00857 MacroType dqdn = NormalDerivative(q,q1,q2);
00858 Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00859 qn(0) = q1(0) - c1oCsSq2*L(0);
00860 qn(1) = q(1) + c1o2cs*L(0)/q(0);
00861 qn(2) = q1(2) - L(1);
00862 }
00863 else if (naux==3) {
00864 q(0) = q1(0);
00865 MacroType dqdn = NormalDerivative(q,q1,q2);
00866 Vector<DataType,2> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2));
00867 qn(0) = q1(0) - c1oCsSq2*L(0);
00868 qn(1) = q(1) + c1o2cs*L(0)/q(0);
00869 qn(2) = q(2) - L(1);
00870 }
00871 f[b+base::idx(mxs,j,Nx)] = Equilibrium(qn);
00872 }
00873 break;
00874 case base::Bottom:
00875 for (register int i=mxs; i<=mxe; i++) {
00876 q1 = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00877 q2 = MacroVariables(f[b+base::idx(i,mye+2,Nx)]);
00878 if (naux==1) {
00879 q(1) = q1(1); q(2) = q1(2);
00880 MacroType dqdn = NormalDerivative(q,q1,q2);
00881 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00882 qn(0) = q(0) - c1oCsSq2*L(0);
00883 qn(1) = q1(1) - L(1);
00884 qn(2) = q1(2) + c1o2cs*L(0)/q(0);
00885 }
00886 else if (naux==2) {
00887 q(0) = q1(0); q(1) = q1(1);
00888 MacroType dqdn = NormalDerivative(q,q1,q2);
00889 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00890 qn(0) = q1(0) - c1oCsSq2*L(0);
00891 qn(1) = q1(1) - L(1);
00892 qn(2) = q(2) + c1o2cs*L(0)/q(0);
00893 }
00894 else if (naux==3) {
00895 q(0) = q1(0);
00896 MacroType dqdn = NormalDerivative(q,q1,q2);
00897 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00898 qn(0) = q1(0) - c1oCsSq2*L(0);
00899 qn(1) = q(1) - L(1);
00900 qn(2) = q(2) + c1o2cs*L(0)/q(0);
00901 }
00902 f[b+base::idx(i,mye,Nx)] = Equilibrium(qn);
00903 }
00904 break;
00905 case base::Top:
00906 for (register int i=mxs; i<=mxe; i++) {
00907 q1 = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00908 q2 = MacroVariables(f[b+base::idx(i,mys-2,Nx)]);
00909 if (naux==1) {
00910 q(1) = q1(1); q(2) = q1(2);
00911 MacroType dqdn = NormalDerivative(q,q1,q2);
00912 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00913 qn(0) = q(0) - c1oCsSq2*L(0);
00914 qn(1) = q1(1) - L(1);
00915 qn(2) = q1(2) + c1o2cs*L(0)/q(0);
00916 }
00917 else if (naux==2) {
00918 q(0) = q1(0); q(1) = q1(1);
00919 MacroType dqdn = NormalDerivative(q,q1,q2);
00920 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00921 qn(0) = q1(0) - c1oCsSq2*L(0);
00922 qn(1) = q1(1) - L(1);
00923 qn(2) = q(2) + c1o2cs*L(0)/q(0);
00924 }
00925 else if (naux==3) {
00926 q(0) = q1(0);
00927 MacroType dqdn = NormalDerivative(q,q1,q2);
00928 Vector<DataType,2> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1));
00929 qn(0) = q1(0) - c1oCsSq2*L(0);
00930 qn(1) = q(1) - L(1);
00931 qn(2) = q(2) + c1o2cs*L(0)/q(0);
00932 }
00933 f[b+base::idx(i,mys,Nx)] = Equilibrium(qn);
00934 }
00935 break;
00936 }
00937 }
00938
00943 else if (type==NoSlipWallEntranceExit) {
00944 DCoords wc;
00945 int lc_indx[2] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1) };
00946 DataType WLB[2], WUB[2];
00947 base::GH().wlb(WLB);
00948 base::GH().wub(WUB);
00949 DataType slip = 0.;
00950 DataType lxs, lxe, lys, lye;
00951 DataType min[2], max[2];
00952 if (naux!=4) assert(0);
00953 lxs = aux[0] - WLB[0]; lxe = WUB[0] - aux[1];
00954 lys = aux[2] - WLB[1]; lye = WUB[1] - aux[3];
00955 min[0] = aux[0]; max[0] = aux[1];
00956 min[1] = aux[2]; max[1] = aux[3];
00957
00958 switch (side) {
00959 case base::Left:
00960 for (register int j=mys; j<=mye; j++) {
00961 lc_indx[1] = j*fvec.stepsize(1);
00962 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
00963 if (wc(1)>=min[1] && wc(1)<=max[1]) {
00964 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00965 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00966 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00967 }
00968 else {
00969 if (wc(1)<WLB[1] || wc(1) >WUB[1])
00970 slip = DataType(0.);
00971 else if (wc(1)<min[1])
00972 slip = (wc(1) - WLB[1])/lys;
00973 else if (wc(1)>max[1])
00974 slip = (WUB[1] - wc(1))/lye;
00975 if (j>mys && j<mye) {
00976 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
00977 + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
00978 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
00979 + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
00980 }
00981 else if (j==mys) {
00982 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](5)
00983 + slip*f[b+base::idx(mxe+1,j+1,Nx)](8);
00984 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,Nx)](8)
00985 + slip*f[b+base::idx(mxe+1,j,Nx)](5);
00986 }
00987 else if (j==mye) {
00988 f[b+base::idx(mxe,j,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,Nx)](5)
00989 + slip*f[b+base::idx(mxe+1,j,Nx)](8);
00990 f[b+base::idx(mxe,j,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,Nx)](8)
00991 + slip*f[b+base::idx(mxe+1,j-1,Nx)](5);
00992 }
00993 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00994 }
00995 }
00996 break;
00997 case base::Right:
00998 for (register int j=mys; j<=mye; j++) {
00999 lc_indx[1] = j*fvec.stepsize(1);
01000 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01001 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01002 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
01003 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
01004 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
01005 }
01006 else {
01007 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01008 slip = DataType(0.);
01009 else if (wc(1)<min[1])
01010 slip = (wc(1) - WLB[1])/lys;
01011 else if (wc(1)>max[1])
01012 slip = (WUB[1] - wc(1))/lye;
01013 if (j>mys && j<mye) {
01014 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
01015 + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
01016 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
01017 + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
01018 }
01019 else if (j==mys) {
01020 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](4)
01021 + slip*f[b+base::idx(mxs-1,j+1,Nx)](7);
01022 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,Nx)](7)
01023 + slip*f[b+base::idx(mxs-1,j,Nx)](4);
01024 }
01025 else if (j==mye) {
01026 f[b+base::idx(mxs,j,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,Nx)](4)
01027 + slip*f[b+base::idx(mxs-1,j,Nx)](7);
01028 f[b+base::idx(mxs,j,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,Nx)](7)
01029 + slip*f[b+base::idx(mxs-1,j-1,Nx)](4);
01030 }
01031 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
01032 }
01033 }
01034 break;
01035 case base::Bottom:
01036 for (register int i=mxs; i<=mxe; i++) {
01037 lc_indx[0] = i*fvec.stepsize(0);
01038 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01039 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01040 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
01041 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
01042 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
01043 }
01044 else {
01045 if (wc(0)<WLB[0] || wc(0)>WUB[0])
01046 slip = DataType(0.);
01047 else if (wc(0)<min[0])
01048 slip = (wc(0) - WLB[0])/lxs;
01049 else if (wc(0)>max[0])
01050 slip = (WUB[0] - wc(0))/lxe;
01051 if (i>mxs && i<mxe) {
01052 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,Nx)](7)
01053 + slip*f[b+base::idx(i+1,mye+1,Nx)](8);
01054 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,Nx)](8)
01055 + slip*f[b+base::idx(i-1,mye+1,Nx)](7);
01056 }
01057 else if (i==mxs) {
01058 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](7)
01059 + slip*f[b+base::idx(i+1,mye+1,Nx)](8);
01060 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,Nx)](8)
01061 + slip*f[b+base::idx(i,mye+1,Nx)](7);
01062 }
01063 else if (i==mxe) {
01064 f[b+base::idx(i,mye,Nx)](4) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](7)
01065 + slip*f[b+base::idx(i,mye+1,Nx)](8);
01066 f[b+base::idx(i,mye,Nx)](5) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,Nx)](8)
01067 + slip*f[b+base::idx(i-1,mye+1,Nx)](7);
01068 }
01069 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
01070 }
01071 }
01072 break;
01073 case base::Top:
01074 for (register int i=mxs; i<=mxe; i++) {
01075 lc_indx[0] = i*fvec.stepsize(0);
01076 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01077 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01078 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
01079 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
01080 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
01081 }
01082 else {
01083 if (wc(0)<WLB[0] || wc(0)>WUB[0])
01084 slip = DataType(0.);
01085 else if (wc(0)<min[0])
01086 slip = (wc(0) - WLB[0])/lxs;
01087 else if (wc(0)>max[0])
01088 slip = (WUB[0] - wc(0))/lxe;
01089 if (i>mxs && i<mxe) {
01090 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,Nx)](4)
01091 + slip*f[b+base::idx(i+1,mys-1,Nx)](5);
01092 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,Nx)](5)
01093 + slip*f[b+base::idx(i-1,mys-1,Nx)](4);
01094 }
01095 else if (i==mxs) {
01096 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,Nx)](4)
01097 + slip*f[b+base::idx(i+1,mys-1,Nx)](5);
01098 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,Nx)](5)
01099 + slip*f[b+base::idx(i,mys-1,Nx)](4);
01100 }
01101 else if (i==mxe) {
01102 f[b+base::idx(i,mys,Nx)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,Nx)](4)
01103 + slip*f[b+base::idx(i,mys-1,Nx)](5);
01104 f[b+base::idx(i,mys,Nx)](7) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,Nx)](5)
01105 + slip*f[b+base::idx(i-1,mys-1,Nx)](4);
01106 }
01107 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
01108 }
01109 }
01110 break;
01111 }
01112 }
01113 }
01114
01115 virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
01116 const MicroType* f, const point_type* xc, const DataType* distance,
01117 const point_type* normal, DataType* aux=0, const int naux=0,
01118 const int scaling=0) const {
01119
01120 if (type==GFMNoSlipWall || type==GFMSlipWall) {
01121 DataType fac = S0;
01122 if (scaling==base::Physical)
01123 fac /= U0;
01124 for (int n=0; n<nc; n++) {
01125 MacroType q=MacroVariables(f[n]);
01126 DataType u=-q(1);
01127 DataType v=-q(2);
01128
01129
01130 if (naux>=2) {
01131 u += fac*aux[naux*n];
01132 v += fac*aux[naux*n+1];
01133 }
01134 if (type==GFMNoSlipWall) {
01135
01136 q(1) += DataType(2.)*u;
01137 q(2) += DataType(2.)*v;
01138 }
01139 else if (type==GFMSlipWall) {
01140
01141 DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
01142 q(1) += vl*normal[n](0);
01143 q(2) += vl*normal[n](1);
01144 }
01145 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01146 }
01147 }
01148
01149 else if (type==GFMExtrapolation) {
01150 for (int n=0; n<nc; n++) {
01151 MacroType q=MacroVariables(f[n]);
01152 DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1);
01153 q(1) = vl*normal[n](0);
01154 q(2) = vl*normal[n](1);
01155 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01156 }
01157 }
01158 }
01159
01160 inline virtual Vector<DataType,2> WaveAmplitudes(const DataType rho, const DataType vn, const DataType drhodn,
01161 const DataType dvndn, const DataType dvtdn) const {
01162 Vector<DataType,2> L;
01163 L(0) = (vn - LatticeSpeedOfSound())*(cs2*drhodn - rho*LatticeSpeedOfSound()*dvndn);
01164 L(1) = vn*dvtdn;
01165 return L;
01166 }
01167 inline virtual const MacroType NormalDerivative(const MacroType qa, const MacroType qb, const MacroType qc) const
01168 { return DataType(-1.5)*qa + DataType(2.)*qb - DataType(0.5)*qc; }
01169
01170 virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01171 const int skip_ghosts=1) const {
01172 int Nx = fvec.extents(0), Ny = fvec.extents(1);
01173 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01174 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01175 MicroType *f = (MicroType *)fvec.databuffer();
01176 DataType *work = (DataType *)workvec.databuffer();
01177 DCoords dx = base::GH().worldStep(fvec.stepsize());
01178 DataType dt = dx(0)*S0/U0;
01179
01180 if ((cnt>=1 && cnt<=9) || (cnt>=16 && cnt<=18) || (cnt>=23 && cnt<=25)) {
01181 for (register int j=mys; j<=mye; j++)
01182 for (register int i=mxs; i<=mxe; i++) {
01183 MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
01184 switch(cnt) {
01185 case 1:
01186 if (method[0]==0) work[base::idx(i,j,Nx)]=q(0)*R0+rhop;
01187 else
01188 work[base::idx(i,j,Nx)]=q(0)*R0;
01189 break;
01190 case 2:
01191 work[base::idx(i,j,Nx)]=q(1)*U0/S0;
01192 break;
01193 case 3:
01194 work[base::idx(i,j,Nx)]=q(2)*U0/S0;
01195 break;
01196 case 5:
01197 work[base::idx(i,j,Nx)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure());
01198 break;
01199 case 6:
01200 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure();
01201 break;
01202 case 8: {
01203 MicroType feq = Equilibrium(q);
01204 work[base::idx(i,j,Nx)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt),
01205 GasSpeedofSound(),dt);
01206 break;
01207 }
01208 case 9:
01209 work[base::idx(i,j,Nx)]=std::sqrt(q(1)*q(1)+q(2)*q(2))*U0/S0;
01210 break;
01211 case 16:
01212 case 17:
01213 case 18: {
01214 MicroType feq = Equilibrium(q);
01215 DataType omega;
01216 if (turbulence == laminar)
01217 omega = Omega(dt);
01218 else if (turbulence == LES_Smagorinsky)
01219 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt);
01220 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*DeviatoricStress(f[base::idx(i,j,Nx)],feq,omega)(cnt-16);
01221 break;
01222 }
01223 case 23:
01224 case 24:
01225 case 25: {
01226 DataType omega;
01227 if (turbulence == laminar)
01228 omega = Omega(dt);
01229 else if (turbulence == LES_Smagorinsky) {
01230 MicroType feq = Equilibrium(q);
01231 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt);
01232 }
01233 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*Stress(f[base::idx(i,j,Nx)],q,omega)(cnt-23);
01234 if (cnt==23 || cnt==25) work[base::idx(i,j,Nx)]-=BasePressure();
01235 break;
01236 }
01237 }
01238 }
01239 }
01240 else if ((cnt>=10 && cnt<=11) || cnt==15 || (cnt>=30 && cnt<=32)) {
01241 macro_grid_data_type qgrid(fvec.bbox());
01242 MacroType *q = (MacroType *)qgrid.databuffer();
01243 for (register int j=mys; j<=mye; j++)
01244 for (register int i=mxs; i<=mxe; i++) {
01245 q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
01246 q[base::idx(i,j,Nx)](0)*=R0;
01247 q[base::idx(i,j,Nx)](1)*=U0/S0;
01248 q[base::idx(i,j,Nx)](2)*=U0/S0;
01249 }
01250
01251 if (skip_ghosts==0) {
01252 switch(cnt) {
01253 case 30:
01254 mxs+=1; mxe-=1;
01255 break;
01256 case 32:
01257 mys+=1; mye-=1;
01258 break;
01259 case 10:
01260 case 11:
01261 case 15:
01262 case 31:
01263 mxs+=1; mxe-=1; mys+=1; mye-=1;
01264 break;
01265 }
01266 }
01267
01268 DataType dxux, dxuy, dyux, dyuy;
01269 for (register int j=mys; j<=mye; j++)
01270 for (register int i=mxs; i<=mxe; i++) {
01271 if (cnt==15 || cnt==30)
01272 dxux=(q[base::idx(i+1,j,Nx)](1)-q[base::idx(i-1,j,Nx)](1))/(2.*dx(0));
01273 if (cnt==15 || cnt==32)
01274 dyuy=(q[base::idx(i,j+1,Nx)](2)-q[base::idx(i,j-1,Nx)](2))/(2.*dx(1));
01275 if (cnt==10 || cnt==11 || cnt==15 || cnt==31) {
01276 dxuy=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(2.*dx(0));
01277 dyux=(q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(2.*dx(1));
01278 }
01279
01280 switch (cnt) {
01281 case 10:
01282 work[base::idx(i,j,Nx)]=dxuy-dyux;
01283 break;
01284 case 11:
01285 work[base::idx(i,j,Nx)]=std::abs(dxuy-dyux);
01286 break;
01287 case 15:
01288 work[base::idx(i,j,Nx)]=DataType(-0.5)*(4.*dxux*dxux+8.*dyux*dxuy+4.*dyuy*dyuy);
01289 break;
01290 case 30:
01291 work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dxux;
01292 break;
01293 case 31:
01294 work[base::idx(i,j,Nx)]=q[base::idx(i,j,Nx)](0)*GasViscosity()*(dxuy+dyux);
01295 break;
01296 case 32:
01297 work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dyuy;
01298 break;
01299 }
01300 }
01301 }
01302 }
01303
01304 virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01305 const int skip_ghosts=1) const {
01306 int Nx = fvec.extents(0), Ny = fvec.extents(1);
01307 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01308 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01309 MicroType *f = (MicroType *)fvec.databuffer();
01310 DataType *work = (DataType *)workvec.databuffer();
01311
01312 for (register int j=mys; j<=mye; j++)
01313 for (register int i=mxs; i<=mxe; i++) {
01314 switch(cnt) {
01315 case 1:
01316 f[base::idx(i,j,Nx)](0)=work[base::idx(i,j,Nx)]/R0;
01317 break;
01318 case 2:
01319 f[base::idx(i,j,Nx)](1)=work[base::idx(i,j,Nx)]/(U0/S0);
01320 break;
01321 case 3:
01322 f[base::idx(i,j,Nx)](2)=work[base::idx(i,j,Nx)]/(U0/S0);
01323 f[base::idx(i,j,Nx)] = Equilibrium((const MacroType &)f[base::idx(i,j,Nx)]);
01324 break;
01325 }
01326 }
01327 }
01328
01329 virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
01330 const int verbose) const {
01331 int Nx = fvec.extents(0);
01332 int b = fvec.bottom();
01333 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01334 BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
01335 int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01336 int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01337 int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01338 int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01339 MicroType *f = (MicroType *)fvec.databuffer();
01340
01341 int result = 1;
01342 for (register int j=mys; j<=mye; j++)
01343 for (register int i=mxs; i<=mxe; i++)
01344 for (register int l=0; l<base::NMicroVar(); l++)
01345 if (!(f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l)>-1.e37)) {
01346 result = 0;
01347 if (verbose)
01348 std::cerr << "D2Q9-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << ")(" << l << ")="
01349 << f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l) << " " << fvec.bbox() << std::endl;
01350 }
01351 return result;
01352 }
01353
01354 inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
01355 TensorType P;
01356 if (method[1]==0) {
01357 P(0) = q(0)*q(1)*q(1)-(f(1)+f(2)+f(4)+f(5)+f(7)+f(8))+cs2*q(0);
01358 P(1) = q(0)*q(1)*q(2)-(f(4)-f(5)-f(7)+f(8));
01359 P(2) = q(0)*q(2)*q(2)-(f(3)+f(4)+f(5)+f(6)+f(7)+f(8))+cs2*q(0);
01360 P *= -(DataType(1.0)-DataType(0.5)*om);
01361 }
01362 else {
01363 MicroType feq = Equilibrium(q);
01364 P = DeviatoricStress(f,feq,om);
01365 }
01366 P(0) -= LatticeBasePressure(q(0));
01367 P(2) -= LatticeBasePressure(q(0));
01368 return P;
01369 }
01370
01371 inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
01372 TensorType Sigma;
01373 if (method[2]==0) {
01374 MicroType fneq = feq - f;
01375 Sigma(0) = fneq(1)+fneq(2)+fneq(4)+fneq(5)+fneq(7)+fneq(8);
01376 Sigma(1) = fneq(4)-fneq(5)-fneq(7)+fneq(8);
01377 Sigma(2) = fneq(3)+fneq(4)+fneq(5)+fneq(6)+fneq(7)+fneq(8);
01378 Sigma *= -(DataType(1.0)-DataType(0.5)*om);
01379 }
01380 else {
01381 MacroType q = MacroVariables(f);
01382 Sigma = Stress(f,q,om);
01383 Sigma(0) += LatticeBasePressure(q(0));
01384 Sigma(2) += LatticeBasePressure(q(0));
01385 }
01386 return Sigma;
01387 }
01388
01389 virtual int NMethodOrder() const { return 2; }
01390
01391 inline const DataType& L0() const { return base::LengthScale(); }
01392 inline const DataType& T0() const { return base::TimeScale(); }
01393 inline void SetDensityScale(const DataType r0) { R0 = r0; }
01394 inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
01395 inline void SetSpeedUp(const DataType s0) { S0 = s0; }
01396 inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
01397
01398 inline const DataType& DensityScale() const { return R0; }
01399 inline const DataType VelocityScale() const { return U0/S0; }
01400
01401 inline const DataType& SpeedUp() const { return S0; }
01402
01403
01404 inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
01405 inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
01406 inline virtual DataType LatticeBasePressure(const DataType rho) const {
01407 return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
01408 }
01409
01410
01411 inline void SetGas(DataType rho, DataType nu, DataType cs) {
01412 rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
01413 SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
01414 }
01415 inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
01416 inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q,
01417 const DataType dt) const {
01418 DataType M2 = 0.0;
01419
01420 M2 += (f(4) - feq(4))*(f(4) - feq(4));
01421 M2 += (f(5) - feq(5))*(f(5) - feq(5));
01422 M2 += (f(7) - feq(7))*(f(7) - feq(7));
01423 M2 += (f(8) - feq(8))*(f(8) - feq(8));
01424 M2 = std::sqrt(M2);
01425
01426 DataType tau = DataType(1.0)/Omega(dt);
01427 DataType om_turb = (DataType(2.0)/( tau
01428 + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
01429
01430 return ( om_turb );
01431 }
01432 inline const int TurbulenceType() const { return turbulence; }
01433 inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
01434 inline const DataType& GasDensity() const { return rhop; }
01435 inline const DataType& GasViscosity() const { return nup; }
01436 inline const DataType& GasSpeedofSound() const { return csp; }
01437 inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
01438 { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
01439
01440 inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
01441 inline virtual DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); }
01442 inline virtual DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); }
01443
01444 protected:
01445 DataType cs2, cs22, cssq;
01446 DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
01447 DataType Cs_Smagorinsky, turbulence;
01448 int method[4], mdx[9], mdy[9];
01449 };
01450
01451
01452 #endif