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 CharacteristicOutlet, CharacteristicInlet, NoSlipWallEntranceExit };
00066 enum GFMPredefined { GFMExtrapolation, GFMSlipWall, GFMNoSlipWall, GFMWallLaw };
00067 enum TurbulenceModel { laminar, LES_Smagorinsky, LES_dynamic };
00068
00069 LBMD3Q19() : base(), R0(1.), U0(1.), S0(1.), rhop(0.), csp(0.), cs2p(0.), nup (0.), gp(0.), Wp(1.), Rp(1.) {
00070 cs2 = DataType(1.)/DataType(3.);
00071 cs22 = DataType(2.)*cs2;
00072 cssq = DataType(2.)/DataType(9.);
00073 method[0] = 2; method[1] = 0; method[2] = 0; method[3] = 0;
00074 turbulence = laminar;
00075 mdx[0]=mdx[3]=mdx[4]=mdx[5]=mdx[6]=mdx[11]=mdx[12]=mdx[13]=mdx[14]=0;
00076 mdx[1]=mdx[7]=mdx[9]=mdx[15]=mdx[17]=1;
00077 mdx[2]=mdx[8]=mdx[10]=mdx[16]=mdx[18]=-1;
00078
00079 mdy[0]=mdy[1]=mdy[2]=mdy[5]=mdy[6]=mdy[15]=mdy[16]=mdy[17]=mdy[18]=0;
00080 mdy[3]=mdy[7]=mdy[10]=mdy[11]=mdy[13]=1;
00081 mdy[4]=mdy[8]=mdy[9]=mdy[12]=mdy[14]=-1;
00082
00083 mdz[0]=mdz[1]=mdz[2]=mdz[3]=mdz[4]=mdz[7]=mdz[8]=mdz[9]=mdz[10]=0;
00084 mdz[5]=mdz[11]=mdz[14]=mdz[15]=mdz[18]=1;
00085 mdz[6]=mdz[12]=mdz[13]=mdz[16]=mdz[17]=-1;
00086
00087 }
00088
00089 virtual ~LBMD3Q19() {}
00090
00091 virtual void register_at(ControlDevice& Ctrl, const std::string& prefix) {
00092 base::LocCtrl = Ctrl.getSubDevice(prefix+"D3Q19");
00093 RegisterAt(base::LocCtrl,"Method(1)",method[0]);
00094 RegisterAt(base::LocCtrl,"Method(2)",method[1]);
00095 RegisterAt(base::LocCtrl,"Method(3)",method[2]);
00096 RegisterAt(base::LocCtrl,"Method(4)",method[3]);
00097 RegisterAt(base::LocCtrl,"Turbulence",turbulence);
00098 RegisterAt(base::LocCtrl,"Cs_Smagorinsky",Cs_Smagorinsky);
00099 RegisterAt(base::LocCtrl,"Gas_rho",rhop);
00100 RegisterAt(base::LocCtrl,"Gas_Cs",csp);
00101 RegisterAt(base::LocCtrl,"Gas_nu",nup);
00102 RegisterAt(base::LocCtrl,"Gas_gamma",gp);
00103 RegisterAt(base::LocCtrl,"Gas_W",Wp);
00104 RegisterAt(base::LocCtrl,"Ideal_R",Rp);
00105 RegisterAt(base::LocCtrl,"SpeedUp",S0);
00106 }
00107
00108 virtual void SetupData(GridHierarchy* gh, const int& ghosts) {
00109 base::SetupData(gh,ghosts);
00110 if (rhop>0.) R0 = rhop;
00111 if (csp>0.) {
00112 cs2p = csp*csp;
00113 SetTimeScale(LatticeSpeedOfSound()/csp*L0());
00114 }
00115 std::cout << "D3Q19: Gas_rho=" << rhop << " Gas_Cs=" << csp
00116 << " Gas_nu=" << nup << std::endl;
00117 std::cout << "D3Q19: L0=" << L0() << " T0=" << T0() << " U0=" << U0
00118 << " R0=" << R0 << std::endl;
00119 if ( TurbulenceType() == laminar ) {
00120 std::cout << "LBM Setup() Laminar" << std::endl;
00121 }
00122 if ( TurbulenceType() == LES_Smagorinsky ) {
00123 std::cout << "LBM Setup() LES_Smagorinsky : Smagorinsky Constant = " << SmagorinskyConstant()
00124 << std::endl;
00125
00126 }
00127 if ( TurbulenceType() == LES_dynamic ) {
00128 std::cout << "LBM Setup() LES_dynamic " << std::endl;
00129 }
00130 WriteInit();
00131 }
00132
00133 virtual void WriteInit() const {
00134 int me = MY_PROC;
00135 if (me == VizServer) {
00136 std::ofstream ofs("chem.dat", std::ios::out);
00137 ofs << "RU " << Rp << std::endl;
00138 ofs << "PA " << BasePressure() << std::endl;
00139 ofs << "Sp(01) Vicosity" << std::endl;
00140 ofs << "W(01) " << Wp << std::endl;
00141 ofs << "Sp(02) |Velocity|" << std::endl;
00142 ofs << "W(02) " << 1. << std::endl;
00143 ofs << "Sp(03) Vorticity-x" << std::endl;
00144 ofs << "W(03) " << 1. << std::endl;
00145 ofs << "Sp(04) Vorticity-y" << std::endl;
00146 ofs << "W(04) " << 1. << std::endl;
00147 ofs << "Sp(05) Vorticity-z" << std::endl;
00148 ofs << "W(05) " << 1. << std::endl;
00149 ofs << "Sp(06) |Vorticity|" << std::endl;
00150 ofs << "W(06) " << 1. << std::endl;
00151 ofs.close();
00152 }
00153 }
00154
00155 inline virtual MacroType MacroVariables(const MicroType &f) const {
00156 MacroType q;
00157 q(0)=f(0)+f(1)+f(2)+f(3)+f(4)+f(5)+f(6)+f(7)+f(8)+f(9)
00158 +f(10)+f(11)+f(12)+f(13)+f(14)+f(15)+f(16)+f(17)+f(18);
00159 q(1)=(f(1)-f(2)+f(7)-f(8)+f(9)-f(10)+f(15)-f(16)+f(17)-f(18))/q(0);
00160 q(2)=(f(3)-f(4)+f(7)-f(8)-f(9)+f(10)+f(11)-f(12)+f(13)-f(14))/q(0);
00161 q(3)=(f(5)-f(6)+f(11)-f(12)-f(13)+f(14)+f(15)-f(16)-f(17)+f(18))/q(0);
00162 return q;
00163 }
00164
00165 inline virtual MicroType Equilibrium(const MacroType &q) const {
00166 MicroType feq;
00167
00168 DataType rho = q(0);
00169 DataType u = q(1);
00170 DataType v = q(2);
00171 DataType w = q(3);
00172
00173 DataType ri, rt0, rt1, rt2;
00174 if (method[0]==0) {
00175 ri = DataType(1.);
00176 rt0 = R0 / DataType(3.);
00177 rt1 = R0 / DataType(18.);
00178 rt2 = R0 / DataType(36.);
00179 }
00180 if (method[0]==1) {
00181 ri = rho;
00182 rt0 = DataType(1.) / DataType(3.);
00183 rt1 = DataType(1.) / DataType(18.);
00184 rt2 = DataType(1.) / DataType(36.);
00185 }
00186 else if (method[0]==2) {
00187 ri = DataType(1.);
00188 rt0 = rho / DataType(3.);
00189 rt1 = rho / DataType(18.);
00190 rt2 = rho / DataType(36.);
00191 }
00192
00193 DataType usq = u * u;
00194 DataType vsq = v * v;
00195 DataType wsq = w * w;
00196 DataType sumsq = (usq + vsq + wsq) / cs22;
00197 DataType sumsqxy = (usq + vsq) / cs22;
00198 DataType sumsqyz = (vsq + wsq) / cs22;
00199 DataType sumsqxz = (usq + wsq) / cs22;
00200 DataType sumsq2xy = sumsqxy * (DataType(1.) - cs2) / cs2;
00201 DataType sumsq2yz = sumsqyz * (DataType(1.) - cs2) / cs2;
00202 DataType sumsq2xz = sumsqxz * (DataType(1.) - cs2) / cs2;
00203 DataType u2 = usq / cssq;
00204 DataType v2 = vsq / cssq;
00205 DataType w2 = wsq / cssq;
00206 DataType u22 = usq / cs22;
00207 DataType v22 = vsq / cs22;
00208 DataType w22 = wsq / cs22;
00209 DataType ui = u / cs2;
00210 DataType vi = v / cs2;
00211 DataType wi = w / cs2;
00212 DataType uv = ui * vi;
00213 DataType uw = ui * wi;
00214 DataType vw = vi * wi;
00215
00216 feq(0) = rt0 * (ri - sumsq);
00217
00218 feq(1) = rt1 * (ri - sumsq + u2 + ui);
00219 feq(2) = rt1 * (ri - sumsq + u2 - ui);
00220
00221 feq(3) = rt1 * (ri - sumsq + v2 + vi);
00222 feq(4) = rt1 * (ri - sumsq + v2 - vi);
00223
00224 feq(5) = rt1 * (ri - sumsq + w2 + wi);
00225 feq(6) = rt1 * (ri - sumsq + w2 - wi);
00226
00227 feq(7) = rt2 * (ri + sumsq2xy -w22 + ui + vi + uv);
00228 feq(8) = rt2 * (ri + sumsq2xy -w22 - ui - vi + uv);
00229 feq(9) = rt2 * (ri + sumsq2xy -w22 + ui - vi - uv);
00230 feq(10) = rt2 * (ri + sumsq2xy -w22 - ui + vi - uv);
00231
00232 feq(11) = rt2 * (ri + sumsq2yz -u22 + vi + wi + vw);
00233 feq(12) = rt2 * (ri + sumsq2yz -u22 - vi - wi + vw);
00234 feq(13) = rt2 * (ri + sumsq2yz -u22 + vi - wi - vw);
00235 feq(14) = rt2 * (ri + sumsq2yz -u22 - vi + wi - vw);
00236
00237 feq(15) = rt2 * (ri + sumsq2xz -v22 + ui + wi + uw);
00238 feq(16) = rt2 * (ri + sumsq2xz -v22 - ui - wi + uw);
00239 feq(17) = rt2 * (ri + sumsq2xz -v22 + ui - wi - uw);
00240 feq(18) = rt2 * (ri + sumsq2xz -v22 - ui + wi - uw);
00241
00242
00243 if (method[3]>0) {
00244 DataType cs4 = DataType(1.)/(DataType(4.)*cs2*cs2);
00245 DataType sumsq3 = (usq + vsq + wsq) * cs4;
00246 DataType sumsqxy3 = (usq + vsq) * cs4;
00247 DataType sumsqyz3 = (vsq + wsq) * cs4;
00248 DataType sumsqxz3 = (usq + wsq) * cs4;
00249 DataType upv = (u+v)*(u+v) * cs4;
00250 DataType umv = (u-v)*(u-v) * cs4;
00251 DataType vpw = (v+w)*(v+w) * cs4;
00252 DataType vmw = (v-w)*(v-w) * cs4;
00253 DataType upw = (u+w)*(u+w) * cs4;
00254 DataType umw = (u-w)*(u-w) * cs4;
00255
00256
00257 feq(1) -= rt1 * u*sumsqyz3;
00258 feq(2) += rt1 * u*sumsqyz3;
00259
00260 feq(3) -= rt1 * v*sumsqxz3;
00261 feq(4) += rt1 * v*sumsqxz3;
00262
00263 feq(5) -= rt1 * w*sumsqxy3;
00264 feq(6) += rt1 * w*sumsqxy3;
00265
00266 feq(7) -= rt2 * (u+v)*(sumsq3-upv);
00267 feq(8) += rt2 * (u+v)*(sumsq3-upv);
00268 feq(9) -= rt2 * (u-v)*(sumsq3-umv);
00269 feq(10) += rt2 * (u-v)*(sumsq3-umv);
00270
00271 feq(11) -= rt2 * (v+w)*(sumsq3-vpw);
00272 feq(12) += rt2 * (v+w)*(sumsq3-vpw);
00273 feq(13) -= rt2 * (v-w)*(sumsq3-vmw);
00274 feq(14) += rt2 * (v-w)*(sumsq3-vmw);
00275
00276 feq(15) -= rt2 * (u+w)*(sumsq3-upw);
00277 feq(16) += rt2 * (u+w)*(sumsq3-upw);
00278 feq(17) -= rt2 * (u-w)*(sumsq3-umw);
00279 feq(18) += rt2 * (u-w)*(sumsq3-umw);
00280 }
00281
00282 return feq;
00283 }
00284
00285 inline virtual void Collision(MicroType &f, const DataType dt) const {
00286 MacroType q = MacroVariables(f);
00287 MicroType feq = Equilibrium(q);
00288
00289 DataType omega;
00290 if (turbulence == laminar)
00291 omega = Omega(dt);
00292 else if (turbulence == LES_Smagorinsky)
00293 omega = Omega_LES_Smagorinsky(f,feq,q,dt);
00294 else if (turbulence == LES_dynamic)
00295 omega = Omega_LES_dynamic(f,feq,q,dt);
00296
00297 f=(DataType(1.)-omega)*f + omega*feq;
00298 }
00299
00300 virtual int IncomingIndices(const int side, int indices[]) const {
00301 switch (side) {
00302 case base::Left:
00303 indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;
00304 break;
00305 case base::Right:
00306 indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;
00307 break;
00308 case base::Bottom:
00309 indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;
00310 break;
00311 case base::Top:
00312 indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;
00313 break;
00314 case base::Front:
00315 indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;
00316 break;
00317 case base::Back:
00318 indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;
00319 break;
00320 }
00321 return 5;
00322 }
00323
00324 virtual int OutgoingIndices(const int side, int indices[]) const {
00325 switch (side) {
00326 case base::Left:
00327 indices[0] = 2; indices[1] = 8; indices[2] = 10; indices[3] = 16; indices[4] = 18;
00328 break;
00329 case base::Right:
00330 indices[0] = 1; indices[1] = 7; indices[2] = 9; indices[3] = 15; indices[4] = 17;
00331 break;
00332 case base::Bottom:
00333 indices[0] = 4; indices[1] = 8; indices[2] = 9; indices[3] = 12; indices[4] = 14;
00334 break;
00335 case base::Top:
00336 indices[0] = 3; indices[1] = 7; indices[2] = 10; indices[3] = 11; indices[4] = 13;
00337 break;
00338 case base::Front:
00339 indices[0] = 6; indices[1] = 12; indices[2] = 13; indices[3] = 16; indices[4] = 17;
00340 break;
00341 case base::Back:
00342 indices[0] = 5; indices[1] = 11; indices[2] = 14; indices[3] = 15; indices[4] = 18;
00343 break;
00344 }
00345 return 5;
00346 }
00347
00348
00349 virtual void ReverseStream(vec_grid_data_type &fvec, const BBox &bb,
00350 const int side) const {
00351 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00352 int b = fvec.bottom();
00353 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00354 && fvec.stepsize(2)==bb.stepsize(2));
00355 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00356 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00357 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00358 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00359 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00360 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00361 MicroType *f = (MicroType *)fvec.databuffer();
00362
00363 #ifdef DEBUG_PRINT_FIXUP
00364 ( comm_service::log() << "Reverse streaming at Side: " << side << " of "
00365 << fvec.bbox() << " using " << bb << "\n").flush();
00366 #endif
00367
00368 switch (side) {
00369 case base::Left:
00370 for (register int k=mzs; k<=mze; k++)
00371 for (register int j=mys; j<=mye; j++) {
00372 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](2);
00373 f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](8);
00374 f[b+base::idx(mxs,j,k,Nx,Ny)](10)= f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](10);
00375 f[b+base::idx(mxs,j,k,Nx,Ny)](16)= f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](16);
00376 f[b+base::idx(mxs,j,k,Nx,Ny)](18)= f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](18);
00377 }
00378 break;
00379 case base::Right:
00380 for (register int k=mzs; k<=mze; k++)
00381 for (register int j=mys; j<=mye; j++) {
00382 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](1);
00383 f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](7);
00384 f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](9);
00385 f[b+base::idx(mxe,j,k,Nx,Ny)](15)= f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](15);
00386 f[b+base::idx(mxe,j,k,Nx,Ny)](17)= f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](17);
00387 }
00388 break;
00389 case base::Bottom:
00390 for (register int k=mzs; k<=mze; k++)
00391 for (register int i=mxs; i<=mxe; i++) {
00392 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](4);
00393 f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](8);
00394 f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i+1,mys-1,k,Nx,Ny)](9);
00395 f[b+base::idx(i,mys,k,Nx,Ny)](12)= f[b+base::idx(i,mys-1,k-1,Nx,Ny)](12);
00396 f[b+base::idx(i,mys,k,Nx,Ny)](14)= f[b+base::idx(i,mys-1,k+1,Nx,Ny)](14);
00397 }
00398 break;
00399 case base::Top:
00400 for (register int k=mzs; k<=mze; k++)
00401 for (register int i=mxs; i<=mxe; i++) {
00402 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](3);
00403 f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i+1,mye+1,k,Nx,Ny)](7);
00404 f[b+base::idx(i,mye,k,Nx,Ny)](10)= f[b+base::idx(i-1,mye+1,k,Nx,Ny)](10);
00405 f[b+base::idx(i,mye,k,Nx,Ny)](11)= f[b+base::idx(i,mye+1,k+1,Nx,Ny)](11);
00406 f[b+base::idx(i,mye,k,Nx,Ny)](13)= f[b+base::idx(i,mye+1,k-1,Nx,Ny)](13);
00407 }
00408 break;
00409 case base::Front:
00410 for (register int j=mys; j<=mye; j++)
00411 for (register int i=mxs; i<=mxe; i++) {
00412 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](6);
00413 f[b+base::idx(i,j,mzs,Nx,Ny)](12)= f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](12);
00414 f[b+base::idx(i,j,mzs,Nx,Ny)](13)= f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](13);
00415 f[b+base::idx(i,j,mzs,Nx,Ny)](16)= f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](16);
00416 f[b+base::idx(i,j,mzs,Nx,Ny)](17)= f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](17);
00417 }
00418 break;
00419 case base::Back:
00420 for (register int j=mys; j<=mye; j++)
00421 for (register int i=mxs; i<=mxe; i++) {
00422 f[b+base::idx(i,j,mzs,Nx,Ny)](5) = f[b+base::idx(i,j,mzs+1,Nx,Ny)](5);
00423 f[b+base::idx(i,j,mzs,Nx,Ny)](11)= f[b+base::idx(i,j+1,mzs+1,Nx,Ny)](11);
00424 f[b+base::idx(i,j,mzs,Nx,Ny)](14)= f[b+base::idx(i,j-1,mzs+1,Nx,Ny)](14);
00425 f[b+base::idx(i,j,mzs,Nx,Ny)](15)= f[b+base::idx(i+1,j,mzs+1,Nx,Ny)](15);
00426 f[b+base::idx(i,j,mzs,Nx,Ny)](18)= f[b+base::idx(i-1,j,mzs+1,Nx,Ny)](18);
00427 }
00428 break;
00429 }
00430 }
00431
00432 virtual void LocalStep(vec_grid_data_type &fvec, vec_grid_data_type &ovec,
00433 const BBox &bb, const double& dt) const {
00434 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00435 int b = fvec.bottom();
00436 MicroType *f = (MicroType *)fvec.databuffer();
00437 MicroType *of = (MicroType *)ovec.databuffer();
00438
00439 #ifdef DEBUG_PRINT_FIXUP
00440 ( comm_service::log() << "Local Step in " << fvec.bbox() << " from " << ovec.bbox()
00441 << " using " << bb << "\n").flush();
00442 #endif
00443
00444 assert (fvec.extents(0)==ovec.extents(0));
00445 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1) && fvec.stepsize(2)==bb.stepsize(2) &&
00446 fvec.stepsize(0)==ovec.stepsize(0) && fvec.stepsize(1)==ovec.stepsize(1) && fvec.stepsize(2)==ovec.stepsize(2));
00447 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00448 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00449 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00450 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00451 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00452 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00453
00454 for (register int k=mzs; k<=mze; k++)
00455 for (register int j=mys; j<=mye; j++)
00456 for (register int i=mxs; i<=mxe; i++) {
00457 for (register int n=0; n<base::NMicroVar(); n++)
00458 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);
00459 Collision(f[b+base::idx(i,j,k,Nx,Ny)],dt);
00460 }
00461 }
00462
00463 virtual double Step(vec_grid_data_type &fvec, vec_grid_data_type &ovec, vec_grid_data_type* Flux[],
00464 const double& t, const double& dt, const int& mpass) const {
00465
00466
00467 DCoords dx = base::GH().worldStep(fvec.stepsize());
00468
00469
00470 if ( (std::fabs(dx(0)/L0()-dx(1)/L0()) > DBL_EPSILON*LB_FACTOR)
00471 || (std::fabs(dx(1)/L0()-dx(2)/L0()) > DBL_EPSILON*LB_FACTOR) ) {
00472 std::cerr << "LBMD3Q19 expects dx(0)=" << dx(0)/L0() << " and dx(1)=" << dx(1)/L0()
00473 << " and dx(2)=" << dx(2)/L0()
00474 << " to be equal." << std::endl << std::flush;
00475 std::exit(-1);
00476 }
00477 if (std::fabs(dt/T0()-dx(0)/L0()) > DBL_EPSILON*LB_FACTOR) {
00478 std::cerr << "LBMD3Q19 expects dt=" << dt/T0() << " and dx(0)=" << dx(0)/L0()
00479 << " to be equal."
00480 << " dt=" << dt << " TO=" << T0()
00481 << std::endl << std::flush;
00482 std::exit(-1);
00483 }
00484
00485 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00486 MicroType *f = (MicroType *)fvec.databuffer();
00487
00488
00489 for (register int k=Nz-1; k>=1; k--)
00490 for (register int j=Ny-1; j>=1; j--)
00491 for (register int i=0; i<=Nx-2; i++) {
00492
00493 f[base::idx(i,j,k,Nx,Ny)](3) = f[base::idx(i,j-1,k,Nx,Ny)](3);
00494 f[base::idx(i,j,k,Nx,Ny)](10) = f[base::idx(i+1,j-1,k,Nx,Ny)](10);
00495
00496
00497
00498 f[base::idx(i,j,k,Nx,Ny)](5) = f[base::idx(i,j,k-1,Nx,Ny)](5);
00499 f[base::idx(i,j,k,Nx,Ny)](18) = f[base::idx(i+1,j,k-1,Nx,Ny)](18);
00500 }
00501
00502 for (register int k=0; k<=Nz-2; k++)
00503 for (register int j=Ny-1; j>=1; j--)
00504 for (register int i=0; i<=Nx-2; i++) {
00505
00506 f[base::idx(i,j,k,Nx,Ny)](13) = f[base::idx(i,j-1,k+1,Nx,Ny)](13);
00507 }
00508 for (register int k=Nz-1; k>=1; k--)
00509 for (register int j=Ny-1; j>=1; j--)
00510 for (register int i=Nx-1; i>=1; i--) {
00511
00512 f[base::idx(i,j,k,Nx,Ny)](1) = f[base::idx(i-1,j,k,Nx,Ny)](1);
00513 f[base::idx(i,j,k,Nx,Ny)](7) = f[base::idx(i-1,j-1,k,Nx,Ny)](7);
00514
00515 f[base::idx(i,j,k,Nx,Ny)](11) = f[base::idx(i,j-1,k-1,Nx,Ny)](11);
00516
00517 f[base::idx(i,j,k,Nx,Ny)](15) = f[base::idx(i-1,j,k-1,Nx,Ny)](15);
00518 }
00519
00520 for (register int k=0; k<=Nz-2; k++)
00521 for (register int j=0; j<=Ny-2; j++)
00522 for (register int i=Nx-1; i>=1; i--) {
00523
00524 f[base::idx(i,j,k,Nx,Ny)](4) = f[base::idx(i,j+1,k,Nx,Ny)](4);
00525 f[base::idx(i,j,k,Nx,Ny)](9) = f[base::idx(i-1,j+1,k,Nx,Ny)](9);
00526
00527
00528
00529 f[base::idx(i,j,k,Nx,Ny)](6) = f[base::idx(i,j,k+1,Nx,Ny)](6);
00530 f[base::idx(i,j,k,Nx,Ny)](17) = f[base::idx(i-1,j,k+1,Nx,Ny)](17);
00531 }
00532 for (register int k=Nz-1; k>=1; k--)
00533 for (register int j=0; j<=Ny-2; j++)
00534 for (register int i=0; i<=Nx-2; i++) {
00535
00536 f[base::idx(i,j,k,Nx,Ny)](14) = f[base::idx(i,j+1,k-1,Nx,Ny)](14);
00537 }
00538 for (register int k=0; k<=Nz-2; k++)
00539 for (register int j=0; j<=Ny-2; j++)
00540 for (register int i=0; i<=Nx-2; i++) {
00541
00542 f[base::idx(i,j,k,Nx,Ny)](2) = f[base::idx(i+1,j,k,Nx,Ny)](2);
00543 f[base::idx(i,j,k,Nx,Ny)](8) = f[base::idx(i+1,j+1,k,Nx,Ny)](8);
00544
00545 f[base::idx(i,j,k,Nx,Ny)](12) = f[base::idx(i,j+1,k+1,Nx,Ny)](12);
00546
00547 f[base::idx(i,j,k,Nx,Ny)](16) = f[base::idx(i+1,j,k+1,Nx,Ny)](16);
00548 }
00549
00550
00551 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00552 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00553 for (register int k=mzs; k<=mze; k++)
00554 for (register int j=mys; j<=mye; j++)
00555 for (register int i=mxs; i<=mxe; i++)
00556 Collision(f[base::idx(i,j,k,Nx,Ny)],dt);
00557
00558 return DataType(1.);
00559 }
00560
00561 virtual void ICStandard(vec_grid_data_type &fvec, const int type,
00562 DataType* aux=0, const int naux=0, const int scaling=0) const {
00563 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
00564 int mxs = base::NGhosts(), mys = base::NGhosts(), mzs = base::NGhosts();
00565 int mxe = Nx-base::NGhosts()-1, mye = Ny-base::NGhosts()-1, mze = Nz-base::NGhosts()-1;
00566 MicroType *f = (MicroType *)fvec.databuffer();
00567
00568
00569
00570 if (type==ConstantMicro) {
00571
00572 MicroType g;
00573 for (register int i=0; i<base::NMicroVar(); i++)
00574 g(i) = aux[i];
00575
00576 for (register int k=mzs; k<=mze; k++)
00577 for (register int j=mys; j<=mye; j++)
00578 for (register int i=mxs; i<=mxe; i++)
00579 f[base::idx(i,j,k,Nx,Ny)] = g;
00580
00581 }
00582
00583 else if (type==ConstantMacro) {
00584
00585 MacroType q;
00586 if (scaling==base::Physical) {
00587 q(0) = aux[0]/R0;
00588 q(1) = aux[1]/U0;
00589 q(2) = aux[2]/U0;
00590 q(3) = aux[3]/U0;
00591 }
00592 else
00593 for (register int i=0; i<base::NMacroVar(); i++)
00594 q(i) = aux[i];
00595 q(1) *= S0; q(2) *= S0; q(3) *= S0;
00596
00597
00598
00599
00600
00601
00602
00603 for (register int k=mzs; k<=mze; k++)
00604 for (register int j=mys; j<=mye; j++)
00605 for (register int i=mxs; i<=mxe; i++)
00606 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00607
00608 }
00609
00610 else if (type==GasAtRest) {
00611
00612 MacroType q;
00613 q(0) = GasDensity()/R0;
00614 q(1) = DataType(0.);
00615 q(2) = DataType(0.);
00616 q(3) = DataType(0.);
00617 for (register int k=mzs; k<=mze; k++)
00618 for (register int j=mys; j<=mye; j++)
00619 for (register int i=mxs; i<=mxe; i++)
00620 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium(q);
00621
00622 }
00623 }
00624
00625
00626 virtual void BCStandard(vec_grid_data_type &fvec, const BBox &bb, const int type, const int side,
00627 DataType* aux=0, const int naux=0, const int scaling=0) const {
00628 int Nx = fvec.extents(0), Ny = fvec.extents(1);
00629 int b = fvec.bottom();
00630 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1)
00631 && fvec.stepsize(2)==bb.stepsize(2));
00632 int mxs = std::max(bb.lower(0)/bb.stepsize(0),fvec.lower(0)/fvec.stepsize(0));
00633 int mys = std::max(bb.lower(1)/bb.stepsize(1),fvec.lower(1)/fvec.stepsize(1));
00634 int mzs = std::max(bb.lower(2)/bb.stepsize(2),fvec.lower(2)/fvec.stepsize(2));
00635 int mxe = std::min(bb.upper(0)/bb.stepsize(0),fvec.upper(0)/fvec.stepsize(0));
00636 int mye = std::min(bb.upper(1)/bb.stepsize(1),fvec.upper(1)/fvec.stepsize(1));
00637 int mze = std::min(bb.upper(2)/bb.stepsize(2),fvec.upper(2)/fvec.stepsize(2));
00638 MicroType *f = (MicroType *)fvec.databuffer();
00639
00640
00641 if (type==Symmetry) {
00642 switch (side) {
00643 case base::Left:
00644 for (register int k=mzs; k<=mze; k++)
00645 for (register int j=mys; j<=mye; j++) {
00646 if (j>mys) f[b+base::idx(mxe,j,k,Nx,Ny)](9) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
00647 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00648 if (j<mye) f[b+base::idx(mxe,j,k,Nx,Ny)](7) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
00649
00650 if (k>mzs) f[b+base::idx(mxe,j,k,Nx,Ny)](15) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
00651 if (k<mze) f[b+base::idx(mxe,j,k,Nx,Ny)](17) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
00652 }
00653 break;
00654 case base::Right:
00655 for (register int k=mzs; k<=mze; k++)
00656 for (register int j=mys; j<=mye; j++) {
00657 if (j>mys) f[b+base::idx(mxs,j,k,Nx,Ny)](10) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
00658 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00659 if (j<mye) f[b+base::idx(mxs,j,k,Nx,Ny)](8) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
00660
00661 if (k>mzs) f[b+base::idx(mxs,j,k,Nx,Ny)](18) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
00662 if (k<mze) f[b+base::idx(mxs,j,k,Nx,Ny)](16) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
00663 }
00664 break;
00665 case base::Bottom:
00666 for (register int k=mzs; k<=mze; k++)
00667 for (register int i=mxs; i<=mxe; i++) {
00668 if (i>mxs) f[b+base::idx(i,mye,k,Nx,Ny)](7) = f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
00669 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00670 if (i<mxe) f[b+base::idx(i,mye,k,Nx,Ny)](10) = f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
00671
00672 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);
00673 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);
00674 }
00675 break;
00676 case base::Top:
00677 for (register int k=mzs; k<=mze; k++)
00678 for (register int i=mxs; i<=mxe; i++) {
00679 if (i>mxs) f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
00680 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00681 if (i<mxe) f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
00682
00683 if (k>mzs) f[b+base::idx(i,mys,k,Nx,Ny)](14) = f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
00684 if (k<mze) f[b+base::idx(i,mys,k,Nx,Ny)](12) = f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
00685 }
00686 break;
00687 case base::Front:
00688
00689 for (register int j=mys; j<=mye; j++)
00690 for (register int i=mxs; i<=mxe; i++) {
00691 if (i>mxs) f[b+base::idx(i,j,mze,Nx,Ny)](18) = f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
00692 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00693 if (i<mxe) f[b+base::idx(i,j,mze,Nx,Ny)](15) = f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
00694
00695 if (j>mys) f[b+base::idx(i,j,mze,Nx,Ny)](14) = f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
00696 if (j<mye) f[b+base::idx(i,j,mze,Nx,Ny)](11) = f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
00697 }
00698 break;
00699 case base::Back:
00700
00701 for (register int j=mys; j<=mye; j++)
00702 for (register int i=mxs; i<=mxe; i++) {
00703
00704
00705
00706
00707 if (i>mxs) f[b+base::idx(i,j,mzs,Nx,Ny)](17) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
00708 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00709 if (i<mxe) f[b+base::idx(i,j,mzs,Nx,Ny)](16) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
00710
00711 if (j>mys) f[b+base::idx(i,j,mzs,Nx,Ny)](13) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
00712 if (j<mye) f[b+base::idx(i,j,mzs,Nx,Ny)](12) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
00713 }
00714 break;
00715 }
00716 }
00717
00718 else if (type==SlipWall) {
00719 switch (side) {
00720 case base::Left:
00721
00722 for (register int k=mzs; k<=mze; k++)
00723 for (register int j=mys; j<=mye; j++) {
00724 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);
00725 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00726 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);
00727
00728 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);
00729 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);
00730 }
00731 break;
00732 case base::Right:
00733
00734 for (register int k=mzs; k<=mze; k++)
00735 for (register int j=mys; j<=mye; j++) {
00736 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);
00737 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00738 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);
00739
00740 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);
00741 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);
00742 }
00743 break;
00744 case base::Bottom:
00745
00746 for (register int k=mzs; k<=mze; k++)
00747 for (register int i=mxs; i<=mxe; i++) {
00748 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);
00749 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00750 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);
00751
00752 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);
00753 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);
00754 }
00755 break;
00756 case base::Top:
00757
00758 for (register int k=mzs; k<=mze; k++)
00759 for (register int i=mxs; i<=mxe; i++) {
00760 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);
00761 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00762 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);
00763
00764 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);
00765 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);
00766 }
00767 break;
00768 case base::Front:
00769
00770 for (register int j=mys; j<=mye; j++)
00771 for (register int i=mxs; i<=mxe; i++) {
00772 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);
00773 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00774 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);
00775
00776 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);
00777 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);
00778 }
00779 break;
00780 case base::Back:
00781
00782 for (register int j=mys; j<=mye; j++)
00783 for (register int i=mxs; i<=mxe; i++) {
00784
00785
00786
00787
00788 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);
00789 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00790 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);
00791
00792 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);
00793 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);
00794 }
00795 break;
00796 }
00797 }
00798
00799 else if (type==NoSlipWall) {
00800 switch (side) {
00801 case base::Left:
00802
00803 for (register int k=mzs; k<=mze; k++)
00804 for (register int j=mys; j<=mye; j++) {
00805 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);
00806 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
00807 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);
00808
00809 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);
00810 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);
00811 }
00812 break;
00813 case base::Right:
00814
00815 for (register int k=mzs; k<=mze; k++)
00816 for (register int j=mys; j<=mye; j++) {
00817 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);
00818 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
00819 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);
00820
00821 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);
00822 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);
00823 }
00824 break;
00825 case base::Bottom:
00826
00827 for (register int k=mzs; k<=mze; k++)
00828 for (register int i=mxs; i<=mxe; i++) {
00829 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);
00830 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
00831 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);
00832
00833 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);
00834 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);
00835 }
00836 break;
00837 case base::Top:
00838
00839 for (register int k=mzs; k<=mze; k++)
00840 for (register int i=mxs; i<=mxe; i++) {
00841 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);
00842 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
00843 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);
00844
00845 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);
00846 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);
00847 }
00848 break;
00849 case base::Front:
00850
00851 for (register int j=mys; j<=mye; j++)
00852 for (register int i=mxs; i<=mxe; i++) {
00853 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);
00854 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
00855 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);
00856
00857 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);
00858 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);
00859 }
00860 break;
00861 case base::Back:
00862
00863 for (register int j=mys; j<=mye; j++)
00864 for (register int i=mxs; i<=mxe; i++) {
00865
00866
00867
00868
00869 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);
00870 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
00871 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);
00872
00873 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);
00874 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);
00875 }
00876 break;
00877 }
00878 }
00879
00880 else if (type==Inlet) {
00881 if (aux==0 || naux<=0) return;
00882 DataType a0=S0*aux[0], a1=0., a2=0.;
00883 if (naux>1) {
00884 a1 = S0*aux[1];
00885 a2 = S0*aux[2];
00886 }
00887 if (scaling==base::Physical) {
00888 a0 /= U0;
00889 a1 /= U0;
00890 a2 /= U0;
00891 }
00892 switch (side) {
00893 case base::Left:
00894 for (register int k=mzs; k<=mze; k++)
00895 for (register int j=mys; j<=mye; j++) {
00896 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00897 q(1) = a0;
00898 if (naux>1) {
00899 q(2) = a1;
00900 q(3) = a2;
00901 }
00902 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
00903 }
00904 break;
00905 case base::Right:
00906 for (register int k=mzs; k<=mze; k++)
00907 for (register int j=mys; j<=mye; j++) {
00908 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00909 q(1) = a0;
00910 if (naux>1) {
00911 q(2) = a1;
00912 q(3) = a2;
00913 }
00914 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
00915 }
00916 break;
00917 case base::Bottom:
00918 for (register int k=mzs; k<=mze; k++)
00919 for (register int i=mxs; i<=mxe; i++) {
00920 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
00921 if (naux==1)
00922 q(2) = a0;
00923 else {
00924 q(1) = a0;
00925 q(2) = a1;
00926 q(3) = a2;
00927 }
00928 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
00929 }
00930 break;
00931 case base::Top:
00932 for (register int k=mzs; k<=mze; k++)
00933 for (register int i=mxs; i<=mxe; i++) {
00934 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
00935 if (naux==1)
00936 q(2) = a0;
00937 else {
00938 q(1) = a0;
00939 q(2) = a1;
00940 q(3) = a2;
00941 }
00942 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
00943 }
00944 break;
00945 case base::Front:
00946 for (register int j=mys; j<=mye; j++)
00947 for (register int i=mxs; i<=mxe; i++) {
00948 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
00949 if (naux==1)
00950 q(3) = a0;
00951 else {
00952 q(1) = a1;
00953 q(2) = a2;
00954 q(3) = a0;
00955 }
00956 f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
00957 }
00958 break;
00959 case base::Back:
00960 for (register int j=mys; j<=mye; j++)
00961 for (register int i=mxs; i<=mxe; i++) {
00962 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
00963 if (naux==1)
00964 q(3) = a0;
00965 else {
00966 q(1) = a1;
00967 q(2) = a2;
00968 q(3) = a0;
00969 }
00970 f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
00971 }
00972 break;
00973 }
00974 }
00975
00976 else if (type==Outlet) {
00977
00978 switch (side) {
00979 case base::Left:
00980
00981 for (register int k=mzs; k<=mze; k++)
00982 for (register int j=mys; j<=mye; j++) {
00983 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
00984
00985 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);
00986 }
00987 break;
00988 case base::Right:
00989
00990 for (register int k=mzs; k<=mze; k++)
00991 for (register int j=mys; j<=mye; j++) {
00992 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
00993
00994 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);
00995 }
00996 break;
00997 case base::Bottom:
00998 for (register int k=mzs; k<=mze; k++)
00999 for (register int i=mxs; i<=mxe; i++) {
01000 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01001 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);
01002 }
01003 break;
01004 case base::Top:
01005 for (register int k=mzs; k<=mze; k++)
01006 for (register int i=mxs; i<=mxe; i++) {
01007 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01008 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);
01009 }
01010 break;
01011 case base::Front:
01012 for (register int j=mys; j<=mye; j++)
01013 for (register int i=mxs; i<=mxe; i++) {
01014 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01015 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);
01016 }
01017 break;
01018 case base::Back:
01019 for (register int j=mys; j<=mye; j++)
01020 for (register int i=mxs; i<=mxe; i++) {
01021 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01022 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);
01023 }
01024 break;
01025 }
01026 }
01027
01028
01029 else if (type==Pressure) {
01030 if (aux==0 || naux<=0) return;
01031 DataType a0=aux[0];
01032 if (scaling==base::Physical)
01033 a0 /= R0;
01034 switch (side) {
01035 case base::Left:
01036 for (register int k=mzs; k<=mze; k++)
01037 for (register int j=mys; j<=mye; j++) {
01038 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01039 q(0) = a0;
01040 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(q);
01041 }
01042 break;
01043 case base::Right:
01044 for (register int k=mzs; k<=mze; k++)
01045 for (register int j=mys; j<=mye; j++) {
01046 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01047 q(0) = a0;
01048 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(q);
01049 }
01050 break;
01051 case base::Bottom:
01052 for (register int k=mzs; k<=mze; k++)
01053 for (register int i=mxs; i<=mxe; i++) {
01054 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01055 q(0) = a0;
01056 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(q);
01057 }
01058 break;
01059 case base::Top:
01060 for (register int k=mzs; k<=mze; k++)
01061 for (register int i=mxs; i<=mxe; i++) {
01062 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01063 q(0) = a0;
01064 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(q);
01065 }
01066 break;
01067 case base::Front:
01068 for (register int j=mys; j<=mye; j++)
01069 for (register int i=mxs; i<=mxe; i++) {
01070 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01071 q(0) = a0;
01072 f[b+base::idx(i,j,mze+1,Nx,Ny)] = Equilibrium(q);
01073 }
01074 break;
01075 case base::Back:
01076 for (register int j=mys; j<=mye; j++)
01077 for (register int i=mxs; i<=mxe; i++) {
01078 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01079 q(0) = a0;
01080 f[b+base::idx(i,j,mzs-1,Nx,Ny)] = Equilibrium(q);
01081 }
01082 break;
01083 }
01084 }
01085
01086 else if (type==SlidingWall) {
01087 if (aux==0 || naux<=0) return;
01088 if (aux==0 || naux<=0) return;
01089 DataType a0=S0*aux[0];
01090 if (scaling==base::Physical)
01091 a0 /= U0;
01092 switch (side) {
01093 case base::Left:
01094 for (register int k=mzs; k<=mze; k++)
01095 for (register int j=mys; j<=mye; j++) {
01096 MacroType q = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01097 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);
01098 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01099 DataType pw = DataType(1.)-qw;
01100 if (j<mye) f[b+base::idx(mxe,j+1,k,Nx,Ny)](9) = pw*rd1+qw*rd2;
01101 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01102 if (j>mys) f[b+base::idx(mxe,j-1,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01103
01104 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);
01105 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);
01106 }
01107 break;
01108 case base::Right:
01109 for (register int k=mzs; k<=mze; k++)
01110 for (register int j=mys; j<=mye; j++) {
01111 MacroType q = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01112 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);
01113 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01114 DataType pw = DataType(1.)-qw;
01115 if (j<mye) f[b+base::idx(mxs,j+1,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01116 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01117 if (j>mys) f[b+base::idx(mxs,j-1,k,Nx,Ny)](10) = qw*rd1+pw*rd2;
01118
01119 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);
01120 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);
01121 }
01122 break;
01123 case base::Bottom:
01124 for (register int k=mzs; k<=mze; k++)
01125 for (register int i=mxs; i<=mxe; i++) {
01126 MacroType q = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01127 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);
01128 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01129 DataType pw = DataType(1.)-qw;
01130 if (i<mxe) f[b+base::idx(i+1,mye,k,Nx,Ny)](10) = pw*rd1+qw*rd2;
01131 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01132 if (i>mxs) f[b+base::idx(i-1,mye,k,Nx,Ny)](7) = qw*rd1+pw*rd2;
01133
01134 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);
01135 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);
01136 }
01137 break;
01138 case base::Top:
01139 for (register int k=mzs; k<=mze; k++)
01140 for (register int i=mxs; i<=mxe; i++) {
01141 MacroType q = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01142 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);
01143 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01144 DataType pw = DataType(1.)-qw;
01145 if (i<mxe) f[b+base::idx(i+1,mys,k,Nx,Ny)](8) = pw*rd1+qw*rd2;
01146 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01147 if (i>mxs) f[b+base::idx(i-1,mys,k,Nx,Ny)](9) = qw*rd1+pw*rd2;
01148
01149 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);
01150 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);
01151 }
01152 break;
01153 case base::Front:
01154 for (register int j=mys; j<=mye; j++)
01155 for (register int i=mxs; i<=mxe; i++) {
01156 MacroType q = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01157 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);
01158 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01159 DataType pw = DataType(1.)-qw;
01160 if (i<mxe) f[b+base::idx(i+1,j,mze,Nx,Ny)](18) = pw*rd1+qw*rd2;
01161 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01162 if (i>mxs) f[b+base::idx(i-1,j,mze,Nx,Ny)](15) = qw*rd1+pw*rd2;
01163
01164 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);
01165 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);
01166 }
01167 break;
01168 case base::Back:
01169 for (register int j=mys; j<=mye; j++)
01170 for (register int i=mxs; i<=mxe; i++) {
01171 MacroType q = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01172 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);
01173 DataType qw = (q(0)/DataType(6.)*a0)/(rd1-rd2);
01174 DataType pw = DataType(1.)-qw;
01175 if (i<mxe) f[b+base::idx(i+1,j,mzs,Nx,Ny)](16) = pw*rd1+qw*rd2;
01176 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
01177 if (i>mxs) f[b+base::idx(i-1,j,mzs,Nx,Ny)](17) = qw*rd1+pw*rd2;
01178
01179 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);
01180 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);
01181 }
01182 break;
01183 }
01184 }
01185
01189 else if (type==CharacteristicOutlet) {
01190 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
01191 DCoords dx = base::GH().worldStep(fvec.stepsize());
01192 DataType dt_temp = dx(0)/VelocityScale();
01193 switch (side) {
01194 case base::Left:
01195 for (register int k=mzs; k<=mze; k++)
01196 for (register int j=mys; j<=mye; j++) {
01197 MicroType ftmp = f[b+base::idx(mxe,j,k,Nx,Ny)];
01198 MacroType q = MacroVariables(ftmp);
01199 Collision(ftmp,dt_temp);
01200 MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01201 MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,k,Nx,Ny)]);
01202 MacroType dqdn = NormalDerivative(q,q1,q2);
01203 Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01204
01205
01206 MacroType qn;
01207 qn(0) = q(0) - c1oCsSq2*L(0);
01208 qn(1) = q(1) + c1o2cs*L(0)/q(0);
01209 qn(2) = q(2) - L(1);
01210 qn(3) = q(3) - L(2);
01211
01212 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(qn);
01213 }
01214 break;
01215 case base::Right:
01216 for (register int k=mzs; k<=mze; k++)
01217 for (register int j=mys; j<=mye; j++) {
01218 MicroType ftmp = f[b+base::idx(mxs,j,k,Nx,Ny)];
01219 MacroType q = MacroVariables(ftmp);
01220 Collision(ftmp,dt_temp);
01221 MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01222 MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,k,Nx,Ny)]);
01223 MacroType dqdn = NormalDerivative(q,q1,q2);
01224 Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01225
01226
01227 MacroType qn;
01228 qn(0) = q(0) - c1oCsSq2*L(0);
01229 qn(1) = q(1) + c1o2cs*L(0)/q(0);
01230 qn(2) = q(2) - L(1);
01231 qn(3) = q(3) - L(2);
01232
01233 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(qn);
01234 }
01235 break;
01236 case base::Bottom:
01237 for (register int k=mzs; k<=mze; k++)
01238 for (register int i=mxs; i<=mxe; i++) {
01239 MicroType ftmp = f[b+base::idx(i,mye,k,Nx,Ny)];
01240 MacroType q = MacroVariables(ftmp);
01241 Collision(ftmp,dt_temp);
01242 MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01243 MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,k,Nx,Ny)]);
01244 MacroType dqdn = NormalDerivative(q,q1,q2);
01245 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01246
01247
01248 MacroType qn;
01249 qn(0) = q(0) - c1oCsSq2*L(0);
01250 qn(1) = q(1) - L(1);
01251 qn(2) = q(2) + c1o2cs*L(0)/q(0);
01252 qn(3) = q(3) - L(2);
01253
01254 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(qn);
01255 }
01256 break;
01257 case base::Top:
01258 for (register int k=mzs; k<=mze; k++)
01259 for (register int i=mxs; i<=mxe; i++) {
01260 MicroType ftmp = f[b+base::idx(i,mys,k,Nx,Ny)];
01261 MacroType q = MacroVariables(ftmp);
01262 Collision(ftmp,dt_temp);
01263 MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01264 MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,k,Nx,Ny)]);
01265 MacroType dqdn = NormalDerivative(q,q1,q2);
01266 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01267
01268
01269 MacroType qn;
01270 qn(0) = q(0) - c1oCsSq2*L(0);
01271 qn(1) = q(1) - L(1);
01272 qn(2) = q(2) + c1o2cs*L(0)/q(0);
01273 qn(3) = q(3) - L(2);
01274
01275 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(qn);
01276 }
01277 break;
01278 case base::Front:
01279 for (register int j=mys; j<=mye; j++)
01280 for (register int i=mxs; i<=mxe; i++) {
01281 MicroType ftmp = f[b+base::idx(i,j,mze,Nx,Ny)];
01282 MacroType q = MacroVariables(ftmp);
01283 Collision(ftmp,dt_temp);
01284 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01285 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mze+2,Nx,Ny)]);
01286 MacroType dqdn = NormalDerivative(q,q1,q2);
01287 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01288
01289
01290 MacroType qn;
01291 qn(0) = q(0) - c1oCsSq2*L(0);
01292 qn(1) = q(1) - L(1);
01293 qn(2) = q(2) - L(2);
01294 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01295
01296 f[b+base::idx(i,j,mze,Nx,Ny)] = Equilibrium(qn);
01297 }
01298 break;
01299 case base::Back:
01300 for (register int j=mys; j<=mye; j++)
01301 for (register int i=mxs; i<=mxe; i++) {
01302 MicroType ftmp = f[b+base::idx(i,j,mzs,Nx,Ny)];
01303 MacroType q = MacroVariables(ftmp);
01304 Collision(ftmp,dt_temp);
01305 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01306 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mzs-2,Nx,Ny)]);
01307 MacroType dqdn = NormalDerivative(q,q1,q2);
01308 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01309
01310
01311 MacroType qn;
01312 qn(0) = q(0) - c1oCsSq2*L(0);
01313 qn(1) = q(1) - L(1);
01314 qn(2) = q(2) - L(2);
01315 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01316
01317 f[b+base::idx(i,j,mzs,Nx,Ny)] = Equilibrium(qn);
01318 }
01319 break;
01320 }
01321 }
01322
01326 else if (type==CharacteristicInlet) {
01327 if (aux==0 || naux<=0) return;
01328 DataType rho=aux[0];
01329 DataType a0=S0*aux[1];
01330 DataType a1=S0*aux[2];
01331 DataType a2=S0*aux[3];
01332 if (scaling==base::Physical) {
01333 rho /= R0;
01334 a0 /= U0;
01335 a1 /= U0;
01336 a2 /= U0;
01337 }
01338 MacroType q, qn;
01339 q(0) = rho;
01340 q(1) = a0;
01341 q(2) = a1;
01342 q(3) = a2;
01343 DataType c1oCsSq2 = DataType(0.5)/cs2, c1o2cs = DataType(0.5)/LatticeSpeedOfSound();
01344 switch (side) {
01345 case base::Left:
01346 for (register int k=mzs; k<=mze; k++)
01347 for (register int j=mys; j<=mye; j++) {
01348 MacroType q1 = MacroVariables(f[b+base::idx(mxe+1,j,k,Nx,Ny)]);
01349 MacroType q2 = MacroVariables(f[b+base::idx(mxe+2,j,k,Nx,Ny)]);
01350 Vector<DataType,3> L;
01351 if (naux==1) {
01352 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01353 MacroType dqdn = NormalDerivative(q,q1,q2);
01354 L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01355 }
01356 else if (naux==2) {
01357 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01358 MacroType dqdn = NormalDerivative(q,q1,q2);
01359 L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01360 }
01361 else if (naux==4) {
01362 q(0) = q1(0);
01363 MacroType dqdn = NormalDerivative(q,q1,q2);
01364 L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01365 }
01366 MacroType qn;
01367 qn(0) = q(0) - c1oCsSq2*L(0);
01368 qn(1) = q(1) + c1o2cs*L(0)/q(0);
01369 qn(2) = q(2) - L(1);
01370 qn(3) = q(3) - L(2);
01371 f[b+base::idx(mxe,j,k,Nx,Ny)] = Equilibrium(qn);
01372 }
01373 break;
01374 case base::Right:
01375 for (register int k=mzs; k<=mze; k++)
01376 for (register int j=mys; j<=mye; j++) {
01377 MacroType q1 = MacroVariables(f[b+base::idx(mxs-1,j,k,Nx,Ny)]);
01378 MacroType q2 = MacroVariables(f[b+base::idx(mxs-2,j,k,Nx,Ny)]);
01379 if (naux==1) {
01380 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01381 MacroType dqdn = NormalDerivative(q,q1,q2);
01382 Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01383 MacroType qn;
01384 qn(0) = q(0) - c1oCsSq2*L(0);
01385 qn(1) = q1(1) + c1o2cs*L(0)/q(0);
01386 qn(2) = q1(2) - L(1);
01387 qn(3) = q1(3) - L(2);
01388 }
01389 else if (naux==2) {
01390 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01391 MacroType dqdn = NormalDerivative(q,q1,q2);
01392 Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01393 MacroType qn;
01394 qn(0) = q1(0) - c1oCsSq2*L(0);
01395 qn(1) = q(1) + c1o2cs*L(0)/q(0);
01396 qn(2) = q1(2) - L(1);
01397 qn(3) = q1(3) - L(2);
01398 }
01399 else if (naux==4) {
01400 q(0) = q1(0);
01401 MacroType dqdn = NormalDerivative(q,q1,q2);
01402 Vector<DataType,3> L = WaveAmplitudes(q(0),q(1),dqdn(0),dqdn(1),dqdn(2),dqdn(3));
01403 MacroType qn;
01404 qn(0) = q1(0) - c1oCsSq2*L(0);
01405 qn(1) = q(1) + c1o2cs*L(0)/q(0);
01406 qn(2) = q(2) - L(1);
01407 qn(3) = q(3) - L(2);
01408 }
01409 f[b+base::idx(mxs,j,k,Nx,Ny)] = Equilibrium(qn);
01410 }
01411 break;
01412 case base::Bottom:
01413 for (register int k=mzs; k<=mze; k++)
01414 for (register int i=mxs; i<=mxe; i++) {
01415 MacroType q1 = MacroVariables(f[b+base::idx(i,mye+1,k,Nx,Ny)]);
01416 MacroType q2 = MacroVariables(f[b+base::idx(i,mye+2,k,Nx,Ny)]);
01417 if (naux==1) {
01418 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01419 MacroType dqdn = NormalDerivative(q,q1,q2);
01420 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01421 MacroType qn;
01422 qn(0) = q(0) - c1oCsSq2*L(0);
01423 qn(1) = q1(1) - L(1);
01424 qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01425 qn(3) = q1(3) - L(2);
01426 }
01427 else if (naux==2) {
01428 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01429 MacroType dqdn = NormalDerivative(q,q1,q2);
01430 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01431 MacroType qn;
01432 qn(0) = q1(0) - c1oCsSq2*L(0);
01433 qn(1) = q(1) - L(1);
01434 qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01435 qn(3) = q1(3) - L(2);
01436 }
01437 else if (naux==4) {
01438 q(0) = q1(0);
01439 MacroType dqdn = NormalDerivative(q,q1,q2);
01440 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01441 MacroType qn;
01442 qn(0) = q1(0) - c1oCsSq2*L(0);
01443 qn(1) = q(1) - L(1);
01444 qn(2) = q(2) + c1o2cs*L(0)/q(0);
01445 qn(3) = q(3) - L(2);
01446 }
01447 f[b+base::idx(i,mye,k,Nx,Ny)] = Equilibrium(qn);
01448 }
01449 break;
01450 case base::Top:
01451 for (register int k=mzs; k<=mze; k++)
01452 for (register int i=mxs; i<=mxe; i++) {
01453 MacroType q1 = MacroVariables(f[b+base::idx(i,mys-1,k,Nx,Ny)]);
01454 MacroType q2 = MacroVariables(f[b+base::idx(i,mys-2,k,Nx,Ny)]);
01455 if (naux==1) {
01456 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01457 MacroType dqdn = NormalDerivative(q,q1,q2);
01458 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01459 MacroType qn;
01460 qn(0) = q(0) - c1oCsSq2*L(0);
01461 qn(1) = q1(1) - L(1);
01462 qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01463 qn(3) = q1(3) - L(2);
01464 }
01465 else if (naux==2) {
01466 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01467 MacroType dqdn = NormalDerivative(q,q1,q2);
01468 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01469 MacroType qn;
01470 qn(0) = q1(0) - c1oCsSq2*L(0);
01471 qn(1) = q(1) - L(1);
01472 qn(2) = q1(2) + c1o2cs*L(0)/q(0);
01473 qn(3) = q1(3) - L(2);
01474 }
01475 else if (naux==4) {
01476 q(0) = q1(0);
01477 MacroType dqdn = NormalDerivative(q,q1,q2);
01478 Vector<DataType,3> L = WaveAmplitudes(q(0),q(2),dqdn(0),dqdn(2),dqdn(1),dqdn(3));
01479 MacroType qn;
01480 qn(0) = q1(0) - c1oCsSq2*L(0);
01481 qn(1) = q(1) - L(1);
01482 qn(2) = q(2) + c1o2cs*L(0)/q(0);
01483 qn(3) = q(3) - L(2);
01484 }
01485 f[b+base::idx(i,mys,k,Nx,Ny)] = Equilibrium(qn);
01486 }
01487 break;
01488 case base::Front:
01489 for (register int j=mys; j<=mye; j++)
01490 for (register int i=mxs; i<=mxe; i++) {
01491 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mze+1,Nx,Ny)]);
01492 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mze+2,Nx,Ny)]);
01493 if (naux==1) {
01494 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01495 MacroType dqdn = NormalDerivative(q,q1,q2);
01496 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01497 MacroType qn;
01498 qn(0) = q(0) - c1oCsSq2*L(0);
01499 qn(1) = q(1) - L(1);
01500 qn(2) = q(2) - L(2);
01501 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01502 }
01503 else if (naux==2) {
01504 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01505 MacroType dqdn = NormalDerivative(q,q1,q2);
01506 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01507 MacroType qn;
01508 qn(0) = q(0) - c1oCsSq2*L(0);
01509 qn(1) = q(1) - L(1);
01510 qn(2) = q(2) - L(2);
01511 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01512 }
01513 else if (naux==4) {
01514 q(0) = q1(0);
01515 MacroType dqdn = NormalDerivative(q,q1,q2);
01516 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01517 MacroType qn;
01518 qn(0) = q(0) - c1oCsSq2*L(0);
01519 qn(1) = q(1) - L(1);
01520 qn(2) = q(2) - L(2);
01521 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01522 }
01523 MacroType dqdn = NormalDerivative(q,q1,q2);
01524 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01525
01526
01527 MacroType qn;
01528 qn(0) = q(0) - c1oCsSq2*L(0);
01529 qn(1) = q(1) - L(1);
01530 qn(2) = q(2) - L(2);
01531 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01532
01533 f[b+base::idx(i,j,mze,Nx,Ny)] = Equilibrium(qn);
01534 }
01535 break;
01536 case base::Back:
01537 for (register int j=mys; j<=mye; j++)
01538 for (register int i=mxs; i<=mxe; i++) {
01539 MacroType q1 = MacroVariables(f[b+base::idx(i,j,mzs-1,Nx,Ny)]);
01540 MacroType q2 = MacroVariables(f[b+base::idx(i,j,mzs-2,Nx,Ny)]);
01541 if (naux==1) {
01542 q(1) = q1(1); q(2) = q1(2); q(3) = q1(3);
01543 MacroType dqdn = NormalDerivative(q,q1,q2);
01544 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01545 MacroType qn;
01546 qn(0) = q(0) - c1oCsSq2*L(0);
01547 qn(1) = q(1) - L(1);
01548 qn(2) = q(2) - L(2);
01549 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01550 }
01551 else if (naux==2) {
01552 q(0) = q1(0); q(2) = q1(2); q(3) = q1(3);
01553 MacroType dqdn = NormalDerivative(q,q1,q2);
01554 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01555 MacroType qn;
01556 qn(0) = q(0) - c1oCsSq2*L(0);
01557 qn(1) = q(1) - L(1);
01558 qn(2) = q(2) - L(2);
01559 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01560 }
01561 else if (naux==4) {
01562 q(0) = q1(0);
01563 MacroType dqdn = NormalDerivative(q,q1,q2);
01564 Vector<DataType,3> L = WaveAmplitudes(q(0),q(3),dqdn(0),dqdn(3),dqdn(1),dqdn(2));
01565 MacroType qn;
01566 qn(0) = q(0) - c1oCsSq2*L(0);
01567 qn(1) = q(1) - L(1);
01568 qn(2) = q(2) - L(2);
01569 qn(3) = q(3) + c1o2cs*L(0)/q(0);
01570 }
01571 f[b+base::idx(i,j,mzs,Nx,Ny)] = Equilibrium(qn);
01572 }
01573 break;
01574 }
01575 }
01576
01577 else if (type==NoSlipWallEntranceExit) {
01578 int lc_indx[3] = { mxs*fvec.stepsize(0), mys*fvec.stepsize(1), mzs*fvec.stepsize(2) };
01579 DataType WLB[3], WUB[3];
01580 base::GH().wlb(WLB);
01581 base::GH().wub(WUB);
01582 DCoords wc;
01583 DataType slip = 0.;
01584 DataType lxs, lxe, lys, lye, lzs, lze;
01585 DataType min[3], max[3];
01586 if (naux!=6) assert(0);
01587 lxs = aux[0] - WLB[0];
01588 lxe = WUB[0] - aux[1];
01589 lys = aux[2] - WLB[1];
01590 lye = WUB[1] - aux[3];
01591 lzs = aux[4] - WLB[2];
01592 lze = WUB[2] - aux[5];
01593 min[0] = aux[0]; max[0] = aux[1];
01594 min[1] = aux[2]; max[1] = aux[3];
01595 min[2] = aux[4]; max[2] = aux[5];
01596
01597 switch (side) {
01598 case base::Left:
01599 lc_indx[0] = mxe*fvec.stepsize(0);
01600 for (register int k=mzs; k<=mze; k++) {
01601 lc_indx[2] = k*fvec.stepsize(2);
01602 for (register int j=mys; j<=mye; j++) {
01603 lc_indx[1] = j*fvec.stepsize(1);
01604 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01605 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01606 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);
01607 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);
01608 }
01609 else {
01610 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01611 slip = DataType(0.);
01612 else if (wc(1)<min[1])
01613 slip = (wc(1) - WLB[1])/lys;
01614 else if (wc(1)>max[1])
01615 slip = (WUB[1] - wc(1))/lye;
01616 if (j>mys && j<mye) {
01617 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)
01618 + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01619 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)
01620 + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01621 }
01622 else if (j==mys) {
01623 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)
01624 + slip*f[b+base::idx(mxe+1,j+1,k,Nx,Ny)](8);
01625 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)
01626 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](10);
01627 }
01628 else if (j==mye) {
01629 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)
01630 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](8);
01631 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)
01632 + slip*f[b+base::idx(mxe+1,j-1,k,Nx,Ny)](10);
01633 }
01634 }
01635 if (wc(2)>=min[2] && wc(2)<=max[2]) {
01636 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);
01637 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);
01638 }
01639 else {
01640 if (wc(2)<WLB[2] || wc(2) >WUB[2])
01641 slip = DataType(0.);
01642 else if (wc(2)<min[2])
01643 slip = (wc(2) - WLB[2])/lzs;
01644 else if (wc(2)>max[2])
01645 slip = (WUB[2] - wc(2))/lze;
01646 if (k>mzs && k<mze) {
01647 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)
01648 + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01649 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)
01650 + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01651 }
01652 else if (k==mzs) {
01653 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)
01654 + slip*f[b+base::idx(mxe+1,j,k+1,Nx,Ny)](16);
01655 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)
01656 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](18);
01657 }
01658 else if (k==mze) {
01659 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)
01660 + slip*f[b+base::idx(mxe+1,j,k,Nx,Ny)](16);
01661 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)
01662 + slip*f[b+base::idx(mxe+1,j,k-1,Nx,Ny)](18);
01663 }
01664 }
01665 f[b+base::idx(mxe,j,k,Nx,Ny)](1) = f[b+base::idx(mxe+1,j,k,Nx,Ny)](2);
01666 }
01667 }
01668 break;
01669 case base::Right:
01670 lc_indx[0] = mxs*fvec.stepsize(0);
01671 for (register int k=mzs; k<=mze; k++) {
01672 lc_indx[2] = k*fvec.stepsize(2);
01673 for (register int j=mys; j<=mye; j++) {
01674 lc_indx[1] = j*fvec.stepsize(1);
01675 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01676 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01677 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);
01678 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);
01679 }
01680
01681 else {
01682 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01683 slip = DataType(0.);
01684 else if (wc(1)<min[1])
01685 slip = (wc(1) - WLB[1])/lys;
01686 else if (wc(1)>max[1])
01687 slip = (WUB[1] - wc(1))/lye;
01688 if (j>mys && j<mye) {
01689 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)
01690 + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
01691 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)
01692 + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
01693 }
01694 else if (j==mys) {
01695 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)
01696 + slip*f[b+base::idx(mxs-1,j+1,k,Nx,Ny)](9);
01697 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)
01698 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](7);
01699 }
01700 else if (j==mye) {
01701 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)
01702 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](9);
01703 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)
01704 + slip*f[b+base::idx(mxs-1,j-1,k,Nx,Ny)](7);
01705 }
01706 }
01707 if (wc(2)>=min[2] && wc(2)<=max[2]) {
01708 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);
01709 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);
01710 }
01711 else {
01712 if (wc(2)<WLB[2] || wc(2) >WUB[2])
01713 slip = DataType(0.);
01714 else if (wc(2)<min[2])
01715 slip = (wc(2) - WLB[2])/lzs;
01716 else if (wc(2)>max[2])
01717 slip = (WUB[2] - wc(2))/lze;
01718 if (k>mzs && k<mze) {
01719 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)
01720 + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01721 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)
01722 + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01723 }
01724 else if (k==mzs) {
01725 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)
01726 + slip*f[b+base::idx(mxs-1,j,k+1,Nx,Ny)](17);
01727 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)
01728 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](15);
01729 }
01730 else if (k==mze) {
01731 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)
01732 + slip*f[b+base::idx(mxs-1,j,k,Nx,Ny)](17);
01733 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)
01734 + slip*f[b+base::idx(mxs-1,j,k-1,Nx,Ny)](15);
01735 }
01736 }
01737 f[b+base::idx(mxs,j,k,Nx,Ny)](2) = f[b+base::idx(mxs-1,j,k,Nx,Ny)](1);
01738 }
01739 }
01740 break;
01741 case base::Bottom:
01742 lc_indx[1] = mye*fvec.stepsize(1);
01743 for (register int k=mzs; k<=mze; k++) {
01744 lc_indx[2] = k*fvec.stepsize(2);
01745 for (register int i=mxs; i<=mxe; i++) {
01746 lc_indx[0] = i*fvec.stepsize(0);
01747 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01748 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01749 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);
01750 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);
01751 }
01752 else {
01753 if (wc(0)<WLB[0] || wc(0)>WUB[0])
01754 slip = DataType(0.);
01755 else if (wc(0)<min[0])
01756 slip = (wc(0) - WLB[0])/lxs;
01757 else if (wc(0)>max[0])
01758 slip = (WUB[0] - wc(0))/lxe;
01759 if (i>mxs && i<mxe) {
01760 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)
01761 + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
01762 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)
01763 + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
01764 }
01765 else if (i==mxs) {
01766 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)
01767 + slip*f[b+base::idx(i+1,mye+1,k,Nx,Ny)](8);
01768 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)
01769 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](9);
01770 }
01771 else if (i==mxe) {
01772 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)
01773 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](8);
01774 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)
01775 + slip*f[b+base::idx(i-1,mye+1,k,Nx,Ny)](9);
01776 }
01777 }
01778 if (wc(2)>=min[2] && wc(2)<=max[2]) {
01779 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);
01780 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);
01781 }
01782 else {
01783 if (wc(2)<WLB[2] || wc(2) >WUB[2])
01784 slip = DataType(0.);
01785 else if (wc(2)<min[2])
01786 slip = (wc(2) - WLB[2])/lzs;
01787 else if (wc(2)>max[2])
01788 slip = (WUB[2] - wc(2))/lze;
01789 if (k>mzs && k<mze) {
01790 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)
01791 + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01792 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)
01793 + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01794 }
01795 else if (k==mzs) {
01796 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)
01797 + slip*f[b+base::idx(i,mye+1,k+1,Nx,Ny)](12);
01798 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)
01799 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](14);
01800 }
01801 else if (k==mze) {
01802 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)
01803 + slip*f[b+base::idx(i,mye+1,k,Nx,Ny)](12);
01804 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)
01805 + slip*f[b+base::idx(i,mye+1,k-1,Nx,Ny)](14);
01806 }
01807 }
01808 f[b+base::idx(i,mye,k,Nx,Ny)](3) = f[b+base::idx(i,mye+1,k,Nx,Ny)](4);
01809 }
01810 }
01811 break;
01812 case base::Top:
01813 lc_indx[1] = mys*fvec.stepsize(1);
01814 for (register int k=mzs; k<=mze; k++) {
01815 lc_indx[2] = k*fvec.stepsize(2);
01816 for (register int i=mxs; i<=mxe; i++) {
01817 lc_indx[0] = i*fvec.stepsize(0);
01818 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01819 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01820 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);
01821 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);
01822 }
01823 else {
01824 if (wc(0)<WLB[0] || wc(0) >WUB[0])
01825 slip = DataType(0.);
01826 else if (wc(0)<min[0])
01827 slip = (wc(0) - WLB[0])/lxs;
01828 else if (wc(0)>max[0])
01829 slip = (WUB[0] - wc(0))/lxe;
01830 if (i>mxs && i<mxe) {
01831 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)
01832 + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
01833 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)
01834 + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
01835 }
01836 else if (i==mxs) {
01837 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)
01838 + slip*f[b+base::idx(i+1,mys-1,k,Nx,Ny)](10);
01839 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)
01840 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](7);
01841 }
01842 else if (i==mxe) {
01843 f[b+base::idx(i,mys,k,Nx,Ny)](8) = f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7)
01844 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](10);
01845 f[b+base::idx(i,mys,k,Nx,Ny)](9) = f[b+base::idx(i,mys-1,k,Nx,Ny)](10)
01846 + slip*f[b+base::idx(i-1,mys-1,k,Nx,Ny)](7);
01847 }
01848 }
01849 if (wc(2)>=min[2] && wc(2)<=max[2]) {
01850 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);
01851 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);
01852 }
01853 else {
01854 if (wc(2)<WLB[2] || wc(2)>WUB[2])
01855 slip = DataType(0.);
01856 else if (wc(2)<min[2])
01857 slip = (wc(2) - WLB[2])/lzs;
01858 else if (wc(2)>max[2])
01859 slip = (WUB[2] - wc(2))/lze;
01860 if (k>mzs && k<mze) {
01861 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)
01862 + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01863 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)
01864 + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01865 }
01866 else if (k==mzs) {
01867 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)
01868 + slip*f[b+base::idx(i,mys-1,k+1,Nx,Ny)](13);
01869 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)
01870 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](11);
01871 }
01872 else if (k==mze) {
01873 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)
01874 + slip*f[b+base::idx(i,mys-1,k,Nx,Ny)](13);
01875 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)
01876 + slip*f[b+base::idx(i,mys-1,k-1,Nx,Ny)](11);
01877 }
01878 }
01879 f[b+base::idx(i,mys,k,Nx,Ny)](4) = f[b+base::idx(i,mys-1,k,Nx,Ny)](3);
01880 }
01881 }
01882 break;
01883 case base::Front:
01884 lc_indx[2] = mze*fvec.stepsize(2);
01885 for (register int j=mys; j<=mye; j++) {
01886 lc_indx[1] = j*fvec.stepsize(1);
01887 for (register int i=mxs; i<=mxe; i++) {
01888 lc_indx[0] = i*fvec.stepsize(0);
01889 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01890 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01891 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);
01892 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);
01893 }
01894 else {
01895 if (wc(0)<WLB[0] || wc(0) >WUB[0])
01896 slip = DataType(0.);
01897 else if (wc(0)<min[0])
01898 slip = (wc(0) - WLB[0])/lxs;
01899 else if (wc(0)>max[0])
01900 slip = (WUB[0] - wc(0))/lxe;
01901 if (i>mxs && i<mxe) {
01902 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)
01903 + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
01904 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)
01905 + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
01906 }
01907 else if (i==mxs) {
01908 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)
01909 + slip*f[b+base::idx(i+1,j,mze+1,Nx,Ny)](16);
01910 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)
01911 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](17);
01912 }
01913 else if (i==mxe) {
01914 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)
01915 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](16);
01916 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)
01917 + slip*f[b+base::idx(i-1,j,mze+1,Nx,Ny)](17);
01918 }
01919 }
01920 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01921 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);
01922 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);
01923 }
01924 else {
01925 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01926 slip = DataType(0.);
01927 else if (wc(1)<min[1])
01928 slip = (wc(1) - WLB[1])/lys;
01929 else if (wc(1)>max[1])
01930 slip = (WUB[1] - wc(1))/lye;
01931 if (j>mys && j<mye) {
01932 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)
01933 + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01934 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)
01935 + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01936 }
01937 else if (j==mys) {
01938 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)
01939 + slip*f[b+base::idx(i,j+1,mze+1,Nx,Ny)](12);
01940 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)
01941 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](13);
01942 }
01943 else if (j==mye) {
01944 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)
01945 + slip*f[b+base::idx(i,j,mze+1,Nx,Ny)](12);
01946 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)
01947 + slip*f[b+base::idx(i,j-1,mze+1,Nx,Ny)](13);
01948 }
01949 }
01950 f[b+base::idx(i,j,mze,Nx,Ny)](5) = f[b+base::idx(i,j,mze+1,Nx,Ny)](6);
01951 }
01952 }
01953 break;
01954 case base::Back:
01955 lc_indx[2] = mzs*fvec.stepsize(2);
01956 for (register int j=mys; j<=mye; j++) {
01957 lc_indx[1] = j*fvec.stepsize(1);
01958 for (register int i=mxs; i<=mxe; i++) {
01959 lc_indx[0] = i*fvec.stepsize(0);
01960 wc = base::GH().worldCoords(lc_indx,fvec.stepsize());
01961 if (wc(0)>=min[0] && wc(0)<=max[0]) {
01962 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);
01963 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);
01964 }
01965 else {
01966 if (wc(0)<WLB[0] || wc(0) >WUB[0])
01967 slip = DataType(0.);
01968 else if (wc(0)<min[0])
01969 slip = (wc(0) - WLB[0])/lxs;
01970 else if (wc(0)>max[0])
01971 slip = (WUB[0] - wc(0))/lxe;
01972 if (i>mxs && i<mxe) {
01973 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)
01974 + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
01975 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)
01976 + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
01977 }
01978 else if (i==mxs) {
01979 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)
01980 + slip*f[b+base::idx(i+1,j,mzs-1,Nx,Ny)](18);
01981 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)
01982 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](15);
01983 }
01984 else if (i==mxe) {
01985 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)
01986 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](18);
01987 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)
01988 + slip*f[b+base::idx(i-1,j,mzs-1,Nx,Ny)](15);
01989 }
01990 }
01991 if (wc(1)>=min[1] && wc(1)<=max[1]) {
01992 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);
01993 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);
01994 }
01995 else {
01996 if (wc(1)<WLB[1] || wc(1) >WUB[1])
01997 slip = DataType(0.);
01998 else if (wc(1)<min[1])
01999 slip = (wc(1) - WLB[1])/lys;
02000 else if (wc(1)>max[1])
02001 slip = (WUB[1] - wc(1))/lye;
02002 if (j>mys && j<mye) {
02003 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)
02004 + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
02005 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)
02006 + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
02007 }
02008 else if (j==mys) {
02009 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)
02010 + slip*f[b+base::idx(i,j+1,mzs-1,Nx,Ny)](14);
02011 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)
02012 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](11);
02013 }
02014 else if (j==mye) {
02015 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)
02016 + slip*f[b+base::idx(i,j,mzs-1,Nx,Ny)](14);
02017 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)
02018 + slip*f[b+base::idx(i,j-1,mzs-1,Nx,Ny)](11);
02019 }
02020 }
02021 f[b+base::idx(i,j,mzs,Nx,Ny)](6) = f[b+base::idx(i,j,mzs-1,Nx,Ny)](5);
02022 }
02023 }
02024 break;
02025 }
02026 }
02027 }
02028
02029 inline virtual const MacroType NormalDerivative(const MacroType qa,const MacroType qb, const MacroType qc) const
02030 { return DataType(-1.5)*qa + DataType(2.)*qb - DataType(0.5)*qc; }
02031 inline virtual Vector<DataType,3> WaveAmplitudes(const DataType rho, const DataType vn, const DataType drhodn,
02032 const DataType dvndn, const DataType dvt0dn, const DataType dvt1dn) const {
02033 Vector<DataType,3> L;
02034 L(0) = (vn - LatticeSpeedOfSound())*(cs2*drhodn - rho*LatticeSpeedOfSound()*dvndn);
02035 L(1) = vn*dvt0dn;
02036 L(2) = vn*dvt1dn;
02037 return L;
02038 }
02039
02040 virtual void GFMBCStandard(vec_grid_data_type &fvec, const int type, const int& nc, const int* idx,
02041 const MicroType* f, const point_type* xc, const DataType* distance,
02042 const point_type* normal, DataType* aux=0, const int naux=0,
02043 const int scaling=0) const {
02044 if (type==GFMNoSlipWall || type==GFMSlipWall) {
02045 DataType fac = S0;
02046 if (scaling==base::Physical)
02047 fac /= U0;
02048 for (int n=0; n<nc; n++) {
02049 MacroType q=MacroVariables(f[n]);
02050 DataType u=-q(1);
02051 DataType v=-q(2);
02052 DataType w=-q(3);
02053
02054
02055 if (naux>=3) {
02056 u += fac*aux[naux*n];
02057 v += fac*aux[naux*n+1];
02058 w += fac*aux[naux*n+2];
02059 }
02060 if (type==GFMNoSlipWall) {
02061
02062 q(1) += DataType(2.)*u;
02063 q(2) += DataType(2.)*v;
02064 q(3) += DataType(2.)*w;
02065 }
02066 else if (type==GFMSlipWall) {
02067
02068 DataType vl = DataType(2.)*(u*normal[n](0)+v*normal[n](1)+w*normal[n](2));
02069 q(1) += vl*normal[n](0);
02070 q(2) += vl*normal[n](1);
02071 q(3) += vl*normal[n](2);
02072 }
02073 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02074 }
02075 }
02076
02077 else if (type==GFMExtrapolation) {
02078 for (int n=0; n<nc; n++) {
02079 MacroType q=MacroVariables(f[n]);
02080 DataType vl = q(1)*normal[n](0)+q(2)*normal[n](1)+q(3)*normal[n](2);
02081 q(1) = vl*normal[n](0);
02082 q(2) = vl*normal[n](1);
02083 q(3) = vl*normal[n](2);
02084 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02085 }
02086 }
02087
02088 else if (type==GFMWallLaw ) {
02089 DCoords dx = base::GH().worldStep(fvec.stepsize());
02090 DataType fac = S0;
02091 if (scaling==base::Physical)
02092 fac /= U0;
02093
02094
02095 DataType yp_ = 0., yp = 0., yMax = 50.;
02096 DataType yPlus[3];
02097 int lc_low[3], lc_up[3], yPlus_lc[3];
02098
02099
02100
02101
02102
02103 DataType shear_vel = 1., karman = 0.41, y0 = 1., law = 1.;
02104
02105 DataType shearStress;
02106 DataType u, v, w, u_, v_, w_, vel0;
02107
02108 MacroType q, q_;
02109 int kMax = 0;
02110
02111 for (int n=0; n<nc; n++) {
02112
02114 q=MacroVariables(f[n]);
02115 q_=MacroVariables(f[n]);
02116
02117 yPlus[0] = (xc[n](0) - DataType(2.)*dx(0)*normal[n](0));
02118 yPlus[1] = (xc[n](1) - DataType(2.)*dx(1)*normal[n](1));
02119 yPlus[2] = (xc[n](2) - DataType(2.)*dx(2)*normal[n](2));
02120
02121 yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
02122 yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
02123 yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
02124
02125 lc_low[0] = fvec.lower(0); lc_low[1] = fvec.lower(1); lc_low[2] = fvec.lower(2);
02126 lc_up[0] = fvec.upper(0); lc_up[1] = fvec.upper(1); lc_up[2] = fvec.upper(2);
02127
02128
02129
02130 if ( yPlus_lc[0] > lc_low[0] ) {
02131 if ( yPlus_lc[1] > lc_low[1] ) {
02132 if ( yPlus_lc[2] > lc_low[2] ) {
02133
02134 if ( yPlus_lc[0] < lc_up[0] ) {
02135 if ( yPlus_lc[1] < lc_up[1] ) {
02136 if ( yPlus_lc[2] < lc_up[2] ) {
02137 q=MacroVariables( fvec( base::GH().localCoords((const DataType *) yPlus) ) );
02138 }
02139 }
02140 }
02141
02142 }
02143 }
02144 }
02145
02146 u_=-q(1);
02147 v_=-q(2);
02148 w_=-q(3);
02149
02150
02151 q=MacroVariables(f[n]);
02152 u=-q(1);
02153 v=-q(2);
02154 w=-q(3);
02155 vel0 = sqrt(u*u+v*v+w*w);
02156
02157 u_+=q(1);
02158 v_+=q(2);
02159 w_+=q(3);
02160 vel0 = ( u_*normal[n](0) + v_*normal[n](1) + w_*normal[n](2) );
02161
02162 u_ = u_ - vel0 * normal[n](0);
02163 v_ = v_ - vel0 * normal[n](1);
02164 w_ = w_ - vel0 * normal[n](2);
02165 vel0 = sqrt(u_*u_ + v_*v_ + w_*w_)*U0/S0;
02166
02167 shearStress = (DataType(0.5)*vel0/dx(0) * q(0)*R0*GasViscosity() );
02168
02169
02170 if (naux>=3) {
02171 u += fac*aux[naux*n];
02172 v += fac*aux[naux*n+1];
02173 w += fac*aux[naux*n+2];
02174 }
02175 q(1) += DataType(2.)*u;
02176 q(2) += DataType(2.)*v;
02177 q(3) += DataType(2.)*w;
02178 fvec(Coords(base::Dim(),idx+base::Dim()*n)) = Equilibrium(q);
02179
02181 shear_vel = std::sqrt( shearStress/(q(0)*R0) );
02182
02184 y0 = karman/30.;
02185
02186
02187
02188 if ( shear_vel > 1.e-5 ) {
02189 kMax = (int) std::ceil( (yMax*GasViscosity())/(dx(0)*shear_vel) );
02190
02191 for (int k=1; k < kMax; k++) {
02192
02193 yPlus[0] = (xc[n](0) - k*dx(0)*normal[n](0));
02194 yPlus[1] = (xc[n](1) - k*dx(1)*normal[n](1));
02195 yPlus[2] = (xc[n](2) - k*dx(2)*normal[n](2));
02196
02197 yPlus_lc[0] = base::GH().localCoords((const DataType *) yPlus )(0);
02198 yPlus_lc[1] = base::GH().localCoords((const DataType *) yPlus )(1);
02199 yPlus_lc[2] = base::GH().localCoords((const DataType *) yPlus )(2);
02200
02201 if ( yPlus_lc[0] > lc_low[0] ) {
02202 if ( yPlus_lc[1] > lc_low[1] ) {
02203 if ( yPlus_lc[2] > lc_low[2] ) {
02204
02205 if ( yPlus_lc[0] < lc_up[0] ) {
02206 if ( yPlus_lc[1] < lc_up[1] ) {
02207 if ( yPlus_lc[2] < lc_up[2] ) {
02208 yp = k*dx(0);
02209 yp_ = yp*shear_vel/GasViscosity();
02210
02213 if ( yp_ <= 5.0 ) {
02214 law = yp_;
02215 }
02217 else if (yp_ > 5 && yp_ < 30.0 ) {
02218
02219 law = std::min( yp_ , shear_vel/karman*std::log( yp/y0 ) );
02220
02221 }
02223 else if (yp_ >= 30.0 ) {
02224 law = shear_vel/karman*std::log( yp/y0 );
02225 }
02226
02227 q_=MacroVariables(f[n]);
02228 u = law*(-q_(1) + fac*aux[naux*n]);
02229 v = law*(-q_(2) + fac*aux[naux*n+1]);
02230 w = law*(-q_(3) + fac*aux[naux*n+2]);
02231
02232 q_(1) += DataType(2.)*u;
02233 q_(2) += DataType(2.)*v;
02234 q_(3) += DataType(2.)*w;
02235 fvec( base::GH().localCoords((const DataType *) yPlus) ) = Equilibrium(q_);
02236 }
02237 }
02238 }
02239
02240
02241 }
02242 }
02243 }
02244
02245 }
02246
02247 }
02248
02249 }
02250
02251 }
02252
02253 }
02254
02255 virtual void Output(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
02256 const int skip_ghosts=1) const {
02257 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
02258 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
02259 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
02260 MicroType *f = (MicroType *)fvec.databuffer();
02261 DataType *work = (DataType *)workvec.databuffer();
02262 DCoords dx = base::GH().worldStep(fvec.stepsize());
02263 DataType dt = dx(0)*S0/U0;
02264
02265 if ((cnt>=1 && cnt<=10) || (cnt>=19 && cnt<=24) || (cnt>=28 && cnt<=33)) {
02266 for (register int k=mzs; k<=mze; k++)
02267 for (register int j=mys; j<=mye; j++)
02268 for (register int i=mxs; i<=mxe; i++) {
02269 MacroType q=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
02270 switch(cnt) {
02271 case 1:
02272 if (method[0]==0) work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0+rhop;
02273 else
02274 work[base::idx(i,j,k,Nx,Ny)]=q(0)*R0;
02275 break;
02276 case 2:
02277 work[base::idx(i,j,k,Nx,Ny)]=q(1)*U0/S0;
02278 break;
02279 case 3:
02280 work[base::idx(i,j,k,Nx,Ny)]=q(2)*U0/S0;
02281 break;
02282 case 4:
02283 work[base::idx(i,j,k,Nx,Ny)]=q(3)*U0/S0;
02284 break;
02285 case 6:
02286 work[base::idx(i,j,k,Nx,Ny)]=TempEquation(U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure());
02287 break;
02288 case 7:
02289 work[base::idx(i,j,k,Nx,Ny)]=U0/S0*U0/S0*R0*LatticeBasePressure(q(0))+BasePressure();
02290 break;
02291 case 9: {
02292 MicroType feq = Equilibrium(q);
02293 work[base::idx(i,j,k,Nx,Ny)]=GasViscosity(Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt),
02294 GasSpeedofSound(),dt);
02295 break;
02296 }
02297 case 10:
02298 work[base::idx(i,j,k,Nx,Ny)]=std::sqrt(q(1)*q(1)+q(2)*q(2)+q(3)*q(3))*U0/S0;
02299 break;
02300 case 19:
02301 case 20:
02302 case 21:
02303 case 22:
02304 case 23:
02305 case 24: {
02306 MicroType feq = Equilibrium(q);
02307 DataType omega;
02308 if (turbulence == laminar)
02309 omega = Omega(dt);
02310 else if (turbulence == LES_Smagorinsky)
02311 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02312 else if (turbulence == LES_dynamic)
02313 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02314 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);
02315 break;
02316 }
02317 case 28:
02318 case 29:
02319 case 30:
02320 case 31:
02321 case 32:
02322 case 33: {
02323 DataType omega;
02324 if (turbulence == laminar)
02325 omega = Omega(dt);
02326 else if (turbulence == LES_Smagorinsky) {
02327 MicroType feq = Equilibrium(q);
02328 omega = Omega_LES_Smagorinsky(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02329 }
02330 else if (turbulence == LES_dynamic) {
02331 MicroType feq = Equilibrium(q);
02332 omega = Omega_LES_dynamic(f[base::idx(i,j,k,Nx,Ny)],feq,q,dt);
02333 }
02334 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);
02335 if (cnt==28 || cnt==30 || cnt==33) work[base::idx(i,j,k,Nx,Ny)]-=BasePressure();
02336 break;
02337 }
02338 }
02339 }
02340 }
02341 else if ((cnt>=11 && cnt<=14) || cnt==18 || (cnt>=40 && cnt<=45)) {
02342 macro_grid_data_type qgrid(fvec.bbox());
02343 MacroType *q = (MacroType *)qgrid.databuffer();
02344 for (register int k=mzs; k<=mze; k++)
02345 for (register int j=mys; j<=mye; j++)
02346 for (register int i=mxs; i<=mxe; i++) {
02347 q[base::idx(i,j,k,Nx,Ny)]=MacroVariables(f[base::idx(i,j,k,Nx,Ny)]);
02348 q[base::idx(i,j,k,Nx,Ny)](0)*=R0;
02349 q[base::idx(i,j,k,Nx,Ny)](1)*=U0/S0;
02350 q[base::idx(i,j,k,Nx,Ny)](2)*=U0/S0;
02351 q[base::idx(i,j,k,Nx,Ny)](3)*=U0/S0;
02352 }
02353
02354 if (skip_ghosts==0) {
02355 switch(cnt) {
02356 case 40:
02357 mxs+=1; mxe-=1;
02358 break;
02359 case 42:
02360 mys+=1; mye-=1;
02361 break;
02362 case 45:
02363 mzs+=1; mze-=1;
02364 break;
02365 case 11:
02366 case 44:
02367 mys+=1; mye-=1; mzs+=1; mze-=1;
02368 break;
02369 case 12:
02370 case 43:
02371 mxs+=1; mxe-=1; mzs+=1; mze-=1;
02372 break;
02373 case 13:
02374 case 41:
02375 mxs+=1; mxe-=1; mys+=1; mye-=1;
02376 break;
02377 case 14:
02378 case 18:
02379 mxs+=1; mxe-=1; mys+=1; mye-=1; mzs+=1; mze-=1;
02380 break;
02381 }
02382 }
02383
02384 DataType dxux, dxuy, dxuz, dyux, dyuy, dyuz, dzux, dzuy, dzuz;
02385 for (register int k=mzs; k<=mze; k++)
02386 for (register int j=mys; j<=mye; j++)
02387 for (register int i=mxs; i<=mxe; i++) {
02388 if (cnt==18 || cnt==40)
02389 dxux=(q[base::idx(i+1,j,k,Nx,Ny)](1)-q[base::idx(i-1,j,k,Nx,Ny)](1))/(2.*dx(0));
02390 if (cnt==18 || cnt==42)
02391 dyuy=(q[base::idx(i,j+1,k,Nx,Ny)](2)-q[base::idx(i,j-1,k,Nx,Ny)](2))/(2.*dx(1));
02392 if (cnt==18 || cnt==45)
02393 dzuz=(q[base::idx(i,j,k+1,Nx,Ny)](3)-q[base::idx(i,j,k-1,Nx,Ny)](3))/(2.*dx(2));
02394 if (cnt==11 || cnt==14 || cnt==18 || cnt==44) {
02395 dyuz=(q[base::idx(i,j+1,k,Nx,Ny)](3)-q[base::idx(i,j-1,k,Nx,Ny)](3))/(2.*dx(1));
02396 dzuy=(q[base::idx(i,j,k+1,Nx,Ny)](2)-q[base::idx(i,j,k-1,Nx,Ny)](2))/(2.*dx(2));
02397 }
02398 if (cnt==12 || cnt==14 || cnt==18 || cnt==43) {
02399 dxuz=(q[base::idx(i+1,j,k,Nx,Ny)](3)-q[base::idx(i-1,j,k,Nx,Ny)](3))/(2.*dx(0));
02400 dzux=(q[base::idx(i,j,k+1,Nx,Ny)](1)-q[base::idx(i,j,k-1,Nx,Ny)](1))/(2.*dx(2));
02401 }
02402 if (cnt==13 || cnt==14 || cnt==18 || cnt==41) {
02403 dxuy=(q[base::idx(i+1,j,k,Nx,Ny)](2)-q[base::idx(i-1,j,k,Nx,Ny)](2))/(2.*dx(0));
02404 dyux=(q[base::idx(i,j+1,k,Nx,Ny)](1)-q[base::idx(i,j-1,k,Nx,Ny)](1))/(2.*dx(1));
02405 }
02406 switch(cnt) {
02407 case 11:
02408 work[base::idx(i,j,k,Nx,Ny)]=dyuz-dzuy;
02409 break;
02410 case 12:
02411 work[base::idx(i,j,k,Nx,Ny)]=dzux-dxuz;
02412 break;
02413 case 13:
02414 work[base::idx(i,j,k,Nx,Ny)]=dxuy-dyux;
02415 break;
02416 case 14: {
02417 DataType ox = dyuz-dzuy, oy = dzux-dxuz, oz = dxuy-dyux;
02418 work[base::idx(i,j,k,Nx,Ny)] = std::sqrt(ox*ox+oy*oy+oz*oz);
02419 break;
02420 }
02421 case 18:
02422 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);
02423 break;
02424 case 40:
02425 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dxux;
02426 break;
02427 case 41:
02428 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuy+dyux);
02429 break;
02430 case 42:
02431 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dyuy;
02432 break;
02433 case 43:
02434 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dxuz+dzux);
02435 break;
02436 case 44:
02437 work[base::idx(i,j,k,Nx,Ny)]=q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*(dyuz+dzuy);
02438 break;
02439 case 45:
02440 work[base::idx(i,j,k,Nx,Ny)]=2.*q[base::idx(i,j,k,Nx,Ny)](0)*GasViscosity()*dzuz;
02441 break;
02442 }
02443 }
02444 }
02445 }
02446
02447 virtual void Input(vec_grid_data_type& fvec, grid_data_type& workvec, const int cnt,
02448 const int skip_ghosts=1) const {
02449 int Nx = fvec.extents(0), Ny = fvec.extents(1), Nz = fvec.extents(2);
02450 int mxs = base::NGhosts()*skip_ghosts, mys = base::NGhosts()*skip_ghosts, mzs = base::NGhosts()*skip_ghosts;
02451 int mxe = Nx-base::NGhosts()*skip_ghosts-1, mye = Ny-base::NGhosts()*skip_ghosts-1, mze = Nz-base::NGhosts()*skip_ghosts-1;
02452 MicroType *f = (MicroType *)fvec.databuffer();
02453 DataType *work = (DataType *)workvec.databuffer();
02454 DCoords dx = base::GH().worldStep(fvec.stepsize());
02455
02456 for (register int k=mzs; k<=mze; k++)
02457 for (register int j=mys; j<=mye; j++)
02458 for (register int i=mxs; i<=mxe; i++) {
02459 switch(cnt) {
02460 case 1:
02461 f[base::idx(i,j,k,Nx,Ny)](0)=work[base::idx(i,j,k,Nx,Ny)]/R0;
02462 break;
02463 case 2:
02464 f[base::idx(i,j,k,Nx,Ny)](1)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02465 break;
02466 case 3:
02467 f[base::idx(i,j,k,Nx,Ny)](2)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02468 break;
02469 case 4:
02470 f[base::idx(i,j,k,Nx,Ny)](3)=work[base::idx(i,j,k,Nx,Ny)]/(U0/S0);
02471 f[base::idx(i,j,k,Nx,Ny)] = Equilibrium((const MacroType &)f[base::idx(i,j,k,Nx,Ny)]);
02472 break;
02473 }
02474 }
02475 }
02476
02477 virtual int Check(vec_grid_data_type& fvec, const BBox &bb, const double& time,
02478 const int verbose) const {
02479 int Nx = fvec.extents(0), Ny = fvec.extents(1);
02480 int b = fvec.bottom();
02481 assert (fvec.stepsize(0)==bb.stepsize(0) && fvec.stepsize(1)==bb.stepsize(1));
02482 BBox bbmax = grow(fvec.bbox(),-base::NGhosts());
02483 int mxs = std::max(bb.lower(0)/bb.stepsize(0),bbmax.lower(0)/bbmax.stepsize(0));
02484 int mys = std::max(bb.lower(1)/bb.stepsize(1),bbmax.lower(1)/bbmax.stepsize(1));
02485 int mzs = std::max(bb.lower(2)/bb.stepsize(2),bbmax.lower(2)/bbmax.stepsize(2));
02486 int mxe = std::min(bb.upper(0)/bb.stepsize(0),bbmax.upper(0)/bbmax.stepsize(0));
02487 int mye = std::min(bb.upper(1)/bb.stepsize(1),bbmax.upper(1)/bbmax.stepsize(1));
02488 int mze = std::min(bb.upper(2)/bb.stepsize(2),bbmax.upper(2)/bbmax.stepsize(2));
02489 MicroType *f = (MicroType *)fvec.databuffer();
02490
02491
02492 int result = 1;
02493 for (register int k=mzs; k<=mze; k++)
02494 for (register int j=mys; j<=mye; j++)
02495 for (register int i=mxs; i<=mxe; i++)
02496 for (register int l=0; l<base::NMicroVar(); l++)
02497 if (!(f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l)>-1.e37)) {
02498 result = 0;
02499 if (verbose) {
02500 DCoords lbcorner = base::GH().worldCoords(fvec.lower(), fvec.stepsize());
02501 DCoords dx = base::GH().worldStep(fvec.stepsize());
02502
02503 std::cerr << "D3Q19-Check(): f(" << i-mdx[l] << "," << j-mdy[l] << "," << k-mdz[l] << ")(" << l << ")="
02504 << f[b+base::idx(i-mdx[l],j-mdy[l],k-mdz[l],Nx,Ny)](l) << " " << fvec.bbox()
02505 << " Coords: (" << (i+0.5)*dx(0)+lbcorner(0) << "," << (j+0.5)*dx(1)+lbcorner(1) << ","
02506 << (k+0.5)*dx(2)+lbcorner(2) << ")" << std::endl;
02507 }
02508 }
02509 return result;
02510 }
02511
02512 inline const TensorType Stress(const MicroType &f, const MacroType &q, const DataType om) const {
02513 TensorType P;
02514 if (method[1]==0) {
02515 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);
02516 P(1) = q(0)*q(1)*q(2)-(f(7)+f(8)-f(9)-f(10));
02517 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);
02518 P(3) = q(0)*q(1)*q(3)-(f(15)+f(16)-f(17)-f(18));
02519 P(4) = q(0)*q(2)*q(3)-(f(11)+f(12)-f(13)-f(14));
02520 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);
02521 P *= -(DataType(1.0)-DataType(0.5)*om);
02522 }
02523 else {
02524 MicroType feq = Equilibrium(q);
02525 P = DeviatoricStress(f,feq,om);
02526 }
02527 P(0) -= LatticeBasePressure(q(0));
02528 P(2) -= LatticeBasePressure(q(0));
02529 P(5) -= LatticeBasePressure(q(0));
02530 return P;
02531 }
02532
02533 inline const TensorType DeviatoricStress(const MicroType &f, const MicroType &feq, const DataType om) const {
02534 TensorType Sigma;
02535 if (method[2]==0) {
02536 MicroType fneq = feq - f;
02537 Sigma(0) = fneq(1)+fneq(2)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
02538 Sigma(1) = fneq(7)+fneq(8)-fneq(9)-fneq(10);
02539 Sigma(2) = fneq(3)+fneq(4)+fneq(7)+fneq(8)+fneq(9)+fneq(10)+fneq(11)+fneq(12)+fneq(13)+fneq(14);
02540 Sigma(3) = fneq(15)+fneq(16)-fneq(17)-fneq(18);
02541 Sigma(4) = fneq(11)+fneq(12)-fneq(13)-fneq(14);
02542 Sigma(5) = fneq(5)+fneq(6)+fneq(11)+fneq(12)+fneq(13)+fneq(14)+fneq(15)+fneq(16)+fneq(17)+fneq(18);
02543 Sigma *= -(DataType(1.0)-DataType(0.5)*om);
02544 }
02545 else {
02546 MacroType q = MacroVariables(f);
02547 Sigma = Stress(f,q,om);
02548 Sigma(0) += LatticeBasePressure(q(0));
02549 Sigma(2) += LatticeBasePressure(q(0));
02550 Sigma(5) += LatticeBasePressure(q(0));
02551 }
02552 return Sigma;
02553 }
02554
02555 virtual int NMethodOrder() const { return 2; }
02556
02557 inline const DataType& L0() const { return base::LengthScale(); }
02558 inline const DataType& T0() const { return base::TimeScale(); }
02559 inline void SetDensityScale(const DataType r0) { R0 = r0; }
02560 inline void SetVelocityScale(const DataType u0) { U0 = u0; base::T0 = L0()/U0; base::T0 *= S0; }
02561 inline void SetSpeedUp(const DataType s0) { S0 = s0; }
02562 inline virtual void SetTimeScale(const DataType t0) { base::T0 = t0; U0 = L0()/T0(); base::T0 *= S0; }
02563
02564 inline const DataType& DensityScale() const { return R0; }
02565 inline const DataType VelocityScale() const { return U0/S0; }
02566
02567 inline const DataType& SpeedUp() const { return S0; }
02568
02569
02570 inline DataType LatticeViscosity(const DataType omega) const { return ((DataType(1.)/omega-DataType(0.5))*cs2); }
02571 inline DataType LatticeSpeedOfSound() const { return (std::sqrt(cs2)); }
02572 inline virtual DataType LatticeBasePressure(const DataType rho) const {
02573 return (method[0]==0 ? cs2*rho/R0 : cs2*(rho-rhop/R0));
02574 }
02575
02576
02577 inline void SetGas(DataType rho, DataType nu, DataType cs) {
02578 rhop = rho; nup = nu; csp = cs; cs2p = cs*cs;
02579 SetTimeScale(LatticeSpeedOfSound()/GasSpeedofSound()*L0());
02580 }
02581 inline virtual const DataType Omega(const DataType dt) const { return (cs2p*dt/S0/(nup*S0+cs2p*dt/S0*DataType(0.5))); }
02582 inline const DataType Omega_LES_Smagorinsky(const MicroType &f, const MicroType &feq, const MacroType& q,
02583 const DataType dt) const {
02584 DataType M2 = 0.0;
02585
02586 M2 += (f(7) - feq(7))*(f(7) - feq(7));
02587 M2 += (f(8) - feq(8))*(f(8) - feq(8));
02588 M2 += (f(9) - feq(9))*(f(9) - feq(9));
02589 M2 += (f(10) - feq(10))*(f(10) - feq(10));
02590
02591 M2 += (f(11) - feq(11))*(f(11) - feq(11));
02592 M2 += (f(12) - feq(12))*(f(12) - feq(12));
02593 M2 += (f(13) - feq(13))*(f(13) - feq(13));
02594 M2 += (f(14) - feq(14))*(f(14) - feq(14));
02595
02596 M2 += (f(15) - feq(15))*(f(15) - feq(15));
02597 M2 += (f(16) - feq(16))*(f(16) - feq(16));
02598 M2 += (f(17) - feq(17))*(f(17) - feq(17));
02599 M2 += (f(18) - feq(18))*(f(18) - feq(18));
02600 M2 = std::sqrt(M2);
02601
02602 DataType tau = DataType(1.0)/Omega(dt);
02603 DataType om_turb = (DataType(2.0)/( tau
02604 + sqrt(tau*tau + (DataType(18.)*Cs_Smagorinsky*Cs_Smagorinsky*M2)/(R0*q(0)*cs2/S0/S0)) ));
02605
02606 return ( om_turb );
02607 }
02608 inline const DataType Omega_LES_dynamic(const MicroType &f, const MicroType &feq, const MacroType& q,
02609 const DataType dt) const {
02610 DataType M2 = 0.0;
02611
02612 M2 += (f(7) - feq(7))*(f(7) - feq(7));
02613 M2 += (f(8) - feq(8))*(f(8) - feq(8));
02614 M2 += (f(9) - feq(9))*(f(9) - feq(9));
02615 M2 += (f(10) - feq(10))*(f(10) - feq(10));
02616
02617 M2 += (f(11) - feq(11))*(f(11) - feq(11));
02618 M2 += (f(12) - feq(12))*(f(12) - feq(12));
02619 M2 += (f(13) - feq(13))*(f(13) - feq(13));
02620 M2 += (f(14) - feq(14))*(f(14) - feq(14));
02621
02622 M2 += (f(15) - feq(15))*(f(15) - feq(15));
02623 M2 += (f(16) - feq(16))*(f(16) - feq(16));
02624 M2 += (f(17) - feq(17))*(f(17) - feq(17));
02625 M2 += (f(18) - feq(18))*(f(18) - feq(18));
02626 M2 = std::sqrt(M2);
02627
02628 DataType tau = DataType(1.0)/Omega(dt);
02629 DataType vel = sqrt( q(1)*q(1)+q(2)*q(2)+q(3)*q(3) );
02630 DataType Cs_SmagDYN = DataType(0.5)*L0()*L0()*( (vel*vel/csp*(M2-vel)) / ((M2-vel)*(M2-vel)) );
02631 DataType om_turb = (DataType(2.0)/( tau
02632 + sqrt(tau*tau + (DataType(18.)*Cs_SmagDYN*M2)/(R0*q(0)*cs2/S0/S0)) ));
02633
02634 return ( om_turb );
02635 }
02636 inline const int TurbulenceType() const { return turbulence; }
02637 inline const DataType& SmagorinskyConstant() { return Cs_Smagorinsky; }
02638 inline const DataType& GasDensity() const { return rhop; }
02639 inline const DataType& GasViscosity() const { return nup; }
02640 inline const DataType& GasSpeedofSound() const { return csp; }
02641 inline const DataType GasViscosity(const DataType omega, const DataType cs, const DataType dt) const
02642 { return ((DataType(1.)/omega-DataType(0.5))*cs*cs*dt/S0/S0); }
02643
02644 inline void SetGasProp(DataType g, DataType W, DataType R) { gp = g; Wp = W; R = Rp; }
02645 inline virtual const DataType BasePressure() const { return (gp > DataType(0.) ? rhop*cs2p/gp : DataType(0.)); }
02646 inline virtual const DataType TempEquation(const DataType p) const { return (p*Wp/(rhop*Rp)); }
02647
02648 protected:
02649 DataType cs2, cs22, cssq;
02650 DataType R0, U0, S0, rhop, csp, cs2p, nup, gp, Wp, Rp;
02651 DataType Cs_Smagorinsky, turbulence;
02652 int method[3], mdx[19], mdy[19], mdz[19];
02653 };
02654
02655
02656 #endif