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