|
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 * Adaptivity functions 00010 * 00011 *****************************************************************************/ 00012 00013 #include "solver_class.h" 00014 #include "solver_relax.h" 00015 #include "particletracer.h" 00016 00017 00018 00019 /*****************************************************************************/ 00021 /*****************************************************************************/ 00022 00023 00024 00025 void LbmFsgrSolver::coarseCalculateFluxareas(int lev) 00026 { 00027 FSGR_FORIJK_BOUNDS(lev) { 00028 if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) { 00029 if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) { 00030 LbmFloat totArea = mFsgrCellArea[0]; // for l=0 00031 for(int l=1; l<this->cDirNum; l++) { 00032 int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l]; 00033 if(RFLAG(lev+1, ni,nj,nk, mLevel[lev+1].setCurr)& 00034 (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused) 00035 ) { 00036 totArea += mFsgrCellArea[l]; 00037 } 00038 } // l 00039 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = totArea; 00040 //continue; 00041 } else 00042 if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & (CFEmpty|CFUnused)) { 00043 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 1.0; 00044 //continue; 00045 } else { 00046 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 0.0; 00047 } 00048 //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) ); 00049 } 00050 } // } TEST DEBUG 00051 if(!this->mSilent){ debMsgStd("coarseCalculateFluxareas",DM_MSG,"level "<<lev<<" calculated", 7); } 00052 } 00053 00054 void LbmFsgrSolver::coarseAdvance(int lev) 00055 { 00056 LbmFloat calcCurrentMass = 0.0; 00057 LbmFloat calcCurrentVolume = 0.0; 00058 00059 LbmFloat *ccel = NULL; 00060 LbmFloat *tcel = NULL; 00061 LbmFloat m[LBM_DFNUM]; 00062 LbmFloat rho, ux, uy, uz, tmp, usqr, lcsmqo; 00063 #if OPT3D==1 00064 LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega; 00065 #endif // OPT3D==true 00066 m[0] = tmp = usqr = 0.0; 00067 00068 coarseCalculateFluxareas(lev); 00069 // copied from fineAdv. 00070 CellFlagType *pFlagSrc = &RFLAG(lev, 1,1,getForZMin1(),SRCS(lev)); 00071 CellFlagType *pFlagDst = &RFLAG(lev, 1,1,getForZMin1(),TSET(lev)); 00072 pFlagSrc -= 1; 00073 pFlagDst -= 1; 00074 ccel = RACPNT(lev, 1,1,getForZMin1() ,SRCS(lev)); // QTEST 00075 ccel -= QCELLSTEP; 00076 tcel = RACPNT(lev, 1,1,getForZMin1() ,TSET(lev)); // QTEST 00077 tcel -= QCELLSTEP; 00078 //if(strstr(this->getName().c_str(),"Debug")){ errMsg("DEBUG","DEBUG!!!!!!!!!!!!!!!!!!!!!!!"); } 00079 00080 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { 00081 for(int j=1;j<mLevel[lev].lSizey-1;++j) { 00082 for(int i=1;i<mLevel[lev].lSizex-1;++i) { 00083 #if FSGR_STRICT_DEBUG==1 00084 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG 00085 #endif 00086 pFlagSrc++; 00087 pFlagDst++; 00088 ccel += QCELLSTEP; 00089 tcel += QCELLSTEP; 00090 00091 // from coarse cells without unused nbs are not necessary...! -> remove 00092 if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 00093 bool invNb = false; 00094 FORDF1 { if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; } } 00095 if(!invNb) { 00096 // WARNING - this modifies source flag array... 00097 *pFlagSrc = CFFluid|CFGrNorm; 00098 #if ELBEEM_PLUGIN!=1 00099 errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK); 00100 #endif // ELBEEM_PLUGIN!=1 00101 // move to perform coarsening? 00102 } 00103 } // */ 00104 00105 #if FSGR_STRICT_DEBUG==1 00106 *pFlagDst = *pFlagSrc; // always set other set... 00107 #else 00108 *pFlagDst = (*pFlagSrc & (~CFGrCoarseInited)); // always set other set... , remove coarse inited flag 00109 #endif 00110 00111 // old INTCFCOARSETEST==1 00112 if((*pFlagSrc) & CFGrFromCoarse) { // interpolateFineFromCoarse test! 00113 if(( this->mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) { 00114 FORDF0 { RAC(tcel,l) = RAC(ccel,l); } 00115 } else { 00116 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false); 00117 this->mNumUsedCells++; 00118 } 00119 continue; // interpolateFineFromCoarse test! 00120 } // interpolateFineFromCoarse test! old INTCFCOARSETEST==1 00121 00122 if( ((*pFlagSrc) & (CFFluid)) ) { 00123 ccel = RACPNT(lev, i,j,k ,SRCS(lev)); 00124 tcel = RACPNT(lev, i,j,k ,TSET(lev)); 00125 00126 if( ((*pFlagSrc) & (CFGrFromFine)) ) { 00127 FORDF0 { RAC(tcel,l) = RAC(ccel,l); } // always copy...? 00128 continue; // comes from fine grid 00129 } 00130 // also ignore CFGrFromCoarse 00131 else if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 00132 FORDF0 { RAC(tcel,l) = RAC(ccel,l); } // always copy...? 00133 continue; 00134 } 00135 00136 OPTIMIZED_STREAMCOLLIDE; 00137 *pFlagDst |= CFNoBndFluid; // test? 00138 calcCurrentVolume += RAC(ccel,dFlux); 00139 calcCurrentMass += RAC(ccel,dFlux)*rho; 00140 //ebugMarkCell(lev+1, 2*i+1,2*j+1,2*k ); 00141 #if FSGR_STRICT_DEBUG==1 00142 if(rho<-1.0){ debugMarkCell(lev, i,j,k ); 00143 errMsg("INVRHOCELL_CHECK"," l"<<lev<<" "<< PRINT_IJK<<" rho:"<<rho ); 00144 CAUSE_PANIC; 00145 } 00146 #endif // FSGR_STRICT_DEBUG==1 00147 this->mNumUsedCells++; 00148 00149 } 00150 } 00151 pFlagSrc+=2; // after x 00152 pFlagDst+=2; // after x 00153 ccel += (QCELLSTEP*2); 00154 tcel += (QCELLSTEP*2); 00155 } 00156 pFlagSrc+= mLevel[lev].lSizex*2; // after y 00157 pFlagDst+= mLevel[lev].lSizex*2; // after y 00158 ccel += (QCELLSTEP*mLevel[lev].lSizex*2); 00159 tcel += (QCELLSTEP*mLevel[lev].lSizex*2); 00160 } // all cell loop k,j,i 00161 00162 00163 //errMsg("coarseAdvance","level "<<lev<<" stepped from "<<mLevel[lev].setCurr<<" to "<<mLevel[lev].setOther); 00164 if(!this->mSilent){ errMsg("coarseAdvance","level "<<lev<<" stepped from "<<SRCS(lev)<<" to "<<TSET(lev)); } 00165 // */ 00166 00167 // update other set 00168 mLevel[lev].setOther = mLevel[lev].setCurr; 00169 mLevel[lev].setCurr ^= 1; 00170 mLevel[lev].lsteps++; 00171 mLevel[lev].lmass = calcCurrentMass * mLevel[lev].lcellfactor; 00172 mLevel[lev].lvolume = calcCurrentVolume * mLevel[lev].lcellfactor; 00173 #if ELBEEM_PLUGIN!=1 00174 debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "mass: lev="<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<" lcf="<< mLevel[lev].lcellfactor, 8 ); 00175 debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "volume: lev="<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<" lcf="<< mLevel[lev].lcellfactor, 8 ); 00176 #endif // ELBEEM_PLUGIN 00177 } 00178 00179 /*****************************************************************************/ 00181 /*****************************************************************************/ 00182 00183 00184 // get dfs from level (lev+1) to (lev) coarse border nodes 00185 void LbmFsgrSolver::coarseRestrictFromFine(int lev) 00186 { 00187 if((lev<0) || ((lev+1)>mMaxRefine)) return; 00188 #if FSGR_STRICT_DEBUG==1 00189 // reset all unused cell values to invalid 00190 int unuCnt = 0; 00191 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { 00192 for(int j=1;j<mLevel[lev].lSizey-1;++j) { 00193 for(int i=1;i<mLevel[lev].lSizex-1;++i) { 00194 CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,mLevel[lev].setCurr); 00195 if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 00196 FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0; } 00197 unuCnt++; 00198 // set here 00199 } else if( ((*pFlagSrc) & (CFFluid|CFGrNorm)) == (CFFluid|CFGrNorm) ) { 00200 // simulated... 00201 } else { 00202 // reset in interpolation 00203 //errMsg("coarseRestrictFromFine"," reset l"<<lev<<" "<<PRINT_IJK); 00204 } 00205 if( ((*pFlagSrc) & (CFEmpty|CFUnused)) ) { // test, also reset? 00206 FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0; } 00207 } // test 00208 } } } 00209 errMsg("coarseRestrictFromFine"," reset l"<<lev<<" fluid|coarseBorder cells: "<<unuCnt); 00210 #endif // FSGR_STRICT_DEBUG==1 00211 const int srcSet = mLevel[lev+1].setCurr; 00212 const int dstSet = mLevel[lev].setCurr; 00213 00214 //restrict 00215 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { 00216 for(int j=1;j<mLevel[lev].lSizey-1;++j) { 00217 for(int i=1;i<mLevel[lev].lSizex-1;++i) { 00218 CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,dstSet); 00219 if((*pFlagSrc) & (CFFluid)) { 00220 if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 00221 // do resctriction 00222 mNumInterdCells++; 00223 coarseRestrictCell(lev, i,j,k,srcSet,dstSet); 00224 00225 this->mNumUsedCells++; 00226 } // from fine & fluid 00227 else { 00228 if(RFLAG(lev+1, 2*i,2*j,2*k,srcSet) & CFGrFromCoarse) { 00229 RFLAG(lev, i,j,k,dstSet) |= CFGrToFine; 00230 } else { 00231 RFLAG(lev, i,j,k,dstSet) &= (~CFGrToFine); 00232 } 00233 } 00234 } // & fluid 00235 }}} 00236 if(!this->mSilent){ errMsg("coarseRestrictFromFine"," from l"<<(lev+1)<<",s"<<mLevel[lev+1].setCurr<<" to l"<<lev<<",s"<<mLevel[lev].setCurr); } 00237 } 00238 00239 bool LbmFsgrSolver::adaptGrid(int lev) { 00240 if((lev<0) || ((lev+1)>mMaxRefine)) return false; 00241 bool change = false; 00242 { // refinement, PASS 1-3 00243 00244 //bool nbsok; 00245 // FIXME remove TIMEINTORDER ? 00246 LbmFloat interTime = 0.0; 00247 // update curr from other, as streaming afterwards works on curr 00248 // thus read only from srcSet, modify other 00249 const int srcSet = mLevel[lev].setOther; 00250 const int dstSet = mLevel[lev].setCurr; 00251 const int srcFineSet = mLevel[lev+1].setCurr; 00252 const bool debugRefinement = false; 00253 00254 // use //template functions for 2D/3D 00255 /*if(strstr(this->getName().c_str(),"Debug")) 00256 if(lev+1==mMaxRefine) { // mixborder 00257 for(int l=0;((l<this->cDirNum) && (!removeFromFine)); l++) { // FARBORD 00258 int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l]; 00259 if(RFLAG(lev+1, ni,nj,nk, srcFineSet)&CFBnd) { // NEWREFT 00260 removeFromFine=true; 00261 } 00262 } 00263 } // FARBORD */ 00264 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { 00265 for(int j=1;j<mLevel[lev].lSizey-1;++j) { 00266 for(int i=1;i<mLevel[lev].lSizex-1;++i) { 00267 00268 if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) { 00269 bool removeFromFine = false; 00270 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine); 00271 CellFlagType reqType = CFGrNorm; 00272 if(lev+1==mMaxRefine) reqType = CFNoBndFluid; 00273 00274 if( (RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & reqType) && 00275 (!(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & (notAllowed)) ) ){ // ok 00276 } else { 00277 removeFromFine=true; 00278 } 00279 00280 if(removeFromFine) { 00281 // dont turn CFGrFromFine above interface cells into CFGrNorm 00282 //errMsg("performRefinement","Removing CFGrFromFine on lev"<<lev<<" " <<PRINT_IJK<<" srcflag:"<<convertCellFlagType2String(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet)) <<" set:"<<dstSet ); 00283 RFLAG(lev, i,j,k, dstSet) = CFEmpty; 00284 #if FSGR_STRICT_DEBUG==1 00285 // for interpolation later on during fine grid fixing 00286 // these cells are still correctly inited 00287 RFLAG(lev, i,j,k, dstSet) |= CFGrCoarseInited; // remove later on? FIXME? 00288 #endif // FSGR_STRICT_DEBUG==1 00289 //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFEmpty; // FLAGTEST 00290 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,i,j,k); 00291 change=true; 00292 mNumFsgrChanges++; 00293 for(int l=1; l<this->cDirNum; l++) { 00294 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 00295 //errMsg("performRefinement","On lev:"<<lev<<" check: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet<<" = "<<convertCellFlagType2String(RFLAG(lev, ni,nj,nk, srcSet)) ); 00296 if( ( RFLAG(lev, ni,nj,nk, srcSet)&CFFluid ) && 00297 (!(RFLAG(lev, ni,nj,nk, srcSet)&CFGrFromFine)) ) { // dont change status of nb. from fine cells 00298 // tag as inited for debugging, cell contains fluid DFs anyway 00299 RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromFine|CFGrCoarseInited; 00300 //errMsg("performRefinement","On lev:"<<lev<<" set to from fine: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet); 00301 //if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 00302 } 00303 } // l 00304 00305 // FIXME fix fine level? 00306 } 00307 00308 // recheck from fine flag 00309 } 00310 }}} // TEST 00311 // PASS 1 */ 00312 00313 00314 /*if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 00315 bool invNb = false; 00316 FORDF1 { 00317 if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; } 00318 } 00319 if(!invNb) { 00320 *pFlagSrc = CFFluid|CFGrNorm; 00321 errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK); 00322 } 00323 } // */ 00324 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST 00325 for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST 00326 for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST 00327 00328 // test from coarseAdvance 00329 // from coarse cells without unused nbs are not necessary...! -> remove 00330 00331 if(RFLAG(lev, i,j,k, srcSet) & CFGrFromCoarse) { 00332 00333 // from coarse cells without unused nbs are not necessary...! -> remove 00334 bool invNb = false; 00335 bool fluidNb = false; 00336 for(int l=1; l<this->cDirNum; l++) { 00337 if(RFLAG_NB(lev, i, j, k, srcSet, l) & CFUnused) { invNb = true; } 00338 if(RFLAG_NB(lev, i, j, k, srcSet, l) & (CFGrNorm)) { fluidNb = true; } 00339 } 00340 if(!invNb) { 00341 // no unused cells around -> calculate normally from now on 00342 RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm; 00343 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 00344 change=true; 00345 mNumFsgrChanges++; 00346 } // from advance 00347 if(!fluidNb) { 00348 // no fluid cells near -> no transfer necessary 00349 RFLAG(lev, i,j,k, dstSet) = CFUnused; 00350 //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFUnused; // FLAGTEST 00351 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 00352 change=true; 00353 mNumFsgrChanges++; 00354 } // from advance 00355 00356 00357 // dont allow double transfer 00358 // this might require fixing the neighborhood 00359 if(RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse)) { 00360 // dont turn CFGrFromFine above interface cells into CFGrNorm 00361 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<lev<<" " <<PRINT_IJK<<" due to finer from coarse cell " ); 00362 RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm; 00363 if(lev>0) RFLAG(lev-1, i/2,j/2,k/2, mLevel[lev-1].setCurr) &= (~CFGrToFine); // TODO add more of these? 00364 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 00365 change=true; 00366 mNumFsgrChanges++; 00367 for(int l=1; l<this->cDirNum; l++) { 00368 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 00369 if(RFLAG(lev, ni,nj,nk, srcSet)&(CFGrNorm)) { //ok 00370 for(int m=1; m<this->cDirNum; m++) { 00371 int mi= ni +this->dfVecX[m], mj= nj +this->dfVecY[m], mk= nk +this->dfVecZ[m]; 00372 if(RFLAG(lev, mi, mj, mk, srcSet)&CFUnused) { 00373 // norm cells in neighborhood with unused nbs have to be new border... 00374 RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromCoarse; 00375 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 00376 } 00377 } 00378 // these alreay have valid values... 00379 } 00380 else if(RFLAG(lev, ni,nj,nk, srcSet)&(CFUnused)) { //ok 00381 // this should work because we have a valid neighborhood here for now 00382 interpolateCellFromCoarse(lev, ni, nj, nk, dstSet, interTime, CFFluid|CFGrFromCoarse, false); 00383 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 00384 mNumFsgrChanges++; 00385 } 00386 } // l 00387 } // double transer 00388 00389 } // from coarse 00390 00391 } } } 00392 // PASS 2 */ 00393 00394 00395 // fix dstSet from fine cells here 00396 // warning - checks CFGrFromFine on dstset changed before! 00397 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST 00398 for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST 00399 for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST 00400 00401 //if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) { 00402 if(RFLAG(lev, i,j,k, dstSet) & CFGrFromFine) { 00403 // modify finer level border 00404 if((RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse))) { 00405 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" from l"<<lev<<" " <<PRINT_IJK ); 00406 CellFlagType setf = CFFluid; 00407 if(lev+1 < mMaxRefine) setf = CFFluid|CFGrNorm; 00408 RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)=setf; 00409 change=true; 00410 mNumFsgrChanges++; 00411 for(int l=1; l<this->cDirNum; l++) { 00412 int bi=(2*i)+this->dfVecX[l], bj=(2*j)+this->dfVecY[l], bk=(2*k)+this->dfVecZ[l]; 00413 if(RFLAG(lev+1, bi, bj, bk, srcFineSet)&(CFGrFromCoarse)) { 00414 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) ); 00415 RFLAG(lev+1, bi, bj, bk, srcFineSet) = setf; 00416 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 00417 } 00418 else if(RFLAG(lev+1, bi, bj, bk, srcFineSet)&(CFUnused )) { 00419 //errMsg("performRefinement","Removing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) ); 00420 interpolateCellFromCoarse(lev+1, bi, bj, bk, srcFineSet, interTime, setf, false); 00421 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 00422 mNumFsgrChanges++; 00423 } 00424 } 00425 for(int l=1; l<this->cDirNum; l++) { 00426 int bi=(2*i)+this->dfVecX[l], bj=(2*j)+this->dfVecY[l], bk=(2*k)+this->dfVecZ[l]; 00427 if( (RFLAG(lev+1, bi, bj, bk, srcFineSet)&CFFluid ) && 00428 (!(RFLAG(lev+1, bi, bj, bk, srcFineSet)&CFGrFromCoarse)) ) { 00429 // all unused nbs now of coarse have to be from coarse 00430 for(int m=1; m<this->cDirNum; m++) { 00431 int mi= bi +this->dfVecX[m], mj= bj +this->dfVecY[m], mk= bk +this->dfVecZ[m]; 00432 if(RFLAG(lev+1, mi, mj, mk, srcFineSet)&CFUnused) { 00433 //errMsg("performRefinement","Changing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(mi,mj,mk) ); 00434 interpolateCellFromCoarse(lev+1, mi, mj, mk, srcFineSet, interTime, CFFluid|CFGrFromCoarse, false); 00435 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,mi,mj,mk); 00436 mNumFsgrChanges++; 00437 } 00438 } 00439 // nbs prepared... 00440 } 00441 } 00442 } 00443 00444 } // convert regions of from fine 00445 }}} // TEST 00446 // PASS 3 */ 00447 00448 if(!this->mSilent){ errMsg("performRefinement"," for l"<<lev<<" done ("<<change<<") " ); } 00449 } // PASS 1-3 00450 // refinement done 00451 00452 //LbmFsgrSolver::performCoarsening(int lev) { 00453 { // PASS 4,5 00454 bool nbsok; 00455 // WARNING 00456 // now work on modified curr set 00457 const int srcSet = mLevel[lev].setCurr; 00458 const int dstlev = lev+1; 00459 const int dstFineSet = mLevel[dstlev].setCurr; 00460 const bool debugCoarsening = false; 00461 00462 // PASS 5 test DEBUG 00463 /*if(this->mInitDone) { 00464 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { 00465 for(int j=1;j<mLevel[lev].lSizey-1;++j) { 00466 for(int i=1;i<mLevel[lev].lSizex-1;++i) { 00467 if(RFLAG(lev, i,j,k, srcSet) & CFEmpty) { 00468 // check empty -> from fine conversion 00469 bool changeToFromFine = false; 00470 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine); 00471 CellFlagType reqType = CFGrNorm; 00472 if(lev+1==mMaxRefine) reqType = CFNoBndFluid; 00473 if( (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) && 00474 (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) ) ){ 00475 changeToFromFine=true; } 00476 if(changeToFromFine) { 00477 change = true; 00478 mNumFsgrChanges++; 00479 RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrFromFine; 00480 if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 00481 // same as restr from fine func! not necessary ?! 00482 // coarseRestrictFromFine part 00483 coarseRestrictCell(lev, i,j,k,srcSet, dstFineSet); 00484 } 00485 } // only check empty cells 00486 }}} // TEST! 00487 } // PASS 5 */ 00488 00489 // use //template functions for 2D/3D 00490 /*if(strstr(this->getName().c_str(),"Debug")) 00491 if((nbsok)&&(lev+1==mMaxRefine)) { // mixborder 00492 for(int l=0;((l<this->cDirNum) && (nbsok)); l++) { // FARBORD 00493 int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l]; 00494 if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT 00495 nbsok=false; 00496 } 00497 } 00498 } // FARBORD */ 00499 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { 00500 for(int j=1;j<mLevel[lev].lSizey-1;++j) { 00501 for(int i=1;i<mLevel[lev].lSizex-1;++i) { 00502 00503 // from coarse cells without unused nbs are not necessary...! -> remove 00504 // perform check from coarseAdvance here? 00505 if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) { 00506 // remove from fine cells now that are completely in fluid 00507 // FIXME? check that new from fine in performRefinement never get deleted here afterwards? 00508 // or more general, one cell never changed more than once? 00509 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine); 00510 //const CellFlagType notNbAllowed = (CFInter|CFBnd|CFGrFromFine); unused 00511 CellFlagType reqType = CFGrNorm; 00512 if(lev+1==mMaxRefine) reqType = CFNoBndFluid; 00513 00514 nbsok = true; 00515 for(int l=0; l<this->cDirNum && nbsok; l++) { 00516 int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l]; 00517 if( (RFLAG(lev+1, ni,nj,nk, dstFineSet) & reqType) && 00518 (!(RFLAG(lev+1, ni,nj,nk, dstFineSet) & (notAllowed)) ) ){ 00519 // ok 00520 } else { 00521 nbsok=false; 00522 } 00523 // FARBORD 00524 } 00525 // dont turn CFGrFromFine above interface cells into CFGrNorm 00526 // now check nbs on same level 00527 for(int l=1; l<this->cDirNum && nbsok; l++) { 00528 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 00529 if(RFLAG(lev, ni,nj,nk, srcSet)&(CFFluid)) { //ok 00530 } else { 00531 nbsok = false; 00532 } 00533 } // l 00534 00535 if(nbsok) { 00536 // conversion to coarse fluid cell 00537 change = true; 00538 mNumFsgrChanges++; 00539 RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrNorm; 00540 // dfs are already ok... 00541 //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine changed to CFGrNorm at lev"<<lev<<" " <<PRINT_IJK ); 00542 if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 00543 00544 // only check complete cubes 00545 for(int dx=-1;dx<=1;dx+=2) { 00546 for(int dy=-1;dy<=1;dy+=2) { 00547 for(int dz=-1*(LBMDIM&1);dz<=1*(LBMDIM&1);dz+=2) { // 2d/3d 00548 // check for norm and from coarse, as the coarse level might just have been refined... 00549 if( 00550 // we now the flag of the current cell! ( RFLAG(lev, i , j , k , srcSet)&(CFGrNorm)) && 00551 ( RFLAG(lev, i+dx, j , k , srcSet)&(CFGrNorm|CFGrFromCoarse)) && 00552 ( RFLAG(lev, i , j+dy, k , srcSet)&(CFGrNorm|CFGrFromCoarse)) && 00553 ( RFLAG(lev, i , j , k+dz, srcSet)&(CFGrNorm|CFGrFromCoarse)) && 00554 00555 ( RFLAG(lev, i+dx, j+dy, k , srcSet)&(CFGrNorm|CFGrFromCoarse)) && 00556 ( RFLAG(lev, i+dx, j , k+dz, srcSet)&(CFGrNorm|CFGrFromCoarse)) && 00557 ( RFLAG(lev, i , j+dy, k+dz, srcSet)&(CFGrNorm|CFGrFromCoarse)) && 00558 ( RFLAG(lev, i+dx, j+dy, k+dz, srcSet)&(CFGrNorm|CFGrFromCoarse)) 00559 ) { 00560 // middle source node on higher level 00561 int dstx = (2*i)+dx; 00562 int dsty = (2*j)+dy; 00563 int dstz = (2*k)+dz; 00564 00565 mNumFsgrChanges++; 00566 RFLAG(dstlev, dstx,dsty,dstz, dstFineSet) = CFUnused; 00567 RFLAG(dstlev, dstx,dsty,dstz, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST 00568 //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init center unused set l"<<dstlev<<" at "<<PRINT_VEC(dstx,dsty,dstz) ); 00569 00570 for(int l=1; l<this->cDirNum; l++) { 00571 int dstni=dstx+this->dfVecX[l], dstnj=dsty+this->dfVecY[l], dstnk=dstz+this->dfVecZ[l]; 00572 if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFFluid)) { 00573 RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse; 00574 } 00575 if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFInter)) { 00576 //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init CHECK Warning - deleting interface cell..."); 00577 this->mFixMass += QCELL( dstlev, dstni,dstnj,dstnk, dstFineSet, dMass); 00578 RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse; 00579 } 00580 } // l 00581 00582 // again check nb flags of all surrounding cells to see if any from coarse 00583 // can be convted to unused 00584 for(int l=1; l<this->cDirNum; l++) { 00585 int dstni=dstx+this->dfVecX[l], dstnj=dsty+this->dfVecY[l], dstnk=dstz+this->dfVecZ[l]; 00586 // have to be at least from coarse here... 00587 //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" "<< convertCellFlagType2String(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)) ); 00588 if(!(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFUnused) )) { 00589 bool delok = true; 00590 // careful long range here... check domain bounds? 00591 for(int m=1; m<this->cDirNum; m++) { 00592 int chkni=dstni+this->dfVecX[m], chknj=dstnj+this->dfVecY[m], chknk=dstnk+this->dfVecZ[m]; 00593 if(RFLAG(dstlev, chkni,chknj,chknk, dstFineSet)&(CFUnused|CFGrFromCoarse)) { 00594 // this nb cell is ok for deletion 00595 } else { 00596 delok=false; // keep it! 00597 } 00598 //errMsg("performCoarsening"," CHECK "<<PRINT_VEC(dstni,dstnj,dstnk)<<" to "<<PRINT_VEC( chkni,chknj,chknk )<<" f:"<< convertCellFlagType2String( RFLAG(dstlev, chkni,chknj,chknk, dstFineSet))<<" nbsok"<<delok ); 00599 } 00600 //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" ok"<<delok ); 00601 if(delok) { 00602 mNumFsgrChanges++; 00603 RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFUnused; 00604 RFLAG(dstlev, dstni,dstnj,dstnk, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST 00605 if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(dstlev,dstni,dstnj,dstnk); 00606 } 00607 } 00608 } // l 00609 // treat subcube 00610 //ebugMarkCell(lev,i+dx,j+dy,k+dz); 00611 //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init, dir:"<<PRINT_VEC(dx,dy,dz) ); 00612 } 00613 } } } 00614 00615 } // ? 00616 } // convert regions of from fine 00617 }}} // TEST! 00618 // PASS 4 */ 00619 00620 // reinit cell area value 00621 /*if( RFLAG(lev, i,j,k,srcSet) & CFFluid) { 00622 if( RFLAG(lev+1, i*2,j*2,k*2,dstFineSet) & CFGrFromCoarse) { 00623 LbmFloat totArea = mFsgrCellArea[0]; // for l=0 00624 for(int l=1; l<this->cDirNum; l++) { 00625 int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l]; 00626 if(RFLAG(lev+1, ni,nj,nk, dstFineSet)& 00627 (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused) 00628 //(CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused) 00629 ) { 00630 //LbmFloat area = 0.25; if(this->dfVecX[l]!=0) area *= 0.5; if(this->dfVecY[l]!=0) area *= 0.5; if(this->dfVecZ[l]!=0) area *= 0.5; 00631 totArea += mFsgrCellArea[l]; 00632 } 00633 } // l 00634 QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 00635 QCELL(lev, i,j,k,srcSet, dFlux) = totArea; 00636 } else { 00637 QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 00638 QCELL(lev, i,j,k,srcSet, dFlux) = 1.0; 00639 } 00640 //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,srcSet, dFlux) ); 00641 } 00642 // */ 00643 00644 00645 // PASS 5 org 00646 /*if(strstr(this->getName().c_str(),"Debug")) 00647 if((changeToFromFine)&&(lev+1==mMaxRefine)) { // mixborder 00648 for(int l=0;((l<this->cDirNum) && (changeToFromFine)); l++) { // FARBORD 00649 int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l]; 00650 if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT 00651 changeToFromFine=false; } 00652 } 00653 }// FARBORD */ 00654 //if(!this->mInitDone) { 00655 for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { 00656 for(int j=1;j<mLevel[lev].lSizey-1;++j) { 00657 for(int i=1;i<mLevel[lev].lSizex-1;++i) { 00658 00659 00660 if(RFLAG(lev, i,j,k, srcSet) & CFEmpty) { 00661 // check empty -> from fine conversion 00662 bool changeToFromFine = false; 00663 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine); 00664 CellFlagType reqType = CFGrNorm; 00665 if(lev+1==mMaxRefine) reqType = CFNoBndFluid; 00666 00667 if( (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) && 00668 (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) ) ){ 00669 // DEBUG 00670 changeToFromFine=true; 00671 } 00672 00673 // FARBORD 00674 00675 if(changeToFromFine) { 00676 change = true; 00677 mNumFsgrChanges++; 00678 RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrFromFine; 00679 if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 00680 // same as restr from fine func! not necessary ?! 00681 // coarseRestrictFromFine part 00682 } 00683 } // only check empty cells 00684 00685 }}} // TEST! 00686 //} // init done 00687 // PASS 5 */ 00688 } // coarsening, PASS 4,5 00689 00690 if(!this->mSilent){ errMsg("adaptGrid"," for l"<<lev<<" done " ); } 00691 return change; 00692 } 00693 00694 /*****************************************************************************/ 00696 /*****************************************************************************/ 00697 00698 void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet) 00699 { 00700 LbmFloat *ccel = RACPNT(lev+1, 2*i,2*j,2*k,srcSet); 00701 LbmFloat *tcel = RACPNT(lev , i,j,k ,dstSet); 00702 00703 LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0; 00704 //LbmFloat *ccel = NULL; 00705 //LbmFloat *tcel = NULL; 00706 #if OPT3D==1 00707 LbmFloat m[LBM_DFNUM]; 00708 // for macro add 00709 LbmFloat usqr; 00710 //LbmFloat *addfcel, *dstcell; 00711 LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM]; 00712 LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale; 00713 #else // OPT3D==true 00714 LbmFloat df[LBM_DFNUM]; 00715 LbmFloat omegaDst, omegaSrc; 00716 LbmFloat feq[LBM_DFNUM]; 00717 LbmFloat dfScale = mDfScaleUp; 00718 #endif // OPT3D==true 00719 00720 # if OPT3D==0 00721 // add up weighted dfs 00722 FORDF0{ df[l] = 0.0;} 00723 for(int n=0;(n<this->cDirNum); n++) { 00724 int ni=2*i+1*this->dfVecX[n], nj=2*j+1*this->dfVecY[n], nk=2*k+1*this->dfVecZ[n]; 00725 ccel = RACPNT(lev+1, ni,nj,nk,srcSet);// CFINTTEST 00726 const LbmFloat weight = mGaussw[n]; 00727 FORDF0{ 00728 LbmFloat cdf = weight * RAC(ccel,l); 00729 # if FSGR_STRICT_DEBUG==1 00730 if( cdf<-1.0 ){ errMsg("INVDFCREST_DFCHECK", PRINT_IJK<<" s"<<dstSet<<" from "<<PRINT_VEC(2*i,2*j,2*k)<<" s"<<srcSet<<" df"<<l<<":"<< df[l]); } 00731 # endif 00732 //errMsg("INVDFCREST_DFCHECK", PRINT_IJK<<" s"<<dstSet<<" from "<<PRINT_VEC(2*i,2*j,2*k)<<" s"<<srcSet<<" df"<<l<<":"<< df[l]<<" = "<<cdf<<" , w"<<weight); 00733 df[l] += cdf; 00734 } 00735 } 00736 00737 // calc rho etc. from weighted dfs 00738 rho = ux = uy = uz = 0.0; 00739 FORDF0{ 00740 LbmFloat cdf = df[l]; 00741 rho += cdf; 00742 ux += (this->dfDvecX[l]*cdf); 00743 uy += (this->dfDvecY[l]*cdf); 00744 uz += (this->dfDvecZ[l]*cdf); 00745 } 00746 00747 FORDF0{ feq[l] = this->getCollideEq(l, rho,ux,uy,uz); } 00748 if(mLevel[lev ].lcsmago>0.0) { 00749 const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feq); 00750 omegaDst = this->getLesOmega(mLevel[lev ].omega,mLevel[lev ].lcsmago,Qo); 00751 omegaSrc = this->getLesOmega(mLevel[lev+1].omega,mLevel[lev+1].lcsmago,Qo); 00752 } else { 00753 omegaDst = mLevel[lev+0].omega; /* NEWSMAGOT*/ 00754 omegaSrc = mLevel[lev+1].omega; 00755 } 00756 dfScale = (mLevel[lev ].timestep/mLevel[lev+1].timestep)* (1.0/omegaDst-1.0)/ (1.0/omegaSrc-1.0); // yu 00757 FORDF0{ 00758 RAC(tcel, l) = feq[l]+ (df[l]-feq[l])*dfScale; 00759 } 00760 # else // OPT3D 00761 // similar to OPTIMIZED_STREAMCOLLIDE_UNUSED 00762 00763 //rho = ux = uy = uz = 0.0; 00764 MSRC_C = CCELG_C(0) ; 00765 MSRC_N = CCELG_N(0) ; 00766 MSRC_S = CCELG_S(0) ; 00767 MSRC_E = CCELG_E(0) ; 00768 MSRC_W = CCELG_W(0) ; 00769 MSRC_T = CCELG_T(0) ; 00770 MSRC_B = CCELG_B(0) ; 00771 MSRC_NE = CCELG_NE(0); 00772 MSRC_NW = CCELG_NW(0); 00773 MSRC_SE = CCELG_SE(0); 00774 MSRC_SW = CCELG_SW(0); 00775 MSRC_NT = CCELG_NT(0); 00776 MSRC_NB = CCELG_NB(0); 00777 MSRC_ST = CCELG_ST(0); 00778 MSRC_SB = CCELG_SB(0); 00779 MSRC_ET = CCELG_ET(0); 00780 MSRC_EB = CCELG_EB(0); 00781 MSRC_WT = CCELG_WT(0); 00782 MSRC_WB = CCELG_WB(0); 00783 for(int n=1;(n<this->cDirNum); n++) { 00784 ccel = RACPNT(lev+1, 2*i+1*this->dfVecX[n], 2*j+1*this->dfVecY[n], 2*k+1*this->dfVecZ[n] ,srcSet); 00785 MSRC_C += CCELG_C(n) ; 00786 MSRC_N += CCELG_N(n) ; 00787 MSRC_S += CCELG_S(n) ; 00788 MSRC_E += CCELG_E(n) ; 00789 MSRC_W += CCELG_W(n) ; 00790 MSRC_T += CCELG_T(n) ; 00791 MSRC_B += CCELG_B(n) ; 00792 MSRC_NE += CCELG_NE(n); 00793 MSRC_NW += CCELG_NW(n); 00794 MSRC_SE += CCELG_SE(n); 00795 MSRC_SW += CCELG_SW(n); 00796 MSRC_NT += CCELG_NT(n); 00797 MSRC_NB += CCELG_NB(n); 00798 MSRC_ST += CCELG_ST(n); 00799 MSRC_SB += CCELG_SB(n); 00800 MSRC_ET += CCELG_ET(n); 00801 MSRC_EB += CCELG_EB(n); 00802 MSRC_WT += CCELG_WT(n); 00803 MSRC_WB += CCELG_WB(n); 00804 } 00805 rho = MSRC_C + MSRC_N + MSRC_S + MSRC_E + MSRC_W + MSRC_T 00806 + MSRC_B + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT 00807 + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; 00808 ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW 00809 + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB; 00810 uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW 00811 + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB; 00812 uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB 00813 + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB; 00814 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \ 00815 \ 00816 lcsmeq[dC] = EQC ; \ 00817 COLL_CALCULATE_DFEQ(lcsmeq); \ 00818 COLL_CALCULATE_NONEQTENSOR(lev+0, MSRC_ )\ 00819 COLL_CALCULATE_CSMOMEGAVAL(lev+0, lcsmDstOmega); \ 00820 COLL_CALCULATE_CSMOMEGAVAL(lev+1, lcsmSrcOmega); \ 00821 \ 00822 lcsmdfscale = (mLevel[lev+0].timestep/mLevel[lev+1].timestep)* (1.0/lcsmDstOmega-1.0)/ (1.0/lcsmSrcOmega-1.0); \ 00823 RAC(tcel, dC ) = (lcsmeq[dC ] + (MSRC_C -lcsmeq[dC ] )*lcsmdfscale); 00824 RAC(tcel, dN ) = (lcsmeq[dN ] + (MSRC_N -lcsmeq[dN ] )*lcsmdfscale); 00825 RAC(tcel, dS ) = (lcsmeq[dS ] + (MSRC_S -lcsmeq[dS ] )*lcsmdfscale); 00826 RAC(tcel, dE ) = (lcsmeq[dE ] + (MSRC_E -lcsmeq[dE ] )*lcsmdfscale); 00827 RAC(tcel, dW ) = (lcsmeq[dW ] + (MSRC_W -lcsmeq[dW ] )*lcsmdfscale); 00828 RAC(tcel, dT ) = (lcsmeq[dT ] + (MSRC_T -lcsmeq[dT ] )*lcsmdfscale); 00829 RAC(tcel, dB ) = (lcsmeq[dB ] + (MSRC_B -lcsmeq[dB ] )*lcsmdfscale); 00830 RAC(tcel, dNE) = (lcsmeq[dNE] + (MSRC_NE-lcsmeq[dNE] )*lcsmdfscale); 00831 RAC(tcel, dNW) = (lcsmeq[dNW] + (MSRC_NW-lcsmeq[dNW] )*lcsmdfscale); 00832 RAC(tcel, dSE) = (lcsmeq[dSE] + (MSRC_SE-lcsmeq[dSE] )*lcsmdfscale); 00833 RAC(tcel, dSW) = (lcsmeq[dSW] + (MSRC_SW-lcsmeq[dSW] )*lcsmdfscale); 00834 RAC(tcel, dNT) = (lcsmeq[dNT] + (MSRC_NT-lcsmeq[dNT] )*lcsmdfscale); 00835 RAC(tcel, dNB) = (lcsmeq[dNB] + (MSRC_NB-lcsmeq[dNB] )*lcsmdfscale); 00836 RAC(tcel, dST) = (lcsmeq[dST] + (MSRC_ST-lcsmeq[dST] )*lcsmdfscale); 00837 RAC(tcel, dSB) = (lcsmeq[dSB] + (MSRC_SB-lcsmeq[dSB] )*lcsmdfscale); 00838 RAC(tcel, dET) = (lcsmeq[dET] + (MSRC_ET-lcsmeq[dET] )*lcsmdfscale); 00839 RAC(tcel, dEB) = (lcsmeq[dEB] + (MSRC_EB-lcsmeq[dEB] )*lcsmdfscale); 00840 RAC(tcel, dWT) = (lcsmeq[dWT] + (MSRC_WT-lcsmeq[dWT] )*lcsmdfscale); 00841 RAC(tcel, dWB) = (lcsmeq[dWB] + (MSRC_WB-lcsmeq[dWB] )*lcsmdfscale); 00842 # endif // OPT3D==0 00843 } 00844 00845 void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet, bool markNbs) { 00846 LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0; 00847 LbmFloat intDf[19] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 }; 00848 00849 #if OPT3D==1 00850 // for macro add 00851 LbmFloat addDfFacT, addVal, usqr; 00852 LbmFloat *addfcel, *dstcell; 00853 LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM]; 00854 LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale; 00855 #endif // OPT3D==true 00856 00857 // SET required nbs to from coarse (this might overwrite flag several times) 00858 // this is not necessary for interpolateFineFromCoarse 00859 if(markNbs) { 00860 FORDF1{ 00861 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l]; 00862 if(RFLAG(lev,ni,nj,nk,dstSet)&CFUnused) { 00863 // parents have to be inited! 00864 interpolateCellFromCoarse(lev, ni, nj, nk, dstSet, t, CFFluid|CFGrFromCoarse, false); 00865 } 00866 } } 00867 00868 // change flag of cell to be interpolated 00869 RFLAG(lev,i,j,k, dstSet) = flagSet; 00870 mNumInterdCells++; 00871 00872 // interpolation lines... 00873 int betx = i&1; 00874 int bety = j&1; 00875 int betz = k&1; 00876 00877 if((!betx) && (!bety) && (!betz)) { 00878 ADD_INT_DFS(lev-1, i/2 ,j/2 ,k/2 , 0.0, 1.0); 00879 } 00880 else if(( betx) && (!bety) && (!betz)) { 00881 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2) , t, WO1D1); 00882 ADD_INT_DFS(lev-1, (i/2)+1,(j/2) ,(k/2) , t, WO1D1); 00883 } 00884 else if((!betx) && ( bety) && (!betz)) { 00885 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2) , t, WO1D1); 00886 ADD_INT_DFS(lev-1, (i/2) ,(j/2)+1,(k/2) , t, WO1D1); 00887 } 00888 else if((!betx) && (!bety) && ( betz)) { 00889 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2) , t, WO1D1); 00890 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2)+1, t, WO1D1); 00891 } 00892 else if(( betx) && ( bety) && (!betz)) { 00893 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2) , t, WO1D2); 00894 ADD_INT_DFS(lev-1, (i/2)+1,(j/2) ,(k/2) , t, WO1D2); 00895 ADD_INT_DFS(lev-1, (i/2) ,(j/2)+1,(k/2) , t, WO1D2); 00896 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2) , t, WO1D2); 00897 } 00898 else if((!betx) && ( bety) && ( betz)) { 00899 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2) , t, WO1D2); 00900 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2)+1, t, WO1D2); 00901 ADD_INT_DFS(lev-1, (i/2) ,(j/2)+1,(k/2) , t, WO1D2); 00902 ADD_INT_DFS(lev-1, (i/2) ,(j/2)+1,(k/2)+1, t, WO1D2); 00903 } 00904 else if(( betx) && (!bety) && ( betz)) { 00905 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2) , t, WO1D2); 00906 ADD_INT_DFS(lev-1, (i/2)+1,(j/2) ,(k/2) , t, WO1D2); 00907 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2)+1, t, WO1D2); 00908 ADD_INT_DFS(lev-1, (i/2)+1,(j/2) ,(k/2)+1, t, WO1D2); 00909 } 00910 else if(( betx) && ( bety) && ( betz)) { 00911 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2) , t, WO1D3); 00912 ADD_INT_DFS(lev-1, (i/2)+1,(j/2) ,(k/2) , t, WO1D3); 00913 ADD_INT_DFS(lev-1, (i/2) ,(j/2) ,(k/2)+1, t, WO1D3); 00914 ADD_INT_DFS(lev-1, (i/2)+1,(j/2) ,(k/2)+1, t, WO1D3); 00915 ADD_INT_DFS(lev-1, (i/2) ,(j/2)+1,(k/2) , t, WO1D3); 00916 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2) , t, WO1D3); 00917 ADD_INT_DFS(lev-1, (i/2) ,(j/2)+1,(k/2)+1, t, WO1D3); 00918 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)+1, t, WO1D3); 00919 } 00920 else { 00921 CAUSE_PANIC; 00922 errFatal("interpolateCellFromCoarse","Invalid!?", SIMWORLD_GENERICERROR); 00923 } 00924 00925 IDF_WRITEBACK; 00926 return; 00927 } 00928 00929 00930 00931 // required globals 00932 extern bool glob_mpactive; 00933 extern int glob_mpnum, glob_mpindex; 00934 #define MPTADAP_INTERV 4 00935 00936 /*****************************************************************************/ 00938 /*****************************************************************************/ 00939 void LbmFsgrSolver::adaptTimestep() { 00940 LbmFloat massTOld=0.0, massTNew=0.0; 00941 LbmFloat volTOld=0.0, volTNew=0.0; 00942 00943 bool rescale = false; // do any rescale at all? 00944 LbmFloat scaleFac = -1.0; // timestep scaling 00945 if(mPanic) return; 00946 00947 LbmFloat levOldOmega[FSGR_MAXNOOFLEVELS]; 00948 LbmFloat levOldStepsize[FSGR_MAXNOOFLEVELS]; 00949 for(int lev=mMaxRefine; lev>=0 ; lev--) { 00950 levOldOmega[lev] = mLevel[lev].omega; 00951 levOldStepsize[lev] = mLevel[lev].timestep; 00952 } 00953 00954 const LbmFloat reduceFac = 0.8; // modify time step by 20%, TODO? do multiple times for large changes? 00955 LbmFloat diffPercent = 0.05; // dont scale if less than 5% 00956 LbmFloat allowMax = mpParam->getTadapMaxSpeed(); // maximum allowed velocity 00957 LbmFloat nextmax = mpParam->getSimulationMaxSpeed() + norm(mLevel[mMaxRefine].gravity); 00958 00959 // sync nextmax 00960 #if LBM_INCLUDE_TESTSOLVERS==1 00961 if(glob_mpactive) { 00962 if(mLevel[mMaxRefine].lsteps % MPTADAP_INTERV != MPTADAP_INTERV-1) { 00963 debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "mpact:"<<glob_mpactive<<","<<glob_mpindex<<"/"<<glob_mpnum<<" step:"<<mLevel[mMaxRefine].lsteps<<" skipping tadap...",8); 00964 return; 00965 } 00966 nextmax = mrInitTadap(nextmax); 00967 } 00968 #endif // LBM_INCLUDE_TESTSOLVERS==1 00969 00970 LbmFloat newdt = mpParam->getTimestep(); // newtr 00971 if(nextmax > allowMax/reduceFac) { 00972 mTimeMaxvelStepCnt++; } 00973 else { mTimeMaxvelStepCnt=0; } 00974 00975 // emergency, or 10 steps with high vel 00976 if((mTimeMaxvelStepCnt>5) || (nextmax> (1.0/3.0)) || (mForceTimeStepReduce) ) { 00977 newdt = mpParam->getTimestep() * reduceFac; 00978 } else { 00979 if(nextmax<allowMax*reduceFac) { 00980 newdt = mpParam->getTimestep() / reduceFac; 00981 } 00982 } // newtr 00983 //errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<reduceFac<<" simmaxv="<< mpParam->getSimulationMaxSpeed() ); 00984 00985 bool minCutoff = false; 00986 LbmFloat desireddt = newdt; 00987 if(newdt>mpParam->getMaxTimestep()){ newdt = mpParam->getMaxTimestep(); } 00988 if(newdt<mpParam->getMinTimestep()){ 00989 newdt = mpParam->getMinTimestep(); 00990 if(nextmax>allowMax/reduceFac){ minCutoff=true; } // only if really large vels... 00991 } 00992 00993 LbmFloat dtdiff = fabs(newdt - mpParam->getTimestep()); 00994 if(!mSilent) { 00995 debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt 00996 <<" max"<<mpParam->getMaxTimestep()<<" min"<<mpParam->getMinTimestep()<<" diff"<<dtdiff 00997 <<" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep)<< 00998 " olddt"<<levOldStepsize[mMaxRefine]<<" redlock"<<mTimestepReduceLock 00999 , 10); } 01000 01001 // in range, and more than X% change? 01002 //if( newdt < mpParam->getTimestep() ) // DEBUG 01003 LbmFloat rhoAvg = mCurrentMass/mCurrentVolume; 01004 if( (newdt<=mpParam->getMaxTimestep()) && (newdt>=mpParam->getMinTimestep()) 01005 && (dtdiff>(mpParam->getTimestep()*diffPercent)) ) { 01006 if((newdt>levOldStepsize[mMaxRefine])&&(mTimestepReduceLock)) { 01007 // wait some more... 01008 //debMsgNnl("LbmFsgrSolver::TAdp",DM_NOTIFY," Delayed... "<<mTimestepReduceLock<<" ",10); 01009 debMsgDirect("D"); 01010 } else { 01011 mpParam->setDesiredTimestep( newdt ); 01012 rescale = true; 01013 if(!mSilent) { 01014 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"\n\n\n\n",10); 01015 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: new="<<newdt<<" old="<<mpParam->getTimestep() 01016 <<" maxSpeed:"<<mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<mStepCnt, 10 ); 01017 //debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: "<< "rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10); 01018 } 01019 } // really change dt 01020 } 01021 01022 if(mTimestepReduceLock>0) mTimestepReduceLock--; 01023 01024 01025 // forced back and forth switchting (for testing) 01026 /*const int tadtogInter = 100; 01027 const double tadtogSwitch = 0.66; 01028 errMsg("TIMESWITCHTOGGLETEST","warning enabled "<< tadtogSwitch<<","<<tadtogSwitch<<" !!!!!!!!!!!!!!!!!!!"); 01029 if( ((mStepCnt% tadtogInter)== (tadtogInter/4*1)-1) || 01030 ((mStepCnt% tadtogInter)== (tadtogInter/4*2)-1) ){ 01031 rescale = true; minCutoff = false; 01032 newdt = tadtogSwitch * mpParam->getTimestep(); 01033 mpParam->setDesiredTimestep( newdt ); 01034 } else 01035 if( ((mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) || 01036 ((mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){ 01037 rescale = true; minCutoff = false; 01038 newdt = mpParam->getTimestep()/tadtogSwitch ; 01039 mpParam->setDesiredTimestep( newdt ); 01040 } else { 01041 rescale = false; minCutoff = false; 01042 } 01043 // */ 01044 01045 // test mass rescale 01046 01047 scaleFac = newdt/mpParam->getTimestep(); 01048 if(rescale) { 01049 // perform acutal rescaling... 01050 mTimeMaxvelStepCnt=0; 01051 mForceTimeStepReduce = false; 01052 01053 // FIXME - approximate by averaging, take gravity direction here? 01054 //mTimestepReduceLock = 4*(mLevel[mMaxRefine].lSizey+mLevel[mMaxRefine].lSizez+mLevel[mMaxRefine].lSizex)/3; 01055 // use z as gravity direction 01056 mTimestepReduceLock = 4*mLevel[mMaxRefine].lSizez; 01057 01058 mTimeSwitchCounts++; 01059 mpParam->calculateAllMissingValues( mSimulationTime, mSilent ); 01060 recalculateObjectSpeeds(); 01061 // calc omega, force for all levels 01062 mLastOmega=1e10; mLastGravity=1e10; 01063 initLevelOmegas(); 01064 if(mpParam->getTimestep()<mMinTimestep) mMinTimestep = mpParam->getTimestep(); 01065 if(mpParam->getTimestep()>mMaxTimestep) mMaxTimestep = mpParam->getTimestep(); 01066 01067 // this might be called during init, before we have any particles 01068 if(mpParticles) { mpParticles->adaptPartTimestep(scaleFac); } 01069 #if LBM_INCLUDE_TESTSOLVERS==1 01070 if((mUseTestdata)&&(mpTest)) { 01071 mpTest->adaptTimestep(scaleFac, mLevel[mMaxRefine].omega, mLevel[mMaxRefine].timestep, vec2L( mpParam->calculateGravity(mSimulationTime)) ); 01072 mpTest->mGrav3d = mLevel[mMaxRefine].gravity; 01073 } 01074 #endif // LBM_INCLUDE_TESTSOLVERS!=1 01075 01076 for(int lev=mMaxRefine; lev>=0 ; lev--) { 01077 LbmFloat newSteptime = mLevel[lev].timestep; 01078 LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]); 01079 01080 if(!mSilent) { 01081 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Level: "<<lev<<" Timestep chlevel: "<< 01082 " scaleFac="<<dfScaleFac<<" newDt="<<newSteptime<<" newOmega="<<mLevel[lev].omega,10); 01083 } 01084 if(lev!=mMaxRefine) coarseCalculateFluxareas(lev); 01085 01086 int wss = 0, wse = 1; 01087 // only change currset (necessary for compressed grids, better for non-compr.gr.) 01088 wss = wse = mLevel[lev].setCurr; 01089 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT 01090 // warning - check sets for higher levels...? 01091 FSGR_FORIJK1(lev) { 01092 if( (RFLAG(lev,i,j,k, workSet) & CFBndMoving) ) { 01093 /* 01094 // paranoid check - shouldnt be necessary! 01095 if(QCELL(lev, i, j, k, workSet, dFlux)!=mSimulationTime) { 01096 errMsg("TTTT","found invalid bnd cell.... removing at "<<PRINT_IJK); 01097 RFLAG(lev,i,j,k, workSet) = CFInter; 01098 // init empty zero vel interface cell... 01099 initVelocityCell(lev, i,j,k, CFInter, 1.0, 0.01, LbmVec(0.) ); 01100 } else {// */ 01101 for(int l=0; l<cDfNum; l++) { 01102 QCELL(lev, i, j, k, workSet, l) = QCELL(lev, i, j, k, workSet, l)* scaleFac; 01103 } 01104 //} // ok 01105 continue; 01106 } 01107 if( 01108 (RFLAG(lev,i,j,k, workSet) & CFFluid) || 01109 (RFLAG(lev,i,j,k, workSet) & CFInter) || 01110 (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 01111 (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 01112 (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 01113 ) { 01114 // these cells have to be scaled... 01115 } else { 01116 continue; 01117 } 01118 01119 // collide on current set 01120 LbmFloat rhoOld; 01121 LbmVec velOld; 01122 LbmFloat rho, ux,uy,uz; 01123 rho=0.0; ux = uy = uz = 0.0; 01124 for(int l=0; l<cDfNum; l++) { 01125 LbmFloat m = QCELL(lev, i, j, k, workSet, l); 01126 rho += m; 01127 ux += (dfDvecX[l]*m); 01128 uy += (dfDvecY[l]*m); 01129 uz += (dfDvecZ[l]*m); 01130 } 01131 rhoOld = rho; 01132 velOld = LbmVec(ux,uy,uz); 01133 01134 LbmFloat rhoNew = (rhoOld-rhoAvg)*scaleFac +rhoAvg; 01135 LbmVec velNew = velOld * scaleFac; 01136 01137 LbmFloat df[LBM_DFNUM]; 01138 LbmFloat feqOld[LBM_DFNUM]; 01139 LbmFloat feqNew[LBM_DFNUM]; 01140 for(int l=0; l<cDfNum; l++) { 01141 feqOld[l] = getCollideEq(l,rhoOld, velOld[0],velOld[1],velOld[2] ); 01142 feqNew[l] = getCollideEq(l,rhoNew, velNew[0],velNew[1],velNew[2] ); 01143 df[l] = QCELL(lev, i,j,k,workSet, l); 01144 } 01145 const LbmFloat Qo = getLesNoneqTensorCoeff(df,feqOld); 01146 const LbmFloat oldOmega = getLesOmega(levOldOmega[lev], mLevel[lev].lcsmago,Qo); 01147 const LbmFloat newOmega = getLesOmega(mLevel[lev].omega,mLevel[lev].lcsmago,Qo); 01148 //newOmega = mLevel[lev].omega; // FIXME debug test 01149 01150 //LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]); 01151 const LbmFloat dfScale = (newSteptime/newOmega)/(levOldStepsize[lev]/oldOmega); 01152 //dfScale = dfScaleFac/newOmega; 01153 01154 for(int l=0; l<cDfNum; l++) { 01155 // org scaling 01156 //df = eqOld + (df-eqOld)*dfScale; df *= (eqNew/eqOld); // non-eq. scaling, important 01157 // new scaling 01158 LbmFloat dfn = feqNew[l] + (df[l]-feqOld[l])*dfScale*feqNew[l]/feqOld[l]; // non-eq. scaling, important 01159 //df = eqNew + (df-eqOld)*dfScale; // modified ig scaling, no real difference? 01160 QCELL(lev, i,j,k,workSet, l) = dfn; 01161 } 01162 01163 if(RFLAG(lev,i,j,k, workSet) & CFInter) { 01164 //if(workSet==mLevel[lev].setCurr) 01165 LbmFloat area = 1.0; 01166 if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux); 01167 massTOld += QCELL(lev, i,j,k,workSet, dMass) * area; 01168 volTOld += QCELL(lev, i,j,k,workSet, dFfrac); 01169 01170 // wrong... QCELL(i,j,k,workSet, dMass] = (QCELL(i,j,k,workSet, dFfrac]*rhoNew); 01171 QCELL(lev, i,j,k,workSet, dMass) = (QCELL(lev, i,j,k,workSet, dMass)/rhoOld*rhoNew); 01172 QCELL(lev, i,j,k,workSet, dFfrac) = (QCELL(lev, i,j,k,workSet, dMass)/rhoNew); 01173 01174 //if(workSet==mLevel[lev].setCurr) 01175 massTNew += QCELL(lev, i,j,k,workSet, dMass); 01176 volTNew += QCELL(lev, i,j,k,workSet, dFfrac); 01177 } 01178 if(RFLAG(lev,i,j,k, workSet) & CFFluid) { // DEBUG 01179 if(RFLAG(lev,i,j,k, workSet) & (CFGrFromFine|CFGrFromCoarse)) { // DEBUG 01180 // dont include 01181 } else { 01182 LbmFloat area = 1.0; 01183 if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux) * mLevel[lev].lcellfactor; 01184 //if(workSet==mLevel[lev].setCurr) 01185 massTOld += rhoOld*area; 01186 //if(workSet==mLevel[lev].setCurr) 01187 massTNew += rhoNew*area; 01188 volTOld += area; 01189 volTNew += area; 01190 } 01191 } 01192 01193 } // IJK 01194 } // workSet 01195 01196 } // lev 01197 01198 if(!mSilent) { 01199 debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<mStepCnt<< 01200 " no"<<mTimeSwitchCounts<<" maxdt"<<mMaxTimestep<< 01201 " mindt"<<mMinTimestep<<" currdt"<<mLevel[mMaxRefine].timestep, 10); 01202 debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE masst:"<<massTNew<<","<<massTOld<<" org:"<<mCurrentMass<<"; "<< 01203 " volt:"<<volTNew<<","<<volTOld<<" org:"<<mCurrentVolume, 10); 01204 } else { 01205 debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep changed by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<<newdt 01206 <<", oldDt:"<<levOldStepsize[mMaxRefine]<<" newOmega:"<<mOmega<<" gStar:"<<mpParam->getCurrentGStar()<<", step:"<<mStepCnt , 10); 01207 } 01208 } // rescale? 01209 //NEWDEBCHECK("tt2"); 01210 01211 //errMsg("adaptTimestep","Warning - brute force rescale off!"); minCutoff = false; // DEBUG 01212 if(minCutoff) { 01213 errMsg("adaptTimestep","Warning - performing Brute-Force rescale... (sim:"<<mName<<" step:"<<mStepCnt<<" newdt="<<desireddt<<" mindt="<<mpParam->getMinTimestep()<<") " ); 01214 //brute force resacle all the time? 01215 01216 for(int lev=mMaxRefine; lev>=0 ; lev--) { 01217 int rescs=0; 01218 int wss = 0, wse = 1; 01219 #if COMPRESSGRIDS==1 01220 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr; 01221 #endif // COMPRESSGRIDS==1 01222 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT 01223 //for(int workSet = 0; workSet<=1; workSet++) { 01224 FSGR_FORIJK1(lev) { 01225 01226 //if( (RFLAG(lev, i,j,k, workSet) & CFFluid) || (RFLAG(lev, i,j,k, workSet) & CFInter) ) { 01227 if( 01228 (RFLAG(lev,i,j,k, workSet) & CFFluid) || 01229 (RFLAG(lev,i,j,k, workSet) & CFInter) || 01230 (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 01231 (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 01232 (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 01233 ) { 01234 // these cells have to be scaled... 01235 } else { 01236 continue; 01237 } 01238 01239 // collide on current set 01240 LbmFloat rho, ux,uy,uz; 01241 rho=0.0; ux = uy = uz = 0.0; 01242 for(int l=0; l<cDfNum; l++) { 01243 LbmFloat m = QCELL(lev, i, j, k, workSet, l); 01244 rho += m; 01245 ux += (dfDvecX[l]*m); 01246 uy += (dfDvecY[l]*m); 01247 uz += (dfDvecZ[l]*m); 01248 } 01249 #ifndef WIN32 01250 if (!finite(rho)) { 01251 errMsg("adaptTimestep","Brute force non-finite rho at"<<PRINT_IJK); // DEBUG! 01252 rho = 1.0; 01253 ux = uy = uz = 0.0; 01254 QCELL(lev, i, j, k, workSet, dMass) = 1.0; 01255 QCELL(lev, i, j, k, workSet, dFfrac) = 1.0; 01256 } 01257 #endif // WIN32 01258 01259 if( (ux*ux+uy*uy+uz*uz)> (allowMax*allowMax) ) { 01260 LbmFloat cfac = allowMax/sqrt(ux*ux+uy*uy+uz*uz); 01261 ux *= cfac; 01262 uy *= cfac; 01263 uz *= cfac; 01264 for(int l=0; l<cDfNum; l++) { 01265 QCELL(lev, i, j, k, workSet, l) = getCollideEq(l, rho, ux,uy,uz); } 01266 rescs++; 01267 debMsgDirect("B"); 01268 } 01269 01270 } } 01271 //if(rescs>0) { errMsg("adaptTimestep","!!!!! Brute force rescaling was necessary !!!!!!!"); } 01272 debMsgStd("adaptTimestep",DM_MSG,"Brute force rescale done. level:"<<lev<<" rescs:"<<rescs, 1); 01273 //TTT mNumProblems += rescs; // add to problem display... 01274 } // lev,set,ijk 01275 01276 } // try brute force rescale? 01277 01278 // time adap done... 01279 } 01280 01281 01282