|
Blender
V2.59
|
00001 00004 /****************************************************************************** 00005 * 00006 * El'Beem - the visual lattice boltzmann freesurface simulator 00007 * All code distributed as part of El'Beem is covered by the version 2 of the 00008 * GNU General Public License. See the file COPYING for details. 00009 * Copyright 2003-2006 Nils Thuerey 00010 * 00011 * Combined 2D/3D Lattice Boltzmann relaxation macros 00012 * 00013 *****************************************************************************/ 00014 00015 #if FSGR_STRICT_DEBUG==1 00016 #define CAUSE_PANIC { this->mPanic=1; /* *((int*)(0x0)) = 1; crash*/ } 00017 #else // FSGR_STRICT_DEBUG==1 00018 #define CAUSE_PANIC { this->mPanic=1; } /*set flag*/ 00019 #endif // FSGR_STRICT_DEBUG==1 00020 00021 /****************************************************************************** 00022 * normal relaxation 00023 *****************************************************************************/ 00024 00025 // standard arrays 00026 #define CSRC_C RAC(ccel , dC ) 00027 #define CSRC_E RAC(ccel + (-1) *(dTotalNum), dE ) 00028 #define CSRC_W RAC(ccel + (+1) *(dTotalNum), dW ) 00029 #define CSRC_N RAC(ccel + (-mLevel[lev].lOffsx) *(dTotalNum), dN ) 00030 #define CSRC_S RAC(ccel + (+mLevel[lev].lOffsx) *(dTotalNum), dS ) 00031 #define CSRC_NE RAC(ccel + (-mLevel[lev].lOffsx-1) *(dTotalNum), dNE) 00032 #define CSRC_NW RAC(ccel + (-mLevel[lev].lOffsx+1) *(dTotalNum), dNW) 00033 #define CSRC_SE RAC(ccel + (+mLevel[lev].lOffsx-1) *(dTotalNum), dSE) 00034 #define CSRC_SW RAC(ccel + (+mLevel[lev].lOffsx+1) *(dTotalNum), dSW) 00035 #define CSRC_T RAC(ccel + (-mLevel[lev].lOffsy) *(dTotalNum), dT ) 00036 #define CSRC_B RAC(ccel + (+mLevel[lev].lOffsy) *(dTotalNum), dB ) 00037 #define CSRC_ET RAC(ccel + (-mLevel[lev].lOffsy-1) *(dTotalNum), dET) 00038 #define CSRC_EB RAC(ccel + (+mLevel[lev].lOffsy-1) *(dTotalNum), dEB) 00039 #define CSRC_WT RAC(ccel + (-mLevel[lev].lOffsy+1) *(dTotalNum), dWT) 00040 #define CSRC_WB RAC(ccel + (+mLevel[lev].lOffsy+1) *(dTotalNum), dWB) 00041 #define CSRC_NT RAC(ccel + (-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNT) 00042 #define CSRC_NB RAC(ccel + (+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNB) 00043 #define CSRC_ST RAC(ccel + (-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dST) 00044 #define CSRC_SB RAC(ccel + (+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dSB) 00045 00046 #define XSRC_C(x) RAC(ccel + (x) *dTotalNum, dC ) 00047 #define XSRC_E(x) RAC(ccel + ((x)-1) *dTotalNum, dE ) 00048 #define XSRC_W(x) RAC(ccel + ((x)+1) *dTotalNum, dW ) 00049 #define XSRC_N(x) RAC(ccel + ((x)-mLevel[lev].lOffsx) *dTotalNum, dN ) 00050 #define XSRC_S(x) RAC(ccel + ((x)+mLevel[lev].lOffsx) *dTotalNum, dS ) 00051 #define XSRC_NE(x) RAC(ccel + ((x)-mLevel[lev].lOffsx-1) *dTotalNum, dNE) 00052 #define XSRC_NW(x) RAC(ccel + ((x)-mLevel[lev].lOffsx+1) *dTotalNum, dNW) 00053 #define XSRC_SE(x) RAC(ccel + ((x)+mLevel[lev].lOffsx-1) *dTotalNum, dSE) 00054 #define XSRC_SW(x) RAC(ccel + ((x)+mLevel[lev].lOffsx+1) *dTotalNum, dSW) 00055 #define XSRC_T(x) RAC(ccel + ((x)-mLevel[lev].lOffsy) *dTotalNum, dT ) 00056 #define XSRC_B(x) RAC(ccel + ((x)+mLevel[lev].lOffsy) *dTotalNum, dB ) 00057 #define XSRC_ET(x) RAC(ccel + ((x)-mLevel[lev].lOffsy-1) *dTotalNum, dET) 00058 #define XSRC_EB(x) RAC(ccel + ((x)+mLevel[lev].lOffsy-1) *dTotalNum, dEB) 00059 #define XSRC_WT(x) RAC(ccel + ((x)-mLevel[lev].lOffsy+1) *dTotalNum, dWT) 00060 #define XSRC_WB(x) RAC(ccel + ((x)+mLevel[lev].lOffsy+1) *dTotalNum, dWB) 00061 #define XSRC_NT(x) RAC(ccel + ((x)-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNT) 00062 #define XSRC_NB(x) RAC(ccel + ((x)+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNB) 00063 #define XSRC_ST(x) RAC(ccel + ((x)-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dST) 00064 #define XSRC_SB(x) RAC(ccel + ((x)+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dSB) 00065 00066 00067 00068 #define OMEGA(l) mLevel[(l)].omega 00069 00070 #define EQC ( DFL1*(rho - usqr)) 00071 #define EQN ( DFL2*(rho + uy*(4.5*uy + 3.0) - usqr)) 00072 #define EQS ( DFL2*(rho + uy*(4.5*uy - 3.0) - usqr)) 00073 #define EQE ( DFL2*(rho + ux*(4.5*ux + 3.0) - usqr)) 00074 #define EQW ( DFL2*(rho + ux*(4.5*ux - 3.0) - usqr)) 00075 #define EQT ( DFL2*(rho + uz*(4.5*uz + 3.0) - usqr)) 00076 #define EQB ( DFL2*(rho + uz*(4.5*uz - 3.0) - usqr)) 00077 00078 #define EQNE ( DFL3*(rho + (+ux+uy)*(4.5*(+ux+uy) + 3.0) - usqr)) 00079 #define EQNW ( DFL3*(rho + (-ux+uy)*(4.5*(-ux+uy) + 3.0) - usqr)) 00080 #define EQSE ( DFL3*(rho + (+ux-uy)*(4.5*(+ux-uy) + 3.0) - usqr)) 00081 #define EQSW ( DFL3*(rho + (-ux-uy)*(4.5*(-ux-uy) + 3.0) - usqr)) 00082 #define EQNT ( DFL3*(rho + (+uy+uz)*(4.5*(+uy+uz) + 3.0) - usqr)) 00083 #define EQNB ( DFL3*(rho + (+uy-uz)*(4.5*(+uy-uz) + 3.0) - usqr)) 00084 #define EQST ( DFL3*(rho + (-uy+uz)*(4.5*(-uy+uz) + 3.0) - usqr)) 00085 #define EQSB ( DFL3*(rho + (-uy-uz)*(4.5*(-uy-uz) + 3.0) - usqr)) 00086 #define EQET ( DFL3*(rho + (+ux+uz)*(4.5*(+ux+uz) + 3.0) - usqr)) 00087 #define EQEB ( DFL3*(rho + (+ux-uz)*(4.5*(+ux-uz) + 3.0) - usqr)) 00088 #define EQWT ( DFL3*(rho + (-ux+uz)*(4.5*(-ux+uz) + 3.0) - usqr)) 00089 #define EQWB ( DFL3*(rho + (-ux-uz)*(4.5*(-ux-uz) + 3.0) - usqr)) 00090 00091 00092 // this is a bit ugly, but necessary for the CSRC_ access... 00093 #define MSRC_C m[dC ] 00094 #define MSRC_N m[dN ] 00095 #define MSRC_S m[dS ] 00096 #define MSRC_E m[dE ] 00097 #define MSRC_W m[dW ] 00098 #define MSRC_T m[dT ] 00099 #define MSRC_B m[dB ] 00100 #define MSRC_NE m[dNE] 00101 #define MSRC_NW m[dNW] 00102 #define MSRC_SE m[dSE] 00103 #define MSRC_SW m[dSW] 00104 #define MSRC_NT m[dNT] 00105 #define MSRC_NB m[dNB] 00106 #define MSRC_ST m[dST] 00107 #define MSRC_SB m[dSB] 00108 #define MSRC_ET m[dET] 00109 #define MSRC_EB m[dEB] 00110 #define MSRC_WT m[dWT] 00111 #define MSRC_WB m[dWB] 00112 00113 // this is a bit ugly, but necessary for the ccel local access... 00114 #define CCEL_C RAC(ccel, dC ) 00115 #define CCEL_N RAC(ccel, dN ) 00116 #define CCEL_S RAC(ccel, dS ) 00117 #define CCEL_E RAC(ccel, dE ) 00118 #define CCEL_W RAC(ccel, dW ) 00119 #define CCEL_T RAC(ccel, dT ) 00120 #define CCEL_B RAC(ccel, dB ) 00121 #define CCEL_NE RAC(ccel, dNE) 00122 #define CCEL_NW RAC(ccel, dNW) 00123 #define CCEL_SE RAC(ccel, dSE) 00124 #define CCEL_SW RAC(ccel, dSW) 00125 #define CCEL_NT RAC(ccel, dNT) 00126 #define CCEL_NB RAC(ccel, dNB) 00127 #define CCEL_ST RAC(ccel, dST) 00128 #define CCEL_SB RAC(ccel, dSB) 00129 #define CCEL_ET RAC(ccel, dET) 00130 #define CCEL_EB RAC(ccel, dEB) 00131 #define CCEL_WT RAC(ccel, dWT) 00132 #define CCEL_WB RAC(ccel, dWB) 00133 // for coarse to fine interpol access 00134 #define CCELG_C(f) (RAC(ccel, dC )*mGaussw[(f)]) 00135 #define CCELG_N(f) (RAC(ccel, dN )*mGaussw[(f)]) 00136 #define CCELG_S(f) (RAC(ccel, dS )*mGaussw[(f)]) 00137 #define CCELG_E(f) (RAC(ccel, dE )*mGaussw[(f)]) 00138 #define CCELG_W(f) (RAC(ccel, dW )*mGaussw[(f)]) 00139 #define CCELG_T(f) (RAC(ccel, dT )*mGaussw[(f)]) 00140 #define CCELG_B(f) (RAC(ccel, dB )*mGaussw[(f)]) 00141 #define CCELG_NE(f) (RAC(ccel, dNE)*mGaussw[(f)]) 00142 #define CCELG_NW(f) (RAC(ccel, dNW)*mGaussw[(f)]) 00143 #define CCELG_SE(f) (RAC(ccel, dSE)*mGaussw[(f)]) 00144 #define CCELG_SW(f) (RAC(ccel, dSW)*mGaussw[(f)]) 00145 #define CCELG_NT(f) (RAC(ccel, dNT)*mGaussw[(f)]) 00146 #define CCELG_NB(f) (RAC(ccel, dNB)*mGaussw[(f)]) 00147 #define CCELG_ST(f) (RAC(ccel, dST)*mGaussw[(f)]) 00148 #define CCELG_SB(f) (RAC(ccel, dSB)*mGaussw[(f)]) 00149 #define CCELG_ET(f) (RAC(ccel, dET)*mGaussw[(f)]) 00150 #define CCELG_EB(f) (RAC(ccel, dEB)*mGaussw[(f)]) 00151 #define CCELG_WT(f) (RAC(ccel, dWT)*mGaussw[(f)]) 00152 #define CCELG_WB(f) (RAC(ccel, dWB)*mGaussw[(f)]) 00153 00154 00155 #if PARALLEL==1 00156 #define CSMOMEGA_STATS(dlev, domega) 00157 #else // PARALLEL==1 00158 #if FSGR_OMEGA_DEBUG==1 00159 #define CSMOMEGA_STATS(dlev, domega) \ 00160 mLevel[dlev].avgOmega += domega; mLevel[dlev].avgOmegaCnt+=1.0; 00161 #else // FSGR_OMEGA_DEBUG==1 00162 #define CSMOMEGA_STATS(dlev, domega) 00163 #endif // FSGR_OMEGA_DEBUG==1 00164 #endif // PARALLEL==1 00165 00166 00167 // used for main loops and grav init 00168 // source set 00169 #define SRCS(l) mLevel[(l)].setCurr 00170 // target set 00171 #define TSET(l) mLevel[(l)].setOther 00172 00173 // handle mov. obj 00174 #if FSGR_STRICT_DEBUG==1 00175 00176 #define LBMDS_ADDMOV(linv,l) \ 00177 \ 00178 if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \ 00179 \ 00180 LbmFloat dte=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)-(mSimulationTime+this->mpParam->getTimestep()); \ 00181 if( ABS(dte)< 1e-15 ) { \ 00182 m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \ 00183 } else { \ 00184 const int sdx = i+this->dfVecX[linv], sdy = j+this->dfVecY[linv], sdz = k+this->dfVecZ[linv]; \ 00185 \ 00186 errMsg("INVALID_MOV_OBJ_TIME"," at "<<PRINT_IJK<<" from l"<<l<<" "<<PRINT_VEC(sdx,sdy,sdz)<<" t="<<(mSimulationTime+this->mpParam->getTimestep())<<" ct="<<QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)<<" dte="<<dte); \ 00187 debugMarkCell(lev,sdx,sdy,sdz); \ 00188 } \ 00189 } \ 00190 00191 00192 00193 #else // FSGR_STRICT_DEBUG==1 00194 00195 #define LBMDS_ADDMOV(linv,l) \ 00196 \ 00197 \ 00198 if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \ 00199 \ 00200 m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \ 00201 } \ 00202 00203 00204 00205 #endif // !FSGR_STRICT_DEBUG==1 00206 00207 // treatment of freeslip reflection 00208 // used both for OPT and nonOPT 00209 #define DEFAULT_STREAM_FREESLIP(l,invl,mnbf) \ 00210 \ 00211 int nb1 = 0, nb2 = 0; \ 00212 LbmFloat newval = 0.0; \ 00213 const int dx = this->dfVecX[invl], dy = this->dfVecY[invl], dz = this->dfVecZ[invl]; \ 00214 \ 00215 \ 00216 \ 00217 const LbmFloat movadd = ( \ 00218 ((nbflag[invl]&CFBndMoving)&&(!(nbflag[l]&CFBnd))) ? \ 00219 (QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l)) : 0.); \ 00220 \ 00221 if(dz==0) { \ 00222 nb1 = !(RFLAG(lev, i, j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \ 00223 nb2 = !(RFLAG(lev, i+dx,j, k, SRCS(lev))&(CFFluid|CFInter)); \ 00224 if((nb1)&&(!nb2)) { \ 00225 \ 00226 newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \ 00227 } else \ 00228 if((!nb1)&&(nb2)) { \ 00229 \ 00230 newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \ 00231 } else { \ 00232 \ 00233 newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \ 00234 } \ 00235 } else \ 00236 if(dy==0) { \ 00237 nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \ 00238 nb2 = !(RFLAG(lev, i+dx,j,k, SRCS(lev))&(CFFluid|CFInter)); \ 00239 if((nb1)&&(!nb2)) { \ 00240 \ 00241 newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \ 00242 } else \ 00243 if((!nb1)&&(nb2)) { \ 00244 \ 00245 newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \ 00246 } else { \ 00247 \ 00248 newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \ 00249 } \ 00250 \ 00251 } else \ 00252 \ 00253 { \ 00254 \ 00255 nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \ 00256 nb2 = !(RFLAG(lev, i,j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \ 00257 if((nb1)&&(!nb2)) { \ 00258 \ 00259 newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \ 00260 } else \ 00261 if((!nb1)&&(nb2)) { \ 00262 \ 00263 newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \ 00264 } else { \ 00265 \ 00266 newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \ 00267 } \ 00268 } \ 00269 \ 00270 if(mnbf & CFBndPartslip) { \ 00271 const LbmFloat partv = mObjectPartslips[(int)(mnbf>>24)]; \ 00272 \ 00273 m[l] = (RAC(ccel, this->dfInv[l] ) +movadd /* d *(1./1.) */ ) * partv + newval * (1.0-partv); \ 00274 } else { \ 00275 m[l] = newval; \ 00276 } \ 00277 \ 00278 00279 00280 00281 00282 // complete default stream&collide, 2d/3d 00283 /* read distribution funtions of adjacent cells = sweep step */ 00284 #if OPT3D==0 00285 00286 #if FSGR_STRICT_DEBUG==1 00287 #define MARKCELLCHECK \ 00288 debugMarkCell(lev,i,j,k); CAUSE_PANIC; 00289 #define STREAMCHECK(id,ni,nj,nk,nl) \ 00290 if((!(m[nl] > -1.0) && (m[nl]<1.0)) ) {\ 00291 errMsg("STREAMCHECK","ID"<<id<<" Invalid streamed DF nl"<<nl<<" value:"<<m[nl]<<" at "<<PRINT_IJK<<" from "<<PRINT_VEC(ni,nj,nk)<<" nl"<<(nl)<<\ 00292 " nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther) ); \ 00293 /*FORDF0{ errMsg("STREAMCHECK"," at "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); } */ \ 00294 MARKCELLCHECK; \ 00295 m[nl] = dfEquil[nl]; /* REPAIR */ \ 00296 } 00297 #define COLLCHECK \ 00298 if( (rho>2.0) || (rho<-1.0) || (ABS(ux)>1.0) || (ABS(uy)>1.0) |(ABS(uz)>1.0) ) {\ 00299 errMsg("COLLCHECK","Invalid collision values r:"<<rho<<" u:"PRINT_VEC(ux,uy,uz)<<" at? "<<PRINT_IJK ); \ 00300 /*FORDF0{ errMsg("COLLCHECK"," at? "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); }*/ \ 00301 rho=ux=uy=uz= 0.; /* REPAIR */ \ 00302 MARKCELLCHECK; \ 00303 } 00304 #else 00305 #define STREAMCHECK(id, ni,nj,nk,nl) 00306 #define COLLCHECK 00307 #endif 00308 00309 // careful ux,uy,uz need to be inited before! 00310 #define DEFAULT_STREAM \ 00311 m[dC] = RAC(ccel,dC); \ 00312 STREAMCHECK(1, i,j,k, dC); \ 00313 FORDF1 { \ 00314 CellFlagType nbf = nbflag[ this->dfInv[l] ]; \ 00315 if(nbf & CFBnd) { \ 00316 if(nbf & CFBndNoslip) { \ 00317 \ 00318 m[l] = RAC(ccel, this->dfInv[l] ); \ 00319 LBMDS_ADDMOV(this->dfInv[l],l); \ 00320 STREAMCHECK(2, i,j,k, l); \ 00321 } else if(nbf & (CFBndFreeslip|CFBndPartslip)) { \ 00322 \ 00323 if(l<=LBMDIM*2) { \ 00324 m[l] = RAC(ccel, this->dfInv[l] ); STREAMCHECK(3, i,j,k, l); \ 00325 LBMDS_ADDMOV(this->dfInv[l],l); \ 00326 } else { \ 00327 const int inv_l = this->dfInv[l]; \ 00328 DEFAULT_STREAM_FREESLIP(l,inv_l,nbf); \ 00329 } \ 00330 \ 00331 } \ 00332 else { \ 00333 errMsg("LbmFsgrSolver","Invalid Bnd type at "<<PRINT_IJK<<" f"<<convertCellFlagType2String(nbf)<<",nbdir"<<this->dfInv[l] ); \ 00334 } \ 00335 } else { \ 00336 m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \ 00337 if(RFLAG(lev, i,j,k, mLevel[lev].setCurr)&CFFluid) { \ 00338 if(!(nbf&(CFFluid|CFInter)) ) { \ 00339 int ni=i+this->dfVecX[this->dfInv[l]], nj=j+this->dfVecY[this->dfInv[l]], nk=k+this->dfVecZ[this->dfInv[l]]; \ 00340 errMsg("STREAMCHECK"," Invalid nbflag, streamed DF l"<<l<<" value:"<<m[l]<<" at "<<PRINT_IJK<<" from "<< \ 00341 PRINT_VEC(ni,nj,nk) <<" l"<<(l)<< \ 00342 " nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther) ); \ 00343 \ 00344 \ 00345 } } \ 00346 STREAMCHECK(4, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \ 00347 } \ 00348 } \ 00349 00350 00351 00352 00353 // careful ux,uy,uz need to be inited before! 00354 #define DEFAULT_COLLIDEG(grav) \ 00355 this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), grav, mLevel[lev].lcsmago, &mDebugOmegaRet, &lcsmqo ); \ 00356 CSMOMEGA_STATS(lev,mDebugOmegaRet); \ 00357 FORDF0 { RAC(tcel,l) = m[l]; } \ 00358 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00359 COLLCHECK; \ 00360 00361 00362 00363 #define OPTIMIZED_STREAMCOLLIDE \ 00364 m[0] = RAC(ccel,0); \ 00365 FORDF1 { \ 00366 \ 00367 if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC; \ 00368 } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \ 00369 STREAMCHECK(8, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \ 00370 } \ 00371 rho=m[0]; \ 00372 DEFAULT_COLLIDEG(mLevel[lev].gravity) \ 00373 00374 00375 00376 #define OPTIMIZED_STREAMCOLLIDE___UNUSED \ 00377 \ 00378 this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo ); \ 00379 CSMOMEGA_STATS(lev,mDebugOmegaRet); \ 00380 FORDF0 { RAC(tcel,l) = m[l]; } \ 00381 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00382 COLLCHECK; \ 00383 00384 00385 00386 #else // 3D, opt OPT3D==true 00387 00388 00389 // default stream opt3d add moving bc val 00390 #define DEFAULT_STREAM \ 00391 m[dC] = RAC(ccel,dC); \ 00392 \ 00393 if((!nbored & CFBnd)) { \ 00394 \ 00395 m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \ 00396 m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \ 00397 m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \ 00398 m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \ 00399 m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \ 00400 m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \ 00401 } else { \ 00402 \ 00403 if(nbflag[dS ]&CFBnd) { m[dN ] = RAC(ccel,dS ); LBMDS_ADDMOV(dS ,dN ); } else { m[dN ] = CSRC_N ; } \ 00404 if(nbflag[dN ]&CFBnd) { m[dS ] = RAC(ccel,dN ); LBMDS_ADDMOV(dN ,dS ); } else { m[dS ] = CSRC_S ; } \ 00405 if(nbflag[dW ]&CFBnd) { m[dE ] = RAC(ccel,dW ); LBMDS_ADDMOV(dW ,dE ); } else { m[dE ] = CSRC_E ; } \ 00406 if(nbflag[dE ]&CFBnd) { m[dW ] = RAC(ccel,dE ); LBMDS_ADDMOV(dE ,dW ); } else { m[dW ] = CSRC_W ; } \ 00407 if(nbflag[dB ]&CFBnd) { m[dT ] = RAC(ccel,dB ); LBMDS_ADDMOV(dB ,dT ); } else { m[dT ] = CSRC_T ; } \ 00408 if(nbflag[dT ]&CFBnd) { m[dB ] = RAC(ccel,dT ); LBMDS_ADDMOV(dT ,dB ); } else { m[dB ] = CSRC_B ; } \ 00409 \ 00410 \ 00411 if(nbflag[dSW]&CFBnd) { if(nbflag[dSW]&CFBndNoslip){ m[dNE] = RAC(ccel,dSW); LBMDS_ADDMOV(dSW,dNE); }else{ DEFAULT_STREAM_FREESLIP(dNE,dSW,nbflag[dSW]);} } else { m[dNE] = CSRC_NE; } \ 00412 if(nbflag[dSE]&CFBnd) { if(nbflag[dSE]&CFBndNoslip){ m[dNW] = RAC(ccel,dSE); LBMDS_ADDMOV(dSE,dNW); }else{ DEFAULT_STREAM_FREESLIP(dNW,dSE,nbflag[dSE]);} } else { m[dNW] = CSRC_NW; } \ 00413 if(nbflag[dNW]&CFBnd) { if(nbflag[dNW]&CFBndNoslip){ m[dSE] = RAC(ccel,dNW); LBMDS_ADDMOV(dNW,dSE); }else{ DEFAULT_STREAM_FREESLIP(dSE,dNW,nbflag[dNW]);} } else { m[dSE] = CSRC_SE; } \ 00414 if(nbflag[dNE]&CFBnd) { if(nbflag[dNE]&CFBndNoslip){ m[dSW] = RAC(ccel,dNE); LBMDS_ADDMOV(dNE,dSW); }else{ DEFAULT_STREAM_FREESLIP(dSW,dNE,nbflag[dNE]);} } else { m[dSW] = CSRC_SW; } \ 00415 if(nbflag[dSB]&CFBnd) { if(nbflag[dSB]&CFBndNoslip){ m[dNT] = RAC(ccel,dSB); LBMDS_ADDMOV(dSB,dNT); }else{ DEFAULT_STREAM_FREESLIP(dNT,dSB,nbflag[dSB]);} } else { m[dNT] = CSRC_NT; } \ 00416 if(nbflag[dST]&CFBnd) { if(nbflag[dST]&CFBndNoslip){ m[dNB] = RAC(ccel,dST); LBMDS_ADDMOV(dST,dNB); }else{ DEFAULT_STREAM_FREESLIP(dNB,dST,nbflag[dST]);} } else { m[dNB] = CSRC_NB; } \ 00417 if(nbflag[dNB]&CFBnd) { if(nbflag[dNB]&CFBndNoslip){ m[dST] = RAC(ccel,dNB); LBMDS_ADDMOV(dNB,dST); }else{ DEFAULT_STREAM_FREESLIP(dST,dNB,nbflag[dNB]);} } else { m[dST] = CSRC_ST; } \ 00418 if(nbflag[dNT]&CFBnd) { if(nbflag[dNT]&CFBndNoslip){ m[dSB] = RAC(ccel,dNT); LBMDS_ADDMOV(dNT,dSB); }else{ DEFAULT_STREAM_FREESLIP(dSB,dNT,nbflag[dNT]);} } else { m[dSB] = CSRC_SB; } \ 00419 if(nbflag[dWB]&CFBnd) { if(nbflag[dWB]&CFBndNoslip){ m[dET] = RAC(ccel,dWB); LBMDS_ADDMOV(dWB,dET); }else{ DEFAULT_STREAM_FREESLIP(dET,dWB,nbflag[dWB]);} } else { m[dET] = CSRC_ET; } \ 00420 if(nbflag[dWT]&CFBnd) { if(nbflag[dWT]&CFBndNoslip){ m[dEB] = RAC(ccel,dWT); LBMDS_ADDMOV(dWT,dEB); }else{ DEFAULT_STREAM_FREESLIP(dEB,dWT,nbflag[dWT]);} } else { m[dEB] = CSRC_EB; } \ 00421 if(nbflag[dEB]&CFBnd) { if(nbflag[dEB]&CFBndNoslip){ m[dWT] = RAC(ccel,dEB); LBMDS_ADDMOV(dEB,dWT); }else{ DEFAULT_STREAM_FREESLIP(dWT,dEB,nbflag[dEB]);} } else { m[dWT] = CSRC_WT; } \ 00422 if(nbflag[dET]&CFBnd) { if(nbflag[dET]&CFBndNoslip){ m[dWB] = RAC(ccel,dET); LBMDS_ADDMOV(dET,dWB); }else{ DEFAULT_STREAM_FREESLIP(dWB,dET,nbflag[dET]);} } else { m[dWB] = CSRC_WB; } \ 00423 } \ 00424 00425 00426 00427 00428 00429 #define COLL_CALCULATE_DFEQ(dstarray) \ 00430 dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \ 00431 dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \ 00432 dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \ 00433 dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \ 00434 dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \ 00435 dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB; \ 00436 00437 00438 00439 #define COLL_CALCULATE_NONEQTENSOR(csolev, srcArray ) \ 00440 lcsmqadd = (srcArray##NE - lcsmeq[ dNE ]); \ 00441 lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \ 00442 lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \ 00443 lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \ 00444 lcsmqo = (lcsmqadd* lcsmqadd); \ 00445 lcsmqadd = (srcArray##ET - lcsmeq[ dET ]); \ 00446 lcsmqadd -= (srcArray##EB - lcsmeq[ dEB ]); \ 00447 lcsmqadd -= (srcArray##WT - lcsmeq[ dWT ]); \ 00448 lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \ 00449 lcsmqo += (lcsmqadd* lcsmqadd); \ 00450 lcsmqadd = (srcArray##NT - lcsmeq[ dNT ]); \ 00451 lcsmqadd -= (srcArray##NB - lcsmeq[ dNB ]); \ 00452 lcsmqadd -= (srcArray##ST - lcsmeq[ dST ]); \ 00453 lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \ 00454 lcsmqo += (lcsmqadd* lcsmqadd); \ 00455 lcsmqo *= 2.0; \ 00456 lcsmqadd = (srcArray##E - lcsmeq[ dE ]); \ 00457 lcsmqadd += (srcArray##W - lcsmeq[ dW ]); \ 00458 lcsmqadd += (srcArray##NE - lcsmeq[ dNE ]); \ 00459 lcsmqadd += (srcArray##NW - lcsmeq[ dNW ]); \ 00460 lcsmqadd += (srcArray##SE - lcsmeq[ dSE ]); \ 00461 lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \ 00462 lcsmqadd += (srcArray##ET - lcsmeq[ dET ]); \ 00463 lcsmqadd += (srcArray##EB - lcsmeq[ dEB ]); \ 00464 lcsmqadd += (srcArray##WT - lcsmeq[ dWT ]); \ 00465 lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \ 00466 lcsmqo += (lcsmqadd* lcsmqadd); \ 00467 lcsmqadd = (srcArray##N - lcsmeq[ dN ]); \ 00468 lcsmqadd += (srcArray##S - lcsmeq[ dS ]); \ 00469 lcsmqadd += (srcArray##NE - lcsmeq[ dNE ]); \ 00470 lcsmqadd += (srcArray##NW - lcsmeq[ dNW ]); \ 00471 lcsmqadd += (srcArray##SE - lcsmeq[ dSE ]); \ 00472 lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \ 00473 lcsmqadd += (srcArray##NT - lcsmeq[ dNT ]); \ 00474 lcsmqadd += (srcArray##NB - lcsmeq[ dNB ]); \ 00475 lcsmqadd += (srcArray##ST - lcsmeq[ dST ]); \ 00476 lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \ 00477 lcsmqo += (lcsmqadd* lcsmqadd); \ 00478 lcsmqadd = (srcArray##T - lcsmeq[ dT ]); \ 00479 lcsmqadd += (srcArray##B - lcsmeq[ dB ]); \ 00480 lcsmqadd += (srcArray##NT - lcsmeq[ dNT ]); \ 00481 lcsmqadd += (srcArray##NB - lcsmeq[ dNB ]); \ 00482 lcsmqadd += (srcArray##ST - lcsmeq[ dST ]); \ 00483 lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \ 00484 lcsmqadd += (srcArray##ET - lcsmeq[ dET ]); \ 00485 lcsmqadd += (srcArray##EB - lcsmeq[ dEB ]); \ 00486 lcsmqadd += (srcArray##WT - lcsmeq[ dWT ]); \ 00487 lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \ 00488 lcsmqo += (lcsmqadd* lcsmqadd); \ 00489 lcsmqo = sqrt(lcsmqo); \ 00490 00491 00492 00493 // COLL_CALCULATE_CSMOMEGAVAL(csolev, lcsmomega); 00494 00495 // careful - need lcsmqo 00496 #define COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega ) \ 00497 dstomega = 1.0/ \ 00498 ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*( \ 00499 -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \ 00500 / (6.0*mLevel[(csolev)].lcsmago_sqr)) \ 00501 ) +0.5 ); \ 00502 00503 00504 00505 #define DEFAULT_COLLIDE_LES(grav) \ 00506 rho = + MSRC_C + MSRC_N \ 00507 + MSRC_S + MSRC_E \ 00508 + MSRC_W + MSRC_T \ 00509 + MSRC_B + MSRC_NE \ 00510 + MSRC_NW + MSRC_SE \ 00511 + MSRC_SW + MSRC_NT \ 00512 + MSRC_NB + MSRC_ST \ 00513 + MSRC_SB + MSRC_ET \ 00514 + MSRC_EB + MSRC_WT \ 00515 + MSRC_WB; \ 00516 \ 00517 ux = MSRC_E - MSRC_W \ 00518 + MSRC_NE - MSRC_NW \ 00519 + MSRC_SE - MSRC_SW \ 00520 + MSRC_ET + MSRC_EB \ 00521 - MSRC_WT - MSRC_WB ; \ 00522 \ 00523 uy = MSRC_N - MSRC_S \ 00524 + MSRC_NE + MSRC_NW \ 00525 - MSRC_SE - MSRC_SW \ 00526 + MSRC_NT + MSRC_NB \ 00527 - MSRC_ST - MSRC_SB ; \ 00528 \ 00529 uz = MSRC_T - MSRC_B \ 00530 + MSRC_NT - MSRC_NB \ 00531 + MSRC_ST - MSRC_SB \ 00532 + MSRC_ET - MSRC_EB \ 00533 + MSRC_WT - MSRC_WB ; \ 00534 PRECOLLIDE_MODS(rho,ux,uy,uz, grav); \ 00535 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00536 COLL_CALCULATE_DFEQ(lcsmeq); \ 00537 COLL_CALCULATE_NONEQTENSOR(lev, MSRC_); \ 00538 COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \ 00539 CSMOMEGA_STATS(lev,lcsmomega); \ 00540 \ 00541 RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C + lcsmomega*EQC ; \ 00542 \ 00543 RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N + lcsmomega*lcsmeq[ dN ]; \ 00544 RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S + lcsmomega*lcsmeq[ dS ]; \ 00545 RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E + lcsmomega*lcsmeq[ dE ]; \ 00546 RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W + lcsmomega*lcsmeq[ dW ]; \ 00547 RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T + lcsmomega*lcsmeq[ dT ]; \ 00548 RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B + lcsmomega*lcsmeq[ dB ]; \ 00549 \ 00550 RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \ 00551 RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \ 00552 RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \ 00553 RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \ 00554 RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \ 00555 RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \ 00556 RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \ 00557 RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \ 00558 RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \ 00559 RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \ 00560 RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \ 00561 RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \ 00562 00563 00564 00565 #define DEFAULT_COLLIDE_NOLES(grav) \ 00566 rho = + MSRC_C + MSRC_N \ 00567 + MSRC_S + MSRC_E \ 00568 + MSRC_W + MSRC_T \ 00569 + MSRC_B + MSRC_NE \ 00570 + MSRC_NW + MSRC_SE \ 00571 + MSRC_SW + MSRC_NT \ 00572 + MSRC_NB + MSRC_ST \ 00573 + MSRC_SB + MSRC_ET \ 00574 + MSRC_EB + MSRC_WT \ 00575 + MSRC_WB; \ 00576 \ 00577 ux = MSRC_E - MSRC_W \ 00578 + MSRC_NE - MSRC_NW \ 00579 + MSRC_SE - MSRC_SW \ 00580 + MSRC_ET + MSRC_EB \ 00581 - MSRC_WT - MSRC_WB ; \ 00582 \ 00583 uy = MSRC_N - MSRC_S \ 00584 + MSRC_NE + MSRC_NW \ 00585 - MSRC_SE - MSRC_SW \ 00586 + MSRC_NT + MSRC_NB \ 00587 - MSRC_ST - MSRC_SB ; \ 00588 \ 00589 uz = MSRC_T - MSRC_B \ 00590 + MSRC_NT - MSRC_NB \ 00591 + MSRC_ST - MSRC_SB \ 00592 + MSRC_ET - MSRC_EB \ 00593 + MSRC_WT - MSRC_WB ; \ 00594 PRECOLLIDE_MODS(rho, ux,uy,uz, grav); \ 00595 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00596 \ 00597 RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C + OMEGA(lev)*EQC ; \ 00598 \ 00599 RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N + OMEGA(lev)*EQN ; \ 00600 RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S + OMEGA(lev)*EQS ; \ 00601 RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E + OMEGA(lev)*EQE ; \ 00602 RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W + OMEGA(lev)*EQW ; \ 00603 RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T + OMEGA(lev)*EQT ; \ 00604 RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B + OMEGA(lev)*EQB ; \ 00605 \ 00606 RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \ 00607 RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \ 00608 RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \ 00609 RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \ 00610 RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \ 00611 RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \ 00612 RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \ 00613 RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \ 00614 RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \ 00615 RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \ 00616 RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \ 00617 RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB; \ 00618 00619 00620 00621 00622 00623 #define OPTIMIZED_STREAMCOLLIDE_LES \ 00624 \ 00625 m[dC ] = CSRC_C ; \ 00626 m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \ 00627 m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \ 00628 m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \ 00629 m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \ 00630 m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \ 00631 m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \ 00632 \ 00633 rho = MSRC_C + MSRC_N + MSRC_S + MSRC_E + MSRC_W + MSRC_T \ 00634 + MSRC_B + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \ 00635 + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \ 00636 ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \ 00637 + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB; \ 00638 uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \ 00639 + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB; \ 00640 uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \ 00641 + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB; \ 00642 PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \ 00643 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00644 COLL_CALCULATE_DFEQ(lcsmeq); \ 00645 COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \ 00646 COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \ 00647 CSMOMEGA_STATS(lev,lcsmomega); \ 00648 \ 00649 RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C + lcsmomega*EQC ; \ 00650 RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N + lcsmomega*lcsmeq[ dN ]; \ 00651 RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S + lcsmomega*lcsmeq[ dS ]; \ 00652 RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E + lcsmomega*lcsmeq[ dE ]; \ 00653 RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W + lcsmomega*lcsmeq[ dW ]; \ 00654 RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T + lcsmomega*lcsmeq[ dT ]; \ 00655 RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B + lcsmomega*lcsmeq[ dB ]; \ 00656 \ 00657 RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \ 00658 RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \ 00659 RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \ 00660 RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \ 00661 \ 00662 RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \ 00663 RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \ 00664 RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \ 00665 RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \ 00666 \ 00667 RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \ 00668 RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \ 00669 RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \ 00670 RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \ 00671 00672 00673 00674 #define OPTIMIZED_STREAMCOLLIDE_UNUSED \ 00675 \ 00676 rho = CSRC_C + CSRC_N + CSRC_S + CSRC_E + CSRC_W + CSRC_T \ 00677 + CSRC_B + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \ 00678 + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \ 00679 ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \ 00680 + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \ 00681 uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \ 00682 + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \ 00683 uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \ 00684 + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \ 00685 PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \ 00686 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00687 COLL_CALCULATE_DFEQ(lcsmeq); \ 00688 COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \ 00689 COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \ 00690 \ 00691 RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C + lcsmomega*EQC ; \ 00692 RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N + lcsmomega*lcsmeq[ dN ]; \ 00693 RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S + lcsmomega*lcsmeq[ dS ]; \ 00694 RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E + lcsmomega*lcsmeq[ dE ]; \ 00695 RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W + lcsmomega*lcsmeq[ dW ]; \ 00696 RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T + lcsmomega*lcsmeq[ dT ]; \ 00697 RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B + lcsmomega*lcsmeq[ dB ]; \ 00698 \ 00699 RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE]; \ 00700 RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW]; \ 00701 RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE]; \ 00702 RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \ 00703 \ 00704 RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT]; \ 00705 RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB]; \ 00706 RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST]; \ 00707 RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \ 00708 \ 00709 RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET]; \ 00710 RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \ 00711 RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT]; \ 00712 RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB]; \ 00713 00714 00715 00716 #define OPTIMIZED_STREAMCOLLIDE_NOLES \ 00717 \ 00718 rho = CSRC_C + CSRC_N + CSRC_S + CSRC_E + CSRC_W + CSRC_T \ 00719 + CSRC_B + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \ 00720 + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \ 00721 ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \ 00722 + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \ 00723 uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \ 00724 + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \ 00725 uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \ 00726 + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \ 00727 PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \ 00728 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00729 RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C + OMEGA(lev)*EQC ; \ 00730 RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N + OMEGA(lev)*EQN ; \ 00731 RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S + OMEGA(lev)*EQS ; \ 00732 RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E + OMEGA(lev)*EQE ; \ 00733 RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W + OMEGA(lev)*EQW ; \ 00734 RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T + OMEGA(lev)*EQT ; \ 00735 RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B + OMEGA(lev)*EQB ; \ 00736 \ 00737 RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE; \ 00738 RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW; \ 00739 RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE; \ 00740 RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \ 00741 \ 00742 RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT; \ 00743 RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB; \ 00744 RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST; \ 00745 RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \ 00746 \ 00747 RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET; \ 00748 RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \ 00749 RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT; \ 00750 RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB; \ 00751 00752 00753 00754 00755 00756 // LES switching for OPT3D 00757 #if USE_LES==1 00758 #define DEFAULT_COLLIDEG(grav) DEFAULT_COLLIDE_LES(grav) 00759 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_LES 00760 #else 00761 #define DEFAULT_COLLIDEG(grav) DEFAULT_COLLIDE_NOLES(grav) 00762 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_NOLES 00763 #endif 00764 00765 #endif // 3D, opt OPT3D==true 00766 00767 #define USQRMAXCHECK(Cusqr,Cux,Cuy,Cuz, CmMaxVlen,CmMxvx,CmMxvy,CmMxvz) \ 00768 if(Cusqr>CmMaxVlen) { \ 00769 CmMxvx = Cux; CmMxvy = Cuy; CmMxvz = Cuz; CmMaxVlen = Cusqr; \ 00770 } /* stats */ 00771 00772 00773 00774 /****************************************************************************** 00775 * interpolateCellFromCoarse macros 00776 *****************************************************************************/ 00777 00778 00779 // WOXDY_N = Weight Order X Dimension Y _ number N 00780 #define WO1D1 ( 1.0/ 2.0) 00781 #define WO1D2 ( 1.0/ 4.0) 00782 #define WO1D3 ( 1.0/ 8.0) 00783 00784 #define WO2D1_1 (-1.0/16.0) 00785 #define WO2D1_9 ( 9.0/16.0) 00786 00787 #define WO2D2_11 (WO2D1_1 * WO2D1_1) 00788 #define WO2D2_19 (WO2D1_9 * WO2D1_1) 00789 #define WO2D2_91 (WO2D1_9 * WO2D1_1) 00790 #define WO2D2_99 (WO2D1_9 * WO2D1_9) 00791 00792 #define WO2D3_111 (WO2D1_1 * WO2D1_1 * WO2D1_1) 00793 #define WO2D3_191 (WO2D1_9 * WO2D1_1 * WO2D1_1) 00794 #define WO2D3_911 (WO2D1_9 * WO2D1_1 * WO2D1_1) 00795 #define WO2D3_991 (WO2D1_9 * WO2D1_9 * WO2D1_1) 00796 #define WO2D3_119 (WO2D1_1 * WO2D1_1 * WO2D1_9) 00797 #define WO2D3_199 (WO2D1_9 * WO2D1_1 * WO2D1_9) 00798 #define WO2D3_919 (WO2D1_9 * WO2D1_1 * WO2D1_9) 00799 #define WO2D3_999 (WO2D1_9 * WO2D1_9 * WO2D1_9) 00800 00801 #if FSGR_STRICT_DEBUG==1 00802 #define ADD_INT_DFSCHECK(alev, ai,aj,ak, at, afac, l) \ 00803 if( (((1.0-(at))>0.0) && (!(QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr , l) > -1.0 ))) || \ 00804 ((( (at))>0.0) && (!(QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setOther, l) > -1.0 ))) ){ \ 00805 errMsg("INVDFSCHECK", " l"<<(alev)<<" "<<PRINT_VEC((ai),(aj),(ak))<<" fc:"<<RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr )<<" fo:"<<RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setOther )<<" dfl"<<l ); \ 00806 debugMarkCell((alev), (ai),(aj),(ak));\ 00807 CAUSE_PANIC; \ 00808 } 00809 // end ADD_INT_DFSCHECK 00810 #define ADD_INT_FLAGCHECK(alev, ai,aj,ak, at, afac) \ 00811 if( (((1.0-(at))>0.0) && (!(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr )&(CFInter|CFFluid|CFGrCoarseInited) ))) || \ 00812 ((( (at))>0.0) && (!(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setOther)&(CFInter|CFFluid|CFGrCoarseInited) ))) ){ \ 00813 errMsg("INVFLAGCINTCHECK", " l"<<(alev)<<" at:"<<(at)<<" "<<PRINT_VEC((ai),(aj),(ak))<<\ 00814 " fc:"<< convertCellFlagType2String(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr )) <<\ 00815 " fold:"<< convertCellFlagType2String(RFLAG((alev), (ai),(aj),(ak),mLevel[(alev)].setOther )) ); \ 00816 debugMarkCell((alev), (ai),(aj),(ak));\ 00817 CAUSE_PANIC; \ 00818 } 00819 // end ADD_INT_DFSCHECK 00820 00821 #define INTUNUTCHECK(ix,iy,iz) \ 00822 if( (RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) != (CFFluid|CFGrFromCoarse)) ){\ 00823 errMsg("INTFLAGUNU_CHECK", PRINT_VEC(i,j,k)<<" child not unused at l"<<(lev+1)<<" "<<PRINT_VEC((ix),(iy),(iz))<<" flag: "<< RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) ); \ 00824 debugMarkCell((lev+1), (ix),(iy),(iz));\ 00825 CAUSE_PANIC; \ 00826 }\ 00827 RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) |= CFGrCoarseInited; \ 00828 // INTUNUTCHECK 00829 #define INTSTRICTCHECK(ix,iy,iz,caseId) \ 00830 if( QCELL(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr, l) <= 0.0 ){\ 00831 errMsg("INVDFCCELLCHECK", "caseId:"<<caseId<<" "<<PRINT_VEC(i,j,k)<<" child inter at "<<PRINT_VEC((ix),(iy),(iz))<<" invalid df "<<l<<" = "<< QCELL(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr, l) ); \ 00832 debugMarkCell((lev+1), (ix),(iy),(iz));\ 00833 CAUSE_PANIC; \ 00834 }\ 00835 // INTSTRICTCHECK 00836 00837 #else// FSGR_STRICT_DEBUG==1 00838 #define ADD_INT_FLAGCHECK(alev, ai,aj,ak, at, afac) 00839 #define ADD_INT_DFSCHECK(alev, ai,aj,ak, at, afac, l) 00840 #define INTSTRICTCHECK(x,y,z,caseId) 00841 #define INTUNUTCHECK(ix,iy,iz) 00842 #endif// FSGR_STRICT_DEBUG==1 00843 00844 00845 #if FSGR_STRICT_DEBUG==1 00846 #define INTDEBOUT \ 00847 { /*LbmFloat rho,ux,uy,uz;*/ \ 00848 rho = ux=uy=uz=0.0; \ 00849 FORDF0{ LbmFloat m = QCELL(lev,i,j,k, dstSet, l); \ 00850 rho += m; ux += (this->dfDvecX[l]*m); uy += (this->dfDvecY[l]*m); uz += (this->dfDvecZ[l]*m); \ 00851 if(ABS(m)>1.0) { errMsg("interpolateCellFromCoarse", "ICFC_DFCHECK cell "<<PRINT_IJK<<" m"<<l<<":"<< m );CAUSE_PANIC;}\ 00852 /*errMsg("interpolateCellFromCoarse", " cell "<<PRINT_IJK<<" df"<<l<<":"<<m );*/ \ 00853 } \ 00854 /*if(this->mPanic) { errMsg("interpolateCellFromCoarse", "ICFC_DFOUT cell "<<PRINT_IJK<<" rho:"<<rho<<" u:"<<PRINT_VEC(ux,uy,uz)<<" b"<<PRINT_VEC(betx,bety,betz) ); }*/ \ 00855 if(markNbs) errMsg("interpolateCellFromCoarse", " cell "<<PRINT_IJK<<" rho:"<<rho<<" u:"<<PRINT_VEC(ux,uy,uz)<<" b"<<PRINT_VEC(betx,bety,betz) ); \ 00856 /*errMsg("interpolateCellFromCoarse", "ICFC_DFDEBUG cell "<<PRINT_IJK<<" rho:"<<rho<<" u:"<<PRINT_VEC(ux,uy,uz)<<" b"<<PRINT_VEC(betx,bety,betz) ); */\ 00857 } \ 00858 /* both cases are ok to interpolate */ \ 00859 if( (!(RFLAG(lev,i,j,k, dstSet) & CFGrFromCoarse)) && \ 00860 (!(RFLAG(lev,i,j,k, dstSet) & CFUnused)) ) { \ 00861 /* might also have CFGrCoarseInited (shouldnt be a problem here)*/ \ 00862 errMsg("interpolateCellFromCoarse", "CHECK cell not CFGrFromCoarse? "<<PRINT_IJK<<" flag:"<< RFLAG(lev,i,j,k, dstSet)<<" fstr:"<<convertCellFlagType2String( RFLAG(lev,i,j,k, dstSet) )); \ 00863 /* FIXME check this warning...? return; this can happen !? */ \ 00864 /*CAUSE_PANIC;*/ \ 00865 } \ 00866 // end INTDEBOUT 00867 #else // FSGR_STRICT_DEBUG==1 00868 #define INTDEBOUT 00869 #endif // FSGR_STRICT_DEBUG==1 00870 00871 00872 // t=0.0 -> only current 00873 // t=0.5 -> mix 00874 // t=1.0 -> only other 00875 #if OPT3D==0 00876 #define ADD_INT_DFS(alev, ai,aj,ak, at, afac) \ 00877 ADD_INT_FLAGCHECK(alev, ai,aj,ak, at, afac); \ 00878 FORDF0{ \ 00879 LbmFloat df = ( \ 00880 QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr , l)*(1.0-(at)) + \ 00881 QCELL((alev), (ai),(aj),(ak),mLevel[(alev)].setOther, l)*( (at)) \ 00882 ) ; \ 00883 ADD_INT_DFSCHECK(alev, ai,aj,ak, at, afac, l); \ 00884 df *= (afac); \ 00885 rho += df; \ 00886 ux += (this->dfDvecX[l]*df); \ 00887 uy += (this->dfDvecY[l]*df); \ 00888 uz += (this->dfDvecZ[l]*df); \ 00889 intDf[l] += df; \ 00890 } 00891 // write interpolated dfs back to cell (correct non-eq. parts) 00892 #define IDF_WRITEBACK_ \ 00893 FORDF0{ \ 00894 LbmFloat eq = getCollideEq(l, rho,ux,uy,uz);\ 00895 QCELL(lev,i,j,k, dstSet, l) = (eq+ (intDf[l]-eq)*mDfScaleDown);\ 00896 } \ 00897 /* check that all values are ok */ \ 00898 INTDEBOUT 00899 #define IDF_WRITEBACK \ 00900 LbmFloat omegaDst, omegaSrc;\ 00901 /* smago new */ \ 00902 LbmFloat feq[LBM_DFNUM]; \ 00903 LbmFloat dfScale = mDfScaleDown; \ 00904 FORDF0{ \ 00905 feq[l] = getCollideEq(l, rho,ux,uy,uz); \ 00906 } \ 00907 if(mLevel[lev ].lcsmago>0.0) {\ 00908 LbmFloat Qo = this->getLesNoneqTensorCoeff(intDf,feq); \ 00909 omegaDst = this->getLesOmega(mLevel[lev+0].omega,mLevel[lev+0].lcsmago,Qo); \ 00910 omegaSrc = this->getLesOmega(mLevel[lev-1].omega,mLevel[lev-1].lcsmago,Qo); \ 00911 } else {\ 00912 omegaDst = mLevel[lev+0].omega; \ 00913 omegaSrc = mLevel[lev-1].omega;\ 00914 } \ 00915 \ 00916 dfScale = (mLevel[lev+0].timestep/mLevel[lev-1].timestep)* (1.0/omegaDst-1.0)/ (1.0/omegaSrc-1.0); \ 00917 FORDF0{ \ 00918 /*errMsg("SMAGO"," org"<<mDfScaleDown<<" n"<<dfScale<<" qc"<< QCELL(lev,i,j,k, dstSet, l)<<" idf"<<intDf[l]<<" eq"<<feq[l] ); */ \ 00919 QCELL(lev,i,j,k, dstSet, l) = (feq[l]+ (intDf[l]-feq[l])*dfScale);\ 00920 } \ 00921 /* check that all values are ok */ \ 00922 INTDEBOUT 00923 00924 #else //OPT3D==0 00925 00926 #define ADDALLVALS \ 00927 addVal = addDfFacT * RAC(addfcel , dC ); \ 00928 intDf[dC ] += addVal; rho += addVal; \ 00929 addVal = addDfFacT * RAC(addfcel , dN ); \ 00930 uy+=addVal; intDf[dN ] += addVal; rho += addVal; \ 00931 addVal = addDfFacT * RAC(addfcel , dS ); \ 00932 uy-=addVal; intDf[dS ] += addVal; rho += addVal; \ 00933 addVal = addDfFacT * RAC(addfcel , dE ); \ 00934 ux+=addVal; intDf[dE ] += addVal; rho += addVal; \ 00935 addVal = addDfFacT * RAC(addfcel , dW ); \ 00936 ux-=addVal; intDf[dW ] += addVal; rho += addVal; \ 00937 addVal = addDfFacT * RAC(addfcel , dT ); \ 00938 uz+=addVal; intDf[dT ] += addVal; rho += addVal; \ 00939 addVal = addDfFacT * RAC(addfcel , dB ); \ 00940 uz-=addVal; intDf[dB ] += addVal; rho += addVal; \ 00941 addVal = addDfFacT * RAC(addfcel , dNE); \ 00942 ux+=addVal; uy+=addVal; intDf[dNE] += addVal; rho += addVal; \ 00943 addVal = addDfFacT * RAC(addfcel , dNW); \ 00944 ux-=addVal; uy+=addVal; intDf[dNW] += addVal; rho += addVal; \ 00945 addVal = addDfFacT * RAC(addfcel , dSE); \ 00946 ux+=addVal; uy-=addVal; intDf[dSE] += addVal; rho += addVal; \ 00947 addVal = addDfFacT * RAC(addfcel , dSW); \ 00948 ux-=addVal; uy-=addVal; intDf[dSW] += addVal; rho += addVal; \ 00949 addVal = addDfFacT * RAC(addfcel , dNT); \ 00950 uy+=addVal; uz+=addVal; intDf[dNT] += addVal; rho += addVal; \ 00951 addVal = addDfFacT * RAC(addfcel , dNB); \ 00952 uy+=addVal; uz-=addVal; intDf[dNB] += addVal; rho += addVal; \ 00953 addVal = addDfFacT * RAC(addfcel , dST); \ 00954 uy-=addVal; uz+=addVal; intDf[dST] += addVal; rho += addVal; \ 00955 addVal = addDfFacT * RAC(addfcel , dSB); \ 00956 uy-=addVal; uz-=addVal; intDf[dSB] += addVal; rho += addVal; \ 00957 addVal = addDfFacT * RAC(addfcel , dET); \ 00958 ux+=addVal; uz+=addVal; intDf[dET] += addVal; rho += addVal; \ 00959 addVal = addDfFacT * RAC(addfcel , dEB); \ 00960 ux+=addVal; uz-=addVal; intDf[dEB] += addVal; rho += addVal; \ 00961 addVal = addDfFacT * RAC(addfcel , dWT); \ 00962 ux-=addVal; uz+=addVal; intDf[dWT] += addVal; rho += addVal; \ 00963 addVal = addDfFacT * RAC(addfcel , dWB); \ 00964 ux-=addVal; uz-=addVal; intDf[dWB] += addVal; rho += addVal; 00965 00966 #define ADD_INT_DFS(alev, ai,aj,ak, at, afac) \ 00967 addDfFacT = at*afac; \ 00968 addfcel = RACPNT((alev), (ai),(aj),(ak),mLevel[(alev)].setOther); \ 00969 ADDALLVALS\ 00970 addDfFacT = (1.0-at)*afac; \ 00971 addfcel = RACPNT((alev), (ai),(aj),(ak),mLevel[(alev)].setCurr); \ 00972 ADDALLVALS 00973 00974 // also ugly... 00975 #define INTDF_C intDf[dC ] 00976 #define INTDF_N intDf[dN ] 00977 #define INTDF_S intDf[dS ] 00978 #define INTDF_E intDf[dE ] 00979 #define INTDF_W intDf[dW ] 00980 #define INTDF_T intDf[dT ] 00981 #define INTDF_B intDf[dB ] 00982 #define INTDF_NE intDf[dNE] 00983 #define INTDF_NW intDf[dNW] 00984 #define INTDF_SE intDf[dSE] 00985 #define INTDF_SW intDf[dSW] 00986 #define INTDF_NT intDf[dNT] 00987 #define INTDF_NB intDf[dNB] 00988 #define INTDF_ST intDf[dST] 00989 #define INTDF_SB intDf[dSB] 00990 #define INTDF_ET intDf[dET] 00991 #define INTDF_EB intDf[dEB] 00992 #define INTDF_WT intDf[dWT] 00993 #define INTDF_WB intDf[dWB] 00994 00995 00996 // write interpolated dfs back to cell (correct non-eq. parts) 00997 #define IDF_WRITEBACK_LES \ 00998 dstcell = RACPNT(lev, i,j,k,dstSet); \ 00999 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 01000 \ 01001 lcsmeq[dC] = EQC ; \ 01002 COLL_CALCULATE_DFEQ(lcsmeq); \ 01003 COLL_CALCULATE_NONEQTENSOR(lev, INTDF_ )\ 01004 COLL_CALCULATE_CSMOMEGAVAL(lev+0, lcsmDstOmega); \ 01005 COLL_CALCULATE_CSMOMEGAVAL(lev-1, lcsmSrcOmega); \ 01006 \ 01007 lcsmdfscale = (mLevel[lev+0].timestep/mLevel[lev-1].timestep)* (1.0/lcsmDstOmega-1.0)/ (1.0/lcsmSrcOmega-1.0); \ 01008 RAC(dstcell, dC ) = (lcsmeq[dC ] + (intDf[dC ]-lcsmeq[dC ] )*lcsmdfscale);\ 01009 RAC(dstcell, dN ) = (lcsmeq[dN ] + (intDf[dN ]-lcsmeq[dN ] )*lcsmdfscale);\ 01010 RAC(dstcell, dS ) = (lcsmeq[dS ] + (intDf[dS ]-lcsmeq[dS ] )*lcsmdfscale);\ 01011 RAC(dstcell, dE ) = (lcsmeq[dE ] + (intDf[dE ]-lcsmeq[dE ] )*lcsmdfscale);\ 01012 RAC(dstcell, dW ) = (lcsmeq[dW ] + (intDf[dW ]-lcsmeq[dW ] )*lcsmdfscale);\ 01013 RAC(dstcell, dT ) = (lcsmeq[dT ] + (intDf[dT ]-lcsmeq[dT ] )*lcsmdfscale);\ 01014 RAC(dstcell, dB ) = (lcsmeq[dB ] + (intDf[dB ]-lcsmeq[dB ] )*lcsmdfscale);\ 01015 RAC(dstcell, dNE) = (lcsmeq[dNE] + (intDf[dNE]-lcsmeq[dNE] )*lcsmdfscale);\ 01016 RAC(dstcell, dNW) = (lcsmeq[dNW] + (intDf[dNW]-lcsmeq[dNW] )*lcsmdfscale);\ 01017 RAC(dstcell, dSE) = (lcsmeq[dSE] + (intDf[dSE]-lcsmeq[dSE] )*lcsmdfscale);\ 01018 RAC(dstcell, dSW) = (lcsmeq[dSW] + (intDf[dSW]-lcsmeq[dSW] )*lcsmdfscale);\ 01019 RAC(dstcell, dNT) = (lcsmeq[dNT] + (intDf[dNT]-lcsmeq[dNT] )*lcsmdfscale);\ 01020 RAC(dstcell, dNB) = (lcsmeq[dNB] + (intDf[dNB]-lcsmeq[dNB] )*lcsmdfscale);\ 01021 RAC(dstcell, dST) = (lcsmeq[dST] + (intDf[dST]-lcsmeq[dST] )*lcsmdfscale);\ 01022 RAC(dstcell, dSB) = (lcsmeq[dSB] + (intDf[dSB]-lcsmeq[dSB] )*lcsmdfscale);\ 01023 RAC(dstcell, dET) = (lcsmeq[dET] + (intDf[dET]-lcsmeq[dET] )*lcsmdfscale);\ 01024 RAC(dstcell, dEB) = (lcsmeq[dEB] + (intDf[dEB]-lcsmeq[dEB] )*lcsmdfscale);\ 01025 RAC(dstcell, dWT) = (lcsmeq[dWT] + (intDf[dWT]-lcsmeq[dWT] )*lcsmdfscale);\ 01026 RAC(dstcell, dWB) = (lcsmeq[dWB] + (intDf[dWB]-lcsmeq[dWB] )*lcsmdfscale);\ 01027 /* IDF_WRITEBACK optimized */ 01028 01029 #define IDF_WRITEBACK_NOLES \ 01030 dstcell = RACPNT(lev, i,j,k,dstSet); \ 01031 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 01032 \ 01033 RAC(dstcell, dC ) = (EQC + (intDf[dC ]-EQC )*mDfScaleDown);\ 01034 RAC(dstcell, dN ) = (EQN + (intDf[dN ]-EQN )*mDfScaleDown);\ 01035 RAC(dstcell, dS ) = (EQS + (intDf[dS ]-EQS )*mDfScaleDown);\ 01036 /*old*/ RAC(dstcell, dE ) = (EQE + (intDf[dE ]-EQE )*mDfScaleDown);\ 01037 RAC(dstcell, dW ) = (EQW + (intDf[dW ]-EQW )*mDfScaleDown);\ 01038 RAC(dstcell, dT ) = (EQT + (intDf[dT ]-EQT )*mDfScaleDown);\ 01039 RAC(dstcell, dB ) = (EQB + (intDf[dB ]-EQB )*mDfScaleDown);\ 01040 /*old*/ RAC(dstcell, dNE) = (EQNE + (intDf[dNE]-EQNE )*mDfScaleDown);\ 01041 RAC(dstcell, dNW) = (EQNW + (intDf[dNW]-EQNW )*mDfScaleDown);\ 01042 RAC(dstcell, dSE) = (EQSE + (intDf[dSE]-EQSE )*mDfScaleDown);\ 01043 RAC(dstcell, dSW) = (EQSW + (intDf[dSW]-EQSW )*mDfScaleDown);\ 01044 RAC(dstcell, dNT) = (EQNT + (intDf[dNT]-EQNT )*mDfScaleDown);\ 01045 RAC(dstcell, dNB) = (EQNB + (intDf[dNB]-EQNB )*mDfScaleDown);\ 01046 RAC(dstcell, dST) = (EQST + (intDf[dST]-EQST )*mDfScaleDown);\ 01047 RAC(dstcell, dSB) = (EQSB + (intDf[dSB]-EQSB )*mDfScaleDown);\ 01048 RAC(dstcell, dET) = (EQET + (intDf[dET]-EQET )*mDfScaleDown);\ 01049 /*old*/ RAC(dstcell, dEB) = (EQEB + (intDf[dEB]-EQEB )*mDfScaleDown);\ 01050 RAC(dstcell, dWT) = (EQWT + (intDf[dWT]-EQWT )*mDfScaleDown);\ 01051 RAC(dstcell, dWB) = (EQWB + (intDf[dWB]-EQWB )*mDfScaleDown);\ 01052 /* IDF_WRITEBACK optimized */ 01053 01054 #if USE_LES==1 01055 #define IDF_WRITEBACK IDF_WRITEBACK_LES 01056 #else 01057 #define IDF_WRITEBACK IDF_WRITEBACK_NOLES 01058 #endif 01059 01060 #endif// OPT3D==0 01061 01062 01063 01064 /******************************************************************************/ 01066 /******************************************************************************/ 01067 01068 01069 inline LbmFloat LbmFsgrSolver::getLesNoneqTensorCoeff( 01070 LbmFloat df[], 01071 LbmFloat feq[] ) { 01072 LbmFloat Qo = 0.0; 01073 for(int m=0; m< ((LBMDIM*LBMDIM)-LBMDIM)/2 ; m++) { 01074 LbmFloat qadd = 0.0; 01075 for(int l=1; l<this->cDfNum; l++) { 01076 if(this->lesCoeffOffdiag[m][l]==0.0) continue; 01077 qadd += this->lesCoeffOffdiag[m][l]*(df[l]-feq[l]); 01078 } 01079 Qo += (qadd*qadd); 01080 } 01081 Qo *= 2.0; // off diag twice 01082 for(int m=0; m<LBMDIM; m++) { 01083 LbmFloat qadd = 0.0; 01084 for(int l=1; l<this->cDfNum; l++) { 01085 if(this->lesCoeffDiag[m][l]==0.0) continue; 01086 qadd += this->lesCoeffDiag[m][l]*(df[l]-feq[l]); 01087 } 01088 Qo += (qadd*qadd); 01089 } 01090 Qo = sqrt(Qo); 01091 return Qo; 01092 }; 01093 01094 inline LbmFloat LbmFsgrSolver::getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo) { 01095 const LbmFloat tau = 1.0/omega; 01096 const LbmFloat nu = (2.0*tau-1.0) * (1.0/6.0); 01097 const LbmFloat C = csmago; 01098 const LbmFloat Csqr = C*C; 01099 LbmFloat S = -nu + sqrt( nu*nu + 18.0*Csqr*Qo ) / (6.0*Csqr); 01100 return( 1.0/( 3.0*( nu+Csqr*S ) +0.5 ) ); 01101 } 01102 01103 #define DEBUG_CALCPRINTCELL(str,df) {\ 01104 LbmFloat prho=df[0], pux=0., puy=0., puz=0.; \ 01105 for(int dfl=1; dfl<this->cDfNum; dfl++) { \ 01106 prho += df[dfl]; \ 01107 pux += (this->dfDvecX[dfl]*df[dfl]); \ 01108 puy += (this->dfDvecY[dfl]*df[dfl]); \ 01109 puz += (this->dfDvecZ[dfl]*df[dfl]); \ 01110 } \ 01111 errMsg("DEBUG_CALCPRINTCELL",">"<<str<<" rho="<<prho<<" vel="<<ntlVec3Gfx(pux,puy,puz) ); \ 01112 } /* END DEBUG_CALCPRINTCELL */ 01113 01114 // "normal" collision 01115 inline void LbmFsgrSolver::collideArrays( 01116 int lev, int i, int j, int k, // position - more for debugging 01117 LbmFloat df[], 01118 LbmFloat &outrho, // out only! 01119 // velocity modifiers (returns actual velocity!) 01120 LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 01121 LbmFloat omega, 01122 LbmVec gravity, 01123 LbmFloat csmago, 01124 LbmFloat *newOmegaRet, LbmFloat *newQoRet 01125 ) { 01126 int l; 01127 LbmFloat rho=df[0]; 01128 LbmFloat ux = 0; //mux; 01129 LbmFloat uy = 0; //muy; 01130 LbmFloat uz = 0; //muz; 01131 LbmFloat feq[19]; 01132 LbmFloat omegaNew; 01133 LbmFloat Qo = 0.0; 01134 01135 for(l=1; l<this->cDfNum; l++) { 01136 rho += df[l]; 01137 ux += (this->dfDvecX[l]*df[l]); 01138 uy += (this->dfDvecY[l]*df[l]); 01139 uz += (this->dfDvecZ[l]*df[l]); 01140 } 01141 01142 01143 PRECOLLIDE_MODS(rho,ux,uy,uz, gravity); 01144 for(l=0; l<this->cDfNum; l++) { 01145 feq[l] = getCollideEq(l,rho,ux,uy,uz); 01146 } 01147 01148 if(csmago>0.0) { 01149 Qo = getLesNoneqTensorCoeff(df,feq); 01150 omegaNew = getLesOmega(omega,csmago,Qo); 01151 } else { 01152 omegaNew = omega; // smago off... 01153 } 01154 if(newOmegaRet) *newOmegaRet = omegaNew; // return value for stats 01155 if(newQoRet) *newQoRet = Qo; // return value of non-eq. stress tensor 01156 01157 for(l=0; l<this->cDfNum; l++) { 01158 df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; 01159 } 01160 //if((i==16)&&(j==10)) DEBUG_CALCPRINTCELL( "2dcoll "<<PRINT_IJK, df); 01161 01162 mux = ux; 01163 muy = uy; 01164 muz = uz; 01165 outrho = rho; 01166 01167 lev=i=j=k; // debug, remove warnings 01168 }; 01169