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