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
00050 template <class DataType>
00051 class LBMD3Q19 : public LBMBase<Vector<DataType,19>, Vector<DataType,4>, 3> {
00052 typedef LBMBase<Vector<DataType,19>, Vector<DataType,4>, 3> base;
00053 public:
00054 typedef typename base::vec_grid_data_type vec_grid_data_type;
00055 typedef typename base::grid_data_type grid_data_type;
00056 typedef typename base::MicroType MicroType;
00057 typedef typename base::MacroType MacroType;
00058 typedef GridData<MacroType,3> macro_grid_data_type;
00059 typedef typename base::SideName SideName;
00060 typedef typename base::point_type point_type;
00061 typedef Vector<DataType,6> TensorType;
00062
00063 enum ICPredefined { GasAtRest, ConstantMacro, ConstantMicro };
00064 enum BCPredefined { Symmetry, SlipWall, NoSlipWall, Inlet, Outlet, Pressure, SlidingWall };
00065 enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall, GFMWallLaw, GFMBounceBack };
00066 enum TurbulenceModel { laminar, LES_Smagorinsky, LES_dynamic };
00067
00068 LBMD3Q19() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.) {
00069 cs2 = DataType(1.)/DataType(3.);
00070 cs22 = DataType(2.)*cs2;
00071 cssq = DataType(2.)/DataType(9.);
00072 method[0] = 2; method[1] = 0; method[2] = 0;
00073 turbulence = laminar;
00074 mdx[0]=mdx[3]=mdx[4]=mdx[5]=mdx[6]=mdx[11]=mdx[12]=mdx[13]=mdx[14]=0;
00075 mdx[1]=mdx[7]=mdx[9]=mdx[15]=mdx[17]=1;
00076 mdx[2]=mdx[8]=mdx[10]=mdx[16]=mdx[18]=-1;
00077
00078 mdy[0]=mdy[1]=mdy[2]=mdy[5]=mdy[6]=mdy[15]=mdy[16]=mdy[17]=mdy[18]=0;
00079 mdy[3]=mdy[7]=mdy[10]=mdy[11]=mdy[13]=1;
00080 mdy[4]=mdy[8]=mdy[9]=mdy[12]=mdy[14]=-1;
00081
00082 mdz[0]=mdz[1]=mdz[2]=mdz[3]=mdz[4]=mdz[7]=mdz[8]=mdz[9]=mdz[10]=0;
00083 mdz[5]=mdz[11]=mdz[14]=mdz[15]=mdz[18]=1;
00084 mdz[6]=mdz[12]=mdz[13]=mdz[16]=mdz[17]=-1;
00085
00086 }
00087
00088 virtual ~LBMD3Q19() {}
00089
00090 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00091 base::LocCtrl = Ctrl.getSubDevice(prefix+"D3Q19");
00092 RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00093 RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00094 RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00095 RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00096 RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00097 RegisterAt(base::LocCtrl,"Gas_nu",nup);
00098 RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00099 RegisterAt(base::LocCtrl,"Gas_W",Wp);
00100 RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00101 RegisterAt(base::LocCtrl,"SpeedUp",S0);
00102 }
00103
00104 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00105 base::SetupData(gh,ghosts);
00106 if (rhop>0.) R0 = rhop;
00107 if (csp>0.) {
00108 cs2p = csp*csp;
00109 SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00110 }
00111 std::cout << "D3Q19: Gas_rho=" << rhop << " Gas_Cs=" << csp
00112 << " Gas_nu=" << nup << std::endl;
00113 std::cout << "D3Q19: L0=" << L0() << " T0=" << T0() << " U0=" << U0
00114 << " R0=" << R0 << std::endl;
00115 if ( TurbulenceType() == laminar ) {
00116 std::cout << "LBM Setup() Laminar" << std::endl;
00117 }
00118 if ( TurbulenceType() == LES_Smagorinsky ) {
00119 std::cout << "LBM Setup() LES_Smagorinsky : Smagorinsky Constant = " << SmagorinskyConstant()
00120 << std::endl;
00121
00122 }
00123 if ( TurbulenceType() == LES_dynamic ) {
00124 std::cout << "LBM Setup() LES_dynamic " << std::endl;
00125 }
00126 WriteInit();
00127 }
00128
00129 virtual void WriteInit() const {
00130 int me = MY_PROC;
00131 if (me == VizServer) {
00132 std::ofstream ofs("chem.dat", std::ios::out);
00133 ofs << "RU " << Rp << std::endl;
00134 ofs << "PA " << BasePressure() << std::endl;
00135 ofs << "Sp(01) Vicosity" << std::endl;
00136 ofs << "W(01) " << Wp << std::endl;
00137 ofs << "Sp(02) |Velocity|" << std::endl;
00138 ofs << "W(02) " << 1. << std::endl;
00139 ofs << "Sp(03) Vorticity-x" << std::endl;
00140 ofs << "W(03) " << 1. << std::endl;
00141 ofs << "Sp(04) Vorticity-y" << std::endl;
00142 ofs << "W(04) " << 1. << std::endl;
00143 ofs << "Sp(05) Vorticity-z" << std::endl;
00144 ofs << "W(05) " << 1. << std::endl;
00145 ofs << "Sp(06) |Vorticity|" << std::endl;
00146 ofs << "W(06) " << 1. << std::endl;
00147 ofs.close();
00148 }
00149 }
00150
00151 inline virtual MacroType MacroVariables(const MicroType &f) const {
00152 MacroType q;
00153 q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8)+f(9)
00154 +f(10)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18);
00155 q(1)=(f(1)-f(2)+f(7)-f(8)+f(9)-f(10)+f(15)-f(16)+f(17)-f(18))/q(0);
00156 q(2)=(f(3)-f(4)+f(7)-f(8)-f(9)+f(10)+f(11)-f(12)+f(13)-f(14))/q(0);
00157 q(3)=(f(5)-f(6)+f(11)-f(12)-f(13)+f(14)+f(15)-f(16)-f(17)+f(18))/q(0);
00158 return q;
00159 }
00160
00161 inline virtual MicroType Equilibrium(const MacroType &q) const {
00162 MicroType feq;
00163
00164 DataType rho = q(0);
00165 DataType u = q(1);
00166 DataType v = q(2);
00167 DataType w = q(3);
00168
00169 DataType ri, rt0, rt1, rt2;
00170 if (method[0]==0) {
00171 ri = DataType(1.);
00172 rt0 = R0 / DataType(3.);
00173 rt1 = R0 / DataType(18.);
00174 rt2 = R0 / DataType(36.);
00175 }
00176 if (method[0]==1) {
00177 ri = rho;
00178 rt0 = DataType(1.) / DataType(3.);
00179 rt1 = DataType(1.) / DataType(18.);
00180 rt2 = DataType(1.) / DataType(36.);
00181 }
00182 else if (method[0]==2) {
00183 ri = DataType(1.);
00184 rt0 = rho / DataType(3.);
00185 rt1 = rho / DataType(18.);
00186 rt2 = rho / DataType(36.);
00187 }
00188
00189 DataType usq = u * u;
00190 DataType vsq = v * v;
00191 DataType wsq = w * w;
00192 DataType sumsq = (usq + vsq + wsq) / cs22;
00193 DataType sumsqxy = (usq + vsq) / cs22;
00194 DataType sumsqyz = (vsq + wsq) / cs22;
00195 DataType sumsqxz = (usq + wsq) / cs22;
00196 DataType sumsq2xy = sumsqxy * (DataType(1.) - cs2) / cs2;
00197 DataType sumsq2yz = sumsqyz * (DataType(1.) - cs2) / cs2;
00198 DataType sumsq2xz = sumsqxz * (DataType(1.) - cs2) / cs2;
00199 DataType u2 = usq / cssq;
00200 DataType v2 = vsq / cssq;
00201 DataType w2 = wsq / cssq;
00202 DataType u22 = usq / cs22;
00203 DataType v22 = vsq / cs22;
00204 DataType w22 = wsq / cs22;
00205 DataType ui = u / cs2;
00206 DataType vi = v / cs2;
00207 DataType wi = w / cs2;
00208 DataType uv = ui * vi;
00209 DataType uw = ui * wi;
00210 DataType vw = vi * wi;
00211
00212 feq(0) = rt0 * (ri - sumsq);
00213
00214 feq(1) = rt1 * (ri - sumsq + u2 + ui);
00215 feq(2) = rt1 * (ri - sumsq + u2 - ui);
00216
00217 feq(3) = rt1 * (ri - sumsq + v2 + vi);
00218 feq(4) = rt1 * (ri - sumsq + v2 - vi);
00219
00220 feq(5) = rt1 * (ri - sumsq + w2 + wi);
00221 feq(6) = rt1 * (ri - sumsq + w2 - wi);
00222
00223 feq(7) = rt2 * (ri + sumsq2xy -w22 + ui + vi + uv);
00224 feq(8) = rt2 * (ri + sumsq2xy -w22 - ui - vi + uv);
00225 feq(9) = rt2 * (ri + sumsq2xy -w22 + ui - vi - uv);
00226 feq(10) = rt2 * (ri + sumsq2xy -w22 - ui + vi - uv);
00227
00228 feq(11) = rt2 * (ri + sumsq2yz -u22 + vi + wi + vw);
00229 feq(12) = rt2 * (ri + sumsq2yz -u22 - vi - wi + vw);
00230 feq(13) = rt2 * (ri + sumsq2yz -u22 + vi - wi - vw);
00231 feq(14) = rt2 * (ri + sumsq2yz -u22 - vi + wi - vw);
00232
00233 feq(15) = rt2 * (ri + sumsq2xz -v22 + ui + wi + uw);
00234 feq(16) = rt2 * (ri + sumsq2xz -v22 - ui - wi + uw);
00235 feq(17) = rt2 * (ri + sumsq2xz -v22 + ui - wi - uw);
00236 feq(18) = rt2 * (ri + sumsq2xz -v22 - ui + wi - uw);
00237
00238 return feq;
00239 }
00240
00241 inline virtual void Collision(MicroType &f, const DataType dt) const {
00242 MacroType q = MacroVariables(f);
00243 MicroType feq = Equilibrium(q);
00244
00245 DataType omega;
00246 if (turbulence == laminar)
00247 omega = Omega(dt);
00248 else if (turbulence == LES_Smagorinsky)
00249 omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00250 else if (turbulence == LES_dynamic)
00251 omega = Omega_LES_dynamic(f,feq,q,dt);
00252
00253 f=(DataType(1.)-omega)*f + omega*feq;
00254 }
00255
00256 virtual int IncomingIndices(const int side, int indices[]) const {
00257 switch (side) {
00258 case base::Left:
00259 indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;
00260 break;
00261 case base::Right:
00262 indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;
00263 break;
00264 case base::Bottom:
00265 indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;
00266 break;
00267 case base::Top:
00268 indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;
00269 break;
00270 case base::Front:
00271 indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;
00272 break;
00273 case base::Back:
00274 indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;
00275 break;
00276 }
00277 return 5;
00278 }
00279
00280 virtual int OutgoingIndices(const int side, int indices[]) const {
00281 switch (side) {
00282 case base::Left:
00283 indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;
00284 break;
00285 case base::Right:
00286 indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;
00287 break;
00288 case base::Bottom:
00289 indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;
00290 break;
00291 case base::Top:
00292 indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;
00293 break;
00294 case base::Front:
00295 indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;
00296 break;
00297 case base::Back:
00298 indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;
00299 break;
00300 }
00301 return 5;
00302 }
00303
00304
00305 virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb,
00306 const int side) const {
00307 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00308 int b = fvec.bottom();
00309 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00310 && fvec.stepsize(2)==bb.stepsize(2));
00311 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00312 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00313 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00314 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00315 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00316 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00317 MicroType *f = (MicroType *)fvec.databuffer();
00318
00319 #ifdef DEBUG_PRINT_FIXUP
00320 ( comm_service::log() << "Reverse streaming at Side: " << side << " of "
00321 << fvec.bbox() << " using " << bb << "\n").flush();
00322 #endif
00323
00324 switch (side) {
00325 case base::Left:
00326 for (register int k=mzs; k<=mze; k++)
00327 for (register int j=mys; j<=mye; j++) {
00328 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](2);
00329 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](8);
00330 f[b+base::idx(mxs,j,k,Nx,Ny)](10)= f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](10);
00331 f[b+base::idx(mxs,j,k,Nx,Ny)](16)= f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](16);
00332 f[b+base::idx(mxs,j,k,Nx,Ny)](18)= f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](18);
00333 }
00334 break;
00335 case base::Right:
00336 for (register int k=mzs; k<=mze; k++)
00337 for (register int j=mys; j<=mye; j++) {
00338 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](1);
00339 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](7);
00340 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](9);
00341 f[b+base::idx(mxe,j,k,Nx,Ny)](15)= f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](15);
00342 f[b+base::idx(mxe,j,k,Nx,Ny)](17)= f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](17);
00343 }
00344 break;
00345 case base::Bottom:
00346 for (register int k=mzs; k<=mze; k++)
00347 for (register int i=mxs; i<=mxe; i++) {
00348 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](4);
00349 f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](8);
00350 f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](9);
00351 f[b+base::idx(i,mys,k,Nx,Ny)](12)= f[b+base::idx(i,mys-1,k-1,Nx,Ny)](12);
00352 f[b+base::idx(i,mys,k,Nx,Ny)](14)= f[b+base::idx(i,mys-1,k+1,Nx,Ny)](14);
00353 }
00354 break;
00355 case base::Top:
00356 for (register int k=mzs; k<=mze; k++)
00357 for (register int i=mxs; i<=mxe; i++) {
00358 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](3);
00359 f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](7);
00360 f[b+base::idx(i,mye,k,Nx,Ny)](10)= f[b+base::idx(i-1,mye+1,k,Nx,Ny)](10);
00361 f[b+base::idx(i,mye,k,Nx,Ny)](11)= f[b+base::idx(i,mye+1,k+1,Nx,Ny)](11);
00362 f[b+base::idx(i,mye,k,Nx,Ny)](13)= f[b+base::idx(i,mye+1,k-1,Nx,Ny)](13);
00363 }
00364 break;
00365 case base::Front:
00366 for (register int j=mys; j<=mye; j++)
00367 for (register int i=mxs; i<=mxe; i++) {
00368 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](6);
00369 f[b+base::idx(i,j,mzs,Nx,Ny)](12)= f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](12);
00370 f[b+base::idx(i,j,mzs,Nx,Ny)](13)= f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](13);
00371 f[b+base::idx(i,j,mzs,Nx,Ny)](16)= f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](16);
00372 f[b+base::idx(i,j,mzs,Nx,Ny)](17)= f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](17);
00373 }
00374 break;
00375 case base::Back:
00376 for (register int j=mys; j<=mye; j++)
00377 for (register int i=mxs; i<=mxe; i++) {
00378 f[b+base::idx(i,j,mzs,Nx,Ny)](5) = f[b+base::idx(i,j,mzs+1,Nx,Ny)](5);
00379 f[b+base::idx(i,j,mzs,Nx,Ny)](11)= f[b+base::idx(i,j+1,mzs+1,Nx,Ny)](11);
00380 f[b+base::idx(i,j,mzs,Nx,Ny)](14)= f[b+base::idx(i,j-1,mzs+1,Nx,Ny)](14);
00381 f[b+base::idx(i,j,mzs,Nx,Ny)](15)= f[b+base::idx(i+1,j,mzs+1,Nx,Ny)](15);
00382 f[b+base::idx(i,j,mzs,Nx,Ny)](18)= f[b+base::idx(i-1,j,mzs+1,Nx,Ny)](18);
00383 }
00384 break;
00385 }
00386 }
00387
00388 virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec,
00389 const BBox &bb, const double& dt) const {
00390 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00391 int b = fvec.bottom();
00392 MicroType *f = (MicroType *)fvec.databuffer();
00393 MicroType *of = (MicroType *)ovec.databuffer();
00394
00395 #ifdef DEBUG_PRINT_FIXUP
00396 ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox()
00397 << " using " << bb << "\n").flush();
00398 #endif
00399
00400 assert (fvec.extents(0)==ovec.extents(0));
00401 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) && fvec.stepsize(2)==bb.stepsize(2) &&
00402 fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1) && fvec.stepsize(2)==ovec.stepsize(2));
00403 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00404 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00405 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00406 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00407 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00408 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00409
00410 for (register int k=mzs; k<=mze; k++)
00411 for (register int j=mys; j<=mye; j++)
00412 for (register int i=mxs; i<=mxe; i++) {
00413 for (register int n=0; n<base::NMicroVar(); n++)
00414 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);
00415 Collision(f[b+base::idx(i,j,k,Nx,Ny)],dt);
00416 }
00417 }
00418
00419 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00420 const double& t, const double& dt, const int& mpass) const {
00421
00422
00423 DCoords dx = base::GH().worldStep(fvec.stepsize());
00424
00425
00426 if ( (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR)
00427 || (std::fabs(dx(1)/L0()-dx(2)/L0()) > DBL_EPSILON*LB_FACTOR) ) {
00428 std::cerr << "LBMD3Q19 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00429 << " and dx(2)=" << dx(2)/L0()
00430 << " to be equal." << std::endl << std::flush;
00431 std::exit(-1);
00432 }
00433 if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR) {
00434 std::cerr << "LBMD3Q19 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00435 << " to be equal."
00436 << " dt=" << dt << " TO=" << T0()
00437 << std::endl << std::flush;
00438 std::exit(-1);
00439 }
00440
00441 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00442 MicroType *f = (MicroType *)fvec.databuffer();
00443
00444
00445 for (register int k=Nz-1; k>=1; k--)
00446 for (register int j=Ny-1; j>=1; j--)
00447 for (register int i=0; i<=Nx-2; i++) {
00448
00449 f[base::idx(i,j,k,Nx,Ny)](3) = f[base::idx(i,j-1,k,Nx,Ny)](3);
00450 f[base::idx(i,j,k,Nx,Ny)](10) = f[base::idx(i+1,j-1,k,Nx,Ny)](10);
00451
00452
00453
00454 f[base::idx(i,j,k,Nx,Ny)](5) = f[base::idx(i,j,k-1,Nx,Ny)](5);
00455 f[base::idx(i,j,k,Nx,Ny)](18) = f[base::idx(i+1,j,k-1,Nx,Ny)](18);
00456 }
00457
00458 for (register int k=0; k<=Nz-2; k++)
00459 for (register int j=Ny-1; j>=1; j--)
00460 for (register int i=0; i<=Nx-2; i++) {
00461
00462 f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);
00463 }
00464 for (register int k=Nz-1; k>=1; k--)
00465 for (register int j=Ny-1; j>=1; j--)
00466 for (register int i=Nx-1; i>=1; i--) {
00467
00468 f[base::idx(i,j,k,Nx,Ny)](1) = f[base::idx(i-1,j,k,Nx,Ny)](1);
00469 f[base::idx(i,j,k,Nx,Ny)](7) = f[base::idx(i-1,j-1,k,Nx,Ny)](7);
00470
00471 f[base::idx(i,j,k,Nx,Ny)](11) = f[base::idx(i,j-1,k-1,Nx,Ny)](11);
00472
00473 f[base::idx(i,j,k,Nx,Ny)](15) = f[base::idx(i-1,j,k-1,Nx,Ny)](15);
00474 }
00475
00476 for (register int k=0; k<=Nz-2; k++)
00477 for (register int j=0; j<=Ny-2; j++)
00478 for (register int i=Nx-1; i>=1; i--) {
00479
00480 f[base::idx(i,j,k,Nx,Ny)](4) = f[base::idx(i,j+1,k,Nx,Ny)](4);
00481 f[base::idx(i,j,k,Nx,Ny)](9) = f[base::idx(i-1,j+1,k,Nx,Ny)](9);
00482
00483
00484
00485 f[base::idx(i,j,k,Nx,Ny)](6) = f[base::idx(i,j,k+1,Nx,Ny)](6);
00486 f[base::idx(i,j,k,Nx,Ny)](17) = f[base::idx(i-1,j,k+1,Nx,Ny)](17);
00487 }
00488 for (register int k=Nz-1; k>=1; k--)
00489 for (register int j=0; j<=Ny-2; j++)
00490 for (register int i=0; i<=Nx-2; i++) {
00491
00492 f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);
00493 }
00494 for (register int k=0; k<=Nz-2; k++)
00495 for (register int j=0; j<=Ny-2; j++)
00496 for (register int i=0; i<=Nx-2; i++) {
00497
00498 f[base::idx(i,j,k,Nx,Ny)](2) = f[base::idx(i+1,j,k,Nx,Ny)](2);
00499 f[base::idx(i,j,k,Nx,Ny)](8) = f[base::idx(i+1,j+1,k,Nx,Ny)](8);
00500
00501 f[base::idx(i,j,k,Nx,Ny)](12) = f[base::idx(i,j+1,k+1,Nx,Ny)](12);
00502
00503 f[base::idx(i,j,k,Nx,Ny)](16) = f[base::idx(i+1,j,k+1,Nx,Ny)](16);
00504 }
00505
00506
00507 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00508 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00509 for (register int k=mzs; k<=mze; k++)
00510 for (register int j=mys; j<=mye; j++)
00511 for (register int i=mxs; i<=mxe; i++)
00512 Collision(f[base::idx(i,j,k,Nx,Ny)],dt);
00513
00514 return DataType(1.);
00515 }
00516
00517 virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00518 DataType* aux=0, const int naux=0, const int scaling=0) const {
00519 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00520 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00521 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00522 MicroType *f = (MicroType *)fvec.databuffer();
00523
00524
00525
00526 if (type==ConstantMicro) {
00527
00528 MicroType g;
00529 for (register int i=0; i<base::NMicroVar(); i++)
00530 g(i) = aux[i];
00531
00532 for (register int k=mzs; k<=mze; k++)
00533 for (register int j=mys; j<=mye; j++)
00534 for (register int i=mxs; i<=mxe; i++)
00535 f[base::idx(i,j,k,Nx,Ny)] = g;
00536
00537 }
00538
00539 else if (type==ConstantMacro) {
00540
00541 MacroType q;
00542 if (scaling==base::Physical) {
00543 q(0) = aux[0]/R0;
00544 q(1) = aux[1]/U0;
00545 q(2) = aux[2]/U0;
00546 q(3) = aux[3]/U0;
00547 }
00548 else
00549 for (register int i=0; i<base::NMacroVar(); i++)
00550 q(i) = aux[i];
00551 q(1) *= S0; q(2) *= S0; q(3) *= S0;
00552
00553
00554
00555
00556
00557
00558
00559 for (register int k=mzs; k<=mze; k++)
00560 for (register int j=mys; j<=mye; j++)
00561 for (register int i=mxs; i<=mxe; i++)
00562 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00563
00564 }
00565
00566 else if (type==GasAtRest) {
00567
00568 MacroType q;
00569 q(0) = GasDensity()/R0;
00570 q(1) = DataType(0.);
00571 q(2) = DataType(0.);
00572 q(3) = DataType(0.);
00573 for (register int k=mzs; k<=mze; k++)
00574 for (register int j=mys; j<=mye; j++)
00575 for (register int i=mxs; i<=mxe; i++)
00576 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00577
00578 }
00579 }
00580
00581
00582 virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00583 DataType* aux=0, const int naux=0, const int scaling=0) const {
00584 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00585 int b = fvec.bottom();
00586 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00587 && fvec.stepsize(2)==bb.stepsize(2));
00588 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00589 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00590 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00591 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00592 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00593 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00594 MicroType *f = (MicroType *)fvec.databuffer();
00595
00596
00597 if (type==Symmetry) {
00598 switch (side) {
00599 case base::Left:
00600 for (register int k=mzs; k<=mze; k++)
00601 for (register int j=mys; j<=mye; j++) {
00602 if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
00603 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00604 if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
00605
00606 if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
00607 if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
00608 }
00609 break;
00610 case base::Right:
00611 for (register int k=mzs; k<=mze; k++)
00612 for (register int j=mys; j<=mye; j++) {
00613 if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
00614 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00615 if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
00616
00617 if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
00618 if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
00619 }
00620 break;
00621 case base::Bottom:
00622 for (register int k=mzs; k<=mze; k++)
00623 for (register int i=mxs; i<=mxe; i++) {
00624 if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
00625 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00626 if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
00627
00628 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);
00629 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);
00630 }
00631 break;
00632 case base::Top:
00633 for (register int k=mzs; k<=mze; k++)
00634 for (register int i=mxs; i<=mxe; i++) {
00635 if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
00636 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00637 if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
00638
00639 if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
00640 if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
00641 }
00642 break;
00643 case base::Front:
00644
00645 for (register int j=mys; j<=mye; j++)
00646 for (register int i=mxs; i<=mxe; i++) {
00647 if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
00648 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00649 if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
00650
00651 if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
00652 if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
00653 }
00654 break;
00655 case base::Back:
00656
00657 for (register int j=mys; j<=mye; j++)
00658 for (register int i=mxs; i<=mxe; i++) {
00659
00660
00661
00662
00663 if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
00664 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00665 if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
00666
00667 if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
00668 if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
00669 }
00670 break;
00671 }
00672 }
00673
00674 else if (type==SlipWall) {
00675 switch (side) {
00676 case base::Left:
00677
00678 for (register int k=mzs; k<=mze; k++)
00679 for (register int j=mys; j<=mye; j++) {
00680 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);
00681 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00682 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);
00683
00684 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);
00685 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);
00686 }
00687 break;
00688 case base::Right:
00689
00690 for (register int k=mzs; k<=mze; k++)
00691 for (register int j=mys; j<=mye; j++) {
00692 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);
00693 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00694 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);
00695
00696 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);
00697 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);
00698 }
00699 break;
00700 case base::Bottom:
00701
00702 for (register int k=mzs; k<=mze; k++)
00703 for (register int i=mxs; i<=mxe; i++) {
00704 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);
00705 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00706 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);
00707
00708 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);
00709 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);
00710 }
00711 break;
00712 case base::Top:
00713
00714 for (register int k=mzs; k<=mze; k++)
00715 for (register int i=mxs; i<=mxe; i++) {
00716 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);
00717 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00718 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);
00719
00720 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);
00721 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);
00722 }
00723 break;
00724 case base::Front:
00725
00726 for (register int j=mys; j<=mye; j++)
00727 for (register int i=mxs; i<=mxe; i++) {
00728 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);
00729 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00730 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);
00731
00732 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);
00733 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);
00734 }
00735 break;
00736 case base::Back:
00737
00738 for (register int j=mys; j<=mye; j++)
00739 for (register int i=mxs; i<=mxe; i++) {
00740
00741
00742
00743
00744 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);
00745 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00746 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);
00747
00748 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);
00749 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);
00750 }
00751 break;
00752 }
00753 }
00754
00755 else if (type==NoSlipWall) {
00756 switch (side) {
00757 case base::Left:
00758
00759 for (register int k=mzs; k<=mze; k++)
00760 for (register int j=mys; j<=mye; j++) {
00761 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);
00762 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00763 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);
00764
00765 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);
00766 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);
00767 }
00768 break;
00769 case base::Right:
00770
00771 for (register int k=mzs; k<=mze; k++)
00772 for (register int j=mys; j<=mye; j++) {
00773 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);
00774 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00775 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);
00776
00777 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);
00778 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);
00779 }
00780 break;
00781 case base::Bottom:
00782
00783 for (register int k=mzs; k<=mze; k++)
00784 for (register int i=mxs; i<=mxe; i++) {
00785 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);
00786 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00787 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);
00788
00789 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);
00790 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);
00791 }
00792 break;
00793 case base::Top:
00794
00795 for (register int k=mzs; k<=mze; k++)
00796 for (register int i=mxs; i<=mxe; i++) {
00797 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);
00798 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00799 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);
00800
00801 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);
00802 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);
00803 }
00804 break;
00805 case base::Front:
00806
00807 for (register int j=mys; j<=mye; j++)
00808 for (register int i=mxs; i<=mxe; i++) {
00809 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);
00810 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00811 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);
00812
00813 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);
00814 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);
00815 }
00816 break;
00817 case base::Back:
00818
00819 for (register int j=mys; j<=mye; j++)
00820 for (register int i=mxs; i<=mxe; i++) {
00821
00822
00823
00824
00825 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);
00826 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00827 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);
00828
00829 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);
00830 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);
00831 }
00832 break;
00833 }
00834 }
00835
00836 else if (type==Inlet) {
00837 if (aux==0 || naux<=0) return;
00838 DataType a0=S0*aux[0], a1=0., a2=0.;
00839 if (naux>1) {
00840 a1 = S0*aux[1];
00841 a2 = S0*aux[2];
00842 }
00843 if (scaling==base::Physical) {
00844 a0 /= U0;
00845 a1 /= U0;
00846 a2 /= U0;
00847 }
00848 switch (side) {
00849 case base::Left:
00850 for (register int k=mzs; k<=mze; k++)
00851 for (register int j=mys; j<=mye; j++) {
00852 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00853 q(1) = a0;
00854 if (naux>1) {
00855 q(2) = a1;
00856 q(3) = a2;
00857 }
00858 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
00859 }
00860 break;
00861 case base::Right:
00862 for (register int k=mzs; k<=mze; k++)
00863 for (register int j=mys; j<=mye; j++) {
00864 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00865 q(1) = a0;
00866 if (naux>1) {
00867 q(2) = a1;
00868 q(3) = a2;
00869 }
00870 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
00871 }
00872 break;
00873 case base::Bottom:
00874 for (register int k=mzs; k<=mze; k++)
00875 for (register int i=mxs; i<=mxe; i++) {
00876 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
00877 if (naux==1)
00878 q(2) = a0;
00879 else {
00880 q(1) = a0;
00881 q(2) = a1;
00882 q(3) = a2;
00883 }
00884 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
00885 }
00886 break;
00887 case base::Top:
00888 for (register int k=mzs; k<=mze; k++)
00889 for (register int i=mxs; i<=mxe; i++) {
00890 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
00891 if (naux==1)
00892 q(2) = a0;
00893 else {
00894 q(1) = a0;
00895 q(2) = a1;
00896 q(3) = a2;
00897 }
00898 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
00899 }
00900 break;
00901 case base::Front:
00902 for (register int j=mys; j<=mye; j++)
00903 for (register int i=mxs; i<=mxe; i++) {
00904 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
00905 if (naux==1)
00906 q(3) = a0;
00907 else {
00908 q(1) = a1;
00909 q(2) = a2;
00910 q(3) = a0;
00911 }
00912 f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
00913 }
00914 break;
00915 case base::Back:
00916 for (register int j=mys; j<=mye; j++)
00917 for (register int i=mxs; i<=mxe; i++) {
00918 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
00919 if (naux==1)
00920 q(3) = a0;
00921 else {
00922 q(1) = a1;
00923 q(2) = a2;
00924 q(3) = a0;
00925 }
00926 f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
00927 }
00928 break;
00929 }
00930 }
00931
00932 else if (type==Outlet) {
00933
00934 switch (side) {
00935 case base::Left:
00936
00937 for (register int k=mzs; k<=mze; k++)
00938 for (register int j=mys; j<=mye; j++) {
00939 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00940
00941 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);
00942 }
00943 break;
00944 case base::Right:
00945
00946 for (register int k=mzs; k<=mze; k++)
00947 for (register int j=mys; j<=mye; j++) {
00948 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00949
00950 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);
00951 }
00952 break;
00953 case base::Bottom:
00954 for (register int k=mzs; k<=mze; k++)
00955 for (register int i=mxs; i<=mxe; i++) {
00956 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
00957 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);
00958 }
00959 break;
00960 case base::Top:
00961 for (register int k=mzs; k<=mze; k++)
00962 for (register int i=mxs; i<=mxe; i++) {
00963 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
00964 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);
00965 }
00966 break;
00967 case base::Front:
00968 for (register int j=mys; j<=mye; j++)
00969 for (register int i=mxs; i<=mxe; i++) {
00970 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
00971 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);
00972 }
00973 break;
00974 case base::Back:
00975 for (register int j=mys; j<=mye; j++)
00976 for (register int i=mxs; i<=mxe; i++) {
00977 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
00978 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);
00979 }
00980 break;
00981 }
00982 }
00983
00984
00985 else if (type==Pressure) {
00986 if (aux==0 || naux<=0) return;
00987 DataType a0=aux[0];
00988 if (scaling==base::Physical)
00989 a0 /= R0;
00990 switch (side) {
00991 case base::Left:
00992 for (register int k=mzs; k<=mze; k++)
00993 for (register int j=mys; j<=mye; j++) {
00994 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00995 q(0) = a0;
00996 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
00997 }
00998 break;
00999 case base::Right:
01000 for (register int k=mzs; k<=mze; k++)
01001 for (register int j=mys; j<=mye; j++) {
01002 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01003 q(0) = a0;
01004 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
01005 }
01006 break;
01007 case base::Bottom:
01008 for (register int k=mzs; k<=mze; k++)
01009 for (register int i=mxs; i<=mxe; i++) {
01010 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01011 q(0) = a0;
01012 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
01013 }
01014 break;
01015 case base::Top:
01016 for (register int k=mzs; k<=mze; k++)
01017 for (register int i=mxs; i<=mxe; i++) {
01018 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01019 q(0) = a0;
01020 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
01021 }
01022 break;
01023 case base::Front:
01024 for (register int j=mys; j<=mye; j++)
01025 for (register int i=mxs; i<=mxe; i++) {
01026 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01027 q(0) = a0;
01028 f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
01029 }
01030 break;
01031 case base::Back:
01032 for (register int j=mys; j<=mye; j++)
01033 for (register int i=mxs; i<=mxe; i++) {
01034 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01035 q(0) = a0;
01036 f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
01037 }
01038 break;
01039 }
01040 }
01041
01042 else if (type==SlidingWall) {
01043 if (aux==0 || naux<=0) return;
01044 if (aux==0 || naux<=0) return;
01045 DataType a0=S0*aux[0];
01046 if (scaling==base::Physical)
01047 a0 /= U0;
01048 switch (side) {
01049 case base::Left:
01050 for (register int k=mzs; k<=mze; k++)
01051 for (register int j=mys; j<=mye; j++) {
01052 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01053 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);
01054 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01055 DataType pw = DataType(1.)-qw;
01056 if (j<mye) f[b+base::idx(mxe,j+1,k,Nx,Ny)](9) = pw*rd1+qw*rd2;
01057 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01058 if (j>mys) f[b+base::idx(mxe,j-1,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01059
01060 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);
01061 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);
01062 }
01063 break;
01064 case base::Right:
01065 for (register int k=mzs; k<=mze; k++)
01066 for (register int j=mys; j<=mye; j++) {
01067 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01068 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);
01069 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01070 DataType pw = DataType(1.)-qw;
01071 if (j<mye) f[b+base::idx(mxs,j+1,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01072 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01073 if (j>mys) f[b+base::idx(mxs,j-1,k,Nx,Ny)](10) = qw*rd1+pw*rd2;
01074
01075 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);
01076 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);
01077 }
01078 break;
01079 case base::Bottom:
01080 for (register int k=mzs; k<=mze; k++)
01081 for (register int i=mxs; i<=mxe; i++) {
01082 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01083 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);
01084 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01085 DataType pw = DataType(1.)-qw;
01086 if (i<mxe) f[b+base::idx(i+1,mye,k,Nx,Ny)](10) = pw*rd1+qw*rd2;
01087 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01088 if (i>mxs) f[b+base::idx(i-1,mye,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01089
01090 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);
01091 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);
01092 }
01093 break;
01094 case base::Top:
01095 for (register int k=mzs; k<=mze; k++)
01096 for (register int i=mxs; i<=mxe; i++) {
01097 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01098 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);
01099 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01100 DataType pw = DataType(1.)-qw;
01101 if (i<mxe) f[b+base::idx(i+1,mys,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01102 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01103 if (i>mxs) f[b+base::idx(i-1,mys,k,Nx,Ny)](9) = qw*rd1+pw*rd2;
01104
01105 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);
01106 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);
01107 }
01108 break;
01109 case base::Front:
01110 for (register int j=mys; j<=mye; j++)
01111 for (register int i=mxs; i<=mxe; i++) {
01112 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01113 DataType rd1 = f[b+base::idx(i,j,mze+1,Nx,Ny)](17), rd2 = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
01114 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01115 DataType pw = DataType(1.)-qw;
01116 if (i<mxe) f[b+base::idx(i+1,j,mze,Nx,Ny)](18) = pw*rd1+qw*rd2;
01117 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01118 if (i>mxs) f[b+base::idx(i-1,j,mze,Nx,Ny)](15) = qw*rd1+pw*rd2;
01119
01120 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);
01121 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);
01122 }
01123 break;
01124 case base::Back:
01125 for (register int j=mys; j<=mye; j++)
01126 for (register int i=mxs; i<=mxe; i++) {
01127 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01128 DataType rd1 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15), rd2 = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
01129 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01130 DataType pw = DataType(1.)-qw;
01131 if (i<mxe) f[b+base::idx(i+1,j,mzs,Nx,Ny)](16) = pw*rd1+qw*rd2;
01132 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
01133 if (i>mxs) f[b+base::idx(i-1,j,mzs,Nx,Ny)](17) = qw*rd1+pw*rd2;
01134
01135 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);
01136 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);
01137 }
01138 break;
01139 }
01140 }
01141 }
01142
01143 virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
01144 const MicroType* f, const point_type* xc, const DataType* distance,
01145 const point_type* normal, DataType* aux=0, const int naux=0,
01146 const int scaling=0) const {
01147 if (type==GFMNoSlipWall || type==GFMSlipWall) {
01148 DataType fac = S0;
01149 if (scaling==base::Physical)
01150 fac /= U0;
01151 for (int n=0; n<nc; n++) {
01152 MacroType q=MacroVariables(f[n]);
01153 DataType u=-q(1);
01154 DataType v=-q(2);
01155 DataType w=-q(3);
01156
01157
01158 if (naux>=3) {
01159 u += fac*aux[naux*n];
01160 v += fac*aux[naux*n+1];
01161 w += fac*aux[naux*n+2];
01162 }
01163 if (type==GFMNoSlipWall) {
01164
01165 q(1) += DataType(2.)*u;
01166 q(2) += DataType(2.)*v;
01167 q(3) += DataType(2.)*w;
01168 }
01169 else if (type==GFMSlipWall) {
01170
01171 DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1)+w*normal[n](2));
01172 q(1) += vl*normal[n](0);
01173 q(2) += vl*normal[n](1);
01174 q(3) += vl*normal[n](2);
01175 }
01176 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01177 }
01178 }
01179
01180 else if (type==GFMExtrapolation) {
01181 for (int n=0; n<nc; n++) {
01182 MacroType q=MacroVariables(f[n]);
01183 DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1)+q(3)*normal[n](2);
01184 q(1) = vl*normal[n](0);
01185 q(2) = vl*normal[n](1);
01186 q(3) = vl*normal[n](2);
01187 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01188 }
01189 }
01190
01191 else if (type==GFMWallLaw ) {
01192 DCoords dx = base::GH().worldStep(fvec.stepsize());
01193 DataType fac = S0;
01194 if (scaling==base::Physical)
01195 fac /= U0;
01196
01197
01198 DataType yp_ = 0., yp = 0., yMax = 50.;
01199 DataType yPlus[3];
01200 int lc_low[3], lc_up[3], yPlus_lc[3];
01201
01202
01203
01204
01205
01206 DataType shear_vel = 1., karman = 0.41, y0 = 1., law = 1.;
01207
01208 DataType shearStress;
01209 DataType u, v, w, u_, v_, w_, vel0;
01210
01211 MacroType q, q_;
01212 int kMax = 0;
01213
01214 for (int n=0; n<nc; n++) {
01215
01217 q=MacroVariables(f[n]);
01218 q_=MacroVariables(f[n]);
01219
01220 yPlus[0] = (xc[n](0) - DataType(2.)*dx(0)*normal[n](0));
01221 yPlus[1] = (xc[n](1) - DataType(2.)*dx(1)*normal[n](1));
01222 yPlus[2] = (xc[n](2) - DataType(2.)*dx(2)*normal[n](2));
01223
01224 yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
01225 yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
01226 yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
01227
01228 lc_low[0] = fvec.lower(0); lc_low[1] = fvec.lower(1); lc_low[2] = fvec.lower(2);
01229 lc_up[0] = fvec.upper(0); lc_up[1] = fvec.upper(1); lc_up[2] = fvec.upper(2);
01230
01231
01232
01233 if ( yPlus_lc[0] > lc_low[0] ) {
01234 if ( yPlus_lc[1] > lc_low[1] ) {
01235 if ( yPlus_lc[2] > lc_low[2] ) {
01236
01237 if ( yPlus_lc[0] < lc_up[0] ) {
01238 if ( yPlus_lc[1] < lc_up[1] ) {
01239 if ( yPlus_lc[2] < lc_up[2] ) {
01240 q=MacroVariables( fvec( base::GH().localCoords((const DataType *) yPlus) ) );
01241 }
01242 }
01243 }
01244
01245 }
01246 }
01247 }
01248
01249 u_=-q(1);
01250 v_=-q(2);
01251 w_=-q(3);
01252
01253
01254 q=MacroVariables(f[n]);
01255 u=-q(1);
01256 v=-q(2);
01257 w=-q(3);
01258 vel0 = sqrt(u*u+v*v+w*w);
01259
01260 u_+=q(1);
01261 v_+=q(2);
01262 w_+=q(3);
01263 vel0 = ( u_*normal[n](0) + v_*normal[n](1) + w_*normal[n](2) );
01264
01265 u_ = u_ - vel0 * normal[n](0);
01266 v_ = v_ - vel0 * normal[n](1);
01267 w_ = w_ - vel0 * normal[n](2);
01268 vel0 = sqrt(u_*u_ + v_*v_ + w_*w_)*U0/S0;
01269
01270 shearStress = (DataType(0.5)*vel0/dx(0) * q(0)*R0*GasViscosity() );
01271
01272
01273 if (naux>=3) {
01274 u += fac*aux[naux*n];
01275 v += fac*aux[naux*n+1];
01276 w += fac*aux[naux*n+2];
01277 }
01278 q(1) += DataType(2.)*u;
01279 q(2) += DataType(2.)*v;
01280 q(3) += DataType(2.)*w;
01281 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
01282
01284 shear_vel = std::sqrt( shearStress/(q(0)*R0) );
01285
01287 y0 = karman/30.;
01288
01289
01290
01291 if ( shear_vel > 1.e-5 ) {
01292 kMax = (int) std::ceil( (yMax*GasViscosity())/(dx(0)*shear_vel) );
01293
01294 for (int k=1; k < kMax; k++) {
01295
01296 yPlus[0] = (xc[n](0) - k*dx(0)*normal[n](0));
01297 yPlus[1] = (xc[n](1) - k*dx(1)*normal[n](1));
01298 yPlus[2] = (xc[n](2) - k*dx(2)*normal[n](2));
01299
01300 yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
01301 yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
01302 yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
01303
01304 if ( yPlus_lc[0] > lc_low[0] ) {
01305 if ( yPlus_lc[1] > lc_low[1] ) {
01306 if ( yPlus_lc[2] > lc_low[2] ) {
01307
01308 if ( yPlus_lc[0] < lc_up[0] ) {
01309 if ( yPlus_lc[1] < lc_up[1] ) {
01310 if ( yPlus_lc[2] < lc_up[2] ) {
01311 yp = k*dx(0);
01312 yp_ = yp*shear_vel/GasViscosity();
01313
01316 if ( yp_ <= 5.0 ) {
01317 law = yp_;
01318 }
01320 else if (yp_ > 5 && yp_ < 30.0 ) {
01321
01322 law = std::min( yp_ , shear_vel/karman*std::log( yp/y0 ) );
01323
01324 }
01326 else if (yp_ >= 30.0 ) {
01327 law = shear_vel/karman*std::log( yp/y0 );
01328 }
01329
01330 q_=MacroVariables(f[n]);
01331 u = law*(-q_(1) + fac*aux[naux*n]);
01332 v = law*(-q_(2) + fac*aux[naux*n+1]);
01333 w = law*(-q_(3) + fac*aux[naux*n+2]);
01334
01335 q_(1) += DataType(2.)*u;
01336 q_(2) += DataType(2.)*v;
01337 q_(3) += DataType(2.)*w;
01338 fvec( base::GH().localCoords((const DataType *) yPlus) ) = Equilibrium(q_);
01339 }
01340 }
01341 }
01342
01343
01344 }
01345 }
01346 }
01347
01348 }
01349
01350 }
01351
01352 }
01353
01354 }
01355
01356
01357
01358
01359 else if (type==GFMBounceBack ) {
01360
01361 DCoords dx = base::GH().worldStep(fvec.stepsize());
01362 DataType fac = S0;
01363 if(scaling==base::Physical) fac/=U0;
01364
01365
01366
01367
01368
01369
01370 DataType distB[3], distBScalar, cosa, leng, d_dist, qNorm, mdx2, mdy2, mdz2, zaehler, nenner, lengXi;
01371
01372 for (int n=0; n<nc; n++){
01373 distB[0] = (double &)distance[n]*normal[n](0);
01374 distB[1] = (double &)distance[n]*normal[n](1);
01375 distB[2] = (double &)distance[n]*normal[n](2);
01376
01377
01378
01379
01380
01381
01382 for(int ind=1; ind<base::NMicroVar(); ind++){
01383
01384 mdx2 = mdx[ind]*mdx[ind];
01385 mdy2 = mdy[ind]*mdy[ind];
01386 mdz2 = mdz[ind]*mdz[ind];
01387 distBScalar = std::sqrt(std::pow(distB[0],2)+std::pow(distB[1],2)+std::pow(distB[2],2));
01388 leng = distBScalar*distBScalar*std::sqrt(mdx2+mdy2+mdz2)/(distB[0]*mdx[ind]+distB[1]*mdy[ind]+distB[2]*mdz[ind]);
01389 lengXi = std::sqrt(mdx2*dx(0)*dx(0)+mdy2*dx(1)*dx(1)+mdz2*dx(2)*dx(2));
01390 d_dist = lengXi-leng;
01391 qNorm = d_dist/lengXi;
01392
01393 if(qNorm > 0 && qNorm < 1.){
01394
01395
01396
01397 int idxA, idxB;
01398 if(ind==0.) idxA=0.;
01399 else if(ind % 2 == 0.) idxA = ind-1.;
01400 else idxA = ind+1.;
01401
01402 idxB=idxA;
01403 if(qNorm>=DataType(0.5)) idxB=ind;
01404
01405
01406 DataType cn[3], cnA[3], cnB[3];
01407 cn[0] = xc[n](0);
01408 cn[1] = xc[n](1);
01409 cn[2] = xc[n](2);
01410
01411
01412 cnA[0] = cn[0] + mdx[ind]*dx(0);
01413 cnA[1] = cn[1] + mdy[ind]*dx(1);
01414 cnA[2] = cn[2] + mdz[ind]*dx(2);
01415
01416
01417 cnB[0] = cnA[0];
01418 cnB[1] = cnA[1];
01419 cnB[2] = cnA[2];
01420
01421 if(qNorm<DataType(0.5)){
01422 cnB[0] = cn[0] + DataType(2.)*mdx[ind]*dx(0);
01423 cnB[1] = cn[1] + DataType(2.)*mdy[ind]*dx(1);
01424 cnB[2] = cn[2] + DataType(2.)*mdz[ind]*dx(2);
01425 }
01426
01427
01428
01429
01430
01431
01432
01433
01434
01435 DataType distFuncA, distFuncB;
01436
01437
01438
01439
01440 int cnA_lc[3], cnB_lc[3];
01441 cnA_lc[0] = base::GH().localCoords((const DataType *) cnA)(0);
01442 cnA_lc[1] = base::GH().localCoords((const DataType *) cnA)(1);
01443 cnA_lc[2] = base::GH().localCoords((const DataType *) cnA)(2);
01444 cnB_lc[0] = base::GH().localCoords((const DataType *) cnB)(0);
01445 cnB_lc[1] = base::GH().localCoords((const DataType *) cnB)(1);
01446 cnB_lc[2] = base::GH().localCoords((const DataType *) cnB)(2);
01447
01448 int lc_low[3], lc_up[3];
01449 lc_low[0] = fvec.lower(0); lc_low[1] = fvec.lower(1); lc_low[2] = fvec.lower(2);
01450 lc_up[0] = fvec.upper(0); lc_up[1] = fvec.upper(1); lc_up[2] = fvec.upper(2);
01451
01452 if ( (cnA_lc[0] > lc_low[0]) && (cnB_lc[0] > lc_low[0]) ) {
01453 if ( (cnA_lc[1] > lc_low[1]) && (cnB_lc[1] > lc_low[1]) ) {
01454 if ( (cnA_lc[2] > lc_low[2]) && (cnB_lc[2] > lc_low[2]) ) {
01455
01456 if ( (cnA_lc[0] < lc_up[0]) && (cnB_lc[0] < lc_up[0]) ) {
01457 if ( (cnA_lc[1] < lc_up[1]) && (cnB_lc[1] < lc_up[1]) ) {
01458 if ( (cnA_lc[2] < lc_up[2]) && (cnB_lc[2] < lc_up[2]) ) {
01459
01460 distFuncA = fvec( base::GH().localCoords((const DataType *) cnA) )(idxA);
01461 distFuncB = fvec( base::GH().localCoords((const DataType *) cnB) )(idxB);
01462
01463 DataType intp1,intp2,intp,q2;
01464 q2 = DataType(2.)*qNorm;
01465
01466
01467 if(qNorm<DataType(0.5)){
01468 intp1=q2*distFuncA;
01469 intp2=(DataType(1.)-q2)*distFuncB;
01470 }
01471 else{
01472 if(q2!=0){
01473 intp1 = (DataType(1.)/q2)*distFuncA;
01474 intp2 = ((q2-DataType(1.))/q2)*distFuncB;
01475 }
01476 else{
01477 intp1 = DataType(0.);
01478
01479 }
01480 }
01481
01482 intp = intp1+intp2;
01483
01484 fvec( base::GH().localCoords((const DataType *) cn) )(ind) = intp;
01485
01486 }
01487 }
01488 }
01489
01490 }
01491 }
01492 }
01493
01494
01495
01496
01497 if(naux>=3){
01498 DataType u = fac*aux[naux*n];
01499 DataType v = fac*aux[naux*n+1];
01500 DataType w = fac*aux[naux*n+2];
01501
01502 DataType weight;
01503 weight = DataType(1.)/DataType(12.);
01504 if(ind==0) weight *= 0;
01505 else if(ind>=1 && ind<=6) weight *= DataType(2.);
01506
01507 DataType velForce;
01508
01509 if(qNorm<DataType(0.5)) velForce = DataType(2.)*weight;
01510 else velForce = (DataType(1.)/qNorm)*weight;
01511
01512 velForce *= ((mdx[ind]*u)+(mdy[ind]*v)+(mdz[ind]*w));
01513 fvec( base::GH().localCoords((const DataType*) cn) )(ind) +=velForce;
01514 }
01515
01516 }
01517 }
01518 }
01519
01520 }
01521
01522 }
01523
01524 virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01525 const int skip_ghosts=1) const {
01526 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
01527 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
01528 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
01529 MicroType *f = (MicroType *)fvec.databuffer();
01530 DataType *work = (DataType *)workvec.databuffer();
01531 DCoords dx = base::GH().worldStep(fvec.stepsize());
01532 DataType dt = dx(0)*S0/U0;
01533
01534 if ((cnt>=1 && cnt<=10) || (cnt>=19 && cnt<=24) || (cnt>=28 && cnt<=33)) {
01535 for (register int k=mzs; k<=mze; k++)
01536 for (register int j=mys; j<=mye; j++)
01537 for (register int i=mxs; i<=mxe; i++) {
01538 MacroType q=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
01539 switch(cnt) {
01540 case 1:
01541 if (method[0]==0) work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0+rhop;
01542 else
01543 work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0;
01544 break;
01545 case 2:
01546 work[base::idx(i,j,k,Nx,Ny)]=q(1)*U0/S0;
01547 break;
01548 case 3:
01549 work[base::idx(i,j,k,Nx,Ny)]=q(2)*U0/S0;
01550 break;
01551 case 4:
01552 work[base::idx(i,j,k,Nx,Ny)]=q(3)*U0/S0;
01553 break;
01554 case 6:
01555 work[base::idx(i,j,k,Nx,Ny)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure());
01556 break;
01557 case 7:
01558 work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure();
01559 break;
01560 case 9: {
01561 MicroType feq = Equilibrium(q);
01562 work[base::idx(i,j,k,Nx,Ny)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt),
01563 GasSpeedofSound(),dt);
01564 break;
01565 }
01566 case 10:
01567 work[base::idx(i,j,k,Nx,Ny)]=std::sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3))*U0/S0;
01568 break;
01569 case 19:
01570 case 20:
01571 case 21:
01572 case 22:
01573 case 23:
01574 case 24: {
01575 MicroType feq = Equilibrium(q);
01576 DataType omega;
01577 if (turbulence == laminar)
01578 omega = Omega(dt);
01579 else if (turbulence == LES_Smagorinsky)
01580 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01581 else if (turbulence == LES_dynamic)
01582 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01583 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);
01584 break;
01585 }
01586 case 28:
01587 case 29:
01588 case 30:
01589 case 31:
01590 case 32:
01591 case 33: {
01592 DataType omega;
01593 if (turbulence == laminar)
01594 omega = Omega(dt);
01595 else if (turbulence == LES_Smagorinsky) {
01596 MicroType feq = Equilibrium(q);
01597 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01598 }
01599 else if (turbulence == LES_dynamic) {
01600 MicroType feq = Equilibrium(q);
01601 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
01602 }
01603 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);
01604 if (cnt==28 || cnt==30 || cnt==33) work[base::idx(i,j,k,Nx,Ny)]-=BasePressure();
01605 break;
01606 }
01607 }
01608 }
01609 }
01610 else if ((cnt>=11 && cnt<=14) || cnt==18 || (cnt>=40 && cnt<=45)) {
01611 macro_grid_data_type qgrid(fvec.bbox());
01612 MacroType *q = (MacroType *)qgrid.databuffer();
01613 for (register int k=mzs; k<=mze; k++)
01614 for (register int j=mys; j<=mye; j++)
01615 for (register int i=mxs; i<=mxe; i++) {
01616 q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
01617 q[base::idx(i,j,k,Nx,Ny)](0)*=R0;
01618 q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
01619 q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
01620 q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
01621 }
01622
01623 if (skip_ghosts==0) {
01624 switch(cnt) {
01625 case 40:
01626 mxs+=1; mxe-=1;
01627 break;
01628 case 42:
01629 mys+=1; mye-=1;
01630 break;
01631 case 45:
01632 mzs+=1; mze-=1;
01633 break;
01634 case 11:
01635 case 44:
01636 mys+=1; mye-=1; mzs+=1; mze-=1;
01637 break;
01638 case 12:
01639 case 43:
01640 mxs+=1; mxe-=1; mzs+=1; mze-=1;
01641 break;
01642 case 13:
01643 case 41:
01644 mxs+=1; mxe-=1; mys+=1; mye-=1;
01645 break;
01646 case 14:
01647 case 18:
01648 mxs+=1; mxe-=1; mys+=1; mye-=1; mzs+=1; mze-=1;
01649 break;
01650 }
01651 }
01652
01653 DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz;
01654 for (register int k=mzs; k<=mze; k++)
01655 for (register int j=mys; j<=mye; j++)
01656 for (register int i=mxs; i<=mxe; i++) {
01657 if (cnt==18 || cnt==40)
01658 dxux=(q[base::idx(i+1,j,k,Nx,Ny)](1)-q[base::idx(i-1,j,k,Nx,Ny)](1))/(2.*dx(0));
01659 if (cnt==18 || cnt==42)
01660 dyuy=(q[base::idx(i,j+1,k,Nx,Ny)](2)-q[base::idx(i,j-1,k,Nx,Ny)](2))/(2.*dx(1));
01661 if (cnt==18 || cnt==45)
01662 dzuz=(q[base::idx(i,j,k+1,Nx,Ny)](3)-q[base::idx(i,j,k-1,Nx,Ny)](3))/(2.*dx(2));
01663 if (cnt==11 || cnt==14 || cnt==18 || cnt==44) {
01664 dyuz=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(2.*dx(1));
01665 dzuy=(q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(2.*dx(2));
01666 }
01667 if (cnt==12 || cnt==14 || cnt==18 || cnt==43) {
01668 dxuz=(q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(2.*dx(0));
01669 dzux=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(2.*dx(2));
01670 }
01671 if (cnt==13 || cnt==14 || cnt==18 || cnt==41) {
01672 dxuy=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(2.*dx(0));
01673 dyux=(q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(2.*dx(1));
01674 }
01675 switch(cnt) {
01676 case 11:
01677 work[base::idx(i,j,k,Nx,Ny)]=dyuz-dzuy;
01678 break;
01679 case 12:
01680 work[base::idx(i,j,k,Nx,Ny)]=dzux-dxuz;
01681 break;
01682 case 13:
01683 work[base::idx(i,j,k,Nx,Ny)]=dxuy-dyux;
01684 break;
01685 case 14: {
01686 DataType ox = dyuz-dzuy, oy = dzux-dxuz, oz = dxuy-dyux;
01687 work[base::idx(i,j,k,Nx,Ny)] = std::sqrt(ox*ox+oy*oy+oz*oz);
01688 break;
01689 }
01690 case 18:
01691 work[base::idx(i,j,k,Nx,Ny)] = DataType(-0.5)*(4.*dxux*dxux+4.*dyuy*dyuy+4.*dzuz*dzuz+4.*dxuy*dyux+4.*dxuz*dzux+4.*dyuz*dzuy);
01692 break;
01693 case 40:
01694 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dxux;
01695 break;
01696 case 41:
01697 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuy+dyux);
01698 break;
01699 case 42:
01700 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dyuy;
01701 break;
01702 case 43:
01703 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuz+dzux);
01704 break;
01705 case 44:
01706 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dyuz+dzuy);
01707 break;
01708 case 45:
01709 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dzuz;
01710 break;
01711 }
01712 }
01713 }
01714 }
01715
01716 virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
01717 const int skip_ghosts=1) const {
01718 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
01719 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
01720 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
01721 MicroType *f = (MicroType *)fvec.databuffer();
01722 DataType *work = (DataType *)workvec.databuffer();
01723 DCoords dx = base::GH().worldStep(fvec.stepsize());
01724
01725 for (register int k=mzs; k<=mze; k++)
01726 for (register int j=mys; j<=mye; j++)
01727 for (register int i=mxs; i<=mxe; i++) {
01728 switch(cnt) {
01729 case 1:
01730 f[base::idx(i,j,k,Nx,Ny)](0)=work[base::idx(i,j,k,Nx,Ny)]/R0;
01731 break;
01732 case 2:
01733 f[base::idx(i,j,k,Nx,Ny)](1)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
01734 break;
01735 case 3:
01736 f[base::idx(i,j,k,Nx,Ny)](2)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
01737 break;
01738 case 4:
01739 f[base::idx(i,j,k,Nx,Ny)](3)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
01740 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium((const MacroType &)f[base::idx(i,j,k,Nx,Ny)]);
01741 break;
01742 }
01743 }
01744 }
01745
01746 virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
01747 const int verbose) const {
01748 int Nx = fvec.extents(0), Ny = fvec.extents(1);
01749 int b = fvec.bottom();
01750 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
01751 BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
01752 int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
01753 int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
01754 int mzs = std::max(bb.lower(2)/bb.stepsize(2),bbmax.lower(2)/bbmax.stepsize(2));
01755 int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
01756 int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
01757 int mze = std::min(bb.upper(2)/bb.stepsize(2),bbmax.upper(2)/bbmax.stepsize(2));
01758 MicroType *f = (MicroType *)fvec.databuffer();
01759
01760
01761 int result = 1;
01762 for (register int k=mzs; k<=mze; k++)
01763 for (register int j=mys; j<=mye; j++)
01764 for (register int i=mxs; i<=mxe; i++)
01765 for (register int l=0; l<base::NMicroVar(); l++)
01766 if (!(f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l)>-1.e37)) {
01767 result = 0;
01768 if (verbose) {
01769 DCoords lbcorner = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
01770 DCoords dx = base::GH().worldStep(fvec.stepsize());
01771
01772 std::cerr << "D3Q19-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << "," << k-mdz[l] << ")(" << l << ")="
01773 << f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l) << " " << fvec.bbox()
01774 << " Coords: (" << (i+0.5)*dx(0)+lbcorner(0) << "," << (j+0.5)*dx(1)+lbcorner(1) << ","
01775 << (k+0.5)*dx(2)+lbcorner(2) << ")" << std::endl;
01776 }
01777 }
01778 return result;
01779 }
01780
01781 inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
01782 TensorType P;
01783 if (method[1]==0) {
01784 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);
01785 P(1) = q(0)*q(1)*q(2)-(f(7)+f(8)-f(9)-f(10));
01786 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);
01787 P(3) = q(0)*q(1)*q(3)-(f(15)+f(16)-f(17)-f(18));
01788 P(4) = q(0)*q(2)*q(3)-(f(11)+f(12)-f(13)-f(14));
01789 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);
01790 P *= -(DataType(1.0)-DataType(0.5)*om);
01791 }
01792 else {
01793 MicroType feq = Equilibrium(q);
01794 P = DeviatoricStress(f,feq,om);
01795 }
01796 P(0) -= LatticeBasePressure(q(0));
01797 P(2) -= LatticeBasePressure(q(0));
01798 P(5) -= LatticeBasePressure(q(0));
01799 return P;
01800 }
01801
01802 inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
01803 TensorType Sigma;
01804 if (method[2]==0) {
01805 MicroType fneq = feq - f;
01806 Sigma(0) = fneq(1)+fneq(2)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
01807 Sigma(1) = fneq(7)+fneq(8)-fneq(9)-fneq(10);
01808 Sigma(2) = fneq(3)+fneq(4)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(11)+fneq(12)+fneq(13)+fneq(14);
01809 Sigma(3) = fneq(15)+fneq(16)-fneq(17)-fneq(18);
01810 Sigma(4) = fneq(11)+fneq(12)-fneq(13)-fneq(14);
01811 Sigma(5) = fneq(5)+fneq(6)+fneq(11)+fneq(12)+fneq(13)+fneq(14)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
01812 Sigma *= -(DataType(1.0)-DataType(0.5)*om);
01813 }
01814 else {
01815 MacroType q = MacroVariables(f);
01816 Sigma = Stress(f,q,om);
01817 Sigma(0) += LatticeBasePressure(q(0));
01818 Sigma(2) += LatticeBasePressure(q(0));
01819 Sigma(5) += LatticeBasePressure(q(0));
01820 }
01821 return Sigma;
01822 }
01823
01824 virtual int NMethodOrder() const { return 2; }
01825
01826 inline const DataType& L0() const { return base::LengthScale(); }
01827 inline const DataType& T0() const { return base::TimeScale(); }
01828 inline void SetDensityScale(const DataType r0) { R0 = r0; }
01829 inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
01830 inline void SetSpeedUp(const DataType s0) { S0 = s0; }
01831 inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
01832
01833 inline const DataType& DensityScale() const { return R0; }
01834 inline const DataType VelocityScale() const { return U0/S0; }
01835
01836 inline const DataType& SpeedUp() const { return S0; }
01837
01838
01839 inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
01840 inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
01841 inline virtual DataType LatticeBasePressure(const DataType rho) const {
01842 return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
01843 }
01844
01845
01846 inline void SetGas(DataType rho, DataType nu, DataType cs) {
01847 rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
01848 SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
01849 }
01850 inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
01851 inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q,
01852 const DataType dt) const {
01853 DataType M2 = 0.0;
01854
01855 M2 += (f(7) - feq(7))*(f(7) - feq(7));
01856 M2 += (f(8) - feq(8))*(f(8) - feq(8));
01857 M2 += (f(9) - feq(9))*(f(9) - feq(9));
01858 M2 += (f(10) - feq(10))*(f(10) - feq(10));
01859
01860 M2 += (f(11) - feq(11))*(f(11) - feq(11));
01861 M2 += (f(12) - feq(12))*(f(12) - feq(12));
01862 M2 += (f(13) - feq(13))*(f(13) - feq(13));
01863 M2 += (f(14) - feq(14))*(f(14) - feq(14));
01864
01865 M2 += (f(15) - feq(15))*(f(15) - feq(15));
01866 M2 += (f(16) - feq(16))*(f(16) - feq(16));
01867 M2 += (f(17) - feq(17))*(f(17) - feq(17));
01868 M2 += (f(18) - feq(18))*(f(18) - feq(18));
01869 M2 = std::sqrt(M2);
01870
01871 DataType tau = DataType(1.0)/Omega(dt);
01872 DataType om_turb = (DataType(2.0)/( tau
01873 + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
01874
01875 return ( om_turb );
01876 }
01877 inline const DataType Omega_LES_dynamic(const MicroType &f, const MicroType &feq, const MacroType& q,
01878 const DataType dt) const {
01879 DataType M2 = 0.0;
01880
01881 M2 += (f(7) - feq(7))*(f(7) - feq(7));
01882 M2 += (f(8) - feq(8))*(f(8) - feq(8));
01883 M2 += (f(9) - feq(9))*(f(9) - feq(9));
01884 M2 += (f(10) - feq(10))*(f(10) - feq(10));
01885
01886 M2 += (f(11) - feq(11))*(f(11) - feq(11));
01887 M2 += (f(12) - feq(12))*(f(12) - feq(12));
01888 M2 += (f(13) - feq(13))*(f(13) - feq(13));
01889 M2 += (f(14) - feq(14))*(f(14) - feq(14));
01890
01891 M2 += (f(15) - feq(15))*(f(15) - feq(15));
01892 M2 += (f(16) - feq(16))*(f(16) - feq(16));
01893 M2 += (f(17) - feq(17))*(f(17) - feq(17));
01894 M2 += (f(18) - feq(18))*(f(18) - feq(18));
01895 M2 = std::sqrt(M2);
01896
01897 DataType tau = DataType(1.0)/Omega(dt);
01898 DataType vel = sqrt( q(1)*q(1)+q(2)*q(2)+q(3)*q(3) );
01899 DataType Cs_SmagDYN = DataType(0.5)*L0()*L0()*( (vel*vel/csp*(M2-vel)) / ((M2-vel)*(M2-vel)) );
01900 DataType om_turb = (DataType(2.0)/( tau
01901 + sqrt(tau*tau + (DataType(18.)*Cs_SmagDYN*M2)/(R0*q(0)*cs2/S0/S0)) ));
01902
01903 return ( om_turb );
01904 }
01905 inline const int TurbulenceType() const { return turbulence; }
01906 inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
01907 inline const DataType& GasDensity() const { return rhop; }
01908 inline const DataType& GasViscosity() const { return nup; }
01909 inline const DataType& GasSpeedofSound() const { return csp; }
01910 inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
01911 { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
01912
01913 inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
01914 inline virtual const DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); }
01915 inline virtual const DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); }
01916
01917 protected:
01918 DataType cs2, cs22, cssq;
01919 DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
01920 DataType Cs_Smagorinsky, turbulence;
01921 int method[3], mdx[19], mdy[19], mdz[19];
01922 };
01923
01924
01925 #endif