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 enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall, GFMBounceBack };
00051 enum TurbulenceModel { laminar, LES_Smagorinsky };
00052
00053 LBMD2Q9() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.) {
00054 cs2 = DataType(1.)/DataType(3.);
00055 cs22 = DataType(2.)*cs2;
00056 cssq = DataType(2.)/DataType(9.);
00057 method[0] = 2; method[1] = 0; method[2] = 0;
00058 turbulence = laminar;
00059 mdx[0]=mdx[3]=mdx[6]=0; mdx[1]=mdx[4]=mdx[7]=1; mdx[2]=mdx[5]=mdx[8]=-1;
00060 mdy[0]=mdy[1]=mdy[2]=0; mdy[3]=mdy[4]=mdy[5]=1; mdy[6]=mdy[7]=mdy[8]=-1;
00061 }
00062
00063 virtual ~LBMD2Q9() {}
00064
00065 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00066 base::LocCtrl = Ctrl.getSubDevice(prefix+"D2Q9");
00067 RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00068 RegisterAt(base::LocCtrl,"Method(2)",method[1]);
00069 RegisterAt(base::LocCtrl,"Method(3)",method[2]);
00070 RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00071 RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00072 RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00073 RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00074 RegisterAt(base::LocCtrl,"Gas_nu",nup);
00075 RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00076 RegisterAt(base::LocCtrl,"Gas_W",Wp);
00077 RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00078 RegisterAt(base::LocCtrl,"SpeedUp",S0);
00079 }
00080
00081 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00082 base::SetupData(gh,ghosts);
00083 if (rhop>0.) R0 = rhop;
00084 if (csp>0.) {
00085 cs2p = csp*csp;
00086 SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00087 }
00088 std::cout << "D2Q9: Gas_rho=" << rhop << " Gas_Cs=" << csp
00089 << " Gas_nu=" << nup << std::endl;
00090 std::cout << "D2Q9: L0=" << L0() << " T0=" << T0() << " U0=" << U0
00091 << " R0=" << R0 << std::endl;
00092 WriteInit();
00093 }
00094
00095 virtual void WriteInit() const {
00096 int me = MY_PROC;
00097 if (me == VizServer) {
00098 std::ofstream ofs("chem.dat", std::ios::out);
00099 ofs << "RU " << Rp << std::endl;
00100 ofs << "PA " << BasePressure() << std::endl;
00101 ofs << "Sp(01) Vicosity" << std::endl;
00102 ofs << "W(01) " << Wp << std::endl;
00103 ofs << "Sp(02) |Velocity|" << std::endl;
00104 ofs << "W(02) " << 1.0 << std::endl;
00105 ofs << "Sp(03) Vorticity" << std::endl;
00106 ofs << "W(03) " << 1.0 << std::endl;
00107 ofs << "Sp(04) |Vorticity|" << std::endl;
00108 ofs << "W(04) " << 1.0 << std::endl;
00109 ofs.close();
00110 }
00111 }
00112
00113 inline virtual MacroType MacroVariables(const MicroType &f) const {
00114 MacroType q;
00115 q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8);
00116 q(1)=(f(1)-f(2)+f(4)-f(5)+f(7)-f(8))/q(0);
00117 q(2)=(f(3)+f(4)+f(5)-f(6)-f(7)-f(8))/q(0);
00118 return q;
00119 }
00120
00121 inline virtual MicroType Equilibrium(const MacroType &q) const {
00122 MicroType feq;
00123
00124 DataType rho = q(0);
00125 DataType u = q(1);
00126 DataType v = q(2);
00127
00128 DataType ri, rt0, rt1, rt2;
00129 if (method[0]==0) {
00130 ri = DataType(1.);
00131 rt0 = R0 * DataType(4.) / DataType(9.);
00132 rt1 = R0 / DataType(9.);
00133 rt2 = R0 / DataType(36.);
00134 }
00135 if (method[0]==1) {
00136 ri = rho;
00137 rt0 = DataType(4.) / DataType(9.);
00138 rt1 = DataType(1.) / DataType(9.);
00139 rt2 = DataType(1.) / DataType(36.);
00140 }
00141 else if (method[0]==2) {
00142 ri = DataType(1.);
00143 rt0 = rho * DataType(4.) / DataType(9.);
00144 rt1 = rho / DataType(9.);
00145 rt2 = rho / DataType(36.);
00146 }
00147
00148 DataType usq = u * u;
00149 DataType vsq = v * v;
00150 DataType sumsq = (usq + vsq) / cs22;
00151 DataType sumsq2 = sumsq * (DataType(1.) - cs2) / cs2;
00152 DataType u2 = usq / cssq;
00153 DataType v2 = vsq / cssq;
00154 DataType ui = u / cs2;
00155 DataType vi = v / cs2;
00156 DataType uv = ui * vi;
00157
00158 feq(0) = rt0 * (ri - sumsq);
00159
00160 feq(1) = rt1 * (ri - sumsq + u2 + ui);
00161 feq(2) = rt1 * (ri - sumsq + u2 - ui);
00162 feq(3) = rt1 * (ri - sumsq + v2 + vi);
00163 feq(6) = rt1 * (ri - sumsq + v2 - vi);
00164
00165 feq(4) = rt2 * (ri + sumsq2 + ui + vi + uv);
00166 feq(5) = rt2 * (ri + sumsq2 - ui + vi - uv);
00167 feq(7) = rt2 * (ri + sumsq2 + ui - vi - uv);
00168 feq(8) = rt2 * (ri + sumsq2 - ui - vi + uv);
00169
00170 return feq;
00171 }
00172
00173 inline virtual void Collision(MicroType &f, const DataType dt) const {
00174 MacroType q = MacroVariables(f);
00175 MicroType feq = Equilibrium(q);
00176
00177 DataType omega;
00178 if (turbulence == laminar)
00179 omega = Omega(dt);
00180 else if (turbulence == LES_Smagorinsky)
00181 omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00182
00183 f=(DataType(1.)-omega)*f + omega*feq;
00184 }
00185
00186 virtual int IncomingIndices(const int side, int indices[]) const {
00187 switch (side) {
00188 case base::Left:
00189 indices[0] = 1; indices[1] = 4; indices[2] = 7;
00190 break;
00191 case base::Right:
00192 indices[0] = 2; indices[1] = 5; indices[2] = 8;
00193 break;
00194 case base::Bottom:
00195 indices[0] = 3; indices[1] = 4; indices[2] = 5;
00196 break;
00197 case base::Top:
00198 indices[0] = 6; indices[1] = 7; indices[2] = 8;
00199 break;
00200 }
00201 return 3;
00202 }
00203
00204 virtual int OutgoingIndices(const int side, int indices[]) const {
00205 switch (side) {
00206 case base::Left:
00207 indices[0] = 2; indices[1] = 5; indices[2] = 8;
00208 break;
00209 case base::Right:
00210 indices[0] = 1; indices[1] = 4; indices[2] = 7;
00211 break;
00212 case base::Bottom:
00213 indices[0] = 6; indices[1] = 7; indices[2] = 8;
00214 break;
00215 case base::Top:
00216 indices[0] = 3; indices[1] = 4; indices[2] = 5;
00217 break;
00218 }
00219 return 3;
00220 }
00221
00222
00223 virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb,
00224 const int side) const {
00225 int Nx = fvec.extents(0);
00226 int b = fvec.bottom();
00227 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00228 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00229 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00230 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00231 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00232 MicroType *f = (MicroType *)fvec.databuffer();
00233
00234 #ifdef DEBUG_PRINT_FIXUP
00235 ( comm_service::log() << "Reverse streaming at Side: " << side << " of "
00236 << fvec.bbox() << " using " << bb << "\n").flush();
00237 #endif
00238
00239 switch (side) {
00240 case base::Left:
00241 for (register int j=mys; j<=mye; j++) {
00242 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](2);
00243 f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](8);
00244 f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](5);
00245 }
00246 break;
00247 case base::Right:
00248 for (register int j=mys; j<=mye; j++) {
00249 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](1);
00250 f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](7);
00251 f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](4);
00252 }
00253 break;
00254 case base::Bottom:
00255 for (register int i=mxs; i<=mxe; i++) {
00256 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](6);
00257 f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](8);
00258 f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](7);
00259 }
00260 break;
00261 case base::Top:
00262 for (register int i=mxs; i<=mxe; i++) {
00263 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](3);
00264 f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](5);
00265 f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](4);
00266 }
00267 break;
00268 }
00269 }
00270
00271 virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec,
00272 const BBox &bb, const double& dt) const {
00273 int Nx = fvec.extents(0);
00274 int b = fvec.bottom();
00275 MicroType *f = (MicroType *)fvec.databuffer();
00276 MicroType *of = (MicroType *)ovec.databuffer();
00277
00278 #ifdef DEBUG_PRINT_FIXUP
00279 ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox()
00280 << " using " << bb << "\n").flush();
00281 #endif
00282
00283 assert (fvec.extents(0)==ovec.extents(0));
00284 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) &&
00285 fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1));
00286 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00287 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00288 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00289 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00290
00291 for (register int j=mys; j<=mye; j++)
00292 for (register int i=mxs; i<=mxe; i++) {
00293 for (register int n=0; n<base::NMicroVar(); n++)
00294 f[b+base::idx(i,j,Nx)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],Nx)](n);
00295 Collision(f[b+base::idx(i,j,Nx)],dt);
00296 }
00297 }
00298
00299 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00300 const double& t, const double& dt, const int& mpass) const {
00301
00302
00303 DCoords dx = base::GH().worldStep(fvec.stepsize());
00304
00305 if (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR) {
00306 std::cerr << "LBMD2Q9 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00307 << " to be equal." << std::endl << std::flush;
00308 std::exit(-1);
00309 }
00310 if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR) {
00311 std::cerr << "LBMD2Q9 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00312 << " to be equal."
00313 << " dt=" << dt << " TO=" << T0()
00314 << std::endl << std::flush;
00315 std::exit(-1);
00316 }
00317
00318 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00319 MicroType *f = (MicroType *)fvec.databuffer();
00320
00321
00322 for (register int j=Ny-1; j>=1; j--)
00323 for (register int i=0; i<=Nx-2; i++) {
00324 f[base::idx(i,j,Nx)](3) = f[base::idx(i,j-1,Nx)](3);
00325 f[base::idx(i,j,Nx)](5) = f[base::idx(i+1,j-1,Nx)](5);
00326 }
00327
00328 for (register int j=Ny-1; j>=1; j--)
00329 for (register int i=Nx-1; i>=1; i--) {
00330 f[base::idx(i,j,Nx)](1) = f[base::idx(i-1,j,Nx)](1);
00331 f[base::idx(i,j,Nx)](4) = f[base::idx(i-1,j-1,Nx)](4);
00332 }
00333
00334 for (register int j=0; j<=Ny-2; j++)
00335 for (register int i=Nx-1; i>=1; i--) {
00336 f[base::idx(i,j,Nx)](6) = f[base::idx(i,j+1,Nx)](6);
00337 f[base::idx(i,j,Nx)](7) = f[base::idx(i-1,j+1,Nx)](7);
00338 }
00339
00340 for (register int j=0; j<=Ny-2; j++)
00341 for (register int i=0; i<=Nx-2; i++) {
00342 f[base::idx(i,j,Nx)](2) = f[base::idx(i+1,j,Nx)](2);
00343 f[base::idx(i,j,Nx)](8) = f[base::idx(i+1,j+1,Nx)](8);
00344 }
00345
00346
00347 int mxs = base::NGhosts(), mys = base::NGhosts();
00348 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00349 for (register int j=mys; j<=mye; j++)
00350 for (register int i=mxs; i<=mxe; i++)
00351 Collision(f[base::idx(i,j,Nx)],dt);
00352
00353 return DataType(1.);
00354 }
00355
00356 virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00357 DataType* aux=0, const int naux=0, const int scaling=0) const {
00358 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00359 int mxs = base::NGhosts(), mys = base::NGhosts();
00360 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1;
00361 MicroType *f = (MicroType *)fvec.databuffer();
00362
00363 if (type==ConstantMicro) {
00364 MicroType g;
00365 for (register int i=0; i<base::NMicroVar(); i++)
00366 g(i) = aux[i];
00367
00368 for (register int j=mys; j<=mye; j++)
00369 for (register int i=mxs; i<=mxe; i++)
00370 f[base::idx(i,j,Nx)] = g;
00371 }
00372
00373 else if (type==ConstantMacro) {
00374 MacroType q;
00375 if (scaling==base::Physical) {
00376 q(0) = aux[0]/R0;
00377 q(1) = aux[1]/U0;
00378 q(2) = aux[2]/U0;
00379 }
00380 else
00381 for (register int i=0; i<base::NMacroVar(); i++)
00382 q(i) = aux[i];
00383 q(1) *= S0; q(2) *= S0;
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)] = Equilibrium(q);
00388 }
00389
00390 else if (type==GasAtRest) {
00391 MacroType q;
00392 q(0) = GasDensity()/R0;
00393 q(1) = DataType(0.);
00394 q(2) = DataType(0.);
00395 for (register int j=mys; j<=mye; j++)
00396 for (register int i=mxs; i<=mxe; i++)
00397 f[base::idx(i,j,Nx)] = Equilibrium(q);
00398 }
00399 }
00400
00401
00402 virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00403 DataType* aux=0, const int naux=0, const int scaling=0) const {
00404 int Nx = fvec.extents(0);
00405 int b = fvec.bottom();
00406 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
00407 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00408 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00409 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00410 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00411 MicroType *f = (MicroType *)fvec.databuffer();
00412
00413 if (type==Symmetry) {
00414 switch (side) {
00415 case base::Left:
00416 for (register int j=mys; j<=mye; j++) {
00417 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j,Nx)](8);
00418 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00419 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j,Nx)](5);
00420 }
00421 break;
00422 case base::Right:
00423 for (register int j=mys; j<=mye; j++) {
00424 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j,Nx)](7);
00425 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00426 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j,Nx)](4);
00427 }
00428 break;
00429 case base::Bottom:
00430 for (register int i=mxs; i<=mxe; i++) {
00431 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i,mye+1,Nx)](8);
00432 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00433 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i,mye+1,Nx)](7);
00434 }
00435 break;
00436 case base::Top:
00437 for (register int i=mxs; i<=mxe; i++) {
00438 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i,mys-1,Nx)](5);
00439 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00440 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i,mys-1,Nx)](4);
00441 }
00442 break;
00443 }
00444 }
00445
00446 else if (type==SlipWall) {
00447 switch (side) {
00448 case base::Left:
00449 for (register int j=mys; j<=mye; j++) {
00450 if (j>mys) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00451 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00452 if (j<mye) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00453 }
00454 break;
00455 case base::Right:
00456 for (register int j=mys; j<=mye; j++) {
00457 if (j>mys) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00458 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00459 if (j<mye) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00460 }
00461 break;
00462 case base::Bottom:
00463 for (register int i=mxs; i<=mxe; i++) {
00464 if (i>mxs) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i-1,mye+1,Nx)](7);
00465 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00466 if (i<mxe) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i+1,mye+1,Nx)](8);
00467 }
00468 break;
00469 case base::Top:
00470 for (register int i=mxs; i<=mxe; i++) {
00471 if (i>mxs) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i-1,mys-1,Nx)](4);
00472 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00473 if (i<mxe) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i+1,mys-1,Nx)](5);
00474 }
00475 break;
00476 }
00477 }
00478
00479 else if (type==NoSlipWall) {
00480 switch (side) {
00481 case base::Left:
00482 for (register int j=mys; j<=mye; j++) {
00483 if (j>mys) f[b+base::idx(mxe,j,Nx)](7) = f[b+base::idx(mxe+1,j-1,Nx)](5);
00484 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00485 if (j<mye) f[b+base::idx(mxe,j,Nx)](4) = f[b+base::idx(mxe+1,j+1,Nx)](8);
00486 }
00487 break;
00488 case base::Right:
00489 for (register int j=mys; j<=mye; j++) {
00490 if (j>mys) f[b+base::idx(mxs,j,Nx)](8) = f[b+base::idx(mxs-1,j-1,Nx)](4);
00491 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00492 if (j<mye) f[b+base::idx(mxs,j,Nx)](5) = f[b+base::idx(mxs-1,j+1,Nx)](7);
00493 }
00494 break;
00495 case base::Bottom:
00496 for (register int i=mxs; i<=mxe; i++) {
00497 if (i>mxs) f[b+base::idx(i,mye,Nx)](5) = f[b+base::idx(i-1,mye+1,Nx)](7);
00498 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00499 if (i<mxe) f[b+base::idx(i,mye,Nx)](4) = f[b+base::idx(i+1,mye+1,Nx)](8);
00500 }
00501 break;
00502 case base::Top:
00503 for (register int i=mxs; i<=mxe; i++) {
00504 if (i>mxs) f[b+base::idx(i,mys,Nx)](8) = f[b+base::idx(i-1,mys-1,Nx)](4);
00505 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00506 if (i<mxe) f[b+base::idx(i,mys,Nx)](7) = f[b+base::idx(i+1,mys-1,Nx)](5);
00507 }
00508 break;
00509 }
00510 }
00511
00512 else if (type==Inlet) {
00513 if (aux==0 || naux<=0) return;
00514 DataType a0=S0*aux[0], a1=0.;
00515 if (naux>1) a1 = S0*aux[1];
00516 if (scaling==base::Physical) {
00517 a0 /= U0;
00518 a1 /= U0;
00519 }
00520 switch (side) {
00521 case base::Left:
00522 for (register int j=mys; j<=mye; j++) {
00523 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00524 q(1) = a0;
00525 if (naux>1) q(2) = a1;
00526 f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00527 }
00528 break;
00529 case base::Right:
00530 for (register int j=mys; j<=mye; j++) {
00531 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00532 q(1) = a0;
00533 if (naux>1) q(2) = a1;
00534 f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00535 }
00536 break;
00537 case base::Bottom:
00538 for (register int i=mxs; i<=mxe; i++) {
00539 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00540 if (naux==1)
00541 q(2) = a0;
00542 else {
00543 q(1) = a0;
00544 q(2) = a1;
00545 }
00546 f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00547 }
00548 break;
00549 case base::Top:
00550 for (register int i=mxs; i<=mxe; i++) {
00551 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00552 if (naux==1)
00553 q(2) = a0;
00554 else {
00555 q(1) = a0;
00556 q(2) = a1;
00557 }
00558 f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00559 }
00560 break;
00561 }
00562 }
00563
00564 else if (type==Outlet) {
00565 switch (side) {
00566 case base::Left:
00567 for (register int j=mys; j<=mye; j++) {
00568 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00569 if (q(0)>0) f[b+base::idx(mxe,j,Nx)] = f[b+base::idx(mxe+1,j,Nx)]/q(0);
00570 }
00571 break;
00572 case base::Right:
00573 for (register int j=mys; j<=mye; j++) {
00574 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00575 if (q(0)>0) f[b+base::idx(mxs,j,Nx)] = f[b+base::idx(mxs-1,j,Nx)]/q(0);
00576 }
00577 break;
00578 case base::Bottom:
00579 for (register int i=mxs; i<=mxe; i++) {
00580 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00581 if (q(0)>0) f[b+base::idx(i,mye,Nx)] = f[b+base::idx(i,mye+1,Nx)]/q(0);
00582 }
00583 break;
00584 case base::Top:
00585 for (register int i=mxs; i<=mxe; i++) {
00586 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00587 if (q(0)>0) f[b+base::idx(i,mys,Nx)] = f[b+base::idx(i,mys-1,Nx)]/q(0);
00588 }
00589 break;
00590 }
00591 }
00592
00593
00594 else if (type==Pressure) {
00595 if (aux==0 || naux<=0) return;
00596 DataType a0=aux[0];
00597 if (scaling==base::Physical)
00598 a0 /= R0;
00599 switch (side) {
00600 case base::Left:
00601 for (register int j=mys; j<=mye; j++) {
00602 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00603 q(0) = a0;
00604 f[b+base::idx(mxe,j,Nx)] = Equilibrium(q);
00605 }
00606 break;
00607 case base::Right:
00608 for (register int j=mys; j<=mye; j++) {
00609 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00610 q(0) = a0;
00611 f[b+base::idx(mxs,j,Nx)] = Equilibrium(q);
00612 }
00613 break;
00614 case base::Bottom:
00615 for (register int i=mxs; i<=mxe; i++) {
00616 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00617 q(0) = a0;
00618 f[b+base::idx(i,mye,Nx)] = Equilibrium(q);
00619 }
00620 break;
00621 case base::Top:
00622 for (register int i=mxs; i<=mxe; i++) {
00623 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00624 q(0) = a0;
00625 f[b+base::idx(i,mys,Nx)] = Equilibrium(q);
00626 }
00627 break;
00628 }
00629 }
00630
00631 else if (type==SlidingWall) {
00632 if (aux==0 || naux<=0) return;
00633 if (aux==0 || naux<=0) return;
00634 DataType a0=S0*aux[0];
00635 if (scaling==base::Physical)
00636 a0 /= U0;
00637 switch (side) {
00638 case base::Left:
00639 for (register int j=mys; j<=mye; j++) {
00640 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,Nx)]);
00641 DataType rd1 = f[b+base::idx(mxe+1,j,Nx)](5), rd2 = f[b+base::idx(mxe+1,j,Nx)](8);
00642 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00643 DataType pw = DataType(1.)-qw;
00644 if (j<mye) f[b+base::idx(mxe,j+1,Nx)](7) = pw*rd1+qw*rd2;
00645 f[b+base::idx(mxe,j,Nx)](1) = f[b+base::idx(mxe+1,j,Nx)](2);
00646 if (j>mys) f[b+base::idx(mxe,j-1,Nx)](4) = qw*rd1+pw*rd2;
00647 }
00648 break;
00649 case base::Right:
00650 for (register int j=mys; j<=mye; j++) {
00651 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,Nx)]);
00652 DataType rd1 = f[b+base::idx(mxs-1,j,Nx)](4), rd2 = f[b+base::idx(mxs-1,j,Nx)](7);
00653 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00654 DataType pw = DataType(1.)-qw;
00655 if (j<mye) f[b+base::idx(mxs,j+1,Nx)](8) = pw*rd1+qw*rd2;
00656 f[b+base::idx(mxs,j,Nx)](2) = f[b+base::idx(mxs-1,j,Nx)](1);
00657 if (j>mys) f[b+base::idx(mxs,j-1,Nx)](5) = qw*rd1+pw*rd2;
00658 }
00659 break;
00660 case base::Bottom:
00661 for (register int i=mxs; i<=mxe; i++) {
00662 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,Nx)]);
00663 DataType rd1 = f[b+base::idx(i,mye+1,Nx)](7), rd2 = f[b+base::idx(i,mye+1,Nx)](8);
00664 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00665 DataType pw = DataType(1.)-qw;
00666 if (i<mxe) f[b+base::idx(i+1,mye,Nx)](5) = pw*rd1+qw*rd2;
00667 f[b+base::idx(i,mye,Nx)](3) = f[b+base::idx(i,mye+1,Nx)](6);
00668 if (i>mxs) f[b+base::idx(i-1,mye,Nx)](4) = qw*rd1+pw*rd2;
00669 }
00670 break;
00671 case base::Top:
00672 for (register int i=mxs; i<=mxe; i++) {
00673 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,Nx)]);
00674 DataType rd1 = f[b+base::idx(i,mys-1,Nx)](4), rd2 = f[b+base::idx(i,mys-1,Nx)](5);
00675 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
00676 DataType pw = DataType(1.)-qw;
00677 if (i<mxe) f[b+base::idx(i+1,mys,Nx)](8) = pw*rd1+qw*rd2;
00678 f[b+base::idx(i,mys,Nx)](6) = f[b+base::idx(i,mys-1,Nx)](3);
00679 if (i>mxs) f[b+base::idx(i-1,mys,Nx)](7) = qw*rd1+pw*rd2;
00680 }
00681 break;
00682 }
00683 }
00684 }
00685
00686 virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
00687 const MicroType* f, const point_type* xc, const DataType* distance,
00688 const point_type* normal, DataType* aux=0, const int naux=0,
00689 const int scaling=0) const {
00690
00691 if (type==GFMNoSlipWall || type==GFMSlipWall) {
00692 DataType fac = S0;
00693 if (scaling==base::Physical)
00694 fac /= U0;
00695 for (int n=0; n<nc; n++) {
00696 MacroType q=MacroVariables(f[n]);
00697 DataType u=-q(1);
00698 DataType v=-q(2);
00699
00700
00701 if (naux>=2) {
00702 u += fac*aux[naux*n];
00703 v += fac*aux[naux*n+1];
00704 }
00705 if (type==GFMNoSlipWall) {
00706
00707 q(1) += DataType(2.)*u;
00708 q(2) += DataType(2.)*v;
00709 }
00710 else if (type==GFMSlipWall) {
00711
00712 DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1));
00713 q(1) += vl*normal[n](0);
00714 q(2) += vl*normal[n](1);
00715 }
00716 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
00717 }
00718 }
00719
00720 else if (type==GFMExtrapolation) {
00721 for (int n=0; n<nc; n++) {
00722 MacroType q=MacroVariables(f[n]);
00723 DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1);
00724 q(1) = vl*normal[n](0);
00725 q(2) = vl*normal[n](1);
00726 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
00727 }
00728 }
00729
00730
00731 else if (type==GFMBounceBack) {
00732 DCoords dx = base::GH().worldStep(fvec.stepsize());
00733 DataType fac = S0;
00734 if(scaling==base::Physical) fac/= U0;
00735
00736 DataType dt_temp = dx(0)*S0;
00737
00738
00739
00740
00741
00742
00743 DataType distB[2], distBScalar, cosa, leng, d_dist, qNorm, mdx2, mdy2, zaehler, nenner, lengXi;
00744
00745 for (int n=0; n<nc; n++){
00746 distB[0] = (double &)distance[n]*normal[n](0);
00747 distB[1] = (double &)distance[n]*normal[n](1);
00748
00749
00750
00751
00752
00753
00754 for(int ind=0; ind<base::NMicroVar(); ind++){
00755
00756 mdx2 = mdx[ind]*mdx[ind];
00757 mdy2 = mdy[ind]*mdy[ind];
00758 distBScalar = std::sqrt(std::pow(distB[0],2)+std::pow(distB[1],2));
00759 leng = distBScalar*distBScalar*std::sqrt(mdx2+mdy2)/(distB[0]*mdx[ind]+distB[1]*mdy[ind]);
00760 lengXi = std::sqrt(mdx2*dx(0)*dx(0)+mdy2*dx(1)*dx(1));
00761 d_dist = lengXi-leng;
00762 qNorm = d_dist/lengXi;
00763
00764 if(qNorm>0 && qNorm < 1.){
00765
00766
00767
00768 int idxA, idxB;
00769 if(ind==0) idxA=0;
00770 else if(ind==1) idxA=2;
00771 else if(ind==2) idxA=1;
00772 else if(ind==3) idxA=6;
00773 else if(ind==4) idxA=8;
00774 else if(ind==5) idxA=7;
00775 else if(ind==6) idxA=3;
00776 else if(ind==7) idxA=5;
00777 else if(ind==8) idxA=4;
00778
00779 idxB=idxA;
00780 if(qNorm>=DataType(0.5)) idxB=ind;
00781
00782
00783
00784
00785
00786 DataType cn[2], cnA[2], cnB[2];
00787 cn[0] = xc[n](0);
00788 cn[1] = xc[n](1);
00789
00790
00791 cnA[0] = cn[0]+mdx[ind]*dx(0);
00792 cnA[1] = cn[1]+mdy[ind]*dx(1);
00793
00794
00795 cnB[0] = cnA[0];
00796 cnB[1] = cnA[1];
00797
00798 int k = 2;
00799 if(qNorm<DataType(0.5)){
00800 cnB[0] = cn[0] + k*mdx[ind]*dx(0);
00801 cnB[1] = cn[1] + k*mdy[ind]*dx(1);
00802 }
00803
00804
00805 DataType distFuncA, distFuncB;
00806
00807 int cnA_lc[2], cnB_lc[2];
00808 cnA_lc[0] = base::GH().localCoords((const DataType *) cnA)(0);
00809 cnA_lc[1] = base::GH().localCoords((const DataType *) cnA)(1);
00810 cnB_lc[0] = base::GH().localCoords((const DataType *) cnB)(0);
00811 cnB_lc[1] = base::GH().localCoords((const DataType *) cnB)(1);
00812
00813 int lc_low[2], lc_up[2];
00814 lc_low[0] = fvec.lower(0); lc_low[1] = fvec.lower(1);
00815 lc_up[0] = fvec.upper(0); lc_up[1] = fvec.upper(1);
00816
00817 if ( (cnA_lc[0] > lc_low[0]) && (cnB_lc[0] > lc_low[0]) ) {
00818 if ( (cnA_lc[1] > lc_low[1]) && (cnB_lc[1] > lc_low[1]) ) {
00819
00820 if ( (cnA_lc[0] < lc_up[0]) && (cnB_lc[0] < lc_up[0]) ) {
00821 if ( (cnA_lc[1] < lc_up[1]) && (cnB_lc[1] < lc_up[1]) ) {
00822
00823
00824 distFuncA = fvec(base::GH().localCoords((const DataType*) cnA))(idxA);
00825 distFuncB = fvec(base::GH().localCoords((const DataType*) cnB))(idxB);
00826
00827 DataType intp1,intp2,intp,q2;
00828 q2 = DataType(2.)*qNorm;
00829
00830
00831 if(qNorm<DataType(0.5)){
00832 intp1=q2*distFuncA;
00833 intp2=(DataType(1.)-q2)*distFuncB;
00834 }
00835
00836 else{
00837 if(q2!=0){
00838 intp1 = (DataType(1.)/q2)*distFuncA;
00839 intp2 = ((q2-DataType(1.))/q2)*distFuncB;
00840 }
00841 else{
00842 intp1 = DataType(0.);
00843 intp2 = intp1;
00844 }
00845 }
00846
00847
00848 intp = intp1+intp2;
00849
00850 fvec(base::GH().localCoords((const DataType*) cn))(ind) = intp;
00851
00852 }
00853 }
00854 }
00855 }
00856
00857
00858
00859
00860 if(naux>=2){
00861 DataType u = fac*aux[naux*n];
00862 DataType v = fac*aux[naux*n+1];
00863
00864 DataType weight;
00865 weight = DataType(1.)/DataType(3);
00866 if(ind==0) weight *= 0;
00867 else if(ind==4 || ind ==5 || ind==7 || ind==8) weight *= DataType(1.)/DataType(4.);
00868
00869 DataType velForce;
00870
00871 if(qNorm<DataType(0.5)) velForce = DataType(2.)*weight;
00872 else velForce = (DataType(1.)/qNorm)*weight;
00873
00874 velForce *= ((mdx[ind]*u)+(mdy[ind]*v));
00875 fvec(base::GH().localCoords((const DataType*) cn))(ind) +=velForce;
00876 }
00877 }
00878 }
00879 }
00880 }
00881 }
00882
00883 virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
00884 const int skip_ghosts=1) const {
00885 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00886 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
00887 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
00888 MicroType *f = (MicroType *)fvec.databuffer();
00889 DataType *work = (DataType *)workvec.databuffer();
00890 DCoords dx = base::GH().worldStep(fvec.stepsize());
00891 DataType dt = dx(0)*S0/U0;
00892
00893 if ((cnt>=1 && cnt<=9) || (cnt>=16 && cnt<=18) || (cnt>=23 && cnt<=25)) {
00894 for (register int j=mys; j<=mye; j++)
00895 for (register int i=mxs; i<=mxe; i++) {
00896 MacroType q=MacroVariables(f[base::idx(i,j,Nx)]);
00897 switch(cnt) {
00898 case 1:
00899 if (method[0]==0) work[base::idx(i,j,Nx)]=q(0)*R0+rhop;
00900 else
00901 work[base::idx(i,j,Nx)]=q(0)*R0;
00902 break;
00903 case 2:
00904 work[base::idx(i,j,Nx)]=q(1)*U0/S0;
00905 break;
00906 case 3:
00907 work[base::idx(i,j,Nx)]=q(2)*U0/S0;
00908 break;
00909 case 5:
00910 work[base::idx(i,j,Nx)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure());
00911 break;
00912 case 6:
00913 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure();
00914 break;
00915 case 8: {
00916 MicroType feq = Equilibrium(q);
00917 work[base::idx(i,j,Nx)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt),
00918 GasSpeedofSound(),dt);
00919 break;
00920 }
00921 case 9:
00922 work[base::idx(i,j,Nx)]=std::sqrt(q(1)*q(1)+q(2)*q(2))*U0/S0;
00923 break;
00924 case 16:
00925 case 17:
00926 case 18: {
00927 MicroType feq = Equilibrium(q);
00928 DataType omega;
00929 if (turbulence == laminar)
00930 omega = Omega(dt);
00931 else if (turbulence == LES_Smagorinsky)
00932 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt);
00933 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*DeviatoricStress(f[base::idx(i,j,Nx)],feq,omega)(cnt-16);
00934 break;
00935 }
00936 case 23:
00937 case 24:
00938 case 25: {
00939 DataType omega;
00940 if (turbulence == laminar)
00941 omega = Omega(dt);
00942 else if (turbulence == LES_Smagorinsky) {
00943 MicroType feq = Equilibrium(q);
00944 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,Nx)],feq,q,dt);
00945 }
00946 work[base::idx(i,j,Nx)]=U0/S0*U0/S0*R0*Stress(f[base::idx(i,j,Nx)],q,omega)(cnt-23);
00947 if (cnt==23 || cnt==25) work[base::idx(i,j,Nx)]-=BasePressure();
00948 break;
00949 }
00950 }
00951 }
00952 }
00953 else if ((cnt>=10 && cnt<=11) || cnt==15 || (cnt>=30 && cnt<=32)) {
00954 macro_grid_data_type qgrid(fvec.bbox());
00955 MacroType *q = (MacroType *)qgrid.databuffer();
00956 for (register int j=mys; j<=mye; j++)
00957 for (register int i=mxs; i<=mxe; i++) {
00958 q[base::idx(i,j,Nx)]=MacroVariables(f[base::idx(i,j,Nx)]);
00959 q[base::idx(i,j,Nx)](0)*=R0;
00960 q[base::idx(i,j,Nx)](1)*=U0/S0;
00961 q[base::idx(i,j,Nx)](2)*=U0/S0;
00962 }
00963
00964 if (skip_ghosts==0) {
00965 switch(cnt) {
00966 case 30:
00967 mxs+=1; mxe-=1;
00968 break;
00969 case 32:
00970 mys+=1; mye-=1;
00971 break;
00972 case 10:
00973 case 11:
00974 case 15:
00975 case 31:
00976 mxs+=1; mxe-=1; mys+=1; mye-=1;
00977 break;
00978 }
00979 }
00980
00981 DataType dxux, dxuy, dyux, dyuy;
00982 for (register int j=mys; j<=mye; j++)
00983 for (register int i=mxs; i<=mxe; i++) {
00984 if (cnt==15 || cnt==30)
00985 dxux=(q[base::idx(i+1,j,Nx)](1)-q[base::idx(i-1,j,Nx)](1))/(2.*dx(0));
00986 if (cnt==15 || cnt==32)
00987 dyuy=(q[base::idx(i,j+1,Nx)](2)-q[base::idx(i,j-1,Nx)](2))/(2.*dx(1));
00988 if (cnt==10 || cnt==11 || cnt==15 || cnt==31) {
00989 dxuy=(q[base::idx(i+1,j,Nx)](2)-q[base::idx(i-1,j,Nx)](2))/(2.*dx(0));
00990 dyux=(q[base::idx(i,j+1,Nx)](1)-q[base::idx(i,j-1,Nx)](1))/(2.*dx(1));
00991 }
00992
00993 switch (cnt) {
00994 case 10:
00995 work[base::idx(i,j,Nx)]=dxuy-dyux;
00996 break;
00997 case 11:
00998 work[base::idx(i,j,Nx)]=std::abs(dxuy-dyux);
00999 break;
01000 case 15:
01001 work[base::idx(i,j,Nx)]=DataType(-0.5)*(4.*dxux*dxux+8.*dyux*dxuy+4.*dyuy*dyuy);
01002 break;
01003 case 30:
01004 work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dxux;
01005 break;
01006 case 31:
01007 work[base::idx(i,j,Nx)]=q[base::idx(i,j,Nx)](0)*GasViscosity()*(dxuy+dyux);
01008 break;
01009 case 32:
01010 work[base::idx(i,j,Nx)]=2.*q[base::idx(i,j,Nx)](0)*GasViscosity()*dyuy;
01011 break;
01012 }
01013 }
01014 }
01015 }
01016
01017 virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01018 const int skip_ghosts=1) const {
01019 int Nx = fvec.extents(0), Ny = fvec.extents(1);
01020 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts;
01021 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1;
01022 MicroType *f = (MicroType *)fvec.databuffer();
01023 DataType *work = (DataType *)workvec.databuffer();
01024
01025 for (register int j=mys; j<=mye; j++)
01026 for (register int i=mxs; i<=mxe; i++) {
01027 switch(cnt) {
01028 case 1:
01029 f[base::idx(i,j,Nx)](0)=work[base::idx(i,j,Nx)]/R0;
01030 break;
01031 case 2:
01032 f[base::idx(i,j,Nx)](1)=work[base::idx(i,j,Nx)]/(U0/S0);
01033 break;
01034 case 3:
01035 f[base::idx(i,j,Nx)](2)=work[base::idx(i,j,Nx)]/(U0/S0);
01036 f[base::idx(i,j,Nx)] = Equilibrium((const MacroType &)f[base::idx(i,j,Nx)]);
01037 break;
01038 }
01039 }
01040 }
01041
01042 virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
01043 const int verbose) const {
01044 int Nx = fvec.extents(0);
01045 int b = fvec.bottom();
01046 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01047 BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
01048 int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01049 int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01050 int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01051 int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01052 MicroType *f = (MicroType *)fvec.databuffer();
01053
01054 int result = 1;
01055 for (register int j=mys; j<=mye; j++)
01056 for (register int i=mxs; i<=mxe; i++)
01057 for (register int l=0; l<base::NMicroVar(); l++)
01058 if (!(f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l)>-1.e37)) {
01059 result = 0;
01060 if (verbose)
01061 std::cerr << "D2Q9-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << ")(" << l << ")="
01062 << f[b+base::idx(i-mdx[l],j-mdy[l],Nx)](l) << " " << fvec.bbox() << std::endl;
01063 }
01064 return result;
01065 }
01066
01067 inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
01068 TensorType P;
01069 if (method[1]==0) {
01070 P(0) = q(0)*q(1)*q(1)-(f(1)+f(2)+f(4)+f(5)+f(7)+f(8))+cs2*q(0);
01071 P(1) = q(0)*q(1)*q(2)-(f(4)-f(5)-f(7)+f(8));
01072 P(2) = q(0)*q(2)*q(2)-(f(3)+f(4)+f(5)+f(6)+f(7)+f(8))+cs2*q(0);
01073 P *= -(DataType(1.0)-DataType(0.5)*om);
01074 }
01075 else {
01076 MicroType feq = Equilibrium(q);
01077 P = DeviatoricStress(f,feq,om);
01078 }
01079 P(0) -= LatticeBasePressure(q(0));
01080 P(2) -= LatticeBasePressure(q(0));
01081 return P;
01082 }
01083
01084 inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
01085 TensorType Sigma;
01086 if (method[2]==0) {
01087 MicroType fneq = feq - f;
01088 Sigma(0) = fneq(1)+fneq(2)+fneq(4)+fneq(5)+fneq(7)+fneq(8);
01089 Sigma(1) = fneq(4)-fneq(5)-fneq(7)+fneq(8);
01090 Sigma(2) = fneq(3)+fneq(4)+fneq(5)+fneq(6)+fneq(7)+fneq(8);
01091 Sigma *= -(DataType(1.0)-DataType(0.5)*om);
01092 }
01093 else {
01094 MacroType q = MacroVariables(f);
01095 Sigma = Stress(f,q,om);
01096 Sigma(0) += LatticeBasePressure(q(0));
01097 Sigma(2) += LatticeBasePressure(q(0));
01098 }
01099 return Sigma;
01100 }
01101
01102 virtual int NMethodOrder() const { return 2; }
01103
01104 inline const DataType& L0() const { return base::LengthScale(); }
01105 inline const DataType& T0() const { return base::TimeScale(); }
01106 inline void SetDensityScale(const DataType r0) { R0 = r0; }
01107 inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
01108 inline void SetSpeedUp(const DataType s0) { S0 = s0; }
01109 inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
01110
01111 inline const DataType& DensityScale() const { return R0; }
01112 inline const DataType VelocityScale() const { return U0/S0; }
01113
01114 inline const DataType& SpeedUp() const { return S0; }
01115
01116
01117 inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
01118 inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
01119 inline virtual DataType LatticeBasePressure(const DataType rho) const {
01120 return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
01121 }
01122
01123
01124 inline void SetGas(DataType rho, DataType nu, DataType cs) {
01125 rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
01126 SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
01127 }
01128 inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
01129 inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q,
01130 const DataType dt) const {
01131 DataType M2 = 0.0;
01132
01133 M2 += (f(4) - feq(4))*(f(4) - feq(4));
01134 M2 += (f(5) - feq(5))*(f(5) - feq(5));
01135 M2 += (f(7) - feq(7))*(f(7) - feq(7));
01136 M2 += (f(8) - feq(8))*(f(8) - feq(8));
01137 M2 = std::sqrt(M2);
01138
01139 DataType tau = DataType(1.0)/Omega(dt);
01140 DataType om_turb = (DataType(2.0)/( tau
01141 + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
01142
01143 return ( om_turb );
01144 }
01145 inline const int TurbulenceType() const { return turbulence; }
01146 inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
01147 inline const DataType& GasDensity() const { return rhop; }
01148 inline const DataType& GasViscosity() const { return nup; }
01149 inline const DataType& GasSpeedofSound() const { return csp; }
01150 inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
01151 { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
01152
01153 inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
01154 inline virtual DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); }
01155 inline virtual DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); }
01156
01157 protected:
01158 DataType cs2, cs22, cssq;
01159 DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
01160 DataType Cs_Smagorinsky, turbulence;
01161 int method[3], mdx[9], mdy[9];
01162 };
01163
01164
01165 #endif