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