|
Blender
V2.59
|
00001 00004 /****************************************************************************** 00005 * 00006 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method 00007 * Copyright 2003-2006 Nils Thuerey 00008 * 00009 * Standard LBM Factory implementation 00010 * 00011 *****************************************************************************/ 00012 00013 #include "solver_class.h" 00014 #include "solver_relax.h" 00015 #include "particletracer.h" 00016 #include "loop_tools.h" 00017 #include <stdlib.h> 00018 00019 /*****************************************************************************/ 00021 /*****************************************************************************/ 00022 00023 double globdfcnt; 00024 double globdfavg[19]; 00025 double globdfmax[19]; 00026 double globdfmin[19]; 00027 extern int glob_mpindex,glob_mpnum; 00028 extern bool globOutstrForce; 00029 00030 // simulation object interface 00031 void LbmFsgrSolver::step() { 00032 stepMain(); 00033 } 00034 00035 // lbm main step 00036 void messageOutputForce(string from); 00037 void LbmFsgrSolver::stepMain() { 00038 myTime_t timestart = getTime(); 00039 00040 initLevelOmegas(); 00041 markedClearList(); // DMC clearMarkedCellsList 00042 00043 // safety check, counter reset 00044 mNumUsedCells = 0; 00045 mNumInterdCells = 0; 00046 mNumInvIfCells = 0; 00047 00048 //debugOutNnl("LbmFsgrSolver::step : "<<mStepCnt, 10); 00049 if(!mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime, 10); } 00050 //debMsgDirect( "LbmFsgrSolver::step : "<<mStepCnt<<" "); 00051 //myTime_t timestart = 0; 00052 //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 00053 00054 // time adapt 00055 mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0; 00056 00057 // init moving bc's, can change mMaxVlen 00058 initMovingObstacles(false); 00059 00060 // handle fluid control 00061 handleCpdata(); 00062 00063 // important - keep for tadap 00064 LbmFloat lastMass = mCurrentMass; 00065 mCurrentMass = mFixMass; // reset here for next step 00066 mCurrentVolume = 0.0; 00067 00068 //change to single step advance! 00069 int levsteps = 0; 00070 int dsbits = mStepCnt ^ (mStepCnt-1); 00071 //errMsg("S"," step:"<<mStepCnt<<" s-1:"<<(mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits)); 00072 for(int lev=0; lev<=mMaxRefine; lev++) { 00073 //if(! (mStepCnt&(1<<lev)) ) { 00074 if( dsbits & (1<<(mMaxRefine-lev)) ) { 00075 //errMsg("S"," l"<<lev); 00076 00077 if(lev==mMaxRefine) { 00078 // always advance fine level... 00079 fineAdvance(); 00080 } else { 00081 adaptGrid(lev); 00082 coarseRestrictFromFine(lev); 00083 coarseAdvance(lev); 00084 } 00085 #if FSGR_OMEGA_DEBUG==1 00086 errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) ); 00087 mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0; 00088 #endif // FSGR_OMEGA_DEBUG==1 00089 levsteps++; 00090 } 00091 mCurrentMass += mLevel[lev].lmass; 00092 mCurrentVolume += mLevel[lev].lvolume; 00093 } 00094 00095 // prepare next step 00096 mStepCnt++; 00097 00098 00099 // some dbugging output follows 00100 // calculate MLSUPS 00101 myTime_t timeend = getTime(); 00102 00103 mNumUsedCells += mNumInterdCells; // count both types for MLSUPS 00104 mAvgNumUsedCells += mNumUsedCells; 00105 mMLSUPS = ((double)mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000.0); 00106 if(mMLSUPS>10000){ mMLSUPS = -1; } 00107 //else { mAvgMLSUPS += mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups 00108 00109 LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000); 00110 if(totMLSUPS>10000) totMLSUPS = -1; 00111 mNumInvIfTotal += mNumInvIfCells; // debug 00112 00113 // do some formatting 00114 if(!mSilent){ 00115 int avgcls = (int)(mAvgNumUsedCells/(LONGINT)mStepCnt); 00116 debMsgStd("LbmFsgrSolver::step", DM_MSG, mName<<" cnt:"<<mStepCnt<<" t:"<<mSimulationTime<< 00117 " cur-mlsups:"<<mMLSUPS<< //" avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< 00118 " totcls:"<<mNumUsedCells<< " avgcls:"<< avgcls<< 00119 " intd:"<<mNumInterdCells<< " invif:"<<mNumInvIfCells<< 00120 " invift:"<<mNumInvIfTotal<< " fsgrcs:"<<mNumFsgrChanges<< 00121 " filled:"<<mNumFilledCells<<", emptied:"<<mNumEmptiedCells<< 00122 " mMxv:"<<PRINT_VEC(mMxvx,mMxvy,mMxvz)<<", tscnts:"<<mTimeSwitchCounts<< 00123 //" RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */ 00124 " probs:"<<mNumProblems<< " simt:"<<mSimulationTime<< 00125 " took:"<< getTimeString(timeend-timestart)<< 00126 " for '"<<mName<<"' " , 10); 00127 } else { debMsgDirect("."); } 00128 00129 if(mStepCnt==1) { 00130 mMinNoCells = mMaxNoCells = mNumUsedCells; 00131 } else { 00132 if(mNumUsedCells>mMaxNoCells) mMaxNoCells = mNumUsedCells; 00133 if(mNumUsedCells<mMinNoCells) mMinNoCells = mNumUsedCells; 00134 } 00135 00136 // mass scale test 00137 if((mMaxRefine>0)&&(mInitialMass>0.0)) { 00138 LbmFloat mscale = mInitialMass/mCurrentMass; 00139 00140 mscale = 1.0; 00141 const LbmFloat dchh = 0.001; 00142 if(mCurrentMass<mInitialMass) mscale = 1.0+dchh; 00143 if(mCurrentMass>mInitialMass) mscale = 1.0-dchh; 00144 00145 // use mass rescaling? 00146 // with float precision this seems to be nonsense... 00147 const bool MREnable = false; 00148 00149 const int MSInter = 2; 00150 static int mscount = 0; 00151 if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){ 00152 // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37 00153 // mass rescale MASS RESCALE check 00154 errMsg("MDTDD","\n\n"); 00155 errMsg("MDTDD","FORCE RESCALE MASS! " 00156 <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass) 00157 <<" step:"<<mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" " 00158 ); 00159 errMsg("MDTDD","\n\n"); 00160 00161 mscount++; 00162 for(int lev=mMaxRefine; lev>=0 ; lev--) { 00163 //for(int workSet = 0; workSet<=1; workSet++) { 00164 int wss = 0; 00165 int wse = 1; 00166 #if COMPRESSGRIDS==1 00167 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr; 00168 #endif // COMPRESSGRIDS==1 00169 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT 00170 00171 FSGR_FORIJK1(lev) { 00172 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 00173 ) { 00174 00175 FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; } 00176 QCELL(lev, i,j,k,workSet, dMass) *= mscale; 00177 QCELL(lev, i,j,k,workSet, dFfrac) *= mscale; 00178 00179 } else { 00180 continue; 00181 } 00182 } 00183 } 00184 mLevel[lev].lmass *= mscale; 00185 } 00186 } 00187 00188 mCurrentMass *= mscale; 00189 }// if mass scale test */ 00190 else { 00191 // use current mass after full step for initial setting 00192 if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) { 00193 mInitialMass = mCurrentMass; 00194 debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2); 00195 } 00196 } 00197 00198 #if LBM_INCLUDE_TESTSOLVERS==1 00199 if((mUseTestdata)&&(mInitDone)) { handleTestdata(); } 00200 mrExchange(); 00201 #endif 00202 00203 // advance positions with current grid 00204 advanceParticles(); 00205 if(mpParticles) { 00206 mpParticles->checkDumpTextPositions(mSimulationTime); 00207 mpParticles->checkTrails(mSimulationTime); 00208 } 00209 00210 // one of the last things to do - adapt timestep 00211 // was in fineAdvance before... 00212 if(mTimeAdap) { 00213 adaptTimestep(); 00214 } // time adaptivity 00215 00216 00217 #ifndef WIN32 00218 // good indicator for instabilities 00219 if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; } 00220 if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; } 00221 #endif // WIN32 00222 00223 // output total step time 00224 myTime_t timeend2 = getTime(); 00225 double mdelta = (lastMass-mCurrentMass); 00226 if(ABS(mdelta)<1e-12) mdelta=0.; 00227 double effMLSUPS = ((double)mNumUsedCells / ((timeend2-timestart)/(double)1000.0) ) / (1000000.0); 00228 if(mInitDone) { 00229 if(effMLSUPS>10000){ effMLSUPS = -1; } 00230 else { mAvgMLSUPS += effMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups 00231 } 00232 00233 debMsgStd("LbmFsgrSolver::stepMain", DM_MSG, "mmpid:"<<glob_mpindex<<" step:"<<mStepCnt 00234 <<" dccd="<< mCurrentMass 00235 //<<",d"<<mdelta 00236 //<<",ds"<<(mCurrentMass-mObjectMassMovnd[1]) 00237 //<<"/"<<mCurrentVolume<<"(fix="<<mFixMass<<",ini="<<mInitialMass<<"), " 00238 <<" effMLSUPS=("<< effMLSUPS 00239 <<",avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), " 00240 <<" took totst:"<< getTimeString(timeend2-timestart), 3); 00241 // nicer output 00242 //debMsgDirect(std::endl); 00243 // */ 00244 00245 messageOutputForce(""); 00246 //#endif // ELBEEM_PLUGIN!=1 00247 } 00248 00249 #define NEWDEBCHECK(str) \ 00250 if(!this->mPanic){ FSGR_FORIJK_BOUNDS(mMaxRefine) { \ 00251 if(RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr)&(CFFluid|CFInter)) { \ 00252 for(int l=0;l<dTotalNum;l++) { \ 00253 if(!finite(QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr,l))) { errMsg("NNOFIN"," "<<str<<" at "<<PRINT_IJK<<" l"<<l<<" "); }\ 00254 }/*for*/ \ 00255 }/*if*/ \ 00256 } } 00257 00258 void LbmFsgrSolver::fineAdvance() 00259 { 00260 // do the real thing... 00261 //NEWDEBCHECK("t1"); 00262 mainLoop( mMaxRefine ); 00263 if(mUpdateFVHeight) { 00264 // warning assume -Y gravity... 00265 mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez)); 00266 if(mFVHeight<1.0) mFVHeight = 1.0; 00267 mpParam->setFluidVolumeHeight(mFVHeight); 00268 } 00269 00270 // advance time before timestep change 00271 mSimulationTime += mpParam->getTimestep(); 00272 // time adaptivity 00273 mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) ); 00274 //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 00275 if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); } 00276 00277 // update other set 00278 mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr; 00279 mLevel[mMaxRefine].setCurr ^= 1; 00280 mLevel[mMaxRefine].lsteps++; 00281 00282 // flag init... (work on current set, to simplify flag checks) 00283 reinitFlags( mLevel[mMaxRefine].setCurr ); 00284 if(!mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); } 00285 00286 // DEBUG VEL CHECK 00287 if(0) { 00288 int lev = mMaxRefine; 00289 int workSet = mLevel[lev].setCurr; 00290 int mi=0,mj=0,mk=0; 00291 LbmFloat compRho=0.; 00292 LbmFloat compUx=0.; 00293 LbmFloat compUy=0.; 00294 LbmFloat compUz=0.; 00295 LbmFloat maxUlen=0.; 00296 LbmVec maxU(0.); 00297 LbmFloat maxRho=-100.; 00298 int ri=0,rj=0,rk=0; 00299 00300 FSGR_FORIJK1(lev) { 00301 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) ) { 00302 compRho=QCELL(lev, i,j,k,workSet, dC); 00303 compUx = compUy = compUz = 0.0; 00304 for(int l=1; l<this->cDfNum; l++) { 00305 LbmFloat df = QCELL(lev, i,j,k,workSet, l); 00306 compRho += df; 00307 compUx += (this->dfDvecX[l]*df); 00308 compUy += (this->dfDvecY[l]*df); 00309 compUz += (this->dfDvecZ[l]*df); 00310 } 00311 LbmVec u(compUx,compUy,compUz); 00312 LbmFloat nu = norm(u); 00313 if(nu>maxUlen) { 00314 maxU = u; 00315 maxUlen = nu; 00316 mi=i; mj=j; mk=k; 00317 } 00318 if(compRho>maxRho) { 00319 maxRho=compRho; 00320 ri=i; rj=j; rk=k; 00321 } 00322 } else { 00323 continue; 00324 } 00325 } 00326 00327 errMsg("MAXVELCHECK"," at "<<PRINT_VEC(mi,mj,mk)<<" norm:"<<maxUlen<<" u:"<<maxU); 00328 errMsg("MAXRHOCHECK"," at "<<PRINT_VEC(ri,rj,rk)<<" rho:"<<maxRho); 00329 printLbmCell(lev, 30,36,23, -1); 00330 } // DEBUG VEL CHECK 00331 00332 } 00333 00334 00335 00336 // fine step defines 00337 00338 // access to own dfs during step (may be changed to local array) 00339 #define MYDF(l) RAC(ccel, l) 00340 00341 // drop model definitions 00342 #define RWVEL_THRESH 1.5 00343 #define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5) 00344 00345 #if LBMDIM==3 00346 // normal 00347 #define SLOWDOWNREGION (mSizez/4) 00348 #else // LBMDIM==2 00349 // off 00350 #define SLOWDOWNREGION 10 00351 #endif // LBMDIM==2 00352 #define P_LCSMQO 0.01 00353 00354 /*****************************************************************************/ 00356 /*****************************************************************************/ 00357 void 00358 LbmFsgrSolver::mainLoop(int lev) 00359 { 00360 // loops over _only inner_ cells ----------------------------------------------------------------------------------- 00361 00362 // slow surf regions smooth (if below) 00363 const LbmFloat smoothStrength = 0.0; //0.01; 00364 const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03; 00365 const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit; 00366 00367 const int cutMin = 1; 00368 const int cutConst = mCutoff+2; 00369 00370 00371 # if LBM_INCLUDE_TESTSOLVERS==1 00372 // 3d region off... quit 00373 if((mUseTestdata)&&(mpTest->mFarfMode>0)) { return; } 00374 # endif // ELBEEM_PLUGIN!=1 00375 00376 // main loop region 00377 const bool doReduce = true; 00378 const int gridLoopBound=1; 00379 GRID_REGION_INIT(); 00380 #if PARALLEL==1 00381 #pragma omp parallel default(shared) \ 00382 reduction(+: \ 00383 calcCurrentMass,calcCurrentVolume, \ 00384 calcCellsFilled,calcCellsEmptied, \ 00385 calcNumUsedCells ) 00386 GRID_REGION_START(); 00387 #else // PARALLEL==1 00388 GRID_REGION_START(); 00389 #endif // PARALLEL==1 00390 00391 // local to main 00392 CellFlagType nbflag[LBM_DFNUM]; 00393 int oldFlag, newFlag, nbored; 00394 LbmFloat m[LBM_DFNUM]; 00395 LbmFloat rho, ux, uy, uz, tmp, usqr; 00396 00397 // smago vars 00398 LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega; 00399 00400 // ifempty cell conversion flags 00401 bool iffilled, ifemptied; 00402 LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors 00403 int recons[LBM_DFNUM]; // reconstruct this DF? 00404 int numRecons; // how many are reconstructed? 00405 00406 LbmFloat mass=0., change=0., lcsmqo=0.; 00407 rho= ux= uy= uz= usqr= tmp= 0.; 00408 lcsmqadd = lcsmomega = 0.; 00409 FORDF0{ lcsmeq[l] = 0.; } 00410 00411 // --- 00412 // now stream etc. 00413 // use //template functions for 2D/3D 00414 00415 GRID_LOOP_START(); 00416 // loop start 00417 // stream from current set to other, then collide and store 00418 //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id); 00419 00420 # if FSGR_STRICT_DEBUG==1 00421 // safety pointer check 00422 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG 00423 if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 00424 (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) { 00425 errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<< 00426 RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<< 00427 (*pFlagSrc)<<","<<(*pFlagDst) <<", pointers "<< 00428 (long)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(long)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<< 00429 (long)(pFlagSrc)<<","<<(long)(pFlagDst)<<" " 00430 ); 00431 CAUSE_PANIC; 00432 } 00433 if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 00434 (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) { 00435 errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<< 00436 (long)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(long)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<< 00437 (long)(ccel)<<","<<(long)(tcel)<<" " 00438 ); 00439 CAUSE_PANIC; 00440 } 00441 # endif 00442 oldFlag = *pFlagSrc; 00443 00444 // old INTCFCOARSETEST==1 00445 if( (oldFlag & (CFGrFromCoarse)) ) { 00446 if(( mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) { 00447 FORDF0 { RAC(tcel,l) = RAC(ccel,l); } 00448 } else { 00449 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false); 00450 calcNumUsedCells++; 00451 } 00452 continue; // interpolateFineFromCoarse test! 00453 } // interpolateFineFromCoarse test! 00454 00455 if(oldFlag & (CFMbndInflow)) { 00456 // fluid & if are ok, fill if later on 00457 int isValid = oldFlag & (CFFluid|CFInter); 00458 const LbmFloat iniRho = 1.0; 00459 const int OId = oldFlag>>24; 00460 if(!isValid) { 00461 // make new if cell 00462 const LbmVec vel(mObjectSpeeds[OId]); 00463 // TODO add OPT3D treatment 00464 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); } 00465 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho; 00466 RAC(tcel, dFlux) = FLUX_INIT; 00467 changeFlag(lev, i,j,k, TSET(lev), CFInter); 00468 calcCurrentMass += iniRho; 00469 calcCurrentVolume += 1.0; 00470 calcNumUsedCells++; 00471 mInitialMass += iniRho; 00472 // dont treat cell until next step 00473 continue; 00474 } 00475 } 00476 else // these are exclusive 00477 if(oldFlag & (CFMbndOutflow)) { 00478 int isnotValid = oldFlag & (CFFluid); 00479 if(isnotValid) { 00480 // remove fluid cells, shouldnt be here anyway 00481 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; } 00482 mInitialMass -= fluidRho; 00483 const LbmFloat iniRho = 0.0; 00484 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho; 00485 RAC(tcel, dFlux) = FLUX_INIT; 00486 changeFlag(lev, i,j,k, TSET(lev), CFInter); 00487 00488 // same as ifemptied for if below 00489 LbmPoint oemptyp; oemptyp.flag = 0; 00490 oemptyp.x = i; oemptyp.y = j; oemptyp.z = k; 00491 LIST_EMPTY(oemptyp); 00492 calcCellsEmptied++; 00493 continue; 00494 } 00495 } 00496 00497 if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 00498 *pFlagDst = oldFlag; 00499 continue; 00500 } 00501 /*if( oldFlag & CFNoBndFluid ) { // TEST ME FASTER? 00502 OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK; 00503 RAC(tcel,dFfrac) = 1.0; 00504 *pFlagDst = (CellFlagType)oldFlag; // newFlag; 00505 calcCurrentMass += rho; calcCurrentVolume += 1.0; 00506 calcNumUsedCells++; 00507 continue; 00508 }// TEST ME FASTER? */ 00509 00510 // only neighbor flags! not own flag 00511 nbored = 0; 00512 00513 #if OPT3D==0 00514 FORDF1 { 00515 nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l); 00516 nbored |= nbflag[l]; 00517 } 00518 #else 00519 nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB]; 00520 nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB]; 00521 nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB]; 00522 nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB]; 00523 nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB]; 00524 00525 nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW]; 00526 nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS]; 00527 nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE]; 00528 00529 nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW]; 00530 nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE]; 00531 00532 nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW]; 00533 nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN]; 00534 nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE]; 00535 00536 nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST]; 00537 nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT]; 00538 nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT]; 00539 nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET]; 00540 nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT]; 00541 // */ 00542 #endif 00543 00544 // pointer to destination cell 00545 calcNumUsedCells++; 00546 00547 // FLUID cells 00548 if( oldFlag & CFFluid ) { 00549 // only standard fluid cells (with nothing except fluid as nbs 00550 00551 if(oldFlag&CFMbndInflow) { 00552 // force velocity for inflow, necessary to have constant direction of flow 00553 // FIXME , test also set interface cells? 00554 const int OId = oldFlag>>24; 00555 //? DEFAULT_STREAM; 00556 //const LbmFloat fluidRho = 1.0; 00557 // for submerged inflows, streaming would have to be performed... 00558 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; } 00559 const LbmVec vel(mObjectSpeeds[OId]); 00560 ux=vel[0], uy=vel[1], uz=vel[2]; 00561 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); 00562 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } 00563 rho = fluidRho; 00564 //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho); 00565 } else { 00566 if(nbored&CFBnd) { 00567 DEFAULT_STREAM; 00568 //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 00569 DEFAULT_COLLIDEG(mLevel[lev].gravity); 00570 oldFlag &= (~CFNoBndFluid); 00571 } else { 00572 // do standard stream/collide 00573 OPTIMIZED_STREAMCOLLIDE; 00574 oldFlag |= CFNoBndFluid; 00575 } 00576 } 00577 00578 PERFORM_USQRMAXCHECK; 00579 // "normal" fluid cells 00580 RAC(tcel,dFfrac) = 1.0; 00581 *pFlagDst = (CellFlagType)oldFlag; // newFlag; 00582 calcCurrentMass += rho; 00583 calcCurrentVolume += 1.0; 00584 continue; 00585 } 00586 00587 newFlag = oldFlag; 00588 // make sure here: always check which flags to really unset...! 00589 newFlag = newFlag & (~( 00590 CFNoNbFluid|CFNoNbEmpty| CFNoDelete 00591 | CFNoInterpolSrc 00592 | CFNoBndFluid 00593 )); 00594 if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS 00595 newFlag |= CFNoBndFluid; 00596 } 00597 /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS 00598 // explicitly check for noslip neighbors 00599 bool hasnoslipnb = false; 00600 FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; } 00601 if(!hasnoslipnb) newFlag |= CFNoBndFluid; 00602 } // */ 00603 00604 // store own dfs and mass 00605 mass = RAC(ccel,dMass); 00606 00607 // WARNING - only interface cells arrive here! 00608 // read distribution funtions of adjacent cells = stream step 00609 DEFAULT_STREAM; 00610 00611 if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; } 00612 if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; } 00613 00614 // calculate mass exchange for interface cells 00615 LbmFloat myfrac = RAC(ccel,dFfrac); 00616 if(myfrac<0.) myfrac=0.; // NEWSURFT 00617 # define nbdf(l) m[ this->dfInv[(l)] ] 00618 00619 // update mass 00620 // only do boundaries for fluid cells, and interface cells without 00621 // any fluid neighbors (assume that interface cells _with_ fluid 00622 // neighbors are affected enough by these) 00623 // which Df's have to be reconstructed? 00624 // for fluid cells - just the f_i difference from streaming to empty cells ---- 00625 numRecons = 0; 00626 bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip)); 00627 //onlyBndnb = false; // DEBUG test off 00628 00629 FORDF1 { // dfl loop 00630 recons[l] = 0; 00631 nbfracs[l] = 0.0; 00632 // finally, "normal" interface cells ---- 00633 if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!! 00634 change = nbdf(l) - MYDF(l); 00635 } 00636 // interface cells - distuingish cells that shouldn't fill/empty 00637 else if( nbflag[l] & CFInter ) { 00638 00639 LbmFloat mynbfac,nbnbfac; 00640 // NEW TEST t1 00641 // t2 -> off 00642 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) { 00643 mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux); 00644 nbnbfac = 1.0/mynbfac; 00645 onlyBndnb = false; 00646 } else { 00647 mynbfac = nbnbfac = 1.0; // switch calc flux off 00648 goto changeDefault; // NEWSURFT 00649 //change = 0.; goto changeDone; // NEWSURFT 00650 } 00651 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3 00652 00653 // perform interface case handling 00654 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) { 00655 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) { 00656 case 0: 00657 // we are a normal cell so... 00658 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) { 00659 case CFNoNbFluid: 00660 // just fill current cell = empty neighbor 00661 change = nbnbfac*nbdf(l) ; goto changeDone; 00662 case CFNoNbEmpty: 00663 // just empty current cell = fill neighbor 00664 change = - mynbfac*MYDF(l) ; goto changeDone; 00665 } 00666 break; 00667 00668 case CFNoNbFluid: 00669 // we dont have fluid nb's so... 00670 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) { 00671 case 0: 00672 case CFNoNbEmpty: 00673 // we have no fluid nb's -> just empty 00674 change = - mynbfac*MYDF(l) ; goto changeDone; 00675 } 00676 break; 00677 00678 case CFNoNbEmpty: 00679 // we dont have empty nb's so... 00680 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) { 00681 case 0: 00682 case CFNoNbFluid: 00683 // we have no empty nb's -> just fill 00684 change = nbnbfac*nbdf(l); goto changeDone; 00685 } 00686 break; 00687 }} // inter-inter exchange 00688 00689 changeDefault: ; 00690 // just do normal mass exchange... 00691 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ; 00692 changeDone: ; 00693 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac); 00694 if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT 00695 change *= (myfrac + nbfracs[l]) * 0.5; 00696 } // the other cell is interface 00697 00698 // last alternative - reconstruction in this direction 00699 else { 00700 // empty + bnd case 00701 recons[l] = 1; 00702 numRecons++; 00703 change = 0.0; 00704 } 00705 00706 // modify mass at SRCS 00707 mass += change; 00708 } // l 00709 // normal interface, no if empty/fluid 00710 00711 // computenormal 00712 LbmFloat surfaceNormal[3]; 00713 computeFluidSurfaceNormal(ccel,pFlagSrc, surfaceNormal); 00714 00715 if( (ABS(surfaceNormal[0])+ABS(surfaceNormal[1])+ABS(surfaceNormal[2])) > LBM_EPSILON) { 00716 // normal ok and usable... 00717 FORDF1 { 00718 if( (this->dfDvecX[l]*surfaceNormal[0] + this->dfDvecY[l]*surfaceNormal[1] + this->dfDvecZ[l]*surfaceNormal[2]) // dot Dvec,norml 00719 > LBM_EPSILON) { 00720 recons[l] = 2; 00721 numRecons++; 00722 } 00723 } 00724 } 00725 00726 // calculate macroscopic cell values 00727 LbmFloat oldUx, oldUy, oldUz; 00728 LbmFloat oldRho; // OLD rho = ccel->rho; 00729 # define REFERENCE_PRESSURE 1.0 // always atmosphere... 00730 # if OPT3D==0 00731 oldRho=RAC(ccel,0); 00732 oldUx = oldUy = oldUz = 0.0; 00733 for(int l=1; l<this->cDfNum; l++) { 00734 oldRho += RAC(ccel,l); 00735 oldUx += (this->dfDvecX[l]*RAC(ccel,l)); 00736 oldUy += (this->dfDvecY[l]*RAC(ccel,l)); 00737 oldUz += (this->dfDvecZ[l]*RAC(ccel,l)); 00738 } 00739 // reconstruct dist funcs from empty cells 00740 FORDF1 { 00741 if(recons[ l ]) { 00742 m[ this->dfInv[l] ] = 00743 this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 00744 this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 00745 - MYDF( l ); 00746 } 00747 } 00748 ux=oldUx, uy=oldUy, uz=oldUz; // no local vars, only for usqr 00749 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on 00750 # else // OPT3D==0 00751 oldRho = + RAC(ccel,dC) + RAC(ccel,dN ) 00752 + RAC(ccel,dS ) + RAC(ccel,dE ) 00753 + RAC(ccel,dW ) + RAC(ccel,dT ) 00754 + RAC(ccel,dB ) + RAC(ccel,dNE) 00755 + RAC(ccel,dNW) + RAC(ccel,dSE) 00756 + RAC(ccel,dSW) + RAC(ccel,dNT) 00757 + RAC(ccel,dNB) + RAC(ccel,dST) 00758 + RAC(ccel,dSB) + RAC(ccel,dET) 00759 + RAC(ccel,dEB) + RAC(ccel,dWT) 00760 + RAC(ccel,dWB); 00761 00762 oldUx = + RAC(ccel,dE) - RAC(ccel,dW) 00763 + RAC(ccel,dNE) - RAC(ccel,dNW) 00764 + RAC(ccel,dSE) - RAC(ccel,dSW) 00765 + RAC(ccel,dET) + RAC(ccel,dEB) 00766 - RAC(ccel,dWT) - RAC(ccel,dWB); 00767 00768 oldUy = + RAC(ccel,dN) - RAC(ccel,dS) 00769 + RAC(ccel,dNE) + RAC(ccel,dNW) 00770 - RAC(ccel,dSE) - RAC(ccel,dSW) 00771 + RAC(ccel,dNT) + RAC(ccel,dNB) 00772 - RAC(ccel,dST) - RAC(ccel,dSB); 00773 00774 oldUz = + RAC(ccel,dT) - RAC(ccel,dB) 00775 + RAC(ccel,dNT) - RAC(ccel,dNB) 00776 + RAC(ccel,dST) - RAC(ccel,dSB) 00777 + RAC(ccel,dET) - RAC(ccel,dEB) 00778 + RAC(ccel,dWT) - RAC(ccel,dWB); 00779 00780 // now reconstruction 00781 ux=oldUx, uy=oldUy, uz=oldUz; // no local vars, only for usqr 00782 rho = REFERENCE_PRESSURE; 00783 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on 00784 if(recons[dN ]) { m[dS ] = EQN + EQS - MYDF(dN ); } 00785 if(recons[dS ]) { m[dN ] = EQS + EQN - MYDF(dS ); } 00786 if(recons[dE ]) { m[dW ] = EQE + EQW - MYDF(dE ); } 00787 if(recons[dW ]) { m[dE ] = EQW + EQE - MYDF(dW ); } 00788 if(recons[dT ]) { m[dB ] = EQT + EQB - MYDF(dT ); } 00789 if(recons[dB ]) { m[dT ] = EQB + EQT - MYDF(dB ); } 00790 if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); } 00791 if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); } 00792 if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); } 00793 if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); } 00794 if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); } 00795 if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); } 00796 if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); } 00797 if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); } 00798 if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); } 00799 if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); } 00800 if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); } 00801 if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); } 00802 # endif 00803 00804 00805 // inflow bc handling 00806 if(oldFlag & (CFMbndInflow)) { 00807 // fill if cells in inflow region 00808 if(myfrac<0.5) { 00809 mass += 0.25; 00810 mInitialMass += 0.25; 00811 } 00812 const int OId = oldFlag>>24; 00813 const LbmVec vel(mObjectSpeeds[OId]); 00814 ux=vel[0], uy=vel[1], uz=vel[2]; 00815 //? usqr = 1.5 * (ux*ux + uy*uy + uz*uz); 00816 //FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } rho = fluidRho; 00817 rho = REFERENCE_PRESSURE; 00818 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); } 00819 //errMsg("INFLOW_DEBUG","if at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho); 00820 } else { 00821 // NEWSURFT, todo optimize! 00822 if(onlyBndnb) { //if(usqr<0.001*0.001) { 00823 rho = ux = uy = uz = 0.; 00824 FORDF0{ 00825 rho += m[l]; 00826 ux += (this->dfDvecX[l]*m[l]); 00827 uy += (this->dfDvecY[l]*m[l]); 00828 uz += (this->dfDvecZ[l]*m[l]); 00829 } 00830 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); } 00831 } else {// NEWSURFT */ 00832 if(usqr>0.3*0.3) { 00833 // force reset! , warning causes distortions... 00834 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,0.,0.,0.); } 00835 } else { 00836 // normal collide 00837 // mass streaming done... do normal collide 00838 LbmVec grav = mLevel[lev].gravity*mass; 00839 DEFAULT_COLLIDEG(grav); 00840 PERFORM_USQRMAXCHECK; } 00841 // rho init from default collide necessary for fill/empty check below 00842 } // test 00843 } 00844 00845 // testing..., particle generation 00846 // also check oldFlag for CFNoNbFluid, test 00847 // for inflow no pargen test // NOBUBBB! 00848 if((mInitDone) 00849 // dont allow new if cells, or submerged ones 00850 && (!((oldFlag|newFlag)& (CFNoDelete|CFNoNbEmpty) )) 00851 // dont try to subtract from empty cells 00852 && (mass>0.) && (mPartGenProb>0.0)) { 00853 bool doAdd = true; 00854 bool bndOk=true; 00855 if( (i<cutMin)||(i>mSizex-cutMin)|| 00856 (j<cutMin)||(j>mSizey-cutMin)|| 00857 (k<cutMin)||(k>mSizez-cutMin) ) { bndOk=false; } 00858 if(!bndOk) doAdd=false; 00859 00860 LbmFloat prob = (rand()/(RAND_MAX+1.0)); 00861 LbmFloat basethresh = mPartGenProb*lcsmqo*(LbmFloat)(mSizez+mSizey+mSizex)*0.5*0.333; 00862 00863 // physical drop model 00864 if(mPartUsePhysModel) { 00865 LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep); 00866 LbmFloat rux = (ux * realWorldFac); 00867 LbmFloat ruy = (uy * realWorldFac); 00868 LbmFloat ruz = (uz * realWorldFac); 00869 LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz)); 00870 basethresh *= rl; 00871 00872 // reduce probability in outer region? 00873 const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst; 00874 const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst; 00875 LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord); 00876 LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord); 00877 if(pifac<0.) pifac=0.; 00878 if(pjfac<0.) pjfac=0.; 00879 00880 //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) { 00881 if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) { 00882 // add 00883 } else { 00884 doAdd = false; // dont... 00885 } 00886 00887 // "wind" disturbance 00888 // use realworld relative velocity here instead? 00889 if( (doAdd && 00890 ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks 00891 ||(k>mSizez-SLOWDOWNREGION) ) { 00892 LbmFloat nuz = uz; 00893 if(k>mSizez-SLOWDOWNREGION) { 00894 // special case 00895 LbmFloat zfac = (LbmFloat)( k-(mSizez-SLOWDOWNREGION) ); 00896 zfac /= (LbmFloat)(SLOWDOWNREGION); 00897 nuz += (1.0) * zfac; // check max speed? OFF? 00898 //errMsg("TOPT"," at "<<PRINT_IJK<<" zfac"<<zfac); 00899 } else { 00900 // normal probability 00901 //? LbmFloat fac = P_LCSMQO-lcsmqo; 00902 //? jdf *= fac; 00903 } 00904 FORDF1 { 00905 const LbmFloat jdf = 0.05 * (rand()/(RAND_MAX+1.0)); 00906 // TODO use wind velocity? 00907 if(jdf>0.025) { 00908 const LbmFloat add = this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf; 00909 RAC(tcel,l) += add; } 00910 } 00911 //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) ); 00912 } // wind disturbance 00913 } // mPartUsePhysModel 00914 else { 00915 // empirical model 00916 //if((prob<basethresh) && (lcsmqo>0.0095)) { // add 00917 if((prob<basethresh) && (lcsmqo>0.012)) { // add 00918 } else { doAdd = false; }// dont... 00919 } 00920 00921 00922 // remove noise 00923 if(usqr<0.0001) doAdd=false; // TODO check!? 00924 00925 // dont try to subtract from empty cells 00926 // ensure cell has enough mass for new drop 00927 LbmFloat newPartsize = 1.0; 00928 if(mPartUsePhysModel) { 00929 // 1-10 00930 newPartsize += 9.0* (rand()/(RAND_MAX+1.0)); 00931 } else { 00932 // 1-5, overall size has to be less than 00933 // .62 (ca. 0.5) to make drops significantly smaller 00934 // than a full cell! 00935 newPartsize += 4.0* (rand()/(RAND_MAX+1.0)); 00936 } 00937 LbmFloat dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize); //PARTMASS(mPartDropMassSub*newPartsize); // mass: 4/3 pi r^3 rho 00938 while(dropmass>mass) { 00939 newPartsize -= 0.2; 00940 dropmass = ParticleObject::getMass(mPartDropMassSub*newPartsize); 00941 } 00942 if(newPartsize<=1.) doAdd=false; 00943 00944 if( (doAdd) ) { // init new particle 00945 for(int s=0; s<1; s++) { // one part! 00946 const LbmFloat posjitter = 0.05; 00947 const LbmFloat posjitteroffs = posjitter*-0.5; 00948 LbmFloat jpx = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0)); 00949 LbmFloat jpy = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0)); 00950 LbmFloat jpz = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0)); 00951 00952 const LbmFloat jitterstr = 1.0; 00953 const LbmFloat jitteroffs = jitterstr*-0.5; 00954 LbmFloat jx = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0)); 00955 LbmFloat jy = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0)); 00956 LbmFloat jz = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0)); 00957 00958 // average normal & velocity 00959 // -> mostly along velocity dir, many into surface 00960 // fluid velocity (not normalized!) 00961 LbmVec flvelVel = LbmVec(ux,uy,uz); 00962 LbmFloat flvelLen = norm(flvelVel); 00963 // surface normal 00964 LbmVec normVel = LbmVec(surfaceNormal[0],surfaceNormal[1],surfaceNormal[2]); 00965 normalize(normVel); 00966 LbmFloat normScale = (0.01+flvelLen); 00967 // jitter vector, 0.2 * flvel 00968 LbmVec jittVel = LbmVec(jx,jy,jz)*(0.05+flvelLen)*0.1; 00969 // weighten velocities 00970 const LbmFloat flvelWeight = 0.9; 00971 LbmVec newpartVel = normVel*normScale*(1.-flvelWeight) + flvelVel*(flvelWeight) + jittVel; 00972 00973 // offset towards surface (hide popping) 00974 jpx += -normVel[0]*0.4; 00975 jpy += -normVel[1]*0.4; 00976 jpz += -normVel[2]*0.4; 00977 00978 LbmFloat srci=i+0.5+jpx, srcj=j+0.5+jpy, srck=k+0.5+jpz; 00979 int type=0; 00980 type = PART_DROP; 00981 00982 # if LBMDIM==2 00983 newpartVel[2]=0.; srck=0.5; 00984 # endif // LBMDIM==2 00985 // subtract drop mass 00986 mass -= dropmass; 00987 // init new particle 00988 { 00989 ParticleObject np( ntlVec3Gfx(srci,srcj,srck) ); 00990 np.setVel(newpartVel[0],newpartVel[1],newpartVel[2]); 00991 np.setStatus(PART_IN); 00992 np.setType(type); 00993 //if((s%3)==2) np.setType(PART_FLOAT); 00994 np.setSize(newPartsize); 00995 //errMsg("NEWPART"," at "<<PRINT_IJK<<" u="<<norm(LbmVec(ux,uy,uz)) <<" add"<<doAdd<<" pvel="<<norm(newpartVel)<<" size="<<newPartsize ); 00996 //errMsg("NEWPT","u="<<newpartVel<<" norm="<<normVel<<" flvel="<<flvelVel<<" jitt="<<jittVel ); 00997 FSGR_ADDPART(np); 00998 } // new part 00999 //errMsg("PIT","NEW pit"<<np.getId()<<" pos:"<< np.getPos()<<" status:"<<convertFlags2String(np.getFlags())<<" vel:"<< np.getVel() ); 01000 } // multiple parts 01001 } // doAdd 01002 } // */ 01003 01004 01005 // interface cell filled or emptied? 01006 iffilled = ifemptied = 0; 01007 // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags 01008 // interface cell if full? 01009 if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; } 01010 // interface cell if empty? 01011 if( (mass) <= (rho * ( -FSGR_MAGICNR)) ) { ifemptied = 1; } 01012 01013 if(oldFlag & (CFMbndOutflow)) { 01014 mInitialMass -= mass; 01015 mass = myfrac = 0.0; 01016 iffilled = 0; ifemptied = 1; 01017 } 01018 01019 // looks much nicer... LISTTRICK 01020 # if FSGR_LISTTRICK==1 01021 //if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; } 01022 if(newFlag&CFNoBndFluid) { // test NEW TEST 01023 if(!iffilled) { 01024 // remove cells independent from amount of change... 01025 if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&& 01026 ( (mass>(rho*FSGR_LISTTTHRESHFULL)) || ((nbored&CFInter)==0) )) { 01027 //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); }; 01028 iffilled = 1; 01029 } 01030 } 01031 if(!ifemptied) { 01032 if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&& 01033 ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0) )) { 01034 //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); }; 01035 ifemptied = 1; 01036 } 01037 } 01038 } // nobndfluid only */ 01039 # endif 01040 //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!! 01041 01042 01043 // now that all dfs are known, handle last changes 01044 if(iffilled) { 01045 LbmPoint filledp; filledp.flag=0; 01046 if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1; // NEWSURFT 01047 filledp.x = i; filledp.y = j; filledp.z = k; 01048 LIST_FULL(filledp); 01049 //mNumFilledCells++; // DEBUG 01050 calcCellsFilled++; 01051 } 01052 else if(ifemptied) { 01053 LbmPoint emptyp; emptyp.flag=0; 01054 if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; // NEWSURFT 01055 emptyp.x = i; emptyp.y = j; emptyp.z = k; 01056 LIST_EMPTY(emptyp); 01057 //mNumEmptiedCells++; // DEBUG 01058 calcCellsEmptied++; 01059 } 01060 // dont cutoff values -> better cell conversions 01061 RAC(tcel,dFfrac) = (mass/rho); 01062 01063 // init new flux value 01064 float flux = FLUX_INIT; // dxqn on 01065 if(newFlag&CFNoBndFluid) { 01066 //flux = 50.0; // extreme on 01067 for(int nn=1; nn<this->cDfNum; nn++) { 01068 if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += this->dfLength[nn]; } 01069 } 01070 // optical hack - smooth slow moving 01071 // surface regions 01072 if(usqr< sssUsqrLimit) { 01073 for(int nn=1; nn<this->cDfNum; nn++) { 01074 if(nbfracs[nn]!=0.0) { 01075 LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv; 01076 if(curSmooth>1.0) curSmooth=1.0; 01077 flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ; 01078 } 01079 } } 01080 // NEW TEST */ 01081 } 01082 // flux = FLUX_INIT; // calc flux off 01083 QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */ 01084 01085 // perform mass exchange with streamed values 01086 QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST 01087 // set new flag 01088 *pFlagDst = (CellFlagType)newFlag; 01089 calcCurrentMass += mass; 01090 calcCurrentVolume += RAC(tcel,dFfrac); 01091 01092 // interface cell handling done... 01093 01094 #if PARALLEL!=1 01095 GRID_LOOPREG_END(); 01096 #else // PARALLEL==1 01097 #include "paraloopend.h" // = GRID_LOOPREG_END(); 01098 #endif // PARALLEL==1 01099 01100 // write vars from computations to class 01101 mLevel[lev].lmass = calcCurrentMass; 01102 mLevel[lev].lvolume = calcCurrentVolume; 01103 mNumFilledCells = calcCellsFilled; 01104 mNumEmptiedCells = calcCellsEmptied; 01105 mNumUsedCells = calcNumUsedCells; 01106 } 01107 01108 01109 01110 void 01111 LbmFsgrSolver::preinitGrids() 01112 { 01113 const int lev = mMaxRefine; 01114 const bool doReduce = false; 01115 const int gridLoopBound=0; 01116 01117 // preinit both grids 01118 for(int s=0; s<2; s++) { 01119 01120 GRID_REGION_INIT(); 01121 #if PARALLEL==1 01122 #pragma omp parallel default(shared) \ 01123 reduction(+: \ 01124 calcCurrentMass,calcCurrentVolume, \ 01125 calcCellsFilled,calcCellsEmptied, \ 01126 calcNumUsedCells ) 01127 #endif // PARALLEL==1 01128 GRID_REGION_START(); 01129 GRID_LOOP_START(); 01130 for(int l=0; l<dTotalNum; l++) { RAC(ccel,l) = 0.; } 01131 *pFlagSrc =0; 01132 *pFlagDst =0; 01133 //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id); 01134 #if PARALLEL!=1 01135 GRID_LOOPREG_END(); 01136 #else // PARALLEL==1 01137 #include "paraloopend.h" // = GRID_LOOPREG_END(); 01138 #endif // PARALLEL==1 01139 01140 /* dummy, remove warnings */ 01141 calcCurrentMass = calcCurrentVolume = 0.; 01142 calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0; 01143 01144 // change grid 01145 mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr; 01146 mLevel[mMaxRefine].setCurr ^= 1; 01147 } 01148 } 01149 01150 void 01151 LbmFsgrSolver::standingFluidPreinit() 01152 { 01153 const int lev = mMaxRefine; 01154 const bool doReduce = false; 01155 const int gridLoopBound=1; 01156 01157 GRID_REGION_INIT(); 01158 #if PARALLEL==1 01159 #pragma omp parallel default(shared) \ 01160 reduction(+: \ 01161 calcCurrentMass,calcCurrentVolume, \ 01162 calcCellsFilled,calcCellsEmptied, \ 01163 calcNumUsedCells ) 01164 #endif // PARALLEL==1 01165 GRID_REGION_START(); 01166 01167 LbmFloat rho, ux,uy,uz, usqr; 01168 CellFlagType nbflag[LBM_DFNUM]; 01169 LbmFloat m[LBM_DFNUM]; 01170 LbmFloat lcsmqo; 01171 # if OPT3D==1 01172 LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega; 01173 CellFlagType nbored=0; 01174 # endif // OPT3D==true 01175 01176 GRID_LOOP_START(); 01177 //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id); 01178 const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet); 01179 if( (currFlag & (CFEmpty|CFBnd)) ) continue; 01180 01181 if( (currFlag & (CFInter)) ) { 01182 // copy all values 01183 for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); } 01184 continue; 01185 } 01186 01187 if( (currFlag & CFNoBndFluid)) { 01188 OPTIMIZED_STREAMCOLLIDE; 01189 } else { 01190 FORDF1 { 01191 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l); 01192 } 01193 DEFAULT_STREAM; 01194 DEFAULT_COLLIDEG(mLevel[lev].gravity); 01195 } 01196 for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); } 01197 #if PARALLEL!=1 01198 GRID_LOOPREG_END(); 01199 #else // PARALLEL==1 01200 #include "paraloopend.h" // = GRID_LOOPREG_END(); 01201 #endif // PARALLEL==1 01202 01203 /* dummy remove warnings */ 01204 calcCurrentMass = calcCurrentVolume = 0.; 01205 calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0; 01206 01207 // change grid 01208 mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr; 01209 mLevel[mMaxRefine].setCurr ^= 1; 01210 } 01211 01212 01213 /****************************************************************************** 01214 * work on lists from updateCellMass to reinit cell flags 01215 *****************************************************************************/ 01216 01217 LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) { 01218 //return 0.0; // test 01219 int level = mMaxRefine; 01220 LbmFloat *ccel = RACPNT(level, i,j,k, workSet); 01221 01222 // computenormal 01223 CellFlagType *cflag = &RFLAG(level, i,j,k, workSet); 01224 LbmFloat n[3]; 01225 computeFluidSurfaceNormal(ccel,cflag, n); 01226 LbmFloat scal = mDvecNrm[l][0]*n[0] + mDvecNrm[l][1]*n[1] + mDvecNrm[l][2]*n[2]; 01227 01228 LbmFloat ret = 1.0; 01229 // forward direction, add mass (for filling cells): 01230 if(dirForw) { 01231 if(scal<LBM_EPSILON) ret = 0.0; 01232 else ret = scal; 01233 } else { 01234 // backward for emptying 01235 if(scal>-LBM_EPSILON) ret = 0.0; 01236 else ret = scal * -1.0; 01237 } 01238 //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB 01239 return ret; 01240 } 01241 01242 // warning - normal compuations are without 01243 // boundary checks & 01244 // normalization 01245 void LbmFsgrSolver::computeFluidSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt,LbmFloat *snret) { 01246 const int level = mMaxRefine; 01247 LbmFloat nx,ny,nz, nv1,nv2; 01248 const CellFlagType flagE = *(cflagpnt+1); 01249 const CellFlagType flagW = *(cflagpnt-1); 01250 if(flagE &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } 01251 else if(flagE &(CFBnd)){ nv1 = 1.; } 01252 else nv1 = 0.0; 01253 if(flagW &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } 01254 else if(flagW &(CFBnd)){ nv2 = 1.; } 01255 else nv2 = 0.0; 01256 nx = 0.5* (nv2-nv1); 01257 01258 const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx); 01259 const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx); 01260 if(flagN &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 01261 else if(flagN &(CFBnd)){ nv1 = 1.; } 01262 else nv1 = 0.0; 01263 if(flagS &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } 01264 else if(flagS &(CFBnd)){ nv2 = 1.; } 01265 else nv2 = 0.0; 01266 ny = 0.5* (nv2-nv1); 01267 01268 #if LBMDIM==3 01269 const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy); 01270 const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy); 01271 if(flagT &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 01272 else if(flagT &(CFBnd)){ nv1 = 1.; } 01273 else nv1 = 0.0; 01274 if(flagB &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } 01275 else if(flagB &(CFBnd)){ nv2 = 1.; } 01276 else nv2 = 0.0; 01277 nz = 0.5* (nv2-nv1); 01278 #else //LBMDIM==3 01279 nz = 0.0; 01280 #endif //LBMDIM==3 01281 01282 // return vals 01283 snret[0]=nx; snret[1]=ny; snret[2]=nz; 01284 } 01285 void LbmFsgrSolver::computeFluidSurfaceNormalAcc(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) { 01286 LbmFloat nx=0.,ny=0.,nz=0.; 01287 ccel = NULL; cflagpnt=NULL; // remove warning 01288 snret[0]=nx; snret[1]=ny; snret[2]=nz; 01289 } 01290 void LbmFsgrSolver::computeObstacleSurfaceNormal(LbmFloat *ccel, CellFlagType *cflagpnt, LbmFloat *snret) { 01291 const int level = mMaxRefine; 01292 LbmFloat nx,ny,nz, nv1,nv2; 01293 ccel = NULL; // remove warning 01294 01295 const CellFlagType flagE = *(cflagpnt+1); 01296 const CellFlagType flagW = *(cflagpnt-1); 01297 if(flagE &(CFBnd)){ nv1 = 1.; } 01298 else nv1 = 0.0; 01299 if(flagW &(CFBnd)){ nv2 = 1.; } 01300 else nv2 = 0.0; 01301 nx = 0.5* (nv2-nv1); 01302 01303 const CellFlagType flagN = *(cflagpnt+mLevel[level].lOffsx); 01304 const CellFlagType flagS = *(cflagpnt-mLevel[level].lOffsx); 01305 if(flagN &(CFBnd)){ nv1 = 1.; } 01306 else nv1 = 0.0; 01307 if(flagS &(CFBnd)){ nv2 = 1.; } 01308 else nv2 = 0.0; 01309 ny = 0.5* (nv2-nv1); 01310 01311 #if LBMDIM==3 01312 const CellFlagType flagT = *(cflagpnt+mLevel[level].lOffsy); 01313 const CellFlagType flagB = *(cflagpnt-mLevel[level].lOffsy); 01314 if(flagT &(CFBnd)){ nv1 = 1.; } 01315 else nv1 = 0.0; 01316 if(flagB &(CFBnd)){ nv2 = 1.; } 01317 else nv2 = 0.0; 01318 nz = 0.5* (nv2-nv1); 01319 #else //LBMDIM==3 01320 nz = 0.0; 01321 #endif //LBMDIM==3 01322 01323 // return vals 01324 snret[0]=nx; snret[1]=ny; snret[2]=nz; 01325 } 01326 void LbmFsgrSolver::computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret) { 01327 bool nonorm = false; 01328 LbmFloat nx=0.,ny=0.,nz=0.; 01329 if(i<=0) { nx = 1.; nonorm = true; } 01330 if(i>=mSizex-1) { nx = -1.; nonorm = true; } 01331 if(j<=0) { ny = 1.; nonorm = true; } 01332 if(j>=mSizey-1) { ny = -1.; nonorm = true; } 01333 # if LBMDIM==3 01334 if(k<=0) { nz = 1.; nonorm = true; } 01335 if(k>=mSizez-1) { nz = -1.; nonorm = true; } 01336 # endif // LBMDIM==3 01337 if(!nonorm) { 01338 // in domain, revert to helper, use setCurr&mMaxRefine 01339 LbmVec bnormal; 01340 CellFlagType *bflag = &RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr); 01341 LbmFloat *bcell = RACPNT(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr); 01342 computeObstacleSurfaceNormal(bcell,bflag, &bnormal[0]); 01343 // TODO check if there is a normal near here? 01344 // use wider range otherwise... 01345 snret[0]=bnormal[0]; snret[1]=bnormal[0]; snret[2]=bnormal[0]; 01346 return; 01347 } 01348 snret[0]=nx; snret[1]=ny; snret[2]=nz; 01349 } 01350 01351 void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) { 01352 #if FSGR_STRICT_DEBUG==10 01353 // dangerous, this can change the simulation... 01354 /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin(); 01355 iter != mListNewInter.end(); iter++ ) { 01356 if(ni!=iter->x) continue; 01357 if(nj!=iter->y) continue; 01358 if(nk!=iter->z) continue; 01359 // all 3 values match... skip point 01360 return; 01361 } */ 01362 #endif // FSGR_STRICT_DEBUG==1 01363 // store point 01364 LbmPoint newinter; newinter.flag = 0; 01365 newinter.x = ni; newinter.y = nj; newinter.z = nk; 01366 mListNewInter.push_back(newinter); 01367 } 01368 01369 void LbmFsgrSolver::reinitFlags( int workSet ) { 01370 // reinitCellFlags OLD mods: 01371 // add all to intel list? 01372 // check ffrac for new cells 01373 // new if cell inits (last loop) 01374 // vweights handling 01375 01376 const int debugFlagreinit = 0; 01377 01378 // some things need to be read/modified on the other set 01379 int otherSet = (workSet^1); 01380 // fixed level on which to perform 01381 int workLev = mMaxRefine; 01382 01383 /* modify interface cells from lists */ 01384 /* mark filled interface cells as fluid, or emptied as empty */ 01385 /* count neighbors and distribute excess mass to interface neighbor cells 01386 * problems arise when there are no interface neighbors anymore 01387 * then just distribute to any fluid neighbors... 01388 */ 01389 01390 // for symmetry, first init all neighbor cells */ 01391 for( vector<LbmPoint>::iterator iter=mListFull.begin(); 01392 iter != mListFull.end(); iter++ ) { 01393 int i=iter->x, j=iter->y, k=iter->z; 01394 if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM 01395 FORDF1 { 01396 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 01397 //if((LBMDIM>2)&&( (ni<=0) || (nj<=0) || (nk<=0) || (ni>=mLevel[workLev].lSizex-1) || (nj>=mLevel[workLev].lSizey-1) || (nk>=mLevel[workLev].lSizez-1) )) { 01398 if( (ni<=0) || (nj<=0) || 01399 (ni>=mLevel[workLev].lSizex-1) || 01400 (nj>=mLevel[workLev].lSizey-1) 01401 # if LBMDIM==3 01402 || (nk<=0) || (nk>=mLevel[workLev].lSizez-1) 01403 # endif // LBMDIM==1 01404 ) { 01405 continue; } // new bc, dont treat cells on boundary NEWBC 01406 if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){ 01407 01408 // preinit speed, get from average surrounding cells 01409 // interpolate from non-workset to workset, sets are handled in function 01410 01411 // new and empty interface cell, dont change old flag here! 01412 addToNewInterList(ni,nj,nk); 01413 01414 LbmFloat avgrho = 0.0; 01415 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0; 01416 interpolateCellValues(workLev,ni,nj,nk,workSet, avgrho,avgux,avguy,avguz); 01417 01418 // careful with l's... 01419 FORDF0M { 01420 QCELL(workLev,ni,nj,nk, workSet, m) = this->getCollideEq( m,avgrho, avgux, avguy, avguz ); 01421 //QCELL(workLev,ni,nj,nk, workSet, l) = avgnbdf[l]; // CHECK FIXME test? 01422 } 01423 //errMsg("FNEW", PRINT_VEC(ni,nj,nk)<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<<avgrho<<" vel"<<PRINT_VEC(avgux,avguy,avguz) ); // DEBUG SYMM 01424 QCELL(workLev,ni,nj,nk, workSet, dMass) = 0.0; //?? new 01425 QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 0.0; //?? new 01426 //RFLAG(workLev,ni,nj,nk,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc); 01427 changeFlag(workLev,ni,nj,nk,workSet, (CFInter|CFNoInterpolSrc)); 01428 if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ni,nj,nk)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " ); 01429 } 01430 /* prevent surrounding interface cells from getting removed as empty cells 01431 * (also cells that are not newly inited) */ 01432 if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) { 01433 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete); 01434 changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete)); 01435 // also add to list... 01436 addToNewInterList(ni,nj,nk); 01437 } // NEW? 01438 } 01439 01440 // NEW? no extra loop... 01441 //RFLAG(workLev,i,j,k, workSet) = CFFluid; 01442 changeFlag(workLev,i,j,k, workSet,CFFluid); 01443 } 01444 01445 /* remove empty interface cells that are not allowed to be removed anyway 01446 * this is important, otherwise the dreaded cell-type-flickering can occur! */ 01447 for(int n=0; n<(int)mListEmpty.size(); n++) { 01448 int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z; 01449 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) { 01450 // treat as "new inter" 01451 addToNewInterList(i,j,k); 01452 // remove entry 01453 if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM 01454 if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 01455 mListEmpty.pop_back(); 01456 n--; 01457 } 01458 } 01459 01460 01461 /* problems arise when adjacent cells empty&fill -> 01462 let fill cells+surrounding interface cells have the higher importance */ 01463 for( vector<LbmPoint>::iterator iter=mListEmpty.begin(); 01464 iter != mListEmpty.end(); iter++ ) { 01465 int i=iter->x, j=iter->y, k=iter->z; 01466 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG 01467 if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); 01468 01469 /* set surrounding fluid cells to interface cells */ 01470 FORDF1 { 01471 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 01472 if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){ 01473 // init fluid->interface 01474 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 01475 changeFlag(workLev,ni,nj,nk, workSet, CFInter); 01476 /* new mass = current density */ 01477 LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC); 01478 for(int rl=1; rl< this->cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); } 01479 QCELL(workLev,ni,nj,nk, workSet, dMass) = nbrho; 01480 QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 1.0; 01481 01482 // store point 01483 addToNewInterList(ni,nj,nk); 01484 } 01485 if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){ 01486 // test, also add to list... 01487 addToNewInterList(ni,nj,nk); 01488 } // NEW? 01489 } 01490 01491 /* for symmetry, set our flag right now */ 01492 changeFlag(workLev,i,j,k, workSet, CFEmpty); 01493 // mark cell not be changed mass... - not necessary, not in list anymore anyway! 01494 } // emptylist 01495 01496 01497 01498 // precompute weights to get rid of order dependancies 01499 vector<lbmFloatSet> vWeights; 01500 vWeights.resize( mListFull.size() + mListEmpty.size() ); 01501 int weightIndex = 0; 01502 int nbCount = 0; 01503 LbmFloat nbWeights[LBM_DFNUM]; 01504 LbmFloat nbTotWeights = 0.0; 01505 for( vector<LbmPoint>::iterator iter=mListFull.begin(); 01506 iter != mListFull.end(); iter++ ) { 01507 int i=iter->x, j=iter->y, k=iter->z; 01508 nbCount = 0; nbTotWeights = 0.0; 01509 FORDF1 { 01510 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 01511 if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) { 01512 nbCount++; 01513 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT 01514 else nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l); // NEWSURFT 01515 nbTotWeights += nbWeights[l]; 01516 } else { 01517 nbWeights[l] = -100.0; // DEBUG; 01518 } 01519 } 01520 if(nbCount>0) { 01521 //errMsg("FF I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights); 01522 vWeights[weightIndex].val[0] = nbTotWeights; 01523 FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; } 01524 vWeights[weightIndex].numNbs = (LbmFloat)nbCount; 01525 } else { 01526 vWeights[weightIndex].numNbs = 0.0; 01527 } 01528 weightIndex++; 01529 } 01530 for( vector<LbmPoint>::iterator iter=mListEmpty.begin(); 01531 iter != mListEmpty.end(); iter++ ) { 01532 int i=iter->x, j=iter->y, k=iter->z; 01533 nbCount = 0; nbTotWeights = 0.0; 01534 FORDF1 { 01535 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 01536 if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) { 01537 nbCount++; 01538 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT 01539 else nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l); // NEWSURFT 01540 nbTotWeights += nbWeights[l]; 01541 } else { 01542 nbWeights[l] = -100.0; // DEBUG; 01543 } 01544 } 01545 if(nbCount>0) { 01546 //errMsg("EE I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights); 01547 vWeights[weightIndex].val[0] = nbTotWeights; 01548 FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; } 01549 vWeights[weightIndex].numNbs = (LbmFloat)nbCount; 01550 } else { 01551 vWeights[weightIndex].numNbs = 0.0; 01552 } 01553 weightIndex++; 01554 } 01555 weightIndex = 0; 01556 01557 01558 /* process full list entries, filled cells are done after this loop */ 01559 for( vector<LbmPoint>::iterator iter=mListFull.begin(); 01560 iter != mListFull.end(); iter++ ) { 01561 int i=iter->x, j=iter->y, k=iter->z; 01562 01563 LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC); 01564 FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho 01565 01566 LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho; 01567 if(vWeights[weightIndex].numNbs>0.0) { 01568 const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0]; 01569 //errMsg("FF I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp); 01570 FORDF1 { 01571 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 01572 if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) { 01573 LbmFloat change = -1.0; 01574 if(nbTotWeightsp>0.0) { 01575 //change = massChange * ( nbWeights[l]/nbTotWeightsp ); 01576 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp ); 01577 } else { 01578 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs); 01579 } 01580 QCELL(workLev,ni,nj,nk, workSet, dMass) += change; 01581 } 01582 } 01583 massChange = 0.0; 01584 } else { 01585 // Problem! no interface neighbors 01586 mFixMass += massChange; 01587 //TTT mNumProblems++; 01588 //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<mFixMass); 01589 } 01590 weightIndex++; 01591 01592 // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid; 01593 QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused? 01594 QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused? 01595 } // fulllist 01596 01597 01598 /* now, finally handle the empty cells - order is important, has to be after 01599 * full cell handling */ 01600 for( vector<LbmPoint>::iterator iter=mListEmpty.begin(); 01601 iter != mListEmpty.end(); iter++ ) { 01602 int i=iter->x, j=iter->y, k=iter->z; 01603 01604 LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass); 01605 if(vWeights[weightIndex].numNbs>0.0) { 01606 const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0]; 01607 //errMsg("EE I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp); 01608 FORDF1 { 01609 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 01610 if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) { 01611 LbmFloat change = -1.0; 01612 if(nbTotWeightsp>0.0) { 01613 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp ); 01614 } else { 01615 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs); 01616 } 01617 QCELL(workLev, ni,nj,nk, workSet, dMass) += change; 01618 } 01619 } 01620 massChange = 0.0; 01621 } else { 01622 // Problem! no interface neighbors 01623 mFixMass += massChange; 01624 //TTT mNumProblems++; 01625 //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<mFixMass); 01626 } 01627 weightIndex++; 01628 01629 // finally... make it empty 01630 // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty; 01631 QCELL(workLev,i,j,k, workSet, dMass) = 0.0; 01632 QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0; 01633 } 01634 for( vector<LbmPoint>::iterator iter=mListEmpty.begin(); 01635 iter != mListEmpty.end(); iter++ ) { 01636 int i=iter->x, j=iter->y, k=iter->z; 01637 changeFlag(workLev,i,j,k, otherSet, CFEmpty); 01638 } 01639 01640 01641 // check if some of the new interface cells can be removed again 01642 // never happens !!! not necessary 01643 // calculate ffrac for new IF cells NEW 01644 01645 // how many are really new interface cells? 01646 int numNewIf = 0; 01647 for( vector<LbmPoint>::iterator iter=mListNewInter.begin(); 01648 iter != mListNewInter.end(); iter++ ) { 01649 int i=iter->x, j=iter->y, k=iter->z; 01650 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 01651 continue; 01652 // FIXME remove from list? 01653 } 01654 numNewIf++; 01655 } 01656 01657 // redistribute mass, reinit flags 01658 if(debugFlagreinit) errMsg("NEWIF", "total:"<<mListNewInter.size()); 01659 float newIfFac = 1.0/(LbmFloat)numNewIf; 01660 for( vector<LbmPoint>::iterator iter=mListNewInter.begin(); 01661 iter != mListNewInter.end(); iter++ ) { 01662 int i=iter->x, j=iter->y, k=iter->z; 01663 if((i<=0) || (j<=0) || 01664 (i>=mLevel[workLev].lSizex-1) || 01665 (j>=mLevel[workLev].lSizey-1) || 01666 ((LBMDIM==3) && ((k<=0) || (k>=mLevel[workLev].lSizez-1) ) ) 01667 ) { 01668 continue; } // new bc, dont treat cells on boundary NEWBC 01669 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 01670 //errMsg("???"," "<<PRINT_IJK); 01671 continue; 01672 } // */ 01673 01674 QCELL(workLev,i,j,k, workSet, dMass) += (mFixMass * newIfFac); 01675 01676 int nbored = 0; 01677 FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); } 01678 if(!(nbored & CFBndNoslip)) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; } 01679 if(!(nbored & CFFluid)) { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; } 01680 if(!(nbored & CFEmpty)) { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; } 01681 01682 if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) { 01683 RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete); 01684 } 01685 if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< convertCellFlagType2String(RFLAG(workLev,i,j,k, workSet))<<" wl"<<workLev ); 01686 } 01687 01688 // reinit fill fraction 01689 for( vector<LbmPoint>::iterator iter=mListNewInter.begin(); 01690 iter != mListNewInter.end(); iter++ ) { 01691 int i=iter->x, j=iter->y, k=iter->z; 01692 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; } 01693 01694 initInterfaceVars(workLev, i,j,k, workSet, false); //int level, int i,int j,int k,int workSet, bool initMass) { 01695 //LbmFloat nrho = 0.0; 01696 //FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); } 01697 //QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho; 01698 //QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT; 01699 } 01700 01701 if(mListNewInter.size()>0){ 01702 //errMsg("FixMassDisted"," fm:"<<mFixMass<<" nif:"<<mListNewInter.size() ); 01703 mFixMass = 0.0; 01704 } 01705 01706 // empty lists for next step 01707 mListFull.clear(); 01708 mListEmpty.clear(); 01709 mListNewInter.clear(); 01710 } // reinitFlags 01711 01712 01713