00001
00002
00003
00004
00005
00006 #ifndef LBM_D3Q19_H
00007 #define LBM_D3Q19_H
00008
00016 #include "LBMBase.h"
00017 #include <cfloat>
00018
00019 #define LB_FACTOR 1.0e5
00020
00051 #ifndef NUMPLUS
00052 #define NUMPLUS 0
00053 #endif
00054 #define NUMMICRODIST (19+NUMPLUS)
00055
00056 template <class DataType>
00057 class LBMD3Q19 : public LBMBase<Vector<DataType,NUMMICRODIST>, Vector<DataType,4>, 3> {
00058 typedef LBMBase<Vector<DataType,NUMMICRODIST>, Vector<DataType,4>, 3> base;
00059 public:
00060 typedef typename base::vec_grid_data_type vec_grid_data_type;
00061 typedef typename base::grid_data_type grid_data_type;
00062 typedef typename base::MicroType MicroType;
00063 typedef typename base::MacroType MacroType;
00064 typedef GridData<MacroType,3> macro_grid_data_type;
00065 typedef typename base::SideName SideName;
00066 typedef typename base::point_type point_type;
00067 typedef Vector<DataType,6> TensorType;
00068
00069 enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
00070 enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall, CharacteristicOutlet, CharacteristicInlet, NoSlipWallEntranceExit, Periodic };
00071 enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall, GFMWallLaw };
00072 enum TurbulenceModel { laminar, LES_Smagorinsky, LES_SmagorinskyMemory, LES_dynamic, LES_dynamicMemory, WALE, CSM };
00073
00074 LBMD3Q19() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.) {
00075 cs2 = DataType(1.)/DataType(3.);
00076 cs22 = DataType(2.)*cs2;
00077 cssq = DataType(2.)/DataType(9.);
00078 method[0] = 2; method[1] = 0; method[2] = 0;
00079 turbulence = laminar;
00080 Cs_Smagorinsky = DataType(0.2);
00081 mdx[0]=mdx[3]=mdx[4]=mdx[5]=mdx[6]=mdx[11]=mdx[12]=mdx[13]=mdx[14]=0;
00082 mdx[1]=mdx[7]=mdx[9]=mdx[15]=mdx[17]=1;
00083 mdx[2]=mdx[8]=mdx[10]=mdx[16]=mdx[18]=-1;
00084
00085 mdy[0]=mdy[1]=mdy[2]=mdy[5]=mdy[6]=mdy[15]=mdy[16]=mdy[17]=mdy[18]=0;
00086 mdy[3]=mdy[7]=mdy[10]=mdy[11]=mdy[13]=1;
00087 mdy[4]=mdy[8]=mdy[9]=mdy[12]=mdy[14]=-1;
00088
00089 mdz[0]=mdz[1]=mdz[2]=mdz[3]=mdz[4]=mdz[7]=mdz[8]=mdz[9]=mdz[10]=0;
00090 mdz[5]=mdz[11]=mdz[14]=mdz[15]=mdz[18]=1;
00091 mdz[6]=mdz[12]=mdz[13]=mdz[16]=mdz[17]=-1;
00092
00093 stressPath=0;
00094 }
00095
00096 virtual ~LBMD3Q19() {}
00097
00098 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00099 base::LocCtrl = Ctrl.getSubDevice(prefix+"D3Q19");
00100 RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00101 RegisterAt(base::LocCtrl,"Method(2)",method[1]);
00102 RegisterAt(base::LocCtrl,"Method(3)",method[2]);
00103 RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00104 RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00105 RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00106 RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00107 RegisterAt(base::LocCtrl,"Gas_nu",nup);
00108 RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00109 RegisterAt(base::LocCtrl,"Gas_W",Wp);
00110 RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00111 RegisterAt(base::LocCtrl,"SpeedUp",S0);
00112 RegisterAt(base::LocCtrl,"StressPath",stressPath);
00113 }
00114
00115 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00116 base::SetupData(gh,ghosts);
00117 if (rhop>0.) R0 = rhop;
00118 if (csp>0.) {
00119 cs2p = csp*csp;
00120 SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00121 if (gp==DataType(0.))
00122 mfp = GasViscosity()*std::sqrt(std::asin(1.)*DataType(1.4)/cs2p);
00123 else
00124 mfp = GasViscosity()*std::sqrt(std::asin(1.)*gp/cs2p);
00125 }
00126 std::cout << "D3Q19: Gas_rho=" << rhop << " Gas_Cs=" << csp
00127 << " Gas_nu=" << nup << std::endl;
00128 std::cout << "D3Q19: L0=" << L0() << " T0=" << T0() << " U0=" << U0
00129 << " R0=" << R0 << std::endl;
00130 if ( TurbulenceType() == laminar ) {
00131 std::cout << "LBM Setup() Laminar" << std::endl;
00132 }
00133 if ( TurbulenceType() == LES_Smagorinsky || TurbulenceType() == LES_SmagorinskyMemory ) {
00134 std::cout << "LBM Setup() LES_Smagorinsky type "<<TurbulenceType()
00135 <<" : Smagorinsky Constant = " << SmagorinskyConstant() << std::endl;
00136 }
00137 if ( TurbulenceType() == LES_dynamic || TurbulenceType() == LES_dynamicMemory ) {
00138 std::cout << "LBM Setup() LES_dynamic type " <<TurbulenceType()<< std::endl;
00139 assert(base::NGhosts()>3);
00140 }
00141 if ( TurbulenceType() == LES_SmagorinskyMemory || TurbulenceType() == LES_dynamicMemory ) {
00142 if (NUMPLUS<3)
00143 std::cerr << "LBM Setup() LES type "<<TurbulenceType()
00144 <<" requires additional data structures. Use LBMD2Q9plus_LES.h or LBMD2Q9plus_DL.h\n";
00145 assert(NUMPLUS>=3);
00146 }
00147 if ( TurbulenceType() == CSM ) {
00148 std::cout << "LBM Setup() CSM" << std::endl;
00149 }
00150 WriteInit();
00151 }
00152
00153 virtual void WriteInit() const {
00154 int me = MY_PROC;
00155 if (me == VizServer) {
00156 std::ofstream ofs("chem.dat", std::ios::out);
00157 ofs << "RU " << Rp << std::endl;
00158 ofs << "PA " << BasePressure() << std::endl;
00159 ofs << "Sp(01) Vicosity" << std::endl;
00160 ofs << "W(01) " << Wp << std::endl;
00161 ofs << "Sp(02) |Velocity|" << std::endl;
00162 ofs << "W(02) " << 1. << std::endl;
00163 ofs << "Sp(03) vx" << std::endl;
00164 ofs << "W(03) " << 1. << std::endl;
00165 ofs << "Sp(04) vy" << std::endl;
00166 ofs << "W(04) " << 1. << std::endl;
00167 ofs << "Sp(05) vz" << std::endl;
00168 ofs << "W(05) " << 1. << std::endl;
00169 ofs << "Sp(06) |v|" << std::endl;
00170 ofs << "W(06) " << 1. << std::endl;
00171 ofs << "Sp(07) |Stress|" << std::endl;
00172 ofs << "W(07) " << 1. << std::endl;
00173 ofs << "Sp(08) |DevStress|" << std::endl;
00174 ofs << "W(08) " << 1. << std::endl;
00175 ofs << "Sp(09) StrainRate" << std::endl;
00176 ofs << "W(09) " << 1. << std::endl;
00177 ofs << "Sp(10) qcrit" << std::endl;
00178 ofs << "W(10) " << 1. << std::endl;
00179 ofs << "Sp(11) dxx" << std::endl;
00180 ofs << "W(11) " << 1. << std::endl;
00181 ofs << "Sp(12) dxy" << std::endl;
00182 ofs << "W(12) " << 1. << std::endl;
00183 ofs << "Sp(13) dyy" << std::endl;
00184 ofs << "W(13) " << 1. << std::endl;
00185 ofs << "Sp(14) dxz" << std::endl;
00186 ofs << "W(14) " << 1. << std::endl;
00187 ofs << "Sp(15) dyx" << std::endl;
00188 ofs << "W(15) " << 1. << std::endl;
00189 ofs << "Sp(16) dzz" << std::endl;
00190 ofs.close();
00191 }
00192 }
00193
00194 inline virtual MacroType MacroVariables(const MicroType &f) const {
00195 MacroType q;
00196 q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8)+f(9)
00197 +f(10)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18);
00198 if (method[0]!=3) {
00199 q(1)=(f(1)-f(2)+f(7)-f(8)+f(9)-f(10)+f(15)-f(16)+f(17)-f(18))/q(0);
00200 q(2)=(f(3)-f(4)+f(7)-f(8)-f(9)+f(10)+f(11)-f(12)+f(13)-f(14))/q(0);
00201 q(3)=(f(5)-f(6)+f(11)-f(12)-f(13)+f(14)+f(15)-f(16)-f(17)+f(18))/q(0);
00202 }
00203 else {
00204 q(1)=(f(1)-f(2)+f(7)-f(8)+f(9)-f(10)+f(15)-f(16)+f(17)-f(18))/(rhop/R0);
00205 q(2)=(f(3)-f(4)+f(7)-f(8)-f(9)+f(10)+f(11)-f(12)+f(13)-f(14))/(rhop/R0);
00206 q(3)=(f(5)-f(6)+f(11)-f(12)-f(13)+f(14)+f(15)-f(16)-f(17)+f(18))/(rhop/R0);
00207 }
00208 return q;
00209 }
00210
00211 inline virtual MicroType Equilibrium(const MacroType &q) const {
00212 MicroType feq;
00213
00214 DataType rho = q(0);
00215 DataType u = q(1);
00216 DataType v = q(2);
00217 DataType w = q(3);
00218
00219 DataType ri=0., rt0=0., rt1=0., rt2=0.;
00220 if (method[0]<3) {
00221 if (method[0]==0) {
00222 ri = DataType(1.);
00223 rt0 = R0 / DataType(3.);
00224 rt1 = R0 / DataType(18.);
00225 rt2 = R0 / DataType(36.);
00226 }
00227 if (method[0]==1) {
00228 ri = rho;
00229 rt0 = DataType(1.) / DataType(3.);
00230 rt1 = DataType(1.) / DataType(18.);
00231 rt2 = DataType(1.) / DataType(36.);
00232 }
00233 else if (method[0]==2) {
00234 ri = DataType(1.);
00235 rt0 = rho / DataType(3.);
00236 rt1 = rho / DataType(18.);
00237 rt2 = rho / DataType(36.);
00238 }
00239
00240 DataType usq = u * u;
00241 DataType vsq = v * v;
00242 DataType wsq = w * w;
00243 DataType sumsq = (usq + vsq + wsq) / cs22;
00244 DataType sumsqxy = (usq + vsq) / cs22;
00245 DataType sumsqyz = (vsq + wsq) / cs22;
00246 DataType sumsqxz = (usq + wsq) / cs22;
00247 DataType sumsq2xy = sumsqxy * (DataType(1.) - cs2) / cs2;
00248 DataType sumsq2yz = sumsqyz * (DataType(1.) - cs2) / cs2;
00249 DataType sumsq2xz = sumsqxz * (DataType(1.) - cs2) / cs2;
00250 DataType u2 = usq / cssq;
00251 DataType v2 = vsq / cssq;
00252 DataType w2 = wsq / cssq;
00253 DataType u22 = usq / cs22;
00254 DataType v22 = vsq / cs22;
00255 DataType w22 = wsq / cs22;
00256 DataType ui = u / cs2;
00257 DataType vi = v / cs2;
00258 DataType wi = w / cs2;
00259 DataType uv = ui * vi;
00260 DataType uw = ui * wi;
00261 DataType vw = vi * wi;
00262
00263 feq(0) = rt0 * (ri - sumsq);
00264
00265 feq(1) = rt1 * (ri - sumsq + u2 + ui);
00266 feq(2) = rt1 * (ri - sumsq + u2 - ui);
00267
00268 feq(3) = rt1 * (ri - sumsq + v2 + vi);
00269 feq(4) = rt1 * (ri - sumsq + v2 - vi);
00270
00271 feq(5) = rt1 * (ri - sumsq + w2 + wi);
00272 feq(6) = rt1 * (ri - sumsq + w2 - wi);
00273
00274 feq(7) = rt2 * (ri + sumsq2xy -w22 + ui + vi + uv);
00275 feq(8) = rt2 * (ri + sumsq2xy -w22 - ui - vi + uv);
00276 feq(9) = rt2 * (ri + sumsq2xy -w22 + ui - vi - uv);
00277 feq(10) = rt2 * (ri + sumsq2xy -w22 - ui + vi - uv);
00278
00279 feq(11) = rt2 * (ri + sumsq2yz -u22 + vi + wi + vw);
00280 feq(12) = rt2 * (ri + sumsq2yz -u22 - vi - wi + vw);
00281 feq(13) = rt2 * (ri + sumsq2yz -u22 + vi - wi - vw);
00282 feq(14) = rt2 * (ri + sumsq2yz -u22 - vi + wi - vw);
00283
00284 feq(15) = rt2 * (ri + sumsq2xz -v22 + ui + wi + uw);
00285 feq(16) = rt2 * (ri + sumsq2xz -v22 - ui - wi + uw);
00286 feq(17) = rt2 * (ri + sumsq2xz -v22 + ui - wi - uw);
00287 feq(18) = rt2 * (ri + sumsq2xz -v22 - ui + wi - uw);
00288 }
00289 else if (method[0]==3) {
00290 ri = rhop/R0;
00291 rt0 = DataType(4.) / DataType(9.);
00292 rt1 = DataType(1.) / DataType(9.);
00293 rt2 = DataType(1.) / DataType(36.);
00294 DataType usq = u * u;
00295 DataType vsq = v * v;
00296 DataType wsq = w * w;
00297 DataType sumsq = (usq + vsq + wsq) / cs22;
00298 DataType sumsqxy = (usq + vsq) / cs22;
00299 DataType sumsqyz = (vsq + wsq) / cs22;
00300 DataType sumsqxz = (usq + wsq) / cs22;
00301 DataType sumsq2xy = sumsqxy * (DataType(1.) - cs2) / cs2;
00302 DataType sumsq2yz = sumsqyz * (DataType(1.) - cs2) / cs2;
00303 DataType sumsq2xz = sumsqxz * (DataType(1.) - cs2) / cs2;
00304 DataType u2 = usq / cssq;
00305 DataType v2 = vsq / cssq;
00306 DataType w2 = wsq / cssq;
00307 DataType u22 = usq / cs22;
00308 DataType v22 = vsq / cs22;
00309 DataType w22 = wsq / cs22;
00310 DataType ui = u / cs2;
00311 DataType vi = v / cs2;
00312 DataType wi = w / cs2;
00313 DataType uv = ui * vi;
00314 DataType uw = ui * wi;
00315 DataType vw = vi * wi;
00316
00317 feq(0) = rt0 * ( rho + ri*( - sumsq));
00318
00319 feq(1) = rt1 * ( rho + ri*( - sumsq + u2 + ui));
00320 feq(2) = rt1 * ( rho + ri*( - sumsq + u2 - ui));
00321
00322 feq(3) = rt1 * ( rho + ri*( - sumsq + v2 + vi));
00323 feq(4) = rt1 * ( rho + ri*( - sumsq + v2 - vi));
00324
00325 feq(5) = rt1 * ( rho + ri*( - sumsq + w2 + wi));
00326 feq(6) = rt1 * ( rho + ri*( - sumsq + w2 - wi));
00327
00328 feq(7) = rt2 * ( rho + ri*( sumsq2xy -w22 + ui + vi + uv));
00329 feq(8) = rt2 * ( rho + ri*( sumsq2xy -w22 - ui - vi + uv));
00330 feq(9) = rt2 * ( rho + ri*( sumsq2xy -w22 + ui - vi - uv));
00331 feq(10) = rt2 * ( rho + ri*( sumsq2xy -w22 - ui + vi - uv));
00332
00333 feq(11) = rt2 * ( rho + ri*( sumsq2yz -u22 + vi + wi + vw));
00334 feq(12) = rt2 * ( rho + ri*( sumsq2yz -u22 - vi - wi + vw));
00335 feq(13) = rt2 * ( rho + ri*( sumsq2yz -u22 + vi - wi - vw));
00336 feq(14) = rt2 * ( rho + ri*( sumsq2yz -u22 - vi + wi - vw));
00337
00338 feq(15) = rt2 * ( rho + ri*( sumsq2xz -v22 + ui + wi + uw));
00339 feq(16) = rt2 * ( rho + ri*( sumsq2xz -v22 - ui - wi + uw));
00340 feq(17) = rt2 * ( rho + ri*( sumsq2xz -v22 + ui - wi - uw));
00341 feq(18) = rt2 * ( rho + ri*( sumsq2xz -v22 - ui + wi - uw));
00342 }
00343 return feq;
00344 }
00345
00346 inline virtual void Collision(MicroType &f, const DataType dt) const {
00347 MacroType q = MacroVariables(f);
00348 MicroType feq = Equilibrium(q);
00349
00350 DataType omega=1.;
00351 if (turbulence == laminar ) {
00352 omega = Omega(dt);
00353 }
00354 else
00355 omega = Omega_LES_Smagorinsky(f,feq,q(0),Omega(dt),dt);
00356
00357 for (register int qi=0; qi<19; qi++)
00358 f(qi)=(DataType(1.)-omega)*f(qi) + omega*feq(qi);
00359 }
00360
00361 virtual int IncomingIndices(const int side, int indices[]) const {
00362 switch (side) {
00363 case base::Left:
00364 indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;
00365 break;
00366 case base::Right:
00367 indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;
00368 break;
00369 case base::Bottom:
00370 indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;
00371 break;
00372 case base::Top:
00373 indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;
00374 break;
00375 case base::Front:
00376 indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;
00377 break;
00378 case base::Back:
00379 indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;
00380 break;
00381 }
00382 return 5;
00383 }
00384
00385 virtual int OutgoingIndices(const int side, int indices[]) const {
00386 switch (side) {
00387 case base::Left:
00388 indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;
00389 break;
00390 case base::Right:
00391 indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;
00392 break;
00393 case base::Bottom:
00394 indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;
00395 break;
00396 case base::Top:
00397 indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;
00398 break;
00399 case base::Front:
00400 indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;
00401 break;
00402 case base::Back:
00403 indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;
00404 break;
00405 }
00406 return 5;
00407 }
00408
00409
00410 virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb,
00411 const int side) const {
00412 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00413 int b = fvec.bottom();
00414 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00415 && fvec.stepsize(2)==bb.stepsize(2));
00416 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00417 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00418 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00419 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00420 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00421 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00422 MicroType *f = (MicroType *)fvec.databuffer();
00423
00424 #ifdef DEBUG_PRINT_FIXUP
00425 ( comm_service::log() << "Reverse streaming at Side: " << side << " of "
00426 << fvec.bbox() << " using " << bb << "\n").flush();
00427 #endif
00428
00429 switch (side) {
00430 case base::Left:
00431 for (register int k=mzs; k<=mze; k++)
00432 for (register int j=mys; j<=mye; j++) {
00433 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](2);
00434 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](8);
00435 f[b+base::idx(mxs,j,k,Nx,Ny)](10)= f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](10);
00436 f[b+base::idx(mxs,j,k,Nx,Ny)](16)= f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](16);
00437 f[b+base::idx(mxs,j,k,Nx,Ny)](18)= f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](18);
00438 }
00439 break;
00440 case base::Right:
00441 for (register int k=mzs; k<=mze; k++)
00442 for (register int j=mys; j<=mye; j++) {
00443 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](1);
00444 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](7);
00445 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](9);
00446 f[b+base::idx(mxe,j,k,Nx,Ny)](15)= f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](15);
00447 f[b+base::idx(mxe,j,k,Nx,Ny)](17)= f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](17);
00448 }
00449 break;
00450 case base::Bottom:
00451 for (register int k=mzs; k<=mze; k++)
00452 for (register int i=mxs; i<=mxe; i++) {
00453 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](4);
00454 f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](8);
00455 f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](9);
00456 f[b+base::idx(i,mys,k,Nx,Ny)](12)= f[b+base::idx(i,mys-1,k-1,Nx,Ny)](12);
00457 f[b+base::idx(i,mys,k,Nx,Ny)](14)= f[b+base::idx(i,mys-1,k+1,Nx,Ny)](14);
00458 }
00459 break;
00460 case base::Top:
00461 for (register int k=mzs; k<=mze; k++)
00462 for (register int i=mxs; i<=mxe; i++) {
00463 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](3);
00464 f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](7);
00465 f[b+base::idx(i,mye,k,Nx,Ny)](10)= f[b+base::idx(i-1,mye+1,k,Nx,Ny)](10);
00466 f[b+base::idx(i,mye,k,Nx,Ny)](11)= f[b+base::idx(i,mye+1,k+1,Nx,Ny)](11);
00467 f[b+base::idx(i,mye,k,Nx,Ny)](13)= f[b+base::idx(i,mye+1,k-1,Nx,Ny)](13);
00468 }
00469 break;
00470 case base::Front:
00471 for (register int j=mys; j<=mye; j++)
00472 for (register int i=mxs; i<=mxe; i++) {
00473 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](6);
00474 f[b+base::idx(i,j,mzs,Nx,Ny)](12)= f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](12);
00475 f[b+base::idx(i,j,mzs,Nx,Ny)](13)= f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](13);
00476 f[b+base::idx(i,j,mzs,Nx,Ny)](16)= f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](16);
00477 f[b+base::idx(i,j,mzs,Nx,Ny)](17)= f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](17);
00478 }
00479 break;
00480 case base::Back:
00481 for (register int j=mys; j<=mye; j++)
00482 for (register int i=mxs; i<=mxe; i++) {
00483 f[b+base::idx(i,j,mzs,Nx,Ny)](5) = f[b+base::idx(i,j,mzs+1,Nx,Ny)](5);
00484 f[b+base::idx(i,j,mzs,Nx,Ny)](11)= f[b+base::idx(i,j+1,mzs+1,Nx,Ny)](11);
00485 f[b+base::idx(i,j,mzs,Nx,Ny)](14)= f[b+base::idx(i,j-1,mzs+1,Nx,Ny)](14);
00486 f[b+base::idx(i,j,mzs,Nx,Ny)](15)= f[b+base::idx(i+1,j,mzs+1,Nx,Ny)](15);
00487 f[b+base::idx(i,j,mzs,Nx,Ny)](18)= f[b+base::idx(i-1,j,mzs+1,Nx,Ny)](18);
00488 }
00489 break;
00490 }
00491 }
00492
00493 virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec,
00494 const BBox &bb, const double& dt) const {
00495 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00496 int b = fvec.bottom();
00497 MicroType *f = (MicroType *)fvec.databuffer();
00498 MicroType *of = (MicroType *)ovec.databuffer();
00499
00500 #ifdef DEBUG_PRINT_FIXUP
00501 ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox()
00502 << " using " << bb << "\n").flush();
00503 #endif
00504
00505 assert (fvec.extents(0)==ovec.extents(0));
00506 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) && fvec.stepsize(2)==bb.stepsize(2) &&
00507 fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1) && fvec.stepsize(2)==ovec.stepsize(2));
00508 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00509 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00510 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00511 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00512 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00513 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00514
00515 if (turbulence != WALE && turbulence != CSM)
00516 for (register int k=mzs; k<=mze; k++)
00517 for (register int j=mys; j<=mye; j++)
00518 for (register int i=mxs; i<=mxe; i++) {
00519 for (register int n=0; n<base::NMicroVar(); n++)
00520 f[b+base::idx(i,j,k,Nx,Ny)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n],Nx,Ny)](n);
00521 Collision(f[b+base::idx(i,j,k,Nx,Ny)],dt);
00522 }
00523 else if (turbulence == WALE) {
00524 DataType nu = LatticeViscosity(Omega(dt))/S0/S0;
00525 DCoords dx = base::GH().worldStep(fvec.stepsize());
00526 MicroType fxp, fxm, fyp, fym, fzp, fzm;
00527 for (register int k=mzs; k<=mze; k++)
00528 for (register int j=mys; j<=mye; j++)
00529 for (register int i=mxs; i<=mxe; i++) {
00530 for (register int n=0; n<base::NMicroVar(); n++) {
00531 f[b+base::idx(i,j,k,Nx,Ny)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n],Nx,Ny)](n);
00532 fxp(n) = of[b+base::idx(i-mdx[n]+1,j-mdy[n],k-mdz[n],Nx,Ny)](n);
00533 fxm(n) = of[b+base::idx(i-mdx[n]-1,j-mdy[n],k-mdz[n],Nx,Ny)](n);
00534 fyp(n) = of[b+base::idx(i-mdx[n],j-mdy[n]+1,k-mdz[n],Nx,Ny)](n);
00535 fym(n) = of[b+base::idx(i-mdx[n],j-mdy[n]-1,k-mdz[n],Nx,Ny)](n);
00536 fzp(n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n]+1,Nx,Ny)](n);
00537 fzm(n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n]-1,Nx,Ny)](n);
00538 }
00539 LocalCollisionWALE(f[b+base::idx(i,j,k,Nx,Ny)],nu,
00540 MacroVariables(fxp),MacroVariables(fxm),
00541 MacroVariables(fyp),MacroVariables(fym),
00542 MacroVariables(fzp),MacroVariables(fzm),dx,dt);
00543 }
00544 }
00545 else if (turbulence == CSM) {
00546 DCoords dx = base::GH().worldStep(fvec.stepsize());
00547 MicroType fxp, fxm, fyp, fym, fzp, fzm;
00548 for (register int k=mzs; k<=mze; k++)
00549 for (register int j=mys; j<=mye; j++)
00550 for (register int i=mxs; i<=mxe; i++) {
00551 for (register int n=0; n<base::NMicroVar(); n++) {
00552 f[b+base::idx(i,j,k,Nx,Ny)](n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n],Nx,Ny)](n);
00553 fxp(n) = of[b+base::idx(i-mdx[n]+1,j-mdy[n],k-mdz[n],Nx,Ny)](n);
00554 fxm(n) = of[b+base::idx(i-mdx[n]-1,j-mdy[n],k-mdz[n],Nx,Ny)](n);
00555 fyp(n) = of[b+base::idx(i-mdx[n],j-mdy[n]+1,k-mdz[n],Nx,Ny)](n);
00556 fym(n) = of[b+base::idx(i-mdx[n],j-mdy[n]-1,k-mdz[n],Nx,Ny)](n);
00557 fzp(n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n]+1,Nx,Ny)](n);
00558 fzm(n) = of[b+base::idx(i-mdx[n],j-mdy[n],k-mdz[n]-1,Nx,Ny)](n);
00559 }
00560 LocalCollisionCSM(f[b+base::idx(i,j,k,Nx,Ny)],
00561 MacroVariables(fxp),MacroVariables(fxm),
00562 MacroVariables(fyp),MacroVariables(fym),
00563 MacroVariables(fzp),MacroVariables(fzm),dx,dt);
00564 }
00565 }
00566 #ifdef DEBUG_PRINT_FIXUP
00567 ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox()
00568 << " using " << bb << "COMPLETE\n").flush();
00569 #endif
00570 }
00571
00572 inline virtual macro_grid_data_type Filter(vec_grid_data_type &fvec) const {
00573 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00574 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00575 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00576 MicroType *fv = (MicroType *)fvec.databuffer();
00577 macro_grid_data_type qgrid(fvec.bbox());
00578 MacroType *qg = (MacroType *)qgrid.databuffer();
00579 for (register int k=mzs-3; k<=mze+3; k++)
00580 for (register int j=mys-3; j<=mye+3; j++)
00581 for (register int i=mxs-3; i<=mxe+3; i++) {
00582 qg[base::idx(i,j,k,Nx,Ny)] = DataType(0.25)*MacroVariables(fv[base::idx(i-1,j,k,Nx,Ny)])
00583 + DataType(0.5)*MacroVariables(fv[base::idx(i,j,k,Nx,Ny)])
00584 + DataType(0.25)*MacroVariables(fv[base::idx(i+1,j,k,Nx,Ny)]);
00585 }
00586 macro_grid_data_type qgrid0(fvec.bbox());
00587 MacroType *qg0 = (MacroType *)qgrid0.databuffer();
00588 for (register int k=mzs-1; k<=mze+1; k++)
00589 for (register int j=mys-1; j<=mye+1; j++)
00590 for (register int i=mxs-1; i<=mxe+1; i++) {
00591 qg0[base::idx(i,j,k,Nx,Ny)] = DataType(0.25)*qg[base::idx(i,j+1,k,Nx,Ny)]
00592 + DataType(0.5)*qg[base::idx(i,j,k,Nx,Ny)]
00593 + DataType(0.25)*qg[base::idx(i,j-1,k,Nx,Ny)];
00594 }
00595 for (register int k=mzs-1; k<=mze+1; k++)
00596 for (register int j=mys-1; j<=mye+1; j++)
00597 for (register int i=mxs-1; i<=mxe+1; i++) {
00598 qg[base::idx(i,j,k,Nx,Ny)] = DataType(0.25)*qg0[base::idx(i,j,k+1,Nx,Ny)]
00599 + DataType(0.5)*qg0[base::idx(i,j,k,Nx,Ny)]
00600 + DataType(0.25)*qg0[base::idx(i,j,k-1,Nx,Ny)];
00601 }
00602 return qgrid;
00603 }
00604
00605 virtual void CollisionDynamicSmagorinskyLES(vec_grid_data_type &fvec, const double& dt) const {
00615 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00616 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00617 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00618 MicroType *fv = (MicroType *)fvec.databuffer();
00619 MacroType q;
00620 MicroType feq;
00621 DataType csdyn;
00622 DataType dx2=DataType(1.), dxf2 = DataType(4.), dxf = DataType(2.);
00623 DataType M2, SF2;
00624 DataType sfxx, sfxy, sfyy, sfxz, sfyz, sfzz;
00625 DataType lxx, lxy, lyy, lxz, lyz, lzz;
00626 DataType mxx, mxy, myy, mxz, myz, mzz;
00627 DataType omega = Omega(dt);
00628 DataType om = omega;
00629 DataType nu_t = DataType(0.);
00630 TensorType Sigma, S_;
00631 DataType S;
00632 DataType clim = DataType(1.0e7);
00633
00634 macro_grid_data_type qgrid1 = Filter(fvec);
00635 MacroType *sv = (MacroType *)qgrid1.databuffer();
00636 for (register int k=mzs; k<=mze; k++)
00637 for (register int j=mys; j<=mye; j++)
00638 for (register int i=mxs; i<=mxe; i++) {
00639 q = MacroVariables(fv[base::idx(i,j,k,Nx,Ny)]);
00640 feq = Equilibrium(q);
00641 Sigma = DeviatoricStress(fv[base::idx(i,j,k,Nx,Ny)],feq,om)/(DataType(1.0)-DataType(0.5)*om);
00642 S_ = StrainComponents(q(0),Sigma/(DataType(1.0)-DataType(0.5)*om),om,Cs_Smagorinsky);
00643 S = Magnitude(S_);
00644
00645 sfxx = ( sv[base::idx(i+1,j,k,Nx,Ny)](1) - sv[base::idx(i-1,j,k,Nx,Ny)](1) )/dxf;
00646 sfxy = DataType(0.5)*( sv[base::idx(i+1,j,k,Nx,Ny)](2) - sv[base::idx(i-1,j,k,Nx,Ny)](2)
00647 + sv[base::idx(i,j+1,k,Nx,Ny)](1) - sv[base::idx(i,j-1,k,Nx,Ny)](1) )/dxf;
00648 sfyy = ( sv[base::idx(i,j+1,k,Nx,Ny)](2) - sv[base::idx(i,j-1,k,Nx,Ny)](2) )/dxf;
00649 sfxz = DataType(0.5)*( sv[base::idx(i+1,j,k,Nx,Ny)](3) - sv[base::idx(i-1,j,k,Nx,Ny)](3)
00650 + sv[base::idx(i,j,k+1,Nx,Ny)](1) - sv[base::idx(i,j,k-1,Nx,Ny)](1) )/dxf;
00651 sfyz = DataType(0.5)*( sv[base::idx(i,j+1,k,Nx,Ny)](3) - sv[base::idx(i,j-1,k,Nx,Ny)](3)
00652 + sv[base::idx(i,j,k+1,Nx,Ny)](2) - sv[base::idx(i,j,k-1,Nx,Ny)](2) )/dxf;
00653 sfzz = ( sv[base::idx(i,j,k+1,Nx,Ny)](3) - sv[base::idx(i,j,k+1,Nx,Ny)](3) )/dxf;
00654 SF2 = std::sqrt( DataType(2.)*(sfxx*sfxx + DataType(2.)*sfxy*sfxy + sfyy*sfyy
00655 + DataType(2.)*sfxz*sfxz + DataType(2.)*sfyz*sfyz
00656 + sfzz*sfzz) );
00657
00658 mxx = (dxf2*SF2*sfxx - dx2*S*S_(0));
00659 mxy = (dxf2*SF2*sfxy - dx2*S*S_(1));
00660 myy = (dxf2*SF2*sfyy - dx2*S*S_(2));
00661 mxz = (dxf2*SF2*sfxz - dx2*S*S_(3));
00662 myz = (dxf2*SF2*sfyz - dx2*S*S_(4));
00663 mzz = (dxf2*SF2*sfzz - dx2*S*S_(5));
00664 M2 = (mxx*mxx + DataType(2.)*mxy*mxy + myy*myy
00665 + DataType(2.)*mxz*mxz + DataType(2.)*myz*myz
00666 + mzz*mzz);
00667
00668 lxx = (q(1)*q(1) - sv[base::idx(i,j,k,Nx,Ny)](1)* sv[base::idx(i,j,k,Nx,Ny)](1) );
00669 lxy = (q(1)*q(2) - sv[base::idx(i,j,k,Nx,Ny)](1)* sv[base::idx(i,j,k,Nx,Ny)](2) );
00670 lyy = (q(2)*q(2) - sv[base::idx(i,j,k,Nx,Ny)](2)* sv[base::idx(i,j,k,Nx,Ny)](2) );
00671 lxz = (q(1)*q(3) - sv[base::idx(i,j,k,Nx,Ny)](1)* sv[base::idx(i,j,k,Nx,Ny)](3) );
00672 lyz = (q(2)*q(3) - sv[base::idx(i,j,k,Nx,Ny)](2)* sv[base::idx(i,j,k,Nx,Ny)](3) );
00673 lzz = (q(3)*q(3) - sv[base::idx(i,j,k,Nx,Ny)](3)* sv[base::idx(i,j,k,Nx,Ny)](3) );
00674
00675 DataType ntmp = (lxx*lxx + DataType(2.)*lxy*lxy + lyy*lyy
00676 + DataType(2.)*lxz*lxz + DataType(2.)*lyz*lyz
00677 + lzz*lzz);
00678 if (std::isnan(ntmp)==1) { ntmp = 0.0; }
00679 if (std::isnan(M2)==1) { csdyn = Cs_Smagorinsky*Cs_Smagorinsky; }
00680 else if (M2 < mfp/U0) { csdyn = DataType(-0.5)*ntmp/(mfp/U0); }
00681 else if (M2 > clim) { csdyn = DataType(-0.5)*ntmp/(clim); }
00682 else { csdyn = DataType(-0.5)*ntmp/(M2); }
00683 if (csdyn < DataType(-20.0)*Cs_Smagorinsky) { csdyn = DataType(-20.0)*Cs_Smagorinsky; }
00684 else if (csdyn > DataType(20.0)*Cs_Smagorinsky) { csdyn = DataType(20.0)*Cs_Smagorinsky; }
00685 nu_t = std::abs(csdyn)*S;
00686 omega = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));;
00687 if (omega < DataType(1.000001)) { omega = DataType(1.00001); }
00688 else if (omega > DataType(1.999999)) { omega = DataType(1.999999); }
00689 for (register int qi=0; qi<19; qi++)
00690 fv[base::idx(i,j,k,Nx,Ny)](qi) = (DataType(1.)-omega)*fv[base::idx(i,j,k,Nx,Ny)](qi) + omega*feq(qi);
00691 }
00692
00693 }
00694
00695 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00696 const double& t, const double& dt, const int& mpass) const {
00697
00698
00699 DCoords dx = base::GH().worldStep(fvec.stepsize());
00700
00701
00702 if ( (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR)
00703 || (std::fabs(dx(1)/L0()-dx(2)/L0()) > DBL_EPSILON*LB_FACTOR) ) {
00704 std::cerr << "LBMD3Q19 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00705 << " and dx(2)=" << dx(2)/L0()
00706 << " to be equal." << std::endl << std::flush;
00707 std::exit(-1);
00708 }
00709 if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR) {
00710 std::cerr << "LBMD3Q19 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00711 << " to be equal."
00712 << " dt=" << dt << " TO=" << T0()
00713 << std::endl << std::flush;
00714 std::exit(-1);
00715 }
00716
00717 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00718 MicroType *f = (MicroType *)fvec.databuffer();
00719
00720
00721 for (register int k=Nz-1; k>=1; k--)
00722 for (register int j=Ny-1; j>=1; j--)
00723 for (register int i=0; i<=Nx-2; i++) {
00724
00725 f[base::idx(i,j,k,Nx,Ny)](3) = f[base::idx(i,j-1,k,Nx,Ny)](3);
00726 f[base::idx(i,j,k,Nx,Ny)](10) = f[base::idx(i+1,j-1,k,Nx,Ny)](10);
00727
00728
00729
00730 f[base::idx(i,j,k,Nx,Ny)](5) = f[base::idx(i,j,k-1,Nx,Ny)](5);
00731 f[base::idx(i,j,k,Nx,Ny)](18) = f[base::idx(i+1,j,k-1,Nx,Ny)](18);
00732 }
00733
00734 for (register int k=0; k<=Nz-2; k++)
00735 for (register int j=Ny-1; j>=1; j--)
00736 for (register int i=0; i<=Nx-2; i++) {
00737
00738 f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);
00739 }
00740 for (register int k=Nz-1; k>=1; k--)
00741 for (register int j=Ny-1; j>=1; j--)
00742 for (register int i=Nx-1; i>=1; i--) {
00743
00744 f[base::idx(i,j,k,Nx,Ny)](1) = f[base::idx(i-1,j,k,Nx,Ny)](1);
00745 f[base::idx(i,j,k,Nx,Ny)](7) = f[base::idx(i-1,j-1,k,Nx,Ny)](7);
00746
00747 f[base::idx(i,j,k,Nx,Ny)](11) = f[base::idx(i,j-1,k-1,Nx,Ny)](11);
00748
00749 f[base::idx(i,j,k,Nx,Ny)](15) = f[base::idx(i-1,j,k-1,Nx,Ny)](15);
00750 }
00751
00752 for (register int k=0; k<=Nz-2; k++)
00753 for (register int j=0; j<=Ny-2; j++)
00754 for (register int i=Nx-1; i>=1; i--) {
00755
00756 f[base::idx(i,j,k,Nx,Ny)](4) = f[base::idx(i,j+1,k,Nx,Ny)](4);
00757 f[base::idx(i,j,k,Nx,Ny)](9) = f[base::idx(i-1,j+1,k,Nx,Ny)](9);
00758
00759
00760
00761 f[base::idx(i,j,k,Nx,Ny)](6) = f[base::idx(i,j,k+1,Nx,Ny)](6);
00762 f[base::idx(i,j,k,Nx,Ny)](17) = f[base::idx(i-1,j,k+1,Nx,Ny)](17);
00763 }
00764 for (register int k=Nz-1; k>=1; k--)
00765 for (register int j=0; j<=Ny-2; j++)
00766 for (register int i=0; i<=Nx-2; i++) {
00767
00768 f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);
00769 }
00770 for (register int k=0; k<=Nz-2; k++)
00771 for (register int j=0; j<=Ny-2; j++)
00772 for (register int i=0; i<=Nx-2; i++) {
00773
00774 f[base::idx(i,j,k,Nx,Ny)](2) = f[base::idx(i+1,j,k,Nx,Ny)](2);
00775 f[base::idx(i,j,k,Nx,Ny)](8) = f[base::idx(i+1,j+1,k,Nx,Ny)](8);
00776
00777 f[base::idx(i,j,k,Nx,Ny)](12) = f[base::idx(i,j+1,k+1,Nx,Ny)](12);
00778
00779 f[base::idx(i,j,k,Nx,Ny)](16) = f[base::idx(i+1,j,k+1,Nx,Ny)](16);
00780 }
00781
00782
00783 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00784 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00785
00786 if (turbulence==laminar || turbulence==LES_Smagorinsky || turbulence==LES_SmagorinskyMemory) {
00787 for (register int k=mzs; k<=mze; k++)
00788 for (register int j=mys; j<=mye; j++)
00789 for (register int i=mxs; i<=mxe; i++)
00790 Collision(f[base::idx(i,j,k,Nx,Ny)],dt);
00791 }
00792 else if ( turbulence == LES_dynamic || turbulence == LES_dynamicMemory) {
00793 CollisionDynamicSmagorinskyLES(fvec,dt);
00794 }
00795 else if ( turbulence == WALE) {
00796 CollisionWALE(fvec,dt);
00797 }
00798 else if ( turbulence == CSM) {
00799 CollisionCSM(fvec,dt);
00800 }
00801
00802 return DataType(1.);
00803 }
00804
00805 virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00806 DataType* aux=0, const int naux=0, const int scaling=0) const {
00807 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00808 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00809 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00810 MicroType *f = (MicroType *)fvec.databuffer();
00811
00812 if (type==ConstantMicro) {
00813 MicroType g;
00814 for (register int i=0; i<base::NMicroVar(); i++)
00815 g(i) = aux[i];
00816
00817 for (register int k=mzs; k<=mze; k++)
00818 for (register int j=mys; j<=mye; j++)
00819 for (register int i=mxs; i<=mxe; i++)
00820 f[base::idx(i,j,k,Nx,Ny)] = g;
00821 }
00822
00823 else if (type==ConstantMacro) {
00824 MacroType q;
00825 if (scaling==base::Physical) {
00826 q(0) = aux[0]/R0;
00827 q(1) = aux[1]/U0;
00828 q(2) = aux[2]/U0;
00829 q(3) = aux[3]/U0;
00830 }
00831 else
00832 for (register int i=0; i<base::NMacroVar(); i++)
00833 q(i) = aux[i];
00834 q(1) *= S0; q(2) *= S0; q(3) *= S0;
00835
00836 for (register int k=mzs; k<=mze; k++)
00837 for (register int j=mys; j<=mye; j++)
00838 for (register int i=mxs; i<=mxe; i++)
00839 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00840 }
00841
00842 else if (type==GasAtRest) {
00843 MacroType q;
00844 q(0) = GasDensity()/R0;
00845 q(1) = DataType(0.);
00846 q(2) = DataType(0.);
00847 q(3) = DataType(0.);
00848 for (register int k=mzs; k<=mze; k++)
00849 for (register int j=mys; j<=mye; j++)
00850 for (register int i=mxs; i<=mxe; i++)
00851 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00852 }
00853 }
00854
00855
00856 virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00857 DataType* aux=0, const int naux=0, const int scaling=0) const {
00858 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00859 int b = fvec.bottom();
00860 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00861 && fvec.stepsize(2)==bb.stepsize(2));
00862 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00863 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00864 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00865 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00866 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00867 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00868 MicroType *f = (MicroType *)fvec.databuffer();
00869
00870 if (type==Symmetry) {
00871 switch (side) {
00872 case base::Left:
00873 for (register int k=mzs; k<=mze; k++)
00874 for (register int j=mys; j<=mye; j++) {
00875 if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
00876 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00877 if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
00878
00879 if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
00880 if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
00881 }
00882 break;
00883 case base::Right:
00884 for (register int k=mzs; k<=mze; k++)
00885 for (register int j=mys; j<=mye; j++) {
00886 if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
00887 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00888 if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
00889
00890 if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
00891 if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
00892 }
00893 break;
00894 case base::Bottom:
00895 for (register int k=mzs; k<=mze; k++)
00896 for (register int i=mxs; i<=mxe; i++) {
00897 if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
00898 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00899 if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
00900
00901 if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
00902 if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
00903 }
00904 break;
00905 case base::Top:
00906 for (register int k=mzs; k<=mze; k++)
00907 for (register int i=mxs; i<=mxe; i++) {
00908 if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
00909 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00910 if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
00911
00912 if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
00913 if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
00914 }
00915 break;
00916 case base::Front:
00917 for (register int j=mys; j<=mye; j++)
00918 for (register int i=mxs; i<=mxe; i++) {
00919 if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
00920 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00921 if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
00922
00923 if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
00924 if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
00925 }
00926 break;
00927 case base::Back:
00928 for (register int j=mys; j<=mye; j++)
00929 for (register int i=mxs; i<=mxe; i++) {
00930 if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
00931 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00932 if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
00933
00934 if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
00935 if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
00936 }
00937 break;
00938 }
00939 }
00940
00941 else if (type==SlipWall) {
00942 switch (side) {
00943 case base::Left:
00944
00945 for (register int k=mzs; k<=mze; k++)
00946 for (register int j=mys; j<=mye; j++) {
00947 if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
00948 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00949 if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
00950
00951 if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
00952 if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
00953 }
00954 break;
00955 case base::Right:
00956 for (register int k=mzs; k<=mze; k++)
00957 for (register int j=mys; j<=mye; j++) {
00958 if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
00959 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00960 if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
00961
00962 if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
00963 if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
00964 }
00965 break;
00966 case base::Bottom:
00967 for (register int k=mzs; k<=mze; k++)
00968 for (register int i=mxs; i<=mxe; i++) {
00969 if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
00970 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00971 if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
00972
00973 if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
00974 if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
00975 }
00976 break;
00977 case base::Top:
00978 for (register int k=mzs; k<=mze; k++)
00979 for (register int i=mxs; i<=mxe; i++) {
00980 if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
00981 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00982 if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
00983
00984 if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
00985 if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
00986 }
00987 break;
00988 case base::Front:
00989 for (register int j=mys; j<=mye; j++)
00990 for (register int i=mxs; i<=mxe; i++) {
00991 if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
00992 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00993 if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
00994
00995 if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
00996 if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
00997 }
00998 break;
00999 case base::Back:
01000 for (register int j=mys; j<=mye; j++)
01001 for (register int i=mxs; i<=mxe; i++) {
01002 if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
01003 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
01004 if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
01005
01006 if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
01007 if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
01008 }
01009 break;
01010 }
01011 }
01012
01013 else if (type==NoSlipWall) {
01014 switch (side) {
01015 case base::Left:
01016 for (register int k=mzs; k<=mze; k++)
01017 for (register int j=mys; j<=mye; j++) {
01018 if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01019 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01020 if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01021
01022 if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01023 if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01024 }
01025 break;
01026 case base::Right:
01027 for (register int k=mzs; k<=mze; k++)
01028 for (register int j=mys; j<=mye; j++) {
01029 if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
01030 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01031 if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
01032
01033 if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01034 if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01035 }
01036 break;
01037 case base::Bottom:
01038 for (register int k=mzs; k<=mze; k++)
01039 for (register int i=mxs; i<=mxe; i++) {
01040 if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
01041 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01042 if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
01043
01044 if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01045 if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01046 }
01047 break;
01048 case base::Top:
01049 for (register int k=mzs; k<=mze; k++)
01050 for (register int i=mxs; i<=mxe; i++) {
01051 if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
01052 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01053 if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
01054
01055 if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01056 if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01057 }
01058 break;
01059 case base::Front:
01060 for (register int j=mys; j<=mye; j++)
01061 for (register int i=mxs; i<=mxe; i++) {
01062 if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
01063 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01064 if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
01065
01066 if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01067 if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01068 }
01069 break;
01070 case base::Back:
01071 for (register int j=mys; j<=mye; j++)
01072 for (register int i=mxs; i<=mxe; i++) {
01073 if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
01074 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
01075 if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
01076
01077 if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
01078 if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
01079 }
01080 break;
01081 }
01082 }
01083
01084 else if (type==Inlet) {
01085 if (aux==0 || naux<=0) return;
01086 DataType a0=S0*aux[0], a1=0., a2=0.;
01087 if (naux>1) {
01088 a1 = S0*aux[1];
01089 a2 = S0*aux[2];
01090 }
01091 if (scaling==base::Physical) {
01092 a0 /= U0;
01093 a1 /= U0;
01094 a2 /= U0;
01095 }
01096 switch (side) {
01097 case base::Left:
01098 for (register int k=mzs; k<=mze; k++)
01099 for (register int j=mys; j<=mye; j++) {
01100 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01101 q(1) = a0;
01102 if (naux>1) {
01103 q(2) = a1;
01104 q(3) = a2;
01105 }
01106 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
01107 }
01108 break;
01109 case base::Right:
01110 for (register int k=mzs; k<=mze; k++)
01111 for (register int j=mys; j<=mye; j++) {
01112 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01113 q(1) = a0;
01114 if (naux>1) {
01115 q(2) = a1;
01116 q(3) = a2;
01117 }
01118 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
01119 }
01120 break;
01121 case base::Bottom:
01122 for (register int k=mzs; k<=mze; k++)
01123 for (register int i=mxs; i<=mxe; i++) {
01124 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01125 if (naux==1)
01126 q(2) = a0;
01127 else {
01128 q(1) = a0;
01129 q(2) = a1;
01130 q(3) = a2;
01131 }
01132 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
01133 }
01134 break;
01135 case base::Top:
01136 for (register int k=mzs; k<=mze; k++)
01137 for (register int i=mxs; i<=mxe; i++) {
01138 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01139 if (naux==1)
01140 q(2) = a0;
01141 else {
01142 q(1) = a0;
01143 q(2) = a1;
01144 q(3) = a2;
01145 }
01146 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
01147 }
01148 break;
01149 case base::Front:
01150 for (register int j=mys; j<=mye; j++)
01151 for (register int i=mxs; i<=mxe; i++) {
01152 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01153 if (naux==1)
01154 q(3) = a0;
01155 else {
01156 q(1) = a0;
01157 q(2) = a1;
01158 q(3) = a2;
01159 }
01160 f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
01161 }
01162 break;
01163 case base::Back:
01164 for (register int j=mys; j<=mye; j++)
01165 for (register int i=mxs; i<=mxe; i++) {
01166 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01167 if (naux==1)
01168 q(3) = a0;
01169 else {
01170 q(1) = a0;
01171 q(2) = a1;
01172 q(3) = a2;
01173 }
01174 f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
01175 }
01176 break;
01177 }
01178 }
01179
01180 else if (type==Outlet) {
01181 switch (side) {
01182 case base::Left:
01183 for (register int k=mzs; k<=mze; k++)
01184 for (register int j=mys; j<=mye; j++) {
01185 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01186 if (q(0)>0) f[b+base::idx(mxe,j,k,Nx,Ny)] = f[b+base::idx(mxe+1,j,k,Nx,Ny)]/q(0);
01187 }
01188 break;
01189 case base::Right:
01190 for (register int k=mzs; k<=mze; k++)
01191 for (register int j=mys; j<=mye; j++) {
01192 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01193 if (q(0)>0) f[b+base::idx(mxs,j,k,Nx,Ny)] = f[b+base::idx(mxs-1,j,k,Nx,Ny)]/q(0);
01194 }
01195 break;
01196 case base::Bottom:
01197 for (register int k=mzs; k<=mze; k++)
01198 for (register int i=mxs; i<=mxe; i++) {
01199 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01200 if (q(0)>0) f[b+base::idx(i,mye,k,Nx,Ny)] = f[b+base::idx(i,mye+1,k,Nx,Ny)]/q(0);
01201 }
01202 break;
01203 case base::Top:
01204 for (register int k=mzs; k<=mze; k++)
01205 for (register int i=mxs; i<=mxe; i++) {
01206 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01207 if (q(0)>0) f[b+base::idx(i,mys,k,Nx,Ny)] = f[b+base::idx(i,mys-1,k,Nx,Ny)]/q(0);
01208 }
01209 break;
01210 case base::Front:
01211 for (register int j=mys; j<=mye; j++)
01212 for (register int i=mxs; i<=mxe; i++) {
01213 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01214 if (q(0)>0) f[b+base::idx(i,j,mze,Nx,Ny)] = f[b+base::idx(i,j,mze+1,Nx,Ny)]/q(0);
01215 }
01216 break;
01217 case base::Back:
01218 for (register int j=mys; j<=mye; j++)
01219 for (register int i=mxs; i<=mxe; i++) {
01220 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01221 if (q(0)>0) f[b+base::idx(i,j,mzs,Nx,Ny)] = f[b+base::idx(i,j,mzs-1,Nx,Ny)]/q(0);
01222 }
01223 break;
01224 }
01225 }
01226
01227
01228 else if (type==Pressure) {
01229 if (aux==0 || naux<=0) return;
01230 DataType a0=aux[0];
01231 if (scaling==base::Physical)
01232 a0 /= R0;
01233 switch (side) {
01234 case base::Left:
01235 for (register int k=mzs; k<=mze; k++)
01236 for (register int j=mys; j<=mye; j++) {
01237 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01238 q(0) = a0;
01239 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
01240 }
01241 break;
01242 case base::Right:
01243 for (register int k=mzs; k<=mze; k++)
01244 for (register int j=mys; j<=mye; j++) {
01245 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01246 q(0) = a0;
01247 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
01248 }
01249 break;
01250 case base::Bottom:
01251 for (register int k=mzs; k<=mze; k++)
01252 for (register int i=mxs; i<=mxe; i++) {
01253 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01254 q(0) = a0;
01255 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
01256 }
01257 break;
01258 case base::Top:
01259 for (register int k=mzs; k<=mze; k++)
01260 for (register int i=mxs; i<=mxe; i++) {
01261 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01262 q(0) = a0;
01263 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
01264 }
01265 break;
01266 case base::Front:
01267 for (register int j=mys; j<=mye; j++)
01268 for (register int i=mxs; i<=mxe; i++) {
01269 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01270 q(0) = a0;
01271 f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
01272 }
01273 break;
01274 case base::Back:
01275 for (register int j=mys; j<=mye; j++)
01276 for (register int i=mxs; i<=mxe; i++) {
01277 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01278 q(0) = a0;
01279 f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
01280 }
01281 break;
01282 }
01283 }
01284
01285 else if (type==SlidingWall) {
01286 if (aux==0 || naux<=0) return;
01287 DataType a0=S0*aux[0];
01288 if (scaling==base::Physical)
01289 a0 /= U0;
01290 switch (side) {
01291 case base::Left:
01292 for (register int k=mzs; k<=mze; k++)
01293 for (register int j=mys; j<=mye; j++) {
01294 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01295 DataType rd1 = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10), rd2 = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
01296 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01297 DataType pw = DataType(1.)-qw;
01298 if (j<mye) f[b+base::idx(mxe,j+1,k,Nx,Ny)](9) = pw*rd1+qw*rd2;
01299 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01300 if (j>mys) f[b+base::idx(mxe,j-1,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01301
01302 if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01303 if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01304 }
01305 break;
01306 case base::Right:
01307 for (register int k=mzs; k<=mze; k++)
01308 for (register int j=mys; j<=mye; j++) {
01309 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01310 DataType rd1 = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7), rd2 = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
01311 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01312 DataType pw = DataType(1.)-qw;
01313 if (j<mye) f[b+base::idx(mxs,j+1,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01314 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01315 if (j>mys) f[b+base::idx(mxs,j-1,k,Nx,Ny)](10) = qw*rd1+pw*rd2;
01316
01317 if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01318 if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01319 }
01320 break;
01321 case base::Bottom:
01322 for (register int k=mzs; k<=mze; k++)
01323 for (register int i=mxs; i<=mxe; i++) {
01324 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01325 DataType rd1 = f[b+base::idx(i,mye+1,k,Nx,Ny)](9), rd2 = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
01326 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01327 DataType pw = DataType(1.)-qw;
01328 if (i<mxe) f[b+base::idx(i+1,mye,k,Nx,Ny)](10) = pw*rd1+qw*rd2;
01329 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01330 if (i>mxs) f[b+base::idx(i-1,mye,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01331
01332 if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01333 if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01334 }
01335 break;
01336 case base::Top:
01337 for (register int k=mzs; k<=mze; k++)
01338 for (register int i=mxs; i<=mxe; i++) {
01339 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01340 DataType rd1 = f[b+base::idx(i,mys-1,k,Nx,Ny)](7), rd2 = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
01341 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01342 DataType pw = DataType(1.)-qw;
01343 if (i<mxe) f[b+base::idx(i+1,mys,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01344 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01345 if (i>mxs) f[b+base::idx(i-1,mys,k,Nx,Ny)](9) = qw*rd1+pw*rd2;
01346
01347 if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01348 if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01349 }
01350 break;
01351 case base::Front:
01352 for (register int j=mys; j<=mye; j++)
01353 for (register int i=mxs; i<=mxe; i++) {
01354 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01355 DataType rd1 = f[b+base::idx(i,j,mze+1,Nx,Ny)](7), rd2 = f[b+base::idx(i,j,mze+1,Nx,Ny)](10);
01356 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01357 DataType pw = DataType(1.)-qw;
01358 if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = pw*rd1+qw*rd2;
01359 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01360 if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = qw*rd1+pw*rd2;
01361
01362 if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01363 if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01364 }
01365 break;
01366 case base::Back:
01367 for (register int j=mys; j<=mye; j++)
01368 for (register int i=mxs; i<=mxe; i++) {
01369 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01370 DataType rd1 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](7), rd2 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](10);
01371 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01372 DataType pw = DataType(1.)-qw;
01373 if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = pw*rd1+qw*rd2;
01374 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
01375 if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = qw*rd1+pw*rd2;
01376
01377 if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
01378 if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
01379 }
01380 break;
01381 }
01382 }
01383
01387 else if (type==CharacteristicOutlet) {
01388 DataType rhod=DataType(1.0);
01389 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
01390 DCoords dx = base::GH().worldStep(fvec.stepsize());
01391 DataType dt_temp = dx(0)/VelocityScale();
01392 switch (side) {
01393 case base::Left:
01394 for (register int k=mzs; k<=mze; k++)
01395 for (register int j=mys; j<=mye; j++) {
01396 MicroType ftmp = f[b+base::idx(mxe,j,k,Nx,Ny)];
01397 MacroType q = MacroVariables(ftmp);
01398 if (turbulence!=CSM)
01399 Collision(ftmp,dt_temp);
01400 else {
01401 MacroType qxp=0., qxm=0., qyp=0., qym=0., qzp=0., qzm=0.;
01402 qxp=MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01403 qxm=MacroVariables(f[b+base::idx(mxe-1,j,k,Nx,Ny)]);
01404 if (j<mye) qyp=MacroVariables(f[b+base::idx(mxe,j+1,k,Nx,Ny)]);
01405 else qyp=q;
01406 if (j>mys) qym=MacroVariables(f[b+base::idx(mxe,j-1,k,Nx,Ny)]);
01407 else qym=q;
01408 if (k<mze) qzp=MacroVariables(f[b+base::idx(mxe,j,k+1,Nx,Ny)]);
01409 else qzp=q;
01410 if (k>mze) qzm=MacroVariables(f[b+base::idx(mxe,j,k-1,Nx,Ny)]);
01411 else qzm=q;
01412 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,qzp,qzm,dx,dt_temp);
01413 }
01414 MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01415 MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,k,Nx,Ny)]);
01416 MacroType dqdn = NormalDerivative(q,q1,q2);
01417 if (EquilibriumType()!=3)
01418 rhod = q(0);
01419 Vector<DataType,3> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01420
01421
01422 MacroType qn;
01423 qn(0) = q(0) - c1oCsSq2*L(0);
01424 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01425 qn(2) = q(2) - L(1);
01426 qn(3) = q(3) - L(2);
01427
01428 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(qn);
01429 }
01430 break;
01431 case base::Right:
01432 for (register int k=mzs; k<=mze; k++)
01433 for (register int j=mys; j<=mye; j++) {
01434 MicroType ftmp = f[b+base::idx(mxs,j,k,Nx,Ny)];
01435 MacroType q = MacroVariables(ftmp);
01436 if (turbulence!=CSM)
01437 Collision(ftmp,dt_temp);
01438 else {
01439 MacroType qxp=0., qxm=0., qyp=0., qym=0., qzp=0., qzm=0.;
01440 qxp=MacroVariables(f[b+base::idx(mxs+1,j,k,Nx,Ny)]);
01441 qxm=MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01442 if (j<mye) qyp=MacroVariables(f[b+base::idx(mxs,j+1,k,Nx,Ny)]);
01443 else qyp=q;
01444 if (j>mys) qym=MacroVariables(f[b+base::idx(mxs,j-1,k,Nx,Ny)]);
01445 else qym=q;
01446 if (k<mze) qzp=MacroVariables(f[b+base::idx(mxs,j,k+1,Nx,Ny)]);
01447 else qzp=q;
01448 if (k>mzs) qzm=MacroVariables(f[b+base::idx(mxs,j,k-1,Nx,Ny)]);
01449 else qzm=q;
01450 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,qzp,qzm,dx,dt_temp);
01451 }
01452 MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01453 MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,k,Nx,Ny)]);
01454 MacroType dqdn = NormalDerivative(q,q1,q2);
01455 if (EquilibriumType()!=3)
01456 rhod = q(0);
01457 Vector<DataType,3> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01458
01459
01460 MacroType qn;
01461 qn(0) = q(0) - c1oCsSq2*L(0);
01462 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01463 qn(2) = q(2) - L(1);
01464 qn(3) = q(3) - L(2);
01465
01466 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(qn);
01467 }
01468 break;
01469 case base::Bottom:
01470 for (register int k=mzs; k<=mze; k++)
01471 for (register int i=mxs; i<=mxe; i++) {
01472 MicroType ftmp = f[b+base::idx(i,mye,k,Nx,Ny)];
01473 MacroType q = MacroVariables(ftmp);
01474 if (turbulence!=CSM)
01475 Collision(ftmp,dt_temp);
01476 else {
01477 MacroType qxp=0., qxm=0., qyp=0., qym=0., qzp=0., qzm=0.;
01478 if (i<mxe) qxp=MacroVariables(f[b+base::idx(i+1,mye,k,Nx,Ny)]);
01479 else qxp=q;
01480 if (i>mxs) qxm=MacroVariables(f[b+base::idx(i-1,mye,k,Nx,Ny)]);
01481 else qxm=q;
01482 qyp=MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01483 qym=MacroVariables(f[b+base::idx(i,mye-1,k,Nx,Ny)]);
01484 if (k<mze) qzp=MacroVariables(f[b+base::idx(i,mye,k+1,Nx,Ny)]);
01485 else qzp=q;
01486 if (k>mzs) qzm=MacroVariables(f[b+base::idx(i,mye,k-1,Nx,Ny)]);
01487 else qzm=q;
01488 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,qzp,qzm,dx,dt_temp);
01489 }
01490 MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01491 MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,k,Nx,Ny)]);
01492 MacroType dqdn = NormalDerivative(q,q1,q2);
01493 if (EquilibriumType()!=3)
01494 rhod = q(0);
01495 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01496
01497
01498 MacroType qn;
01499 qn(0) = q(0) - c1oCsSq2*L(0);
01500 qn(1) = q(1) - L(1);
01501 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01502 qn(3) = q(3) - L(2);
01503
01504 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(qn);
01505 }
01506 break;
01507 case base::Top:
01508 for (register int k=mzs; k<=mze; k++)
01509 for (register int i=mxs; i<=mxe; i++) {
01510 MicroType ftmp = f[b+base::idx(i,mys,k,Nx,Ny)];
01511 MacroType q = MacroVariables(ftmp);
01512 if (turbulence!=CSM)
01513 Collision(ftmp,dt_temp);
01514 else {
01515 MacroType qxp=0., qxm=0., qyp=0., qym=0., qzp=0., qzm=0.;
01516 if (i<mxe) qxp=MacroVariables(f[b+base::idx(i+1,mys,k,Nx,Ny)]);
01517 else qxp=q;
01518 if (i<mxe) qxm=MacroVariables(f[b+base::idx(i-1,mys,k,Nx,Ny)]);
01519 else qxm=q;
01520 qyp=MacroVariables(f[b+base::idx(i,mys+1,k,Nx,Ny)]);
01521 qym=MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01522 if (k<mze) qzp=MacroVariables(f[b+base::idx(i,mys,k+1,Nx,Ny)]);
01523 else qzp=q;
01524 if (k>mxs) qzm=MacroVariables(f[b+base::idx(i,mys,k-1,Nx,Ny)]);
01525 else qzm=q;
01526 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,qzp,qzm,dx,dt_temp);
01527 }
01528 MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01529 MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,k,Nx,Ny)]);
01530 MacroType dqdn = NormalDerivative(q,q1,q2);
01531 if (EquilibriumType()!=3)
01532 rhod = q(0);
01533 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01534
01535
01536 MacroType qn;
01537 qn(0) = q(0) - c1oCsSq2*L(0);
01538 qn(1) = q(1) - L(1);
01539 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01540 qn(3) = q(3) - L(2);
01541
01542 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(qn);
01543 }
01544 break;
01545 case base::Front:
01546 for (register int j=mys; j<=mye; j++)
01547 for (register int i=mxs; i<=mxe; i++) {
01548 MicroType ftmp = f[b+base::idx(i,j,mze,Nx,Ny)];
01549 MacroType q = MacroVariables(ftmp);
01550 if (turbulence!=CSM)
01551 Collision(ftmp,dt_temp);
01552 else {
01553 MacroType qxp=0., qxm=0., qyp=0., qym=0., qzp=0., qzm=0.;
01554 if (i<mxe) qxp=MacroVariables(f[b+base::idx(i+1,j,mze,Nx,Ny)]);
01555 else qxp=q;
01556 if (i>mxs) qxm=MacroVariables(f[b+base::idx(i-1,j,mze,Nx,Ny)]);
01557 else qxm=q;
01558 if (j<mye) qyp=MacroVariables(f[b+base::idx(i,j+1,mze,Nx,Ny)]);
01559 else qyp=q;
01560 if (j>mys) qym=MacroVariables(f[b+base::idx(i,j-1,mze,Nx,Ny)]);
01561 else qxm=q;
01562 qzp=MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01563 qzm=MacroVariables(f[b+base::idx(i,j,mze-1,Nx,Ny)]);
01564 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,qzp,qzm,dx,dt_temp);
01565 }
01566 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01567 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mze+2,Nx,Ny)]);
01568 MacroType dqdn = NormalDerivative(q,q1,q2);
01569 if (EquilibriumType()!=3)
01570 rhod = q(0);
01571 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01572
01573
01574 MacroType qn;
01575 qn(0) = q(0) - c1oCsSq2*L(0);
01576 qn(1) = q(1) - L(1);
01577 qn(2) = q(2) - L(2);
01578 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01579
01580 f[b+base::idx(i,j,mze,Nx,Ny)] = Equilibrium(qn);
01581 }
01582 break;
01583 case base::Back:
01584 for (register int j=mys; j<=mye; j++)
01585 for (register int i=mxs; i<=mxe; i++) {
01586 MicroType ftmp = f[b+base::idx(i,j,mzs,Nx,Ny)];
01587 MacroType q = MacroVariables(ftmp);
01588 if (turbulence!=CSM)
01589 Collision(ftmp,dt_temp);
01590 else {
01591 MacroType qxp=0., qxm=0., qyp=0., qym=0., qzp=0., qzm=0.;
01592 if (i<mxe) qxp=MacroVariables(f[b+base::idx(i+1,j,mzs,Nx,Ny)]);
01593 else qxp=q;
01594 if (i>mxs) qxm=MacroVariables(f[b+base::idx(i-1,j,mzs,Nx,Ny)]);
01595 else qxm=q;
01596 if (j<mye) qyp=MacroVariables(f[b+base::idx(i,j+1,mzs,Nx,Ny)]);
01597 else qyp=q;
01598 if (j>mys) qym=MacroVariables(f[b+base::idx(i,j-1,mzs,Nx,Ny)]);
01599 else qym=q;
01600 qzp=MacroVariables(f[b+base::idx(i,j,mzs+1,Nx,Ny)]);
01601 qzm=MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01602 LocalCollisionCSM(ftmp,qxp,qxm,qyp,qym,qzp,qzm,dx,dt_temp);
01603 }
01604 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01605 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mzs-2,Nx,Ny)]);
01606 MacroType dqdn = NormalDerivative(q,q1,q2);
01607 if (EquilibriumType()!=3)
01608 rhod = q(0);
01609 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01610
01611
01612 MacroType qn;
01613 qn(0) = q(0) - c1oCsSq2*L(0);
01614 qn(1) = q(1) - L(1);
01615 qn(2) = q(2) - L(2);
01616 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01617
01618 f[b+base::idx(i,j,mzs,Nx,Ny)] = Equilibrium(qn);
01619 }
01620 break;
01621 }
01622 }
01623
01627 else if (type==CharacteristicInlet) {
01628 if (aux==0 || naux<=0) return;
01629 DataType rho=aux[0];
01630 DataType a0=S0*aux[1];
01631 DataType a1=S0*aux[2];
01632 DataType a2=S0*aux[3];
01633 if (scaling==base::Physical) {
01634 rho /= R0;
01635 a0 /= U0;
01636 a1 /= U0;
01637 a2 /= U0;
01638 }
01639 MacroType q, qn;
01640 q(0) = rho;
01641 q(1) = a0;
01642 q(2) = a1;
01643 q(3) = a2;
01644 DataType rhod=DataType(1.0);
01645 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
01646 switch (side) {
01647 case base::Left:
01648 for (register int k=mzs; k<=mze; k++)
01649 for (register int j=mys; j<=mye; j++) {
01650 MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01651 MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,k,Nx,Ny)]);
01652 Vector<DataType,3> L;
01653 if (naux==1) {
01654 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01655 MacroType dqdn = NormalDerivative(q,q1,q2);
01656 if (EquilibriumType()!=3)
01657 rhod = q(0);
01658 L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01659 }
01660 else if (naux==2) {
01661 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01662 MacroType dqdn = NormalDerivative(q,q1,q2);
01663 if (EquilibriumType()!=3)
01664 rhod = q(0);
01665 L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01666 }
01667 else if (naux==4) {
01668 q(0) = q1(0);
01669 MacroType dqdn = NormalDerivative(q,q1,q2);
01670 if (EquilibriumType()!=3)
01671 rhod = q(0);
01672 L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01673 }
01674 MacroType qn;
01675 qn(0) = q(0) - c1oCsSq2*L(0);
01676 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01677 qn(2) = q(2) - L(1);
01678 qn(3) = q(3) - L(2);
01679 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(qn);
01680 }
01681 break;
01682 case base::Right:
01683 for (register int k=mzs; k<=mze; k++)
01684 for (register int j=mys; j<=mye; j++) {
01685 MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01686 MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,k,Nx,Ny)]);
01687 if (naux==1) {
01688 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01689 MacroType dqdn = NormalDerivative(q,q1,q2);
01690 if (EquilibriumType()!=3)
01691 rhod = q(0);
01692 Vector<DataType,3> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01693 MacroType qn;
01694 qn(0) = q(0) - c1oCsSq2*L(0);
01695 qn(1) = q1(1) + c1o2cs*L(0)/rhod;
01696 qn(2) = q1(2) - L(1);
01697 qn(3) = q1(3) - L(2);
01698 }
01699 else if (naux==2) {
01700 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01701 MacroType dqdn = NormalDerivative(q,q1,q2);
01702 if (EquilibriumType()!=3)
01703 rhod = q(0);
01704 Vector<DataType,3> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01705 MacroType qn;
01706 qn(0) = q1(0) - c1oCsSq2*L(0);
01707 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01708 qn(2) = q1(2) - L(1);
01709 qn(3) = q1(3) - L(2);
01710 }
01711 else if (naux==4) {
01712 q(0) = q1(0);
01713 MacroType dqdn = NormalDerivative(q,q1,q2);
01714 if (EquilibriumType()!=3)
01715 rhod = q(0);
01716 Vector<DataType,3> L = WaveAmplitudes(rhod,q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01717 MacroType qn;
01718 qn(0) = q1(0) - c1oCsSq2*L(0);
01719 qn(1) = q(1) + c1o2cs*L(0)/rhod;
01720 qn(2) = q(2) - L(1);
01721 qn(3) = q(3) - L(2);
01722 }
01723 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(qn);
01724 }
01725 break;
01726 case base::Bottom:
01727 for (register int k=mzs; k<=mze; k++)
01728 for (register int i=mxs; i<=mxe; i++) {
01729 MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01730 MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,k,Nx,Ny)]);
01731 if (naux==1) {
01732 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01733 MacroType dqdn = NormalDerivative(q,q1,q2);
01734 if (EquilibriumType()!=3)
01735 rhod = q(0);
01736 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01737 MacroType qn;
01738 qn(0) = q(0) - c1oCsSq2*L(0);
01739 qn(1) = q1(1) - L(1);
01740 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
01741 qn(3) = q1(3) - L(2);
01742 }
01743 else if (naux==2) {
01744 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01745 MacroType dqdn = NormalDerivative(q,q1,q2);
01746 if (EquilibriumType()!=3)
01747 rhod = q(0);
01748 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01749 MacroType qn;
01750 qn(0) = q1(0) - c1oCsSq2*L(0);
01751 qn(1) = q(1) - L(1);
01752 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
01753 qn(3) = q1(3) - L(2);
01754 }
01755 else if (naux==4) {
01756 q(0) = q1(0);
01757 MacroType dqdn = NormalDerivative(q,q1,q2);
01758 if (EquilibriumType()!=3)
01759 rhod = q(0);
01760 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01761 MacroType qn;
01762 qn(0) = q1(0) - c1oCsSq2*L(0);
01763 qn(1) = q(1) - L(1);
01764 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01765 qn(3) = q(3) - L(2);
01766 }
01767 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(qn);
01768 }
01769 break;
01770 case base::Top:
01771 for (register int k=mzs; k<=mze; k++)
01772 for (register int i=mxs; i<=mxe; i++) {
01773 MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01774 MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,k,Nx,Ny)]);
01775 if (naux==1) {
01776 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01777 MacroType dqdn = NormalDerivative(q,q1,q2);
01778 if (EquilibriumType()!=3)
01779 rhod = q(0);
01780 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01781 MacroType qn;
01782 qn(0) = q(0) - c1oCsSq2*L(0);
01783 qn(1) = q1(1) - L(1);
01784 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
01785 qn(3) = q1(3) - L(2);
01786 }
01787 else if (naux==2) {
01788 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01789 MacroType dqdn = NormalDerivative(q,q1,q2);
01790 if (EquilibriumType()!=3)
01791 rhod = q(0);
01792 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01793 MacroType qn;
01794 qn(0) = q1(0) - c1oCsSq2*L(0);
01795 qn(1) = q(1) - L(1);
01796 qn(2) = q1(2) + c1o2cs*L(0)/rhod;
01797 qn(3) = q1(3) - L(2);
01798 }
01799 else if (naux==4) {
01800 q(0) = q1(0);
01801 MacroType dqdn = NormalDerivative(q,q1,q2);
01802 if (EquilibriumType()!=3)
01803 rhod = q(0);
01804 Vector<DataType,3> L = WaveAmplitudes(rhod,q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01805 MacroType qn;
01806 qn(0) = q1(0) - c1oCsSq2*L(0);
01807 qn(1) = q(1) - L(1);
01808 qn(2) = q(2) + c1o2cs*L(0)/rhod;
01809 qn(3) = q(3) - L(2);
01810 }
01811 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(qn);
01812 }
01813 break;
01814 case base::Front:
01815 for (register int j=mys; j<=mye; j++)
01816 for (register int i=mxs; i<=mxe; i++) {
01817 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01818 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mze+2,Nx,Ny)]);
01819 if (naux==1) {
01820 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01821 MacroType dqdn = NormalDerivative(q,q1,q2);
01822 if (EquilibriumType()!=3)
01823 rhod = q(0);
01824 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01825 MacroType qn;
01826 qn(0) = q(0) - c1oCsSq2*L(0);
01827 qn(1) = q(1) - L(1);
01828 qn(2) = q(2) - L(2);
01829 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01830 }
01831 else if (naux==2) {
01832 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01833 MacroType dqdn = NormalDerivative(q,q1,q2);
01834 if (EquilibriumType()!=3)
01835 rhod = q(0);
01836 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01837 MacroType qn;
01838 qn(0) = q(0) - c1oCsSq2*L(0);
01839 qn(1) = q(1) - L(1);
01840 qn(2) = q(2) - L(2);
01841 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01842 }
01843 else if (naux==4) {
01844 q(0) = q1(0);
01845 MacroType dqdn = NormalDerivative(q,q1,q2);
01846 if (EquilibriumType()!=3)
01847 rhod = q(0);
01848 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01849 MacroType qn;
01850 qn(0) = q(0) - c1oCsSq2*L(0);
01851 qn(1) = q(1) - L(1);
01852 qn(2) = q(2) - L(2);
01853 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01854 }
01855 MacroType dqdn = NormalDerivative(q,q1,q2);
01856 if (EquilibriumType()!=3)
01857 rhod = q(0);
01858 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01859
01860
01861 MacroType qn;
01862 qn(0) = q(0) - c1oCsSq2*L(0);
01863 qn(1) = q(1) - L(1);
01864 qn(2) = q(2) - L(2);
01865 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01866
01867 f[b+base::idx(i,j,mze,Nx,Ny)] = Equilibrium(qn);
01868 }
01869 break;
01870 case base::Back:
01871 for (register int j=mys; j<=mye; j++)
01872 for (register int i=mxs; i<=mxe; i++) {
01873 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01874 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mzs-2,Nx,Ny)]);
01875 if (naux==1) {
01876 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01877 MacroType dqdn = NormalDerivative(q,q1,q2);
01878 if (EquilibriumType()!=3)
01879 rhod = q(0);
01880 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01881 MacroType qn;
01882 qn(0) = q(0) - c1oCsSq2*L(0);
01883 qn(1) = q(1) - L(1);
01884 qn(2) = q(2) - L(2);
01885 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01886 }
01887 else if (naux==2) {
01888 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01889 MacroType dqdn = NormalDerivative(q,q1,q2);
01890 if (EquilibriumType()!=3)
01891 rhod = q(0);
01892 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01893 MacroType qn;
01894 qn(0) = q(0) - c1oCsSq2*L(0);
01895 qn(1) = q(1) - L(1);
01896 qn(2) = q(2) - L(2);
01897 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01898 }
01899 else if (naux==4) {
01900 q(0) = q1(0);
01901 MacroType dqdn = NormalDerivative(q,q1,q2);
01902 if (EquilibriumType()!=3)
01903 rhod = q(0);
01904 Vector<DataType,3> L = WaveAmplitudes(rhod,q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01905 MacroType qn;
01906 qn(0) = q(0) - c1oCsSq2*L(0);
01907 qn(1) = q(1) - L(1);
01908 qn(2) = q(2) - L(2);
01909 qn(3) = q(3) + c1o2cs*L(0)/rhod;
01910 }
01911 f[b+base::idx(i,j,mzs,Nx,Ny)] = Equilibrium(qn);
01912 }
01913 break;
01914 }
01915 }
01916
01917 else if (type==NoSlipWallEntranceExit) {
01918 int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), mzs*fvec.stepsize(2) };
01919 DataType WLB[3], WUB[3];
01920 base::GH().wlb(WLB);
01921 base::GH().wub(WUB);
01922 DCoords wc;
01923 DataType slip = 0.;
01924 DataType lxs, lxe, lys, lye, lzs, lze;
01925 DataType min[3], max[3];
01926 if (naux!=6) assert(0);
01927 lxs = aux[0] - WLB[0];
01928 lxe = WUB[0] - aux[1];
01929 lys = aux[2] - WLB[1];
01930 lye = WUB[1] - aux[3];
01931 lzs = aux[4] - WLB[2];
01932 lze = WUB[2] - aux[5];
01933 min[0] = aux[0]; max[0] = aux[1];
01934 min[1] = aux[2]; max[1] = aux[3];
01935 min[2] = aux[4]; max[2] = aux[5];
01936
01937 switch (side) {
01938 case base::Left:
01939 lc_indx[0] = mxe*fvec.stepsize(0);
01940 for (register int k=mzs; k<=mze; k++) {
01941 lc_indx[2] = k*fvec.stepsize(2);
01942 for (register int j=mys; j<=mye; j++) {
01943 lc_indx[1] = j*fvec.stepsize(1);
01944 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01945 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01946 if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01947 if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01948 }
01949 else {
01950 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01951 slip = DataType(0.);
01952 else if (wc(1)<min[1])
01953 slip = (wc(1) - WLB[1])/lys;
01954 else if (wc(1)>max[1])
01955 slip = (WUB[1] - wc(1))/lye;
01956 if (j>mys && j<mye) {
01957 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10)
01958 + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01959 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8)
01960 + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01961 }
01962 else if (j==mys) {
01963 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](10)
01964 + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01965 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8)
01966 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
01967 }
01968 else if (j==mye) {
01969 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10)
01970 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
01971 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](8)
01972 + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01973 }
01974 }
01975 if (wc(2)>=min[2] && wc(2)<=max[2]) {
01976 if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01977 if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01978 }
01979 else {
01980 if (wc(2)<WLB[2] || wc(2) >WUB[2])
01981 slip = DataType(0.);
01982 else if (wc(2)<min[2])
01983 slip = (wc(2) - WLB[2])/lzs;
01984 else if (wc(2)>max[2])
01985 slip = (WUB[2] - wc(2))/lze;
01986 if (k>mzs && k<mze) {
01987 f[b+base::idx(mxe,j,k,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18)
01988 + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01989 f[b+base::idx(mxe,j,k,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16)
01990 + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01991 }
01992 else if (k==mzs) {
01993 f[b+base::idx(mxe,j,k,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](18)
01994 + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01995 f[b+base::idx(mxe,j,k,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16)
01996 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
01997 }
01998 else if (k==mze) {
01999 f[b+base::idx(mxe,j,k,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18)
02000 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
02001 f[b+base::idx(mxe,j,k,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(mxe+1,j,k,Nx,Ny)](16)
02002 + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
02003 }
02004 }
02005 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
02006 }
02007 }
02008 break;
02009 case base::Right:
02010 lc_indx[0] = mxs*fvec.stepsize(0);
02011 for (register int k=mzs; k<=mze; k++) {
02012 lc_indx[2] = k*fvec.stepsize(2);
02013 for (register int j=mys; j<=mye; j++) {
02014 lc_indx[1] = j*fvec.stepsize(1);
02015 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
02016 if (wc(1)>=min[1] && wc(1)<=max[1]) {
02017 if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
02018 if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
02019 }
02020
02021 else {
02022 if (wc(1)<WLB[1] || wc(1) >WUB[1])
02023 slip = DataType(0.);
02024 else if (wc(1)<min[1])
02025 slip = (wc(1) - WLB[1])/lys;
02026 else if (wc(1)>max[1])
02027 slip = (WUB[1] - wc(1))/lye;
02028 if (j>mys && j<mye) {
02029 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7)
02030 + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
02031 f[b+base::idx(mxs,j,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9)
02032 + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
02033 }
02034 else if (j==mys) {
02035 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](7)
02036 + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
02037 f[b+base::idx(mxs,j,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9)
02038 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
02039 }
02040 else if (j==mye) {
02041 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7)
02042 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
02043 f[b+base::idx(mxs,j,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](9)
02044 + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
02045 }
02046 }
02047 if (wc(2)>=min[2] && wc(2)<=max[2]) {
02048 if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
02049 if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
02050 }
02051 else {
02052 if (wc(2)<WLB[2] || wc(2) >WUB[2])
02053 slip = DataType(0.);
02054 else if (wc(2)<min[2])
02055 slip = (wc(2) - WLB[2])/lzs;
02056 else if (wc(2)>max[2])
02057 slip = (WUB[2] - wc(2))/lze;
02058 if (k>mzs && k<mze) {
02059 f[b+base::idx(mxs,j,k,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15)
02060 + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
02061 f[b+base::idx(mxs,j,k,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17)
02062 + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
02063 }
02064 else if (k==mzs) {
02065 f[b+base::idx(mxs,j,k,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](15)
02066 + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
02067 f[b+base::idx(mxs,j,k,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17)
02068 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
02069 }
02070 else if (k==mze) {
02071 f[b+base::idx(mxs,j,k,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15)
02072 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
02073 f[b+base::idx(mxs,j,k,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(mxs-1,j,k,Nx,Ny)](17)
02074 + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
02075 }
02076 }
02077 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
02078 }
02079 }
02080 break;
02081 case base::Bottom:
02082 lc_indx[1] = mye*fvec.stepsize(1);
02083 for (register int k=mzs; k<=mze; k++) {
02084 lc_indx[2] = k*fvec.stepsize(2);
02085 for (register int i=mxs; i<=mxe; i++) {
02086 lc_indx[0] = i*fvec.stepsize(0);
02087 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
02088 if (wc(0)>=min[0] && wc(0)<=max[0]) {
02089 if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
02090 if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
02091 }
02092 else {
02093 if (wc(0)<WLB[0] || wc(0)>WUB[0])
02094 slip = DataType(0.);
02095 else if (wc(0)<min[0])
02096 slip = (wc(0) - WLB[0])/lxs;
02097 else if (wc(0)>max[0])
02098 slip = (WUB[0] - wc(0))/lxe;
02099 if (i>mxs && i<mxe) {
02100 f[b+base::idx(i,mye,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9)
02101 + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
02102 f[b+base::idx(i,mye,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8)
02103 + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
02104 }
02105 else if (i==mxs) {
02106 f[b+base::idx(i,mye,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](9)
02107 + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
02108 f[b+base::idx(i,mye,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8)
02109 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
02110 }
02111 else if (i==mxe) {
02112 f[b+base::idx(i,mye,k,Nx,Ny)](10) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9)
02113 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
02114 f[b+base::idx(i,mye,k,Nx,Ny)](7) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](8)
02115 + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
02116 }
02117 }
02118 if (wc(2)>=min[2] && wc(2)<=max[2]) {
02119 if (k>mzs) f[b+base::idx(i,mye,k,Nx,Ny)](13) = f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
02120 if (k<mze) f[b+base::idx(i,mye,k,Nx,Ny)](11) = f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
02121 }
02122 else {
02123 if (wc(2)<WLB[2] || wc(2) >WUB[2])
02124 slip = DataType(0.);
02125 else if (wc(2)<min[2])
02126 slip = (wc(2) - WLB[2])/lzs;
02127 else if (wc(2)>max[2])
02128 slip = (WUB[2] - wc(2))/lze;
02129 if (k>mzs && k<mze) {
02130 f[b+base::idx(i,mye,k,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14)
02131 + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
02132 f[b+base::idx(i,mye,k,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12)
02133 + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
02134 }
02135 else if (k==mzs) {
02136 f[b+base::idx(i,mye,k,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](14)
02137 + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
02138 f[b+base::idx(i,mye,k,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12)
02139 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](14);
02140 }
02141 else if (k==mze) {
02142 f[b+base::idx(i,mye,k,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14)
02143 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](12);
02144 f[b+base::idx(i,mye,k,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,mye+1,k,Nx,Ny)](12)
02145 + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
02146 }
02147 }
02148 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
02149 }
02150 }
02151 break;
02152 case base::Top:
02153 lc_indx[1] = mys*fvec.stepsize(1);
02154 for (register int k=mzs; k<=mze; k++) {
02155 lc_indx[2] = k*fvec.stepsize(2);
02156 for (register int i=mxs; i<=mxe; i++) {
02157 lc_indx[0] = i*fvec.stepsize(0);
02158 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
02159 if (wc(0)>=min[0] && wc(0)<=max[0]) {
02160 if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
02161 if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
02162 }
02163 else {
02164 if (wc(0)<WLB[0] || wc(0) >WUB[0])
02165 slip = DataType(0.);
02166 else if (wc(0)<min[0])
02167 slip = (wc(0) - WLB[0])/lxs;
02168 else if (wc(0)>max[0])
02169 slip = (WUB[0] - wc(0))/lxe;
02170 if (i>mxs && i<mxe) {
02171 f[b+base::idx(i,mys,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7)
02172 + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
02173 f[b+base::idx(i,mys,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10)
02174 + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
02175 }
02176 else if (i==mxs) {
02177 f[b+base::idx(i,mys,k,Nx,Ny)](8) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k,Nx,Ny)](7)
02178 + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
02179 f[b+base::idx(i,mys,k,Nx,Ny)](9) = (DataType(1.0)-slip)*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10)
02180 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
02181 }
02182 else if (i==mxe) {
02183 f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7)
02184 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
02185 f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10)
02186 + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
02187 }
02188 }
02189 if (wc(2)>=min[2] && wc(2)<=max[2]) {
02190 if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
02191 if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
02192 }
02193 else {
02194 if (wc(2)<WLB[2] || wc(2)>WUB[2])
02195 slip = DataType(0.);
02196 else if (wc(2)<min[2])
02197 slip = (wc(2) - WLB[2])/lzs;
02198 else if (wc(2)>max[2])
02199 slip = (WUB[2] - wc(2))/lze;
02200 if (k>mzs && k<mze) {
02201 f[b+base::idx(i,mys,k,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11)
02202 + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
02203 f[b+base::idx(i,mys,k,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13)
02204 + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
02205 }
02206 else if (k==mzs) {
02207 f[b+base::idx(i,mys,k,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k,Nx,Ny)](11)
02208 + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
02209 f[b+base::idx(i,mys,k,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13)
02210 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
02211 }
02212 else if (k==mze) {
02213 f[b+base::idx(i,mys,k,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11)
02214 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
02215 f[b+base::idx(i,mys,k,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,mys-1,k,Nx,Ny)](13)
02216 + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
02217 }
02218 }
02219 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
02220 }
02221 }
02222 break;
02223 case base::Front:
02224 lc_indx[2] = mze*fvec.stepsize(2);
02225 for (register int j=mys; j<=mye; j++) {
02226 lc_indx[1] = j*fvec.stepsize(1);
02227 for (register int i=mxs; i<=mxe; i++) {
02228 lc_indx[0] = i*fvec.stepsize(0);
02229 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
02230 if (wc(0)>=min[0] && wc(0)<=max[0]) {
02231 if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
02232 if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
02233 }
02234 else {
02235 if (wc(0)<WLB[0] || wc(0) >WUB[0])
02236 slip = DataType(0.);
02237 else if (wc(0)<min[0])
02238 slip = (wc(0) - WLB[0])/lxs;
02239 else if (wc(0)>max[0])
02240 slip = (WUB[0] - wc(0))/lxe;
02241 if (i>mxs && i<mxe) {
02242 f[b+base::idx(i,j,mze,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17)
02243 + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
02244 f[b+base::idx(i,j,mze,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16)
02245 + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
02246 }
02247 else if (i==mxs) {
02248 f[b+base::idx(i,j,mze,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](17)
02249 + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
02250 f[b+base::idx(i,j,mze,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16)
02251 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
02252 }
02253 else if (i==mxe) {
02254 f[b+base::idx(i,j,mze,Nx,Ny)](18) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17)
02255 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
02256 f[b+base::idx(i,j,mze,Nx,Ny)](15) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](16)
02257 + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
02258 }
02259 }
02260 if (wc(1)>=min[1] && wc(1)<=max[1]) {
02261 if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
02262 if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
02263 }
02264 else {
02265 if (wc(1)<WLB[1] || wc(1) >WUB[1])
02266 slip = DataType(0.);
02267 else if (wc(1)<min[1])
02268 slip = (wc(1) - WLB[1])/lys;
02269 else if (wc(1)>max[1])
02270 slip = (WUB[1] - wc(1))/lye;
02271 if (j>mys && j<mye) {
02272 f[b+base::idx(i,j,mze,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13)
02273 + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
02274 f[b+base::idx(i,j,mze,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12)
02275 + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
02276 }
02277 else if (j==mys) {
02278 f[b+base::idx(i,j,mze,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](13)
02279 + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
02280 f[b+base::idx(i,j,mze,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12)
02281 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
02282 }
02283 else if (j==mye) {
02284 f[b+base::idx(i,j,mze,Nx,Ny)](14) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13)
02285 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
02286 f[b+base::idx(i,j,mze,Nx,Ny)](11) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mze+1,Nx,Ny)](12)
02287 + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
02288 }
02289 }
02290 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
02291 }
02292 }
02293 break;
02294 case base::Back:
02295 lc_indx[2] = mzs*fvec.stepsize(2);
02296 for (register int j=mys; j<=mye; j++) {
02297 lc_indx[1] = j*fvec.stepsize(1);
02298 for (register int i=mxs; i<=mxe; i++) {
02299 lc_indx[0] = i*fvec.stepsize(0);
02300 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
02301 if (wc(0)>=min[0] && wc(0)<=max[0]) {
02302 if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
02303 if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
02304 }
02305 else {
02306 if (wc(0)<WLB[0] || wc(0) >WUB[0])
02307 slip = DataType(0.);
02308 else if (wc(0)<min[0])
02309 slip = (wc(0) - WLB[0])/lxs;
02310 else if (wc(0)>max[0])
02311 slip = (WUB[0] - wc(0))/lxe;
02312 if (i>mxs && i<mxe) {
02313 f[b+base::idx(i,j,mzs,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15)
02314 + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
02315 f[b+base::idx(i,j,mzs,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18)
02316 + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
02317 }
02318 else if (i==mxs) {
02319 f[b+base::idx(i,j,mzs,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](15)
02320 + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
02321 f[b+base::idx(i,j,mzs,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18)
02322 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
02323 }
02324 else if (i==mxe) {
02325 f[b+base::idx(i,j,mzs,Nx,Ny)](16) = (DataType(1.0)-slip)*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15)
02326 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
02327 f[b+base::idx(i,j,mzs,Nx,Ny)](17) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](18)
02328 + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
02329 }
02330 }
02331 if (wc(1)>=min[1] && wc(1)<=max[1]) {
02332 if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
02333 if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
02334 }
02335 else {
02336 if (wc(1)<WLB[1] || wc(1) >WUB[1])
02337 slip = DataType(0.);
02338 else if (wc(1)<min[1])
02339 slip = (wc(1) - WLB[1])/lys;
02340 else if (wc(1)>max[1])
02341 slip = (WUB[1] - wc(1))/lye;
02342 if (j>mys && j<mye) {
02343 f[b+base::idx(i,j,mzs,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11)
02344 + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
02345 f[b+base::idx(i,j,mzs,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14)
02346 + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
02347 }
02348 else if (j==mys) {
02349 f[b+base::idx(i,j,mzs,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](11)
02350 + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
02351 f[b+base::idx(i,j,mzs,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14)
02352 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
02353 }
02354 else if (j==mye) {
02355 f[b+base::idx(i,j,mzs,Nx,Ny)](12) = (DataType(1.0)-slip)*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11)
02356 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
02357 f[b+base::idx(i,j,mzs,Nx,Ny)](13) = (DataType(1.0)-slip)*f[b+base::idx(i,j,mzs-1,Nx,Ny)](14)
02358 + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
02359 }
02360 }
02361 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
02362 }
02363 }
02364 break;
02365 }
02366 }
02367
02368 else if (type==Periodic) {
02369 switch (side) {
02370 case base::Left:
02371 for (register int k=mzs; k<=mze; k++)
02372 for (register int j=mys; j<=mye; j++) {
02373 f[b+base::idx(mxe,j,k,Nx,Ny)] = f[b+base::idx(mxe+1,j,k,Nx,Ny)];
02374 }
02375 break;
02376 case base::Right:
02377 for (register int k=mzs; k<=mze; k++)
02378 for (register int j=mys; j<=mye; j++) {
02379 f[b+base::idx(mxs,j,k,Nx,Ny)] = f[b+base::idx(mxs-1,j,k,Nx,Ny)];
02380 }
02381 break;
02382 case base::Bottom:
02383 for (register int k=mzs; k<=mze; k++)
02384 for (register int i=mxs; i<=mxe; i++) {
02385 f[b+base::idx(i,mye,k,Nx,Ny)] = f[b+base::idx(i,mye+1,k,Nx,Ny)];
02386 }
02387 break;
02388 case base::Top:
02389 for (register int k=mzs; k<=mze; k++)
02390 for (register int i=mxs; i<=mxe; i++) {
02391 f[b+base::idx(i,mys,k,Nx,Ny)] = f[b+base::idx(i,mys-1,k,Nx,Ny)];
02392 }
02393 break;
02394 case base::Front:
02395 for (register int j=mys; j<=mye; j++)
02396 for (register int i=mxs; i<=mxe; i++) {
02397 f[b+base::idx(i,j,mze,Nx,Ny)] = f[b+base::idx(i,j,mze+1,Nx,Ny)];
02398 }
02399 break;
02400 case base::Back:
02401 for (register int j=mys; j<=mye; j++)
02402 for (register int i=mxs; i<=mxe; i++) {
02403 f[b+base::idx(i,j,mzs,Nx,Ny)] = f[b+base::idx(i,j,mzs-1,Nx,Ny)];
02404 }
02405 break;
02406 }
02407 }
02408
02409 }
02410
02411 virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
02412 const MicroType* f, const point_type* xc, const DataType* distance,
02413 const point_type* normal, DataType* aux=0, const int naux=0,
02414 const int scaling=0) const {
02415 if (type==GFMNoSlipWall || type==GFMSlipWall) {
02416 DataType fac = S0;
02417 if (scaling==base::Physical)
02418 fac /= U0;
02419 for (int n=0; n<nc; n++) {
02420 MacroType q=MacroVariables(f[n]);
02421 DataType u=-q(1);
02422 DataType v=-q(2);
02423 DataType w=-q(3);
02424
02425
02426 if (naux>=3) {
02427 u += fac*aux[naux*n];
02428 v += fac*aux[naux*n+1];
02429 w += fac*aux[naux*n+2];
02430 }
02431 if (type==GFMNoSlipWall) {
02432
02433 q(1) += DataType(2.)*u;
02434 q(2) += DataType(2.)*v;
02435 q(3) += DataType(2.)*w;
02436 }
02437 else if (type==GFMSlipWall) {
02438
02439 DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1)+w*normal[n](2));
02440 q(1) += vl*normal[n](0);
02441 q(2) += vl*normal[n](1);
02442 q(3) += vl*normal[n](2);
02443 }
02444 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02445 }
02446 }
02447
02448 else if (type==GFMExtrapolation) {
02449 for (int n=0; n<nc; n++) {
02450 MacroType q=MacroVariables(f[n]);
02451 DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1)+q(3)*normal[n](2);
02452 q(1) = vl*normal[n](0);
02453 q(2) = vl*normal[n](1);
02454 q(3) = vl*normal[n](2);
02455 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02456 }
02457 }
02458
02459 }
02460
02461 virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
02462 const int skip_ghosts=1) const {
02463 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
02464 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
02465 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
02466 MicroType *f = (MicroType *)fvec.databuffer();
02467 DataType *work = (DataType *)workvec.databuffer();
02468 DCoords dx = base::GH().worldStep(fvec.stepsize());
02469 DataType dt = dx(0)*S0/U0;
02470
02471 if ((cnt>=1 && cnt<=10) || (cnt>=19 && cnt<=24) || (cnt>=28 && cnt<=33)) {
02472 for (register int k=mzs; k<=mze; k++)
02473 for (register int j=mys; j<=mye; j++)
02474 for (register int i=mxs; i<=mxe; i++) {
02475 MacroType q=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
02476 switch(cnt) {
02477 case 1:
02478 if (method[0]==0) work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0+rhop;
02479 else
02480 work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0;
02481 break;
02482 case 2:
02483 work[base::idx(i,j,k,Nx,Ny)]=q(1)*U0/S0;
02484 break;
02485 case 3:
02486 work[base::idx(i,j,k,Nx,Ny)]=q(2)*U0/S0;
02487 break;
02488 case 4:
02489 work[base::idx(i,j,k,Nx,Ny)]=q(3)*U0/S0;
02490 break;
02491 case 6:
02492 work[base::idx(i,j,k,Nx,Ny)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure());
02493 break;
02494 case 7:
02495 work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure();
02496 break;
02497 case 9: {
02498 MicroType feq = Equilibrium(q);
02499 work[base::idx(i,j,k,Nx,Ny)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q(0),Omega(dt),dt),
02500 GasSpeedofSound(),dt);
02501 break;
02502 }
02503 case 10:
02504 work[base::idx(i,j,k,Nx,Ny)]=std::sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3))*U0/S0;
02505 break;
02506 case 19:
02507 case 20:
02508 case 21:
02509 case 22:
02510 case 23:
02511 case 24: {
02512 MicroType feq = Equilibrium(q);
02513 DataType omega;
02514 if (turbulence == laminar)
02515 omega = Omega(dt);
02516 else if (turbulence == LES_Smagorinsky)
02517 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q(0),Omega(dt),dt);
02518 else if (turbulence == WALE) {
02519 DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz;
02520 LocalGradVel(f,i,j,k,Nx,Ny,dx,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz);
02521 omega = Omega_WALE(nup,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
02522 }
02523 else if (turbulence == CSM) {
02524 DataType dxux=0., dxuy=0., dxuz=0., dyux=0., dyuy=0., dyuz=0., dzux=0., dzuy=0., dzuz=0.;
02525 if (skip_ghosts!=1) {
02526 if (i>mxs && i<mxe)
02527 if (j>mys && j<mye)
02528 if (k>mzs && k<mze)
02529 LocalGradVel(f,i,j,k,Nx,Ny,dx,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz);
02530 }
02531 else
02532 LocalGradVel(f,i,j,k,Nx,Ny,dx,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz);
02533 omega = Omega_CSM(dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
02534 }
02535 work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*DeviatoricStress(f[base::idx(i,j,k,Nx,Ny)],feq,omega)(cnt-19);
02536 break;
02537 }
02538 case 28:
02539 case 29:
02540 case 30:
02541 case 31:
02542 case 32:
02543 case 33: {
02544 DataType omega;
02545 if (turbulence == laminar)
02546 omega = Omega(dt);
02547 else if (turbulence == LES_Smagorinsky) {
02548 MicroType feq = Equilibrium(q);
02549 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q(0),Omega(dt),dt);
02550 }
02551 else if (turbulence == WALE) {
02552 DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz;
02553 LocalGradVel(f,i,j,k,Nx,Ny,dx,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz);
02554 omega = Omega_WALE(nup,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
02555 }
02556 else if (turbulence == CSM) {
02557 DataType dxux=0., dxuy=0., dxuz=0., dyux=0., dyuy=0., dyuz=0., dzux=0., dzuy=0., dzuz=0.;
02558 if (skip_ghosts!=1) {
02559 if (i>mxs && i<mxe)
02560 if (j>mys && j<mye)
02561 if (k>mzs && k<mze)
02562 LocalGradVel(f,i,j,k,Nx,Ny,dx,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz);
02563 }
02564 else
02565 LocalGradVel(f,i,j,k,Nx,Ny,dx,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz);
02566 omega = Omega_CSM(dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
02567 }
02568 work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*Stress(f[base::idx(i,j,k,Nx,Ny)],q,omega)(cnt-28);
02569 if (cnt==28 || cnt==30 || cnt==33) work[base::idx(i,j,k,Nx,Ny)]-=BasePressure();
02570 break;
02571 }
02572 }
02573 }
02574 }
02575 else if ((cnt>=11 && cnt<=14) || cnt==18 || (cnt>=40 && cnt<=45)) {
02576 macro_grid_data_type qgrid(fvec.bbox());
02577 MacroType *q = (MacroType *)qgrid.databuffer();
02578 for (register int k=mzs; k<=mze; k++)
02579 for (register int j=mys; j<=mye; j++)
02580 for (register int i=mxs; i<=mxe; i++) {
02581 q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
02582 q[base::idx(i,j,k,Nx,Ny)](0)*=R0;
02583 q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
02584 q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
02585 q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
02586 }
02587
02588 if (skip_ghosts==0) {
02589 switch(cnt) {
02590 case 40:
02591 mxs+=1; mxe-=1;
02592 break;
02593 case 42:
02594 mys+=1; mye-=1;
02595 break;
02596 case 45:
02597 mzs+=1; mze-=1;
02598 break;
02599 case 11:
02600 case 44:
02601 mys+=1; mye-=1; mzs+=1; mze-=1;
02602 break;
02603 case 12:
02604 case 43:
02605 mxs+=1; mxe-=1; mzs+=1; mze-=1;
02606 break;
02607 case 13:
02608 case 41:
02609 mxs+=1; mxe-=1; mys+=1; mye-=1;
02610 break;
02611 case 14:
02612 case 18:
02613 mxs+=1; mxe-=1; mys+=1; mye-=1; mzs+=1; mze-=1;
02614 break;
02615 }
02616 }
02617
02618 DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz;
02619 for (register int k=mzs; k<=mze; k++)
02620 for (register int j=mys; j<=mye; j++)
02621 for (register int i=mxs; i<=mxe; i++) {
02622 if (cnt==18 || cnt==40)
02623 dxux=(q[base::idx(i+1,j,k,Nx,Ny)](1)-q[base::idx(i-1,j,k,Nx,Ny)](1))/(2.*dx(0));
02624 if (cnt==18 || cnt==42)
02625 dyuy=(q[base::idx(i,j+1,k,Nx,Ny)](2)-q[base::idx(i,j-1,k,Nx,Ny)](2))/(2.*dx(1));
02626 if (cnt==18 || cnt==45)
02627 dzuz=(q[base::idx(i,j,k+1,Nx,Ny)](3)-q[base::idx(i,j,k-1,Nx,Ny)](3))/(2.*dx(2));
02628 if (cnt==11 || cnt==14 || cnt==18 || cnt==44) {
02629 dyuz=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(2.*dx(1));
02630 dzuy=(q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(2.*dx(2));
02631 }
02632 if (cnt==12 || cnt==14 || cnt==18 || cnt==43) {
02633 dxuz=(q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(2.*dx(0));
02634 dzux=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(2.*dx(2));
02635 }
02636 if (cnt==13 || cnt==14 || cnt==18 || cnt==41) {
02637 dxuy=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(2.*dx(0));
02638 dyux=(q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(2.*dx(1));
02639 }
02640 switch(cnt) {
02641 case 11:
02642 work[base::idx(i,j,k,Nx,Ny)]=dyuz-dzuy;
02643 break;
02644 case 12:
02645 work[base::idx(i,j,k,Nx,Ny)]=dzux-dxuz;
02646 break;
02647 case 13:
02648 work[base::idx(i,j,k,Nx,Ny)]=dxuy-dyux;
02649 break;
02650 case 14: {
02651 DataType ox = dyuz-dzuy, oy = dzux-dxuz, oz = dxuy-dyux;
02652 work[base::idx(i,j,k,Nx,Ny)] = std::sqrt(ox*ox+oy*oy+oz*oz);
02653 break;
02654 }
02655 case 18:
02656 work[base::idx(i,j,k,Nx,Ny)] = DataType(-0.5)*(4.*dxux*dxux+4.*dyuy*dyuy+4.*dzuz*dzuz+8.*dxuy*dyux+8.*dxuz*dzux+8.*dyuz*dzuy);
02657 break;
02658 case 40:
02659 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dxux;
02660 break;
02661 case 41:
02662 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuy+dyux);
02663 break;
02664 case 42:
02665 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dyuy;
02666 break;
02667 case 43:
02668 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuz+dzux);
02669 break;
02670 case 44:
02671 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dyuz+dzuy);
02672 break;
02673 case 45:
02674 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dzuz;
02675 break;
02676 }
02677 }
02678 }
02679 }
02680
02681 virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
02682 const int skip_ghosts=1) const {
02683 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
02684 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
02685 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
02686 MicroType *f = (MicroType *)fvec.databuffer();
02687 DataType *work = (DataType *)workvec.databuffer();
02688 DCoords dx = base::GH().worldStep(fvec.stepsize());
02689
02690 for (register int k=mzs; k<=mze; k++)
02691 for (register int j=mys; j<=mye; j++)
02692 for (register int i=mxs; i<=mxe; i++) {
02693 switch(cnt) {
02694 case 1:
02695 f[base::idx(i,j,k,Nx,Ny)](0)=work[base::idx(i,j,k,Nx,Ny)]/R0;
02696 break;
02697 case 2:
02698 f[base::idx(i,j,k,Nx,Ny)](1)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02699 break;
02700 case 3:
02701 f[base::idx(i,j,k,Nx,Ny)](2)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02702 break;
02703 case 4:
02704 f[base::idx(i,j,k,Nx,Ny)](3)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02705 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium((const MacroType &)f[base::idx(i,j,k,Nx,Ny)]);
02706 break;
02707 }
02708 }
02709 }
02710
02711 virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
02712 const int verbose) const {
02713 int Nx = fvec.extents(0), Ny = fvec.extents(1);
02714 int b = fvec.bottom();
02715 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
02716 BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
02717 int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
02718 int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
02719 int mzs = std::max(bb.lower(2)/bb.stepsize(2),bbmax.lower(2)/bbmax.stepsize(2));
02720 int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
02721 int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
02722 int mze = std::min(bb.upper(2)/bb.stepsize(2),bbmax.upper(2)/bbmax.stepsize(2));
02723 MicroType *f = (MicroType *)fvec.databuffer();
02724
02725 int result = 1;
02726 for (register int k=mzs; k<=mze; k++)
02727 for (register int j=mys; j<=mye; j++)
02728 for (register int i=mxs; i<=mxe; i++)
02729 for (register int l=0; l<base::NMicroVar(); l++)
02730 if (!(f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l)>-1.e37)) {
02731 result = 0;
02732 if (verbose) {
02733 DataType WLB[3];
02734 base::GH().wlb(WLB);
02735 DCoords dx = base::GH().worldStep(fvec.stepsize());
02736 std::cerr << "D3Q19-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << "," << k-mdz[l] << ")(" << l << ")="
02737 << f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l) << " " << fvec.bbox()
02738 << " Coords: (" << (i+0.5)*dx(0)+WLB[0] << "," << (j+0.5)*dx(1)+WLB[1] << ","
02739 << (k+0.5)*dx(2)+WLB[2] << ")" << std::endl;
02740 }
02741 }
02742 return result;
02743 }
02744
02745 inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
02746 TensorType P;
02747 if (stressPath==0) {
02748 if (method[1]==0) {
02749 P(0) = q(0)*q(1)*q(1)-(f(1)+f(2)+f(7)+f(8)+f(9)+f(10)+f(15)+f(16)+f(17)+f(18))+cs2*q(0);
02750 P(1) = q(0)*q(1)*q(2)-(f(7)+f(8)-f(9)-f(10));
02751 P(2) = q(0)*q(2)*q(2)-(f(3)+f(4)+f(7)+f(8)+f(9)+f(10)+f(11)+f(12)+f(13)+f(14))+cs2*q(0);
02752 P(3) = q(0)*q(1)*q(3)-(f(15)+f(16)-f(17)-f(18));
02753 P(4) = q(0)*q(2)*q(3)-(f(11)+f(12)-f(13)-f(14));
02754 P(5) = q(0)*q(3)*q(3)-(f(5)+f(6)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18))+cs2*q(0);
02755 P *= -(DataType(1.0)-DataType(0.5)*om);
02756 }
02757 else {
02758 MicroType feq = Equilibrium(q);
02759 P = DeviatoricStress(f,feq,om);
02760 }
02761 P(0) -= LatticeBasePressure(q(0));
02762 P(2) -= LatticeBasePressure(q(0));
02763 P(5) -= LatticeBasePressure(q(0));
02764 }
02765 else if (stressPath==1) {
02766 P = Stress_velocitySpace(f,Equilibrium(q),om);
02767 }
02768 return P;
02769 }
02770
02771 inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
02772 TensorType Sigma;
02773 if (stressPath==0) {
02774 if (method[2]==0) {
02775 MicroType fneq = feq - f;
02776 Sigma(0) = fneq(1)+fneq(2)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
02777 Sigma(1) = fneq(7)+fneq(8)-fneq(9)-fneq(10);
02778 Sigma(2) = fneq(3)+fneq(4)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(11)+fneq(12)+fneq(13)+fneq(14);
02779 Sigma(3) = fneq(15)+fneq(16)-fneq(17)-fneq(18);
02780 Sigma(4) = fneq(11)+fneq(12)-fneq(13)-fneq(14);
02781 Sigma(5) = fneq(5)+fneq(6)+fneq(11)+fneq(12)+fneq(13)+fneq(14)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
02782 Sigma *= -(DataType(1.0)-DataType(0.5)*om);
02783 }
02784 else {
02785 MacroType q = MacroVariables(f);
02786 Sigma = Stress(f,q,om);
02787 Sigma(0) += LatticeBasePressure(q(0));
02788 Sigma(2) += LatticeBasePressure(q(0));
02789 Sigma(5) += LatticeBasePressure(q(0));
02790 }
02791 }
02792 else if (stressPath==1) {
02793 Sigma = DeviatoricStress_velocitySpace(f,feq,om);
02794 }
02795 return Sigma;
02796 }
02797
02798 virtual int NMethodOrder() const { return 2; }
02799
02800 inline const DataType& L0() const { return base::LengthScale(); }
02801 inline const DataType& T0() const { return base::TimeScale(); }
02802 inline void SetDensityScale(const DataType r0) { R0 = r0; }
02803 inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
02804 inline void SetSpeedUp(const DataType s0) { S0 = s0; }
02805 inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
02806
02807 inline const DataType& DensityScale() const { return R0; }
02808 inline const DataType VelocityScale() const { return U0/S0; }
02809
02810 inline const DataType& SpeedUp() const { return S0; }
02811
02812
02813 inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
02814 inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
02815 inline virtual DataType LatticeBasePressure(const DataType rho) const {
02816 return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
02817 }
02818
02819
02820 inline void SetGas(DataType rho, DataType nu, DataType cs) {
02821 rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
02822 SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
02823 }
02824 inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
02825 inline const TensorType Stress_velocitySpace(const MicroType &f, const MicroType &feq, const DataType om) const {
02826 if (method[0]==3) {
02827 MicroType fneq;
02828 TensorType P;
02829 MacroType q=MacroVariables(f);
02831 DataType cl = DataType(1.0);
02832 DataType cmu_ = cl - q(1);
02833 DataType cmu2_ = cmu_*cmu_;
02834 DataType cpu_ = cl + q(1);
02835 DataType cpu2_ = cpu_*cpu_;
02836 DataType cmv_ = cl - q(2);
02837 DataType cmv2_ = cmv_*cmv_;
02838 DataType cpv_ = cl + q(2);
02839 DataType cpv2_ = cpv_*cpv_;
02840 DataType cmw_ = cl - q(3);
02841 DataType cmw2_ = cmw_*cmw_;
02842 DataType cpw_ = cl + q(3);
02843 DataType cpw2_ = cpw_*cpw_;
02844 DataType u2 = q(1)*q(1);
02845 DataType v2 = q(2)*q(2);
02846 DataType w2 = q(3)*q(3);
02847
02849 fneq = (f - feq);
02850 P(0) = (fneq(0)+fneq(3)+fneq(4)+fneq(5)+fneq(6)+fneq(11)+fneq(12)+fneq(13)+fneq(14))*u2
02851 + (fneq(1)+fneq(7)+fneq(9)+fneq(15)+fneq(17))*cmu2_
02852 + (fneq(2)+fneq(8)+fneq(10)+fneq(16)+fneq(18))*cpu2_;
02853 P(1) = (fneq(0)+fneq(5)+fneq(6))*q(1)*q(2)
02854 - (fneq(3)+fneq(11)+fneq(13))*q(1)*cmv_ + (fneq(4)+fneq(12)+fneq(14))*q(1)*cpv_
02855 - (fneq(1)+fneq(15)+fneq(17))*q(2)*cmu_ + (fneq(2)+fneq(16)+fneq(18))*q(2)*cpu_
02856 + fneq(7)*cmu_*cmv_ - fneq(9)*cmu_*cpv_
02857 - fneq(10)*cpu_*cmv_ + fneq(8)*cpu_*cpv_;
02858 P(2) = (fneq(0)+fneq(1)+fneq(2)+fneq(5)+fneq(6)+fneq(15)+fneq(16)+fneq(17)+fneq(18))*v2
02859 + (fneq(3)+fneq(7)+fneq(10)+fneq(11)+fneq(13))*cmv2_
02860 + (fneq(4)+fneq(8)+fneq(9)+fneq(12)+fneq(14))*cpv2_;
02861 P(3) = (fneq(0)+fneq(3)+fneq(4))*q(1)*q(3)
02862 - (fneq(5)+fneq(11)+fneq(14))*q(1)*cmw_ + (fneq(6)+fneq(12)+fneq(13))*q(1)*cpw_
02863 - (fneq(1)+fneq(7)+fneq(9))*q(3)*cmu_ + (fneq(2)+fneq(8)+fneq(10))*q(3)*cpu_
02864 + fneq(15)*cmu_*cmw_ - fneq(17)*cmu_*cpw_
02865 - fneq(18)*cpu_*cmw_ + fneq(16)*cpu_*cpw_;
02866 P(4) = (fneq(0)+fneq(1)+fneq(2))*q(2)*q(3)
02867 - (fneq(5)+fneq(15)+fneq(18))*q(2)*cmw_ + (fneq(6)+fneq(16)+fneq(17))*q(2)*cpw_
02868 - (fneq(3)+fneq(7)+fneq(10))*q(3)*cmv_ + (fneq(4)+fneq(8)+fneq(9))*q(3)*cpv_
02869 + fneq(11)*cmv_*cmw_ - fneq(13)*cmv_*cpw_
02870 - fneq(14)*cpv_*cmw_ + fneq(12)*cpv_*cpw_;
02871 P(5) = (fneq(0)+fneq(1)+fneq(2)+fneq(3)+fneq(4)+fneq(7)+fneq(8)+fneq(9)+fneq(10))*w2
02872 + (fneq(5)+fneq(11)+fneq(14)+fneq(15)+fneq(18))*cmw2_
02873 + (fneq(6)+fneq(12)+fneq(13)+fneq(16)+fneq(17))*cpw2_;
02874 return P;
02875 }
02876 else {
02877 TensorType P = DeviatoricStress_velocitySpace(f,feq,om);
02878 MacroType q = MacroVariables(f);
02879 P(0) -= cs2*q(0); P(2) -= cs2*q(0); P(5) -= cs2*q(0);
02880 return P;
02881 }
02882 }
02883 inline const TensorType DeviatoricStress_velocitySpace(const MicroType &f, const MicroType &feq, const DataType om) const {
02884 if (method[0]==3) {
02885 TensorType Sigma = Stress_velocitySpace(f,feq,om);
02886 MacroType q = MacroVariables(f);
02887 Sigma(0) += cs2*q(0); Sigma(2) += cs2*q(0); Sigma(5) += cs2*q(0);
02888 return Sigma;
02889 }
02890 else {
02891 MicroType fneq;
02892 TensorType Sigma;
02893
02895 DataType cl = DataType(1.0);
02896 fneq = (f - feq);
02897 Sigma(0) = (cl*cl-DataType(0.5))*( fneq(1)+fneq(2)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(15)+fneq(16)+fneq(17)+fneq(18) );
02898 Sigma(1) = (cl*cl-DataType(0.5))*( fneq(7)+fneq(8) )
02899 - (cl*cl+DataType(0.5))*( fneq(9)+fneq(10) );
02900 Sigma(2) = (cl*cl-DataType(0.5))*( fneq(3)+fneq(4)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(11)+fneq(12)+fneq(13)+fneq(14) );
02901 Sigma(3) = (cl*cl-DataType(0.5))*( fneq(15)+fneq(16) )
02902 - (cl*cl+DataType(0.5))*( fneq(17)+fneq(18) );
02903 Sigma(4) = (cl*cl-DataType(0.5))*( fneq(11)+fneq(12) );
02904 Sigma(5) = (cl*cl-DataType(0.5))*( fneq(5)+fneq(6)+fneq(11)+fneq(12)+fneq(13)+fneq(14)+fneq(15)+fneq(16)+fneq(17)+fneq(18) );
02905 return (DataType(1.0)-DataType(0.5)*om)*Sigma;
02906 }
02907 }
02908 inline const TensorType StrainComponents(const DataType rho, const TensorType& Sigma, const DataType om, const DataType csmag) const {
02909 DataType omTmp = rho*cs2/om;
02910 DataType rhocspcsmTmp = DataType(2.0)*rho*cs2*cs2*csmag*csmag*csmag*csmag;
02911 DataType sigma2 = Magnitude(Sigma);
02912 DataType stmp = ( (rhocspcsmTmp*sigma2 > (mfp/U0)) ? (-omTmp + std::sqrt( (omTmp)*(omTmp) + rhocspcsmTmp*sigma2 ) )/( rhocspcsmTmp*sigma2 )
02913 : (-omTmp + std::sqrt( (omTmp)*(omTmp) + mfp/U0 ) )/( mfp/U0 ));
02914 TensorType Strain = stmp*Sigma;
02915 return Strain;
02916 }
02917 inline const DataType Strain(const DataType rho, const TensorType& Sigma, const DataType om, const DataType csmag) const {
02918 return Magnitude(StrainComponents(rho,Sigma,om,csmag));
02919 }
02920 inline const DataType StrainLaminar(const DataType rho, const TensorType& Sigma, const DataType om) const {
02921 return om/(DataType(2.)*rho*cs2)*Magnitude(Sigma);
02922 }
02923 inline const DataType Magnitude(const TensorType& A) const {
02924 return std::sqrt(DataType(2.0)*(A(0)*A(0) + DataType(2.0)*A(1)*A(1) + A(2)*A(2)
02925 + DataType(2.0)*A(3)*A(3) + DataType(2.0)*A(4)*A(4) + A(5)*A(5)));
02926 }
02927 inline const DataType Omega_LES_Smagorinsky( MicroType &f, const MicroType &feq, const DataType rho,
02928 const DataType om, const DataType dt) const {
02938 TensorType Sigma = DeviatoricStress(f,feq,om)/(DataType(1.0)-DataType(0.5)*om);
02939 DataType S = Strain(rho,Sigma,om,Cs_Smagorinsky);
02940 DataType nu_t = Cs_Smagorinsky*Cs_Smagorinsky*S/S0;
02941 if (nu_t < 0.) nu_t = 0.;
02942 DataType om_eff = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
02943
02944 return ( om_eff );
02945 }
02946 inline const int EquilibriumType() const { return method[0]; }
02947 inline const int TurbulenceType() const { return turbulence; }
02948 inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
02949 inline const DataType& GasDensity() const { return rhop; }
02950 inline const DataType& GasViscosity() const { return nup; }
02951 inline const DataType& GasSpeedofSound() const { return csp; }
02952 inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
02953 { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
02954 inline virtual const MacroType NormalDerivative(const MacroType qa,const MacroType qb, const MacroType qc) const
02955 { return DataType(-1.5)*qa + DataType(2.)*qb - DataType(0.5)*qc; }
02956 inline virtual Vector<DataType,3> WaveAmplitudes(const DataType rho, const DataType vn, const DataType drhodn, const DataType dvndn, const DataType dvt0dn, const DataType dvt1dn) const {
02957 Vector<DataType,3> L;
02958 L(0) = (vn - LatticeSpeedOfSound())*(cs2*drhodn - rho*LatticeSpeedOfSound()*dvndn);
02959 L(1) = vn*dvt0dn;
02960 L(2) = vn*dvt1dn;
02961 return L;
02962 }
02963
02964
02965 inline virtual DataType Omega_WALE( const DataType nu, const DataType dxux, const DataType dxuy, const DataType dxuz,
02966 const DataType dyux, const DataType dyuy, const DataType dyuz,
02967 const DataType dzux, const DataType dzuy, const DataType dzuz,
02968 const DCoords dx, const DataType dt) const {
02974 DataType S2 = DataType(0.5)*( (dxuy + dyux)*(dxuy + dyux) + (dxuz + dzux)*(dxuz + dzux) + (dyuz + dzuy)*(dyuz + dzuy) ) + dxux*dxux + dyuy*dyuy + dzuz*dzuz;
02975 DataType O2 = DataType(0.5)*( (dxuy - dyux)*(dxuy - dyux) + (dxuz - dzux)*(dxuz - dzux) + (dyuz - dzuy)*(dyuz - dzuy) );
02976 DataType IV = ((dyuz/2. - dzuy/2.)*((dxuy/2. + dyux/2.)*(dyuz/2. + dzuy/2.) + dxux*(dxuz/2. + dzux/2.) + dzuz*(dxuz/2. + dzux/2.)) - (dxuy/2. - dyux/2.)*((dxuy/2. + dyux/2.)*(dxuy/2. + dyux/2.) + (dxuz/2. + dzux/2.)*(dxuz/2. + dzux/2.) + dxux*dxux))*(dxuy/2. - dyux/2.) - ((dyuz/2. - dzuy/2.)*((dxuz/2. + dzux/2.)*(dyuz/2. + dzuy/2.) + dxux*(dxuy/2. + dyux/2.) + dyuy*(dxuy/2. + dyux/2.)) + (dxuz/2. - dzux/2.)*((dxuy/2. + dyux/2.)*(dxuy/2. + dyux/2.) + (dxuz/2. + dzux/2.)*(dxuz/2. + dzux/2.) + dxux*dxux))*(dxuz/2. - dzux/2.) - ((dxuz/2. - dzux/2.)*((dxuz/2. + dzux/2.)*(dyuz/2. + dzuy/2.) + dxux*(dxuy/2. + dyux/2.) + dyuy*(dxuy/2. + dyux/2.)) + (dyuz/2. - dzuy/2.)*((dxuy/2. + dyux/2.)*(dxuy/2. + dyux/2.) + (dyuz/2. + dzuy/2.)*(dyuz/2. + dzuy/2.) + dyuy*dyuy))*(dyuz/2. - dzuy/2.) - ((dxuz/2. - dzux/2.)*((dxuy/2. + dyux/2.)*(dxuz/2. + dzux/2.) + dyuy*(dyuz/2. + dzuy/2.) + dzuz*(dyuz/2. + dzuy/2.)) + (dxuy/2. - dyux/2.)*((dxuy/2. + dyux/2.)*(dxuy/2. + dyux/2.) + (dyuz/2. + dzuy/2.)*(dyuz/2. + dzuy/2.) + dyuy*dyuy))*(dxuy/2. - dyux/2.) - ((dxuy/2. - dyux/2.)*((dxuy/2. + dyux/2.)*(dxuz/2. + dzux/2.) + dyuy*(dyuz/2. + dzuy/2.) + dzuz*(dyuz/2. + dzuy/2.)) + (dxuz/2. - dzux/2.)*((dxuz/2. + dzux/2.)*(dxuz/2. + dzux/2.) + (dyuz/2. + dzuy/2.)*(dyuz/2. + dzuy/2.) + dzuz*dzuz))*(dxuz/2. - dzux/2.) + ((dxuy/2. - dyux/2.)*((dxuy/2. + dyux/2.)*(dyuz/2. + dzuy/2.) + dxux*(dxuz/2. + dzux/2.) + dzuz*(dxuz/2. + dzux/2.)) - (dyuz/2. - dzuy/2.)*((dxuz/2. + dzux/2.)*(dxuz/2. + dzux/2.) + (dyuz/2. + dzuy/2.)*(dyuz/2. + dzuy/2.) + dzuz*dzuz))*(dyuz/2. - dzuy/2.);
02977 DataType Sdij2 = DataType(1./6.)*(S2*S2+O2*O2)+DataType(2./3.)*S2*O2+DataType(2.)*IV;
02978 DataType eddyVisc = std::pow(Sdij2,1.5)/( std::pow(S2,2.5) + std::pow(Sdij2,1.25) );
02979 if (std::isfinite(eddyVisc)==0 ) { eddyVisc = DataType(0.0); }
02980 DataType nu_t = DataType(0.25)*dx(0)*dx(1)*std::sqrt(DataType(2.)*eddyVisc);
02981 DataType om_eff = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
02982 if (om_eff < DataType(1.000001)) { std::cout <<" om="<<om_eff<<" nu="<<nu<<" nu_t="<<nu_t<< " lower limit\t"; om_eff = DataType(1.000001); }
02983 else if (om_eff > DataType(1.999999)) { std::cout <<" om="<<om_eff<<" nu="<<nu<<" nu_t="<<nu_t<< " upper limit\t"; om_eff = DataType(1.999999); }
02984 return om_eff;
02985 }
02986
02987 inline virtual void LocalCollisionWALE( MicroType &f, const DataType nu,
02988 const MacroType qxp, const MacroType qxm,
02989 const MacroType qyp, const MacroType qym,
02990 const MacroType qzp, const MacroType qzm,
02991 const DCoords dx, const DataType dt) const {
02992 DataType dxux=(qxp(1)-qxm(1))/(DataType(2.)*dx(0))*U0/S0;
02993 DataType dxuy=(qxp(2)-qxm(2))/(DataType(2.)*dx(0))*U0/S0;
02994 DataType dxuz=(qxp(3)-qxm(3))/(DataType(2.)*dx(0))*U0/S0;
02995 DataType dyux=(qyp(1)-qym(1))/(DataType(2.)*dx(1))*U0/S0;
02996 DataType dyuy=(qyp(2)-qym(2))/(DataType(2.)*dx(1))*U0/S0;
02997 DataType dyuz=(qyp(3)-qym(3))/(DataType(2.)*dx(1))*U0/S0;
02998 DataType dzux=(qzp(1)-qzm(1))/(DataType(2.)*dx(2))*U0/S0;
02999 DataType dzuy=(qzp(2)-qzm(2))/(DataType(2.)*dx(2))*U0/S0;
03000 DataType dzuz=(qzp(3)-qzm(3))/(DataType(2.)*dx(2))*U0/S0;
03001
03002 DataType omega = Omega_WALE(nu,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
03003 MicroType feq = Equilibrium(MacroVariables(f));
03004 for (register int qi=0; qi<19; qi++)
03005 f(qi) = (DataType(1.)-omega)*f(qi) + omega*feq(qi);
03006 }
03007
03008 virtual void CollisionWALE(vec_grid_data_type &fvec, const double& dt) const {
03009 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
03010 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
03011 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
03012 DCoords dx = base::GH().worldStep(fvec.stepsize());
03013 MicroType *f = (MicroType *)fvec.databuffer();
03014 MicroType feq;
03015 DataType omega = Omega(dt);
03016 DataType nu = LatticeViscosity(omega)/S0/S0;
03017
03018 macro_grid_data_type qgrid(fvec.bbox());
03019 MacroType *q = (MacroType *)qgrid.databuffer();
03020 int x=0, y=0, z=0;
03021 for (register int k=mzs-1; k<=mze+1; k++) {
03022 if (k==mzs-1) z=1;
03023 else if (k==mze+1) z=-1;
03024 else z=0;
03025 for (register int j=mys-1; j<=mye+1; j++) {
03026 if (j==mys-1) y=1;
03027 else if (j==mye+1) y=-1;
03028 else y=0;
03029 for (register int i=mxs-1; i<=mxe+1; i++) {
03030 if (i==mxs-1) x=1;
03031 else if (i==mxe+1) x=-1;
03032 else x=0;
03033 q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i+x,j+y,k+z,Nx,Ny)]);
03034 q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
03035 q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
03036 q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
03037 }
03038 }
03039 }
03040 DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz;
03041 for (register int k=mzs; k<=mze; k++)
03042 for (register int j=mys; j<=mye; j++)
03043 for (register int i=mxs; i<=mxe; i++) {
03044 dxux=(q[base::idx(i+1,j,k,Nx,Ny)](1)-q[base::idx(i-1,j,k,Nx,Ny)](1))/(DataType(2.)*dx(0));
03045 dxuy=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(DataType(2.)*dx(0));
03046 dxuz=(q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(DataType(2.)*dx(0));
03047 dyux=(q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(DataType(2.)*dx(1));
03048 dyuy=(q[base::idx(i,j+1,k,Nx,Ny)](2)-q[base::idx(i,j-1,k,Nx,Ny)](2))/(DataType(2.)*dx(1));
03049 dyuz=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(DataType(2.)*dx(1));
03050 dzux=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(DataType(2.)*dx(2));
03051 dzuy=(q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(DataType(2.)*dx(2));
03052 dzuz=(q[base::idx(i,j,k+1,Nx,Ny)](3)-q[base::idx(i,j,k-1,Nx,Ny)](3))/(DataType(2.)*dx(2));
03053
03054 omega = Omega_WALE(nu,dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
03055 if (omega < DataType(1.000001)) { omega = DataType(1.000001); std::cout << "lower limit\t"; }
03056 else if (omega > DataType(1.999999)) { omega = DataType(1.999999); std::cout << "upper limit\t"; }
03057 feq = Equilibrium(MacroVariables(f[base::idx(i,j,k,Nx,Ny)]));
03058 for (register int qi=0; qi<19; qi++)
03059 f[base::idx(i,j,k,Nx,Ny)](qi) = (DataType(1.)-omega)*f[base::idx(i,j,k,Nx,Ny)](qi) + omega*feq(qi);
03060 }
03061 }
03062
03063 inline virtual DataType Omega_CSM(const DataType dxux, const DataType dxuy, const DataType dxuz,
03064 const DataType dyux, const DataType dyuy, const DataType dyuz,
03065 const DataType dzux, const DataType dzuy, const DataType dzuz,
03066 const DCoords dx, const DataType dt) const {
03073 DataType S2, Q, E, FCS, C, eddyVisc, nu_t;
03074 S2 = DataType(0.5)*( (dxuy + dyux)*(dxuy + dyux) + (dxuz + dzux)*(dxuz + dzux) + (dyuz + dzuy)*(dyuz + dzuy) ) + dxux*dxux + dyuy*dyuy + dzuz*dzuz;
03075 Q = DataType(-2.)*(dxux*dxux+dyuy*dyuy+dzuz*dzuz)-DataType(4.)*(dxuy*dyux-dxuz*dzux-dyuz*dzuy);
03076 E = DataType(0.5)*(dxux*dxux +dyuy*dyuy +dzuz*dzuz +dxuy*dxuy +dxuz*dxuz +dyux*dyux +dyuz*dyuz +dzux*dzux + dzuy*dzuy );
03077 if (std::isfinite(S2)==0 || std::isfinite(Q)==0 || std::isfinite(E)==0 || E==0.) {
03078 nu_t = DataType(0.0);
03079 }
03080 else {
03081 FCS = Q/E;
03082 if (std::isfinite(FCS)==0) {
03083 nu_t = DataType(0.0);
03084 }
03085 else {
03086 if (FCS<DataType(-1.0)) { FCS = DataType(-1.0); }
03087 else if (FCS>DataType(1.0)) { FCS = DataType(1.0); }
03088 C = DataType(1./22.)*std::pow(std::fabs(FCS),DataType(1.5))*(DataType(1.)-FCS);
03089 eddyVisc = dx(0)*dx(1)*std::sqrt(DataType(2.)*S2);
03090 if (std::isfinite(eddyVisc)==0 ) { eddyVisc = DataType(0.0); }
03091 if (std::isfinite(C)==0 ) { C = DataType(0.0); }
03092 nu_t = C*eddyVisc;
03093 }
03094 }
03095 DataType om_eff = (cs2p*dt/S0/( (nup + nu_t)*S0 + cs2p*dt/S0*DataType(0.5) ));
03096 return om_eff;
03097 }
03098
03099 inline virtual void LocalCollisionCSM( MicroType &f, const MacroType qxp, const MacroType qxm,
03100 const MacroType qyp, const MacroType qym,
03101 const MacroType qzp, const MacroType qzm,
03102 const DCoords dx, const DataType dt) const {
03103 DataType dxux=(qxp(1)-qxm(1))/(2.*dx(0))*U0/S0;
03104 DataType dxuy=(qxp(2)-qxm(2))/(2.*dx(0))*U0/S0;
03105 DataType dxuz=(qxp(3)-qxm(3))/(2.*dx(0))*U0/S0;
03106 DataType dyux=(qyp(1)-qym(1))/(2.*dx(1))*U0/S0;
03107 DataType dyuy=(qyp(2)-qym(2))/(2.*dx(1))*U0/S0;
03108 DataType dyuz=(qyp(3)-qym(3))/(2.*dx(1))*U0/S0;
03109 DataType dzux=(qzp(1)-qzm(1))/(2.*dx(2))*U0/S0;
03110 DataType dzuy=(qzp(2)-qzm(2))/(2.*dx(2))*U0/S0;
03111 DataType dzuz=(qzp(3)-qzm(3))/(2.*dx(2))*U0/S0;
03112
03113 DataType omega = Omega_CSM(dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
03114 MicroType feq = Equilibrium(MacroVariables(f));
03115 for (register int qi=0; qi<19; qi++)
03116 f(qi) = (DataType(1.)-omega)*f(qi) + omega*feq(qi);
03117 }
03118
03119 virtual void CollisionCSM(vec_grid_data_type &fvec, const double& dt) const {
03120 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
03121 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
03122 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
03123 DCoords dx = base::GH().worldStep(fvec.stepsize());
03124 MicroType *f = (MicroType *)fvec.databuffer();
03125 MicroType feq;
03126 DataType omega = Omega(dt);
03127
03128 macro_grid_data_type qgrid(fvec.bbox());
03129 MacroType *q = (MacroType *)qgrid.databuffer();
03130 int x=0, y=0, z=0;
03131 for (register int k=mzs-1; k<=mze+1; k++) {
03132 for (register int j=mys-1; j<=mye+1; j++) {
03133 for (register int i=mxs-1; i<=mxe+1; i++) {
03134 q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i+x,j+y,k+z,Nx,Ny)]);
03135 q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
03136 q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
03137 q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
03138 }
03139 }
03140 }
03141 DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz;
03142 for (register int k=mzs; k<=mze; k++)
03143 for (register int j=mys; j<=mye; j++)
03144 for (register int i=mxs; i<=mxe; i++) {
03145 dxux=(q[base::idx(i+1,j,k,Nx,Ny)](1)-q[base::idx(i-1,j,k,Nx,Ny)](1))/(DataType(2.)*dx(0));
03146 dxuy=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(DataType(2.)*dx(0));
03147 dxuz=(q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(DataType(2.)*dx(0));
03148
03149 dyux=(q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(DataType(2.)*dx(1));
03150 dyuy=(q[base::idx(i,j+1,k,Nx,Ny)](2)-q[base::idx(i,j-1,k,Nx,Ny)](2))/(DataType(2.)*dx(1));
03151 dyuz=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(DataType(2.)*dx(1));
03152
03153 dzux=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(DataType(2.)*dx(2));
03154 dzuy=(q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(DataType(2.)*dx(2));
03155 dzuz=(q[base::idx(i,j,k+1,Nx,Ny)](3)-q[base::idx(i,j,k-1,Nx,Ny)](3))/(DataType(2.)*dx(2));
03156
03157 omega = Omega_CSM(dxux,dxuy,dxuz,dyux,dyuy,dyuz,dzux,dzuy,dzuz,dx,dt);
03158 feq = Equilibrium(MacroVariables(f[base::idx(i,j,k,Nx,Ny)]));
03159 for (register int qi=0; qi<19; qi++)
03160 f[base::idx(i,j,k,Nx,Ny)](qi) = (DataType(1.)-omega)*f[base::idx(i,j,k,Nx,Ny)](qi) + omega*feq(qi);
03161 }
03162
03163 }
03164
03165 inline void LocalGradVel(const MicroType *f, const int i, const int j, const int k, const int Nx, const int Ny, const DCoords dx,
03166 DataType &dxux, DataType &dxuy, DataType &dxuz,
03167 DataType &dyux, DataType &dyuy, DataType &dyuz,
03168 DataType &dzux, DataType &dzuy, DataType &dzuz) const {
03169 MacroType fxp=MacroVariables(f[base::idx(i+1,j,k,Nx,Ny)]);
03170 MacroType fxm=MacroVariables(f[base::idx(i-1,j,k,Nx,Ny)]);
03171 MacroType fyp=MacroVariables(f[base::idx(i,j+1,k,Nx,Ny)]);
03172 MacroType fym=MacroVariables(f[base::idx(i,j-1,k,Nx,Ny)]);
03173 MacroType fzp=MacroVariables(f[base::idx(i,j,k+1,Nx,Ny)]);
03174 MacroType fzm=MacroVariables(f[base::idx(i,j,k-1,Nx,Ny)]);
03175 dxux = (fxp(1)-fxm(1))/(DataType(2.)*dx(0));
03176 dxuy = (fxp(2)-fxm(2))/(DataType(2.)*dx(0));
03177 dxuz = (fxp(3)-fxm(3))/(DataType(2.)*dx(0));
03178 dyux = (fyp(1)-fym(1))/(DataType(2.)*dx(1));
03179 dyuy = (fyp(2)-fym(2))/(DataType(2.)*dx(1));
03180 dyuz = (fyp(3)-fym(3))/(DataType(2.)*dx(1));
03181 dzux = (fzp(1)-fzm(1))/(DataType(2.)*dx(2));
03182 dzuy = (fzp(2)-fzm(2))/(DataType(2.)*dx(2));
03183 dzuz = (fzp(3)-fzm(3))/(DataType(2.)*dx(2));
03184 }
03185
03186
03187 inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
03188 inline virtual const DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); }
03189 inline virtual const DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); }
03190
03191 protected:
03192 DataType cs2, cs22, cssq;
03193 DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp, mfp;
03194 DataType Cs_Smagorinsky, turbulence;
03195 int method[3], mdx[NUMMICRODIST], mdy[NUMMICRODIST], mdz[NUMMICRODIST];
03196 int stressPath;
03197 };
03198
03199
03200 #endif