This version now includes the fluid control sources, however the Blender
[blender.git] / intern / elbeem / intern / solver_adap.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2003-2006 Nils Thuerey
5  *
6  * Adaptivity functions
7  *
8  *****************************************************************************/
9
10 #include "solver_class.h"
11 #include "solver_relax.h"
12 #include "particletracer.h"
13
14
15
16 /*****************************************************************************/
17 //! coarse step functions
18 /*****************************************************************************/
19
20
21
22 void LbmFsgrSolver::coarseCalculateFluxareas(int lev)
23 {
24         FSGR_FORIJK_BOUNDS(lev) {
25                 if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
26                         if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
27                                 LbmFloat totArea = mFsgrCellArea[0]; // for l=0
28                                 for(int l=1; l<this->cDirNum; l++) { 
29                                         int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l];
30                                         if(RFLAG(lev+1, ni,nj,nk, mLevel[lev+1].setCurr)&
31                                                         (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
32                                                         ) { 
33                                                 totArea += mFsgrCellArea[l];
34                                         }
35                                 } // l
36                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = totArea;
37                                 //continue;
38                         } else
39                         if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & (CFEmpty|CFUnused)) {
40                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 1.0;
41                                 //continue;
42                         } else {
43                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 0.0;
44                         }
45                 //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) ); 
46                 }
47         } // } TEST DEBUG
48         if(!this->mSilent){ debMsgStd("coarseCalculateFluxareas",DM_MSG,"level "<<lev<<" calculated", 7); }
49 }
50         
51 void LbmFsgrSolver::coarseAdvance(int lev)
52 {
53         LbmFloat calcCurrentMass = 0.0;
54         LbmFloat calcCurrentVolume = 0.0;
55
56         LbmFloat *ccel = NULL;
57         LbmFloat *tcel = NULL;
58         LbmFloat m[LBM_DFNUM];
59         LbmFloat rho, ux, uy, uz, tmp, usqr, lcsmqo;
60 #if OPT3D==1 
61         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
62 #endif // OPT3D==true 
63         m[0] = tmp = usqr = 0.0;
64
65         coarseCalculateFluxareas(lev);
66         // copied from fineAdv.
67         CellFlagType *pFlagSrc = &RFLAG(lev, 1,1,getForZMin1(),SRCS(lev));
68         CellFlagType *pFlagDst = &RFLAG(lev, 1,1,getForZMin1(),TSET(lev));
69         pFlagSrc -= 1;
70         pFlagDst -= 1;
71         ccel = RACPNT(lev, 1,1,getForZMin1() ,SRCS(lev)); // QTEST
72         ccel -= QCELLSTEP;
73         tcel = RACPNT(lev, 1,1,getForZMin1() ,TSET(lev)); // QTEST
74         tcel -= QCELLSTEP;
75         //if(strstr(this->getName().c_str(),"Debug")){ errMsg("DEBUG","DEBUG!!!!!!!!!!!!!!!!!!!!!!!"); }
76
77         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
78   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
79   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
80 #if FSGR_STRICT_DEBUG==1
81                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
82 #endif
83                 pFlagSrc++;
84                 pFlagDst++;
85                 ccel += QCELLSTEP;
86                 tcel += QCELLSTEP;
87
88                 // from coarse cells without unused nbs are not necessary...! -> remove
89                 if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
90                         bool invNb = false;
91                         FORDF1 { if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; } }   
92                         if(!invNb) {
93                                 // WARNING - this modifies source flag array...
94                                 *pFlagSrc = CFFluid|CFGrNorm;
95 #if ELBEEM_PLUGIN!=1
96                                 errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK);
97 #endif // ELBEEM_PLUGIN!=1
98                                 // move to perform coarsening?
99                         }
100                 } // */
101
102 #if FSGR_STRICT_DEBUG==1
103                 *pFlagDst = *pFlagSrc; // always set other set...
104 #else
105                 *pFlagDst = (*pFlagSrc & (~CFGrCoarseInited)); // always set other set... , remove coarse inited flag
106 #endif
107
108                 // old INTCFCOARSETEST==1
109                 if((*pFlagSrc) & CFGrFromCoarse) {  // interpolateFineFromCoarse test!
110                         if(( this->mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
111                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
112                         } else {
113                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
114                                 this->mNumUsedCells++;
115                         }
116                         continue; // interpolateFineFromCoarse test!
117                 } // interpolateFineFromCoarse test! old INTCFCOARSETEST==1
118
119                 if( ((*pFlagSrc) & (CFFluid)) ) { 
120                         ccel = RACPNT(lev, i,j,k ,SRCS(lev)); 
121                         tcel = RACPNT(lev, i,j,k ,TSET(lev));
122
123                         if( ((*pFlagSrc) & (CFGrFromFine)) ) { 
124                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }    // always copy...?
125                                 continue; // comes from fine grid
126                         }
127                         // also ignore CFGrFromCoarse
128                         else if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
129                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }    // always copy...?
130                                 continue; 
131                         }
132
133                         OPTIMIZED_STREAMCOLLIDE;
134                         *pFlagDst |= CFNoBndFluid; // test?
135                         calcCurrentVolume += RAC(ccel,dFlux); 
136                         calcCurrentMass   += RAC(ccel,dFlux)*rho;
137                         //ebugMarkCell(lev+1, 2*i+1,2*j+1,2*k  );
138 #if FSGR_STRICT_DEBUG==1
139                         if(rho<-1.0){ debugMarkCell(lev, i,j,k ); 
140                                 errMsg("INVRHOCELL_CHECK"," l"<<lev<<" "<< PRINT_IJK<<" rho:"<<rho ); 
141                                 CAUSE_PANIC;
142                         }
143 #endif // FSGR_STRICT_DEBUG==1
144                         this->mNumUsedCells++;
145
146                 }
147         } 
148         pFlagSrc+=2; // after x
149         pFlagDst+=2; // after x
150         ccel += (QCELLSTEP*2);
151         tcel += (QCELLSTEP*2);
152         } 
153         pFlagSrc+= mLevel[lev].lSizex*2; // after y
154         pFlagDst+= mLevel[lev].lSizex*2; // after y
155         ccel += (QCELLSTEP*mLevel[lev].lSizex*2);
156         tcel += (QCELLSTEP*mLevel[lev].lSizex*2);
157         } // all cell loop k,j,i
158         
159
160         //errMsg("coarseAdvance","level "<<lev<<" stepped from "<<mLevel[lev].setCurr<<" to "<<mLevel[lev].setOther);
161         if(!this->mSilent){ errMsg("coarseAdvance","level "<<lev<<" stepped from "<<SRCS(lev)<<" to "<<TSET(lev)); }
162         // */
163
164         // update other set
165   mLevel[lev].setOther   = mLevel[lev].setCurr;
166   mLevel[lev].setCurr   ^= 1;
167   mLevel[lev].lsteps++;
168   mLevel[lev].lmass   = calcCurrentMass   * mLevel[lev].lcellfactor;
169   mLevel[lev].lvolume = calcCurrentVolume * mLevel[lev].lcellfactor;
170 #if ELBEEM_PLUGIN!=1
171   debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "mass: lev="<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<"  lcf="<< mLevel[lev].lcellfactor, 8 );
172   debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "volume: lev="<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<"  lcf="<< mLevel[lev].lcellfactor, 8 );
173 #endif // ELBEEM_PLUGIN
174 }
175
176 /*****************************************************************************/
177 //! multi level functions
178 /*****************************************************************************/
179
180
181 // get dfs from level (lev+1) to (lev) coarse border nodes
182 void LbmFsgrSolver::coarseRestrictFromFine(int lev)
183 {
184         if((lev<0) || ((lev+1)>mMaxRefine)) return;
185 #if FSGR_STRICT_DEBUG==1
186         // reset all unused cell values to invalid
187         int unuCnt = 0;
188         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
189         for(int j=1;j<mLevel[lev].lSizey-1;++j) {
190         for(int i=1;i<mLevel[lev].lSizex-1;++i) {
191                 CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,mLevel[lev].setCurr);
192                 if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 
193                         FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0;    }
194                         unuCnt++;
195                         // set here
196                 } else if( ((*pFlagSrc) & (CFFluid|CFGrNorm)) == (CFFluid|CFGrNorm) ) { 
197                         // simulated...
198                 } else {
199                         // reset in interpolation
200                         //errMsg("coarseRestrictFromFine"," reset l"<<lev<<" "<<PRINT_IJK);
201                 }
202                 if( ((*pFlagSrc) & (CFEmpty|CFUnused)) ) {  // test, also reset?
203                         FORDF0{ QCELL(lev, i,j,k,mLevel[lev].setCurr, l) = -10000.0;    }
204                 } // test
205         } } }
206         errMsg("coarseRestrictFromFine"," reset l"<<lev<<" fluid|coarseBorder cells: "<<unuCnt);
207 #endif // FSGR_STRICT_DEBUG==1
208         const int srcSet = mLevel[lev+1].setCurr;
209         const int dstSet = mLevel[lev].setCurr;
210
211         //restrict
212         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
213         for(int j=1;j<mLevel[lev].lSizey-1;++j) {
214         for(int i=1;i<mLevel[lev].lSizex-1;++i) {
215                 CellFlagType *pFlagSrc = &RFLAG(lev, i,j,k,dstSet);
216                 if((*pFlagSrc) & (CFFluid)) { 
217                         if( ((*pFlagSrc) & (CFFluid|CFGrFromFine)) == (CFFluid|CFGrFromFine) ) { 
218                                 // do resctriction
219                                 mNumInterdCells++;
220                                 coarseRestrictCell(lev, i,j,k,srcSet,dstSet);
221
222                                 this->mNumUsedCells++;
223                         } // from fine & fluid
224                         else {
225                                 if(RFLAG(lev+1, 2*i,2*j,2*k,srcSet) & CFGrFromCoarse) {
226                                         RFLAG(lev, i,j,k,dstSet) |= CFGrToFine;
227                                 } else {
228                                         RFLAG(lev, i,j,k,dstSet) &= (~CFGrToFine);
229                                 }
230                         }
231                 } // & fluid
232         }}}
233         if(!this->mSilent){ errMsg("coarseRestrictFromFine"," from l"<<(lev+1)<<",s"<<mLevel[lev+1].setCurr<<" to l"<<lev<<",s"<<mLevel[lev].setCurr); }
234 }
235
236 bool LbmFsgrSolver::adaptGrid(int lev) {
237         if((lev<0) || ((lev+1)>mMaxRefine)) return false;
238         bool change = false;
239         { // refinement, PASS 1-3
240
241         //bool nbsok;
242         // FIXME remove TIMEINTORDER ?
243         LbmFloat interTime = 0.0;
244         // update curr from other, as streaming afterwards works on curr
245         // thus read only from srcSet, modify other
246         const int srcSet = mLevel[lev].setOther;
247         const int dstSet = mLevel[lev].setCurr;
248         const int srcFineSet = mLevel[lev+1].setCurr;
249         const bool debugRefinement = false;
250
251         // use //template functions for 2D/3D
252                         /*if(strstr(this->getName().c_str(),"Debug"))
253                         if(lev+1==mMaxRefine) { // mixborder
254                                 for(int l=0;((l<this->cDirNum) && (!removeFromFine)); l++) {  // FARBORD
255                                         int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l];
256                                         if(RFLAG(lev+1, ni,nj,nk, srcFineSet)&CFBnd) { // NEWREFT
257                                                 removeFromFine=true;
258                                         }
259                                 }
260                         } // FARBORD */
261         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
262   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
263   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
264
265                 if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
266                         bool removeFromFine = false;
267                         const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
268                         CellFlagType reqType = CFGrNorm;
269                         if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
270                         
271                         if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & reqType) &&
272                             (!(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet) & (notAllowed)) )  ){ // ok
273                         } else {
274                                 removeFromFine=true;
275                         }
276
277                         if(removeFromFine) {
278                                 // dont turn CFGrFromFine above interface cells into CFGrNorm
279                                 //errMsg("performRefinement","Removing CFGrFromFine on lev"<<lev<<" " <<PRINT_IJK<<" srcflag:"<<convertCellFlagType2String(RFLAG(lev+1, (2*i),(2*j),(2*k), srcFineSet)) <<" set:"<<dstSet );
280                                 RFLAG(lev, i,j,k, dstSet) = CFEmpty;
281 #if FSGR_STRICT_DEBUG==1
282                                 // for interpolation later on during fine grid fixing
283                                 // these cells are still correctly inited
284                                 RFLAG(lev, i,j,k, dstSet) |= CFGrCoarseInited;  // remove later on? FIXME?
285 #endif // FSGR_STRICT_DEBUG==1
286                                 //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFEmpty; // FLAGTEST
287                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,i,j,k); 
288                                 change=true;
289                                 mNumFsgrChanges++;
290                                 for(int l=1; l<this->cDirNum; l++) { 
291                                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
292                                         //errMsg("performRefinement","On lev:"<<lev<<" check: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet<<" = "<<convertCellFlagType2String(RFLAG(lev, ni,nj,nk, srcSet)) );
293                                         if( (  RFLAG(lev, ni,nj,nk, srcSet)&CFFluid      ) &&
294                                                         (!(RFLAG(lev, ni,nj,nk, srcSet)&CFGrFromFine)) ) { // dont change status of nb. from fine cells
295                                                 // tag as inited for debugging, cell contains fluid DFs anyway
296                                                 RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromFine|CFGrCoarseInited;
297                                                 //errMsg("performRefinement","On lev:"<<lev<<" set to from fine: "<<PRINT_VEC(ni,nj,nk)<<" set:"<<dstSet);
298                                                 //if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
299                                         }
300                                 } // l 
301
302                                 // FIXME fix fine level?
303                         }
304
305                         // recheck from fine flag
306                 }
307         }}} // TEST
308         // PASS 1 */
309
310
311                 /*if( ((*pFlagSrc) & (CFGrFromCoarse)) ) { 
312                         bool invNb = false;
313                         FORDF1 { 
314                                 if(RFLAG_NB(lev, i, j, k, SRCS(lev), l) & CFUnused) { invNb = true; }
315                         }   
316                         if(!invNb) {
317                                 *pFlagSrc = CFFluid|CFGrNorm;
318                                 errMsg("coarseAdvance","FC2NRM_CHECK Converted CFGrFromCoarse to Norm at "<<lev<<" "<<PRINT_IJK);
319                         }
320                 } // */
321         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST
322   for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST
323   for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST
324
325                 // test from coarseAdvance
326                 // from coarse cells without unused nbs are not necessary...! -> remove
327
328                 if(RFLAG(lev, i,j,k, srcSet) & CFGrFromCoarse) {
329
330                         // from coarse cells without unused nbs are not necessary...! -> remove
331                         bool invNb = false;
332                         bool fluidNb = false;
333                         for(int l=1; l<this->cDirNum; l++) { 
334                                 if(RFLAG_NB(lev, i, j, k, srcSet, l) & CFUnused) { invNb = true; }
335                                 if(RFLAG_NB(lev, i, j, k, srcSet, l) & (CFGrNorm)) { fluidNb = true; }
336                         }   
337                         if(!invNb) {
338                                 // no unused cells around -> calculate normally from now on
339                                 RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm;
340                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
341                                 change=true;
342                                 mNumFsgrChanges++;
343                         } // from advance 
344                         if(!fluidNb) {
345                                 // no fluid cells near -> no transfer necessary
346                                 RFLAG(lev, i,j,k, dstSet) = CFUnused;
347                                 //RFLAG(lev, i,j,k, mLevel[lev].setOther) = CFUnused; // FLAGTEST
348                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
349                                 change=true;
350                                 mNumFsgrChanges++;
351                         } // from advance 
352
353
354                         // dont allow double transfer
355                         // this might require fixing the neighborhood
356                         if(RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse)) { 
357                                 // dont turn CFGrFromFine above interface cells into CFGrNorm
358                                 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<lev<<" " <<PRINT_IJK<<" due to finer from coarse cell " );
359                                 RFLAG(lev, i,j,k, dstSet) = CFFluid|CFGrNorm;
360                                 if(lev>0) RFLAG(lev-1, i/2,j/2,k/2, mLevel[lev-1].setCurr) &= (~CFGrToFine); // TODO add more of these?
361                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev, i, j, k); 
362                                 change=true;
363                                 mNumFsgrChanges++;
364                                 for(int l=1; l<this->cDirNum; l++) { 
365                                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
366                                         if(RFLAG(lev, ni,nj,nk, srcSet)&(CFGrNorm)) { //ok
367                                                 for(int m=1; m<this->cDirNum; m++) { 
368                                                         int mi=  ni +this->dfVecX[m], mj=  nj +this->dfVecY[m], mk=  nk +this->dfVecZ[m];
369                                                         if(RFLAG(lev,  mi, mj, mk, srcSet)&CFUnused) {
370                                                                 // norm cells in neighborhood with unused nbs have to be new border...
371                                                                 RFLAG(lev, ni,nj,nk, dstSet) = CFFluid|CFGrFromCoarse;
372                                                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
373                                                         }
374                                                 }
375                                                 // these alreay have valid values...
376                                         }
377                                         else if(RFLAG(lev, ni,nj,nk, srcSet)&(CFUnused)) { //ok
378                                                 // this should work because we have a valid neighborhood here for now
379                                                 interpolateCellFromCoarse(lev,  ni, nj, nk, dstSet, interTime, CFFluid|CFGrFromCoarse, false);
380                                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev,ni,nj,nk); 
381                                                 mNumFsgrChanges++;
382                                         }
383                                 } // l 
384                         } // double transer
385
386                 } // from coarse
387
388         } } }
389         // PASS 2 */
390
391
392         // fix dstSet from fine cells here
393         // warning - checks CFGrFromFine on dstset changed before!
394         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) { // TEST
395   for(int j=1;j<mLevel[lev].lSizey-1;++j) { // TEST
396   for(int i=1;i<mLevel[lev].lSizex-1;++i) { // TEST
397
398                 //if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
399                 if(RFLAG(lev, i,j,k, dstSet) & CFGrFromFine) {
400                         // modify finer level border
401                         if((RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)&(CFGrFromCoarse))) { 
402                                 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" from l"<<lev<<" " <<PRINT_IJK );
403                                 CellFlagType setf = CFFluid;
404                                 if(lev+1 < mMaxRefine) setf = CFFluid|CFGrNorm;
405                                 RFLAG(lev+1, 2*i,2*j,2*k, srcFineSet)=setf;
406                                 change=true;
407                                 mNumFsgrChanges++;
408                                 for(int l=1; l<this->cDirNum; l++) { 
409                                         int bi=(2*i)+this->dfVecX[l], bj=(2*j)+this->dfVecY[l], bk=(2*k)+this->dfVecZ[l];
410                                         if(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&(CFGrFromCoarse)) {
411                                                 //errMsg("performRefinement","Removing CFGrFromCoarse on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) );
412                                                 RFLAG(lev+1,  bi, bj, bk, srcFineSet) = setf;
413                                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 
414                                         }
415                                         else if(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&(CFUnused      )) { 
416                                                 //errMsg("performRefinement","Removing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(bi,bj,bk) );
417                                                 interpolateCellFromCoarse(lev+1,  bi, bj, bk, srcFineSet, interTime, setf, false);
418                                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,bi,bj,bk); 
419                                                 mNumFsgrChanges++;
420                                         }
421                                 }
422                                 for(int l=1; l<this->cDirNum; l++) { 
423                                         int bi=(2*i)+this->dfVecX[l], bj=(2*j)+this->dfVecY[l], bk=(2*k)+this->dfVecZ[l];
424                                         if(   (RFLAG(lev+1,  bi, bj, bk, srcFineSet)&CFFluid       ) &&
425                                                         (!(RFLAG(lev+1,  bi, bj, bk, srcFineSet)&CFGrFromCoarse)) ) {
426                                                 // all unused nbs now of coarse have to be from coarse
427                                                 for(int m=1; m<this->cDirNum; m++) { 
428                                                         int mi=  bi +this->dfVecX[m], mj=  bj +this->dfVecY[m], mk=  bk +this->dfVecZ[m];
429                                                         if(RFLAG(lev+1,  mi, mj, mk, srcFineSet)&CFUnused) {
430                                                                 //errMsg("performRefinement","Changing CFUnused on lev"<<(lev+1)<<" "<<PRINT_VEC(mi,mj,mk) );
431                                                                 interpolateCellFromCoarse(lev+1,  mi, mj, mk, srcFineSet, interTime, CFFluid|CFGrFromCoarse, false);
432                                                                 if((LBMDIM==2)&&(debugRefinement)) debugMarkCell(lev+1,mi,mj,mk); 
433                                                                 mNumFsgrChanges++;
434                                                         }
435                                                 }
436                                                 // nbs prepared...
437                                         }
438                                 }
439                         }
440                         
441                 } // convert regions of from fine
442         }}} // TEST
443         // PASS 3 */
444
445         if(!this->mSilent){ errMsg("performRefinement"," for l"<<lev<<" done ("<<change<<") " ); }
446         } // PASS 1-3
447         // refinement done
448
449         //LbmFsgrSolver::performCoarsening(int lev) {
450         { // PASS 4,5
451         bool nbsok;
452         // WARNING
453         // now work on modified curr set
454         const int srcSet = mLevel[lev].setCurr;
455         const int dstlev = lev+1;
456         const int dstFineSet = mLevel[dstlev].setCurr;
457         const bool debugCoarsening = false;
458
459         // PASS 5 test DEBUG
460         /*if(this->mInitDone) {
461         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
462   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
463   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
464                         if(RFLAG(lev, i,j,k, srcSet) & CFEmpty) {
465                                 // check empty -> from fine conversion
466                                 bool changeToFromFine = false;
467                                 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
468                                 CellFlagType reqType = CFGrNorm;
469                                 if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
470                                 if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) &&
471                                     (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
472                                         changeToFromFine=true; }
473                                 if(changeToFromFine) {
474                                         change = true;
475                                         mNumFsgrChanges++;
476                                         RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrFromFine;
477                                         if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
478                                         // same as restr from fine func! not necessary ?!
479                                         // coarseRestrictFromFine part 
480                                         coarseRestrictCell(lev, i,j,k,srcSet, dstFineSet);
481                                 }
482                         } // only check empty cells
483         }}} // TEST!
484         } // PASS 5 */
485
486         // use //template functions for 2D/3D
487                                         /*if(strstr(this->getName().c_str(),"Debug"))
488                                         if((nbsok)&&(lev+1==mMaxRefine)) { // mixborder
489                                                 for(int l=0;((l<this->cDirNum) && (nbsok)); l++) {  // FARBORD
490                                                         int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l];
491                                                         if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT
492                                                                 nbsok=false;
493                                                         }
494                                                 }
495                                         } // FARBORD */
496         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
497   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
498   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
499
500                         // from coarse cells without unused nbs are not necessary...! -> remove
501                         // perform check from coarseAdvance here?
502                         if(RFLAG(lev, i,j,k, srcSet) & CFGrFromFine) {
503                                 // remove from fine cells now that are completely in fluid
504                                 // FIXME? check that new from fine in performRefinement never get deleted here afterwards?
505                                 // or more general, one cell never changed more than once?
506                                 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
507                                 //const CellFlagType notNbAllowed = (CFInter|CFBnd|CFGrFromFine); unused
508                                 CellFlagType reqType = CFGrNorm;
509                                 if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
510
511                                 nbsok = true;
512                                 for(int l=0; l<this->cDirNum && nbsok; l++) { 
513                                         int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l];
514                                         if(   (RFLAG(lev+1, ni,nj,nk, dstFineSet) & reqType) &&
515                                                         (!(RFLAG(lev+1, ni,nj,nk, dstFineSet) & (notAllowed)) )  ){
516                                                 // ok
517                                         } else {
518                                                 nbsok=false;
519                                         }
520                                         // FARBORD
521                                 }
522                                 // dont turn CFGrFromFine above interface cells into CFGrNorm
523                                 // now check nbs on same level
524                                 for(int l=1; l<this->cDirNum && nbsok; l++) { 
525                                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
526                                         if(RFLAG(lev, ni,nj,nk, srcSet)&(CFFluid)) { //ok
527                                         } else {
528                                                 nbsok = false;
529                                         }
530                                 } // l
531
532                                 if(nbsok) {
533                                         // conversion to coarse fluid cell
534                                         change = true;
535                                         mNumFsgrChanges++;
536                                         RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrNorm;
537                                         // dfs are already ok...
538                                         //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine changed to CFGrNorm at lev"<<lev<<" " <<PRINT_IJK );
539                                         if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
540
541                                         // only check complete cubes
542                                         for(int dx=-1;dx<=1;dx+=2) {
543                                         for(int dy=-1;dy<=1;dy+=2) {
544                                         for(int dz=-1*(LBMDIM&1);dz<=1*(LBMDIM&1);dz+=2) { // 2d/3d
545                                                 // check for norm and from coarse, as the coarse level might just have been refined...
546                                                 if( 
547                                                                 // we now the flag of the current cell! ( RFLAG(lev, i   , j   , k   ,  srcSet)&(CFGrNorm)) &&
548                                                                 ( RFLAG(lev, i+dx, j   , k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
549                                                                 ( RFLAG(lev, i   , j+dy, k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
550                                                                 ( RFLAG(lev, i   , j   , k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
551
552                                                                 ( RFLAG(lev, i+dx, j+dy, k   ,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
553                                                                 ( RFLAG(lev, i+dx, j   , k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
554                                                                 ( RFLAG(lev, i   , j+dy, k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) &&
555                                                                 ( RFLAG(lev, i+dx, j+dy, k+dz,  srcSet)&(CFGrNorm|CFGrFromCoarse)) 
556                                                         ) {
557                                                         // middle source node on higher level
558                                                         int dstx = (2*i)+dx;
559                                                         int dsty = (2*j)+dy;
560                                                         int dstz = (2*k)+dz;
561
562                                                         mNumFsgrChanges++;
563                                                         RFLAG(dstlev, dstx,dsty,dstz, dstFineSet) = CFUnused;
564                                                         RFLAG(dstlev, dstx,dsty,dstz, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST
565                                                         //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init center unused set l"<<dstlev<<" at "<<PRINT_VEC(dstx,dsty,dstz) );
566
567                                                         for(int l=1; l<this->cDirNum; l++) { 
568                                                                 int dstni=dstx+this->dfVecX[l], dstnj=dsty+this->dfVecY[l], dstnk=dstz+this->dfVecZ[l];
569                                                                 if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFFluid)) { 
570                                                                         RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse;
571                                                                 }
572                                                                 if(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFInter)) { 
573                                                                         //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init CHECK Warning - deleting interface cell...");
574                                                                         this->mFixMass += QCELL( dstlev, dstni,dstnj,dstnk, dstFineSet, dMass);
575                                                                         RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFFluid|CFGrFromCoarse;
576                                                                 }
577                                                         } // l
578
579                                                         // again check nb flags of all surrounding cells to see if any from coarse
580                                                         // can be convted to unused
581                                                         for(int l=1; l<this->cDirNum; l++) { 
582                                                                 int dstni=dstx+this->dfVecX[l], dstnj=dsty+this->dfVecY[l], dstnk=dstz+this->dfVecZ[l];
583                                                                 // have to be at least from coarse here...
584                                                                 //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" "<< convertCellFlagType2String(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)) );
585                                                                 if(!(RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet)&(CFUnused) )) { 
586                                                                         bool delok = true;
587                                                                         // careful long range here... check domain bounds?
588                                                                         for(int m=1; m<this->cDirNum; m++) {                                                                            
589                                                                                 int chkni=dstni+this->dfVecX[m], chknj=dstnj+this->dfVecY[m], chknk=dstnk+this->dfVecZ[m];
590                                                                                 if(RFLAG(dstlev, chkni,chknj,chknk, dstFineSet)&(CFUnused|CFGrFromCoarse)) { 
591                                                                                         // this nb cell is ok for deletion
592                                                                                 } else { 
593                                                                                         delok=false; // keep it!
594                                                                                 }
595                                                                                 //errMsg("performCoarsening"," CHECK "<<PRINT_VEC(dstni,dstnj,dstnk)<<" to "<<PRINT_VEC( chkni,chknj,chknk )<<" f:"<< convertCellFlagType2String( RFLAG(dstlev, chkni,chknj,chknk, dstFineSet))<<" nbsok"<<delok );
596                                                                         }
597                                                                         //errMsg("performCoarsening","CFGrFromFine subcube init unused check l"<<dstlev<<" at "<<PRINT_VEC(dstni,dstnj,dstnk)<<" ok"<<delok );
598                                                                         if(delok) {
599                                                                                 mNumFsgrChanges++;
600                                                                                 RFLAG(dstlev, dstni,dstnj,dstnk, dstFineSet) = CFUnused;
601                                                                                 RFLAG(dstlev, dstni,dstnj,dstnk, mLevel[dstlev].setOther) = CFUnused; // FLAGTEST
602                                                                                 if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(dstlev,dstni,dstnj,dstnk); 
603                                                                         }
604                                                                 }
605                                                         } // l
606                                                         // treat subcube
607                                                         //ebugMarkCell(lev,i+dx,j+dy,k+dz); 
608                                                         //if(this->mInitDone) errMsg("performCoarsening","CFGrFromFine subcube init, dir:"<<PRINT_VEC(dx,dy,dz) );
609                                                 }
610                                         } } }
611
612                                 }   // ?
613                         } // convert regions of from fine
614         }}} // TEST!
615         // PASS 4 */
616
617                 // reinit cell area value
618                 /*if( RFLAG(lev, i,j,k,srcSet) & CFFluid) {
619                         if( RFLAG(lev+1, i*2,j*2,k*2,dstFineSet) & CFGrFromCoarse) {
620                                 LbmFloat totArea = mFsgrCellArea[0]; // for l=0
621                                 for(int l=1; l<this->cDirNum; l++) { 
622                                         int ni=(2*i)+this->dfVecX[l], nj=(2*j)+this->dfVecY[l], nk=(2*k)+this->dfVecZ[l];
623                                         if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&
624                                                         (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
625                                                         //(CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
626                                                         ) { 
627                                                 //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;
628                                                 totArea += mFsgrCellArea[l];
629                                         }
630                                 } // l
631                                 QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 
632                                 QCELL(lev, i,j,k,srcSet, dFlux) = totArea;
633                         } else {
634                                 QCELL(lev, i,j,k,mLevel[lev].setOther, dFlux) = 
635                                 QCELL(lev, i,j,k,srcSet, dFlux) = 1.0;
636                         }
637                         //errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,srcSet, dFlux) );
638                 }
639         // */
640
641
642         // PASS 5 org
643         /*if(strstr(this->getName().c_str(),"Debug"))
644         if((changeToFromFine)&&(lev+1==mMaxRefine)) { // mixborder
645                 for(int l=0;((l<this->cDirNum) && (changeToFromFine)); l++) {  // FARBORD
646                         int ni=2*i+2*this->dfVecX[l], nj=2*j+2*this->dfVecY[l], nk=2*k+2*this->dfVecZ[l];
647                         if(RFLAG(lev+1, ni,nj,nk, dstFineSet)&CFBnd) { // NEWREFT
648                                 changeToFromFine=false; }
649                 } 
650         }// FARBORD */
651         //if(!this->mInitDone) {
652         for(int k= getForZMin1(); k< getForZMax1(lev); ++k) {
653   for(int j=1;j<mLevel[lev].lSizey-1;++j) {
654   for(int i=1;i<mLevel[lev].lSizex-1;++i) {
655
656
657                         if(RFLAG(lev, i,j,k, srcSet) & CFEmpty) {
658                                 // check empty -> from fine conversion
659                                 bool changeToFromFine = false;
660                                 const CellFlagType notAllowed = (CFInter|CFGrFromFine|CFGrToFine);
661                                 CellFlagType reqType = CFGrNorm;
662                                 if(lev+1==mMaxRefine) reqType = CFNoBndFluid;
663
664                                 if(   (RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & reqType) &&
665                                     (!(RFLAG(lev+1, (2*i),(2*j),(2*k), dstFineSet) & (notAllowed)) )  ){
666                                         // DEBUG 
667                                         changeToFromFine=true;
668                                 }
669
670                                 // FARBORD
671
672                                 if(changeToFromFine) {
673                                         change = true;
674                                         mNumFsgrChanges++;
675                                         RFLAG(lev, i,j,k, srcSet) = CFFluid|CFGrFromFine;
676                                         if((LBMDIM==2)&&(debugCoarsening)) debugMarkCell(lev,i,j,k); 
677                                         // same as restr from fine func! not necessary ?!
678                                         // coarseRestrictFromFine part 
679                                 }
680                         } // only check empty cells
681
682         }}} // TEST!
683         //} // init done
684         // PASS 5 */
685         } // coarsening, PASS 4,5
686
687         if(!this->mSilent){ errMsg("adaptGrid"," for l"<<lev<<" done " ); }
688         return change;
689 }
690
691 /*****************************************************************************/
692 //! cell restriction and prolongation
693 /*****************************************************************************/
694
695 void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet)
696 {
697         LbmFloat *ccel = RACPNT(lev+1, 2*i,2*j,2*k,srcSet);
698         LbmFloat *tcel = RACPNT(lev  , i,j,k      ,dstSet);
699
700         LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;                       
701         //LbmFloat *ccel = NULL;
702         //LbmFloat *tcel = NULL;
703 #if OPT3D==1 
704         LbmFloat m[LBM_DFNUM];
705         // for macro add
706         LbmFloat usqr;
707         //LbmFloat *addfcel, *dstcell;
708         LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM];
709         LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale;
710 #else // OPT3D==true 
711         LbmFloat df[LBM_DFNUM];
712         LbmFloat omegaDst, omegaSrc;
713         LbmFloat feq[LBM_DFNUM];
714         LbmFloat dfScale = mDfScaleUp;
715 #endif // OPT3D==true 
716
717 #                               if OPT3D==0
718         // add up weighted dfs
719         FORDF0{ df[l] = 0.0;}
720         for(int n=0;(n<this->cDirNum); n++) { 
721                 int ni=2*i+1*this->dfVecX[n], nj=2*j+1*this->dfVecY[n], nk=2*k+1*this->dfVecZ[n];
722                 ccel = RACPNT(lev+1, ni,nj,nk,srcSet);// CFINTTEST
723                 const LbmFloat weight = mGaussw[n];
724                 FORDF0{
725                         LbmFloat cdf = weight * RAC(ccel,l);
726 #                                               if FSGR_STRICT_DEBUG==1
727                         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]); }
728 #                                               endif
729                         //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); 
730                         df[l] += cdf;
731                 }
732         }
733
734         // calc rho etc. from weighted dfs
735         rho = ux  = uy  = uz  = 0.0;
736         FORDF0{
737                 LbmFloat cdf = df[l];
738                 rho += cdf; 
739                 ux  += (this->dfDvecX[l]*cdf); 
740                 uy  += (this->dfDvecY[l]*cdf);  
741                 uz  += (this->dfDvecZ[l]*cdf);  
742         }
743
744         FORDF0{ feq[l] = this->getCollideEq(l, rho,ux,uy,uz); }
745         if(mLevel[lev  ].lcsmago>0.0) {
746                 const LbmFloat Qo = this->getLesNoneqTensorCoeff(df,feq);
747                 omegaDst  = this->getLesOmega(mLevel[lev  ].omega,mLevel[lev  ].lcsmago,Qo);
748                 omegaSrc = this->getLesOmega(mLevel[lev+1].omega,mLevel[lev+1].lcsmago,Qo);
749         } else {
750                 omegaDst = mLevel[lev+0].omega; /* NEWSMAGOT*/ 
751                 omegaSrc = mLevel[lev+1].omega;
752         }
753         dfScale   = (mLevel[lev  ].timestep/mLevel[lev+1].timestep)* (1.0/omegaDst-1.0)/ (1.0/omegaSrc-1.0); // yu
754         FORDF0{
755                 RAC(tcel, l) = feq[l]+ (df[l]-feq[l])*dfScale;
756         } 
757 #                               else // OPT3D
758         // similar to OPTIMIZED_STREAMCOLLIDE_UNUSED
759                                                                 
760         //rho = ux = uy = uz = 0.0;
761         MSRC_C  = CCELG_C(0) ;
762         MSRC_N  = CCELG_N(0) ;
763         MSRC_S  = CCELG_S(0) ;
764         MSRC_E  = CCELG_E(0) ;
765         MSRC_W  = CCELG_W(0) ;
766         MSRC_T  = CCELG_T(0) ;
767         MSRC_B  = CCELG_B(0) ;
768         MSRC_NE = CCELG_NE(0);
769         MSRC_NW = CCELG_NW(0);
770         MSRC_SE = CCELG_SE(0);
771         MSRC_SW = CCELG_SW(0);
772         MSRC_NT = CCELG_NT(0);
773         MSRC_NB = CCELG_NB(0);
774         MSRC_ST = CCELG_ST(0);
775         MSRC_SB = CCELG_SB(0);
776         MSRC_ET = CCELG_ET(0);
777         MSRC_EB = CCELG_EB(0);
778         MSRC_WT = CCELG_WT(0);
779         MSRC_WB = CCELG_WB(0);
780         for(int n=1;(n<this->cDirNum); n++) { 
781                 ccel = RACPNT(lev+1,  2*i+1*this->dfVecX[n], 2*j+1*this->dfVecY[n], 2*k+1*this->dfVecZ[n]  ,srcSet);
782                 MSRC_C  += CCELG_C(n) ;
783                 MSRC_N  += CCELG_N(n) ;
784                 MSRC_S  += CCELG_S(n) ;
785                 MSRC_E  += CCELG_E(n) ;
786                 MSRC_W  += CCELG_W(n) ;
787                 MSRC_T  += CCELG_T(n) ;
788                 MSRC_B  += CCELG_B(n) ;
789                 MSRC_NE += CCELG_NE(n);
790                 MSRC_NW += CCELG_NW(n);
791                 MSRC_SE += CCELG_SE(n);
792                 MSRC_SW += CCELG_SW(n);
793                 MSRC_NT += CCELG_NT(n);
794                 MSRC_NB += CCELG_NB(n);
795                 MSRC_ST += CCELG_ST(n);
796                 MSRC_SB += CCELG_SB(n);
797                 MSRC_ET += CCELG_ET(n);
798                 MSRC_EB += CCELG_EB(n);
799                 MSRC_WT += CCELG_WT(n);
800                 MSRC_WB += CCELG_WB(n);
801         }
802         rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T  
803                 + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT 
804                 + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; 
805         ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW 
806                 + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB;  
807         uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW 
808                 + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB;  
809         uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB 
810                 + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB;  
811         usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
812         \
813         lcsmeq[dC] = EQC ; \
814         COLL_CALCULATE_DFEQ(lcsmeq); \
815         COLL_CALCULATE_NONEQTENSOR(lev+0, MSRC_ )\
816         COLL_CALCULATE_CSMOMEGAVAL(lev+0, lcsmDstOmega); \
817         COLL_CALCULATE_CSMOMEGAVAL(lev+1, lcsmSrcOmega); \
818         \
819         lcsmdfscale   = (mLevel[lev+0].timestep/mLevel[lev+1].timestep)* (1.0/lcsmDstOmega-1.0)/ (1.0/lcsmSrcOmega-1.0);  \
820         RAC(tcel, dC ) = (lcsmeq[dC ] + (MSRC_C -lcsmeq[dC ] )*lcsmdfscale);
821         RAC(tcel, dN ) = (lcsmeq[dN ] + (MSRC_N -lcsmeq[dN ] )*lcsmdfscale);
822         RAC(tcel, dS ) = (lcsmeq[dS ] + (MSRC_S -lcsmeq[dS ] )*lcsmdfscale);
823         RAC(tcel, dE ) = (lcsmeq[dE ] + (MSRC_E -lcsmeq[dE ] )*lcsmdfscale);
824         RAC(tcel, dW ) = (lcsmeq[dW ] + (MSRC_W -lcsmeq[dW ] )*lcsmdfscale);
825         RAC(tcel, dT ) = (lcsmeq[dT ] + (MSRC_T -lcsmeq[dT ] )*lcsmdfscale);
826         RAC(tcel, dB ) = (lcsmeq[dB ] + (MSRC_B -lcsmeq[dB ] )*lcsmdfscale);
827         RAC(tcel, dNE) = (lcsmeq[dNE] + (MSRC_NE-lcsmeq[dNE] )*lcsmdfscale);
828         RAC(tcel, dNW) = (lcsmeq[dNW] + (MSRC_NW-lcsmeq[dNW] )*lcsmdfscale);
829         RAC(tcel, dSE) = (lcsmeq[dSE] + (MSRC_SE-lcsmeq[dSE] )*lcsmdfscale);
830         RAC(tcel, dSW) = (lcsmeq[dSW] + (MSRC_SW-lcsmeq[dSW] )*lcsmdfscale);
831         RAC(tcel, dNT) = (lcsmeq[dNT] + (MSRC_NT-lcsmeq[dNT] )*lcsmdfscale);
832         RAC(tcel, dNB) = (lcsmeq[dNB] + (MSRC_NB-lcsmeq[dNB] )*lcsmdfscale);
833         RAC(tcel, dST) = (lcsmeq[dST] + (MSRC_ST-lcsmeq[dST] )*lcsmdfscale);
834         RAC(tcel, dSB) = (lcsmeq[dSB] + (MSRC_SB-lcsmeq[dSB] )*lcsmdfscale);
835         RAC(tcel, dET) = (lcsmeq[dET] + (MSRC_ET-lcsmeq[dET] )*lcsmdfscale);
836         RAC(tcel, dEB) = (lcsmeq[dEB] + (MSRC_EB-lcsmeq[dEB] )*lcsmdfscale);
837         RAC(tcel, dWT) = (lcsmeq[dWT] + (MSRC_WT-lcsmeq[dWT] )*lcsmdfscale);
838         RAC(tcel, dWB) = (lcsmeq[dWB] + (MSRC_WB-lcsmeq[dWB] )*lcsmdfscale);
839 #                               endif // OPT3D==0
840 }
841
842 void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet, bool markNbs) {
843         LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;
844         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 };
845
846 #if OPT3D==1 
847         // for macro add
848         LbmFloat addDfFacT, addVal, usqr;
849         LbmFloat *addfcel, *dstcell;
850         LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM];
851         LbmFloat lcsmDstOmega, lcsmSrcOmega, lcsmdfscale;
852 #endif // OPT3D==true 
853
854         // SET required nbs to from coarse (this might overwrite flag several times)
855         // this is not necessary for interpolateFineFromCoarse
856         if(markNbs) {
857         FORDF1{ 
858                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
859                 if(RFLAG(lev,ni,nj,nk,dstSet)&CFUnused) {
860                         // parents have to be inited!
861                         interpolateCellFromCoarse(lev, ni, nj, nk, dstSet, t, CFFluid|CFGrFromCoarse, false);
862                 }
863         } }
864
865         // change flag of cell to be interpolated
866         RFLAG(lev,i,j,k, dstSet) = flagSet;
867         mNumInterdCells++;
868
869         // interpolation lines...
870         int betx = i&1;
871         int bety = j&1;
872         int betz = k&1;
873         
874         if((!betx) && (!bety) && (!betz)) {
875                 ADD_INT_DFS(lev-1, i/2  ,j/2  ,k/2  , 0.0, 1.0);
876         }
877         else if(( betx) && (!bety) && (!betz)) {
878                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
879                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D1);
880         }
881         else if((!betx) && ( bety) && (!betz)) {
882                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
883                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D1);
884         }
885         else if((!betx) && (!bety) && ( betz)) {
886                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D1);
887                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D1);
888         }
889         else if(( betx) && ( bety) && (!betz)) {
890                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
891                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D2);
892                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D2);
893                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)  , t, WO1D2);
894         }
895         else if((!betx) && ( bety) && ( betz)) {
896                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
897                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D2);
898                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D2);
899                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)+1, t, WO1D2);
900         }
901         else if(( betx) && (!bety) && ( betz)) {
902                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D2);
903                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D2);
904                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D2);
905                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)+1, t, WO1D2);
906         }
907         else if(( betx) && ( bety) && ( betz)) {
908                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)  , t, WO1D3);
909                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)  , t, WO1D3);
910                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)  ,(k/2)+1, t, WO1D3);
911                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)  ,(k/2)+1, t, WO1D3);
912                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)  , t, WO1D3);
913                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)  , t, WO1D3);
914                 ADD_INT_DFS(lev-1, (i/2)  ,(j/2)+1,(k/2)+1, t, WO1D3);
915                 ADD_INT_DFS(lev-1, (i/2)+1,(j/2)+1,(k/2)+1, t, WO1D3);
916         }
917         else {
918                 CAUSE_PANIC;
919                 errFatal("interpolateCellFromCoarse","Invalid!?", SIMWORLD_GENERICERROR);       
920         }
921
922         IDF_WRITEBACK;
923         return;
924 }
925
926
927
928 // required globals
929 extern bool glob_mpactive;
930 extern int glob_mpnum, glob_mpindex;
931 #define MPTADAP_INTERV 4
932
933 /*****************************************************************************/
934 /*! change the  size of the LBM time step */
935 /*****************************************************************************/
936 void LbmFsgrSolver::adaptTimestep() {
937         LbmFloat massTOld=0.0, massTNew=0.0;
938         LbmFloat volTOld=0.0, volTNew=0.0;
939
940         bool rescale = false;  // do any rescale at all?
941         LbmFloat scaleFac = -1.0; // timestep scaling
942         if(mPanic) return;
943
944         LbmFloat levOldOmega[FSGR_MAXNOOFLEVELS];
945         LbmFloat levOldStepsize[FSGR_MAXNOOFLEVELS];
946         for(int lev=mMaxRefine; lev>=0 ; lev--) {
947                 levOldOmega[lev] = mLevel[lev].omega;
948                 levOldStepsize[lev] = mLevel[lev].timestep;
949         }
950
951         const LbmFloat reduceFac = 0.8;          // modify time step by 20%, TODO? do multiple times for large changes?
952         LbmFloat diffPercent = 0.05; // dont scale if less than 5%
953         LbmFloat allowMax = mpParam->getTadapMaxSpeed();  // maximum allowed velocity
954         LbmFloat nextmax = mpParam->getSimulationMaxSpeed() + norm(mLevel[mMaxRefine].gravity);
955
956         // sync nextmax
957 #if LBM_INCLUDE_TESTSOLVERS==1
958         if(glob_mpactive) {
959                 if(mLevel[mMaxRefine].lsteps % MPTADAP_INTERV != MPTADAP_INTERV-1) {
960                         debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "mpact:"<<glob_mpactive<<","<<glob_mpindex<<"/"<<glob_mpnum<<" step:"<<mLevel[mMaxRefine].lsteps<<" skipping tadap...",8);
961                         return;
962                 }
963                 nextmax = mrInitTadap(nextmax);
964         }
965 #endif // LBM_INCLUDE_TESTSOLVERS==1
966
967         LbmFloat newdt = mpParam->getTimestep(); // newtr
968         if(nextmax > allowMax/reduceFac) {
969                 mTimeMaxvelStepCnt++; }
970         else { mTimeMaxvelStepCnt=0; }
971         
972         // emergency, or 10 steps with high vel
973         if((mTimeMaxvelStepCnt>5) || (nextmax> (1.0/3.0)) || (mForceTimeStepReduce) ) {
974                 newdt = mpParam->getTimestep() * reduceFac;
975         } else {
976                 if(nextmax<allowMax*reduceFac) {
977                         newdt = mpParam->getTimestep() / reduceFac;
978                 }
979         } // newtr
980         //errMsg("LbmFsgrSolver::adaptTimestep","nextmax="<<nextmax<<" allowMax="<<allowMax<<" fac="<<reduceFac<<" simmaxv="<< mpParam->getSimulationMaxSpeed() );
981
982         bool minCutoff = false;
983         LbmFloat desireddt = newdt;
984         if(newdt>mpParam->getMaxTimestep()){ newdt = mpParam->getMaxTimestep(); }
985         if(newdt<mpParam->getMinTimestep()){ 
986                 newdt = mpParam->getMinTimestep(); 
987                 if(nextmax>allowMax/reduceFac){ minCutoff=true; } // only if really large vels...
988         }
989
990         LbmFloat dtdiff = fabs(newdt - mpParam->getTimestep());
991         if(!mSilent) {
992                 debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt
993                         <<" max"<<mpParam->getMaxTimestep()<<" min"<<mpParam->getMinTimestep()<<" diff"<<dtdiff
994                         <<" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep)<<
995                         " olddt"<<levOldStepsize[mMaxRefine]<<" redlock"<<mTimestepReduceLock 
996                         , 10); }
997
998         // in range, and more than X% change?
999         //if( newdt <  mpParam->getTimestep() ) // DEBUG
1000         LbmFloat rhoAvg = mCurrentMass/mCurrentVolume;
1001         if( (newdt<=mpParam->getMaxTimestep()) && (newdt>=mpParam->getMinTimestep()) 
1002                         && (dtdiff>(mpParam->getTimestep()*diffPercent)) ) {
1003                 if((newdt>levOldStepsize[mMaxRefine])&&(mTimestepReduceLock)) {
1004                         // wait some more...
1005                         //debMsgNnl("LbmFsgrSolver::TAdp",DM_NOTIFY," Delayed... "<<mTimestepReduceLock<<" ",10);
1006                         debMsgDirect("D");
1007                 } else {
1008                         mpParam->setDesiredTimestep( newdt );
1009                         rescale = true;
1010                         if(!mSilent) {
1011                                 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"\n\n\n\n",10);
1012                                 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: new="<<newdt<<" old="<<mpParam->getTimestep()
1013                                                 <<" maxSpeed:"<<mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<mStepCnt, 10 );
1014                                 //debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep changing: "<< "rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
1015                         }
1016                 } // really change dt
1017         }
1018
1019         if(mTimestepReduceLock>0) mTimestepReduceLock--;
1020
1021         
1022         // forced back and forth switchting (for testing)
1023         /*const int tadtogInter = 100;
1024         const double tadtogSwitch = 0.66;
1025         errMsg("TIMESWITCHTOGGLETEST","warning enabled "<< tadtogSwitch<<","<<tadtogSwitch<<" !!!!!!!!!!!!!!!!!!!");
1026         if( ((mStepCnt% tadtogInter)== (tadtogInter/4*1)-1) ||
1027             ((mStepCnt% tadtogInter)== (tadtogInter/4*2)-1) ){
1028                 rescale = true; minCutoff = false;
1029                 newdt = tadtogSwitch * mpParam->getTimestep();
1030                 mpParam->setDesiredTimestep( newdt );
1031         } else 
1032         if( ((mStepCnt% tadtogInter)== (tadtogInter/4*3)-1) ||
1033             ((mStepCnt% tadtogInter)== (tadtogInter/4*4)-1) ){
1034                 rescale = true; minCutoff = false;
1035                 newdt = mpParam->getTimestep()/tadtogSwitch ;
1036                 mpParam->setDesiredTimestep( newdt );
1037         } else {
1038                 rescale = false; minCutoff = false;
1039         }
1040         // */
1041
1042         // test mass rescale
1043
1044         scaleFac = newdt/mpParam->getTimestep();
1045         if(rescale) {
1046                 // perform acutal rescaling...
1047                 mTimeMaxvelStepCnt=0; 
1048                 mForceTimeStepReduce = false;
1049
1050                 // FIXME - approximate by averaging, take gravity direction here?
1051                 //mTimestepReduceLock = 4*(mLevel[mMaxRefine].lSizey+mLevel[mMaxRefine].lSizez+mLevel[mMaxRefine].lSizex)/3;
1052                 // use z as gravity direction
1053                 mTimestepReduceLock = 4*mLevel[mMaxRefine].lSizez;
1054
1055                 mTimeSwitchCounts++;
1056                 mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
1057                 recalculateObjectSpeeds();
1058                 // calc omega, force for all levels
1059                 mLastOmega=1e10; mLastGravity=1e10;
1060                 initLevelOmegas();
1061                 if(mpParam->getTimestep()<mMinTimestep) mMinTimestep = mpParam->getTimestep();
1062                 if(mpParam->getTimestep()>mMaxTimestep) mMaxTimestep = mpParam->getTimestep();
1063
1064                 // this might be called during init, before we have any particles
1065                 if(mpParticles) { mpParticles->adaptPartTimestep(scaleFac); }
1066 #if LBM_INCLUDE_TESTSOLVERS==1
1067                 if((mUseTestdata)&&(mpTest)) { 
1068                         mpTest->adaptTimestep(scaleFac, mLevel[mMaxRefine].omega, mLevel[mMaxRefine].timestep, vec2L( mpParam->calculateGravity(mSimulationTime)) ); 
1069                         mpTest->mGrav3d = mLevel[mMaxRefine].gravity;
1070                 }
1071 #endif // LBM_INCLUDE_TESTSOLVERS!=1
1072         
1073                 for(int lev=mMaxRefine; lev>=0 ; lev--) {
1074                         LbmFloat newSteptime = mLevel[lev].timestep;
1075                         LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]);
1076
1077                         if(!mSilent) {
1078                                 debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Level: "<<lev<<" Timestep chlevel: "<<
1079                                                 " scaleFac="<<dfScaleFac<<" newDt="<<newSteptime<<" newOmega="<<mLevel[lev].omega,10);
1080                         }
1081                         if(lev!=mMaxRefine) coarseCalculateFluxareas(lev);
1082
1083                         int wss = 0, wse = 1;
1084                         // only change currset (necessary for compressed grids, better for non-compr.gr.)
1085                         wss = wse = mLevel[lev].setCurr;
1086                         for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
1087                                         // warning - check sets for higher levels...?
1088                                 FSGR_FORIJK1(lev) {
1089                                         if( (RFLAG(lev,i,j,k, workSet) & CFBndMoving) ) {
1090                                                 /*
1091                                                 // paranoid check - shouldnt be necessary!
1092                                                 if(QCELL(lev, i, j, k, workSet, dFlux)!=mSimulationTime) {
1093                                                         errMsg("TTTT","found invalid bnd cell.... removing at "<<PRINT_IJK);
1094                                                         RFLAG(lev,i,j,k, workSet) = CFInter;
1095                                                         // init empty zero vel  interface cell...
1096                                                         initVelocityCell(lev, i,j,k, CFInter, 1.0, 0.01, LbmVec(0.) );
1097                                                 } else {// */
1098                                                         for(int l=0; l<cDfNum; l++) {
1099                                                                 QCELL(lev, i, j, k, workSet, l) = QCELL(lev, i, j, k, workSet, l)* scaleFac; 
1100                                                         }
1101                                                 //} //  ok
1102                                                 continue;
1103                                         }
1104                                         if( 
1105                                                         (RFLAG(lev,i,j,k, workSet) & CFFluid) || 
1106                                                         (RFLAG(lev,i,j,k, workSet) & CFInter) ||
1107                                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 
1108                                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 
1109                                                         (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 
1110                                                         ) {
1111                                                 // these cells have to be scaled...
1112                                         } else {
1113                                                 continue;
1114                                         }
1115
1116                                         // collide on current set
1117                                         LbmFloat rhoOld;
1118                                         LbmVec velOld;
1119                                         LbmFloat rho, ux,uy,uz;
1120                                         rho=0.0; ux =  uy = uz = 0.0;
1121                                         for(int l=0; l<cDfNum; l++) {
1122                                                 LbmFloat m = QCELL(lev, i, j, k, workSet, l); 
1123                                                 rho += m;
1124                                                 ux  += (dfDvecX[l]*m);
1125                                                 uy  += (dfDvecY[l]*m); 
1126                                                 uz  += (dfDvecZ[l]*m); 
1127                                         } 
1128                                         rhoOld = rho;
1129                                         velOld = LbmVec(ux,uy,uz);
1130
1131                                         LbmFloat rhoNew = (rhoOld-rhoAvg)*scaleFac +rhoAvg;
1132                                         LbmVec velNew = velOld * scaleFac;
1133
1134                                         LbmFloat df[LBM_DFNUM];
1135                                         LbmFloat feqOld[LBM_DFNUM];
1136                                         LbmFloat feqNew[LBM_DFNUM];
1137                                         for(int l=0; l<cDfNum; l++) {
1138                                                 feqOld[l] = getCollideEq(l,rhoOld, velOld[0],velOld[1],velOld[2] );
1139                                                 feqNew[l] = getCollideEq(l,rhoNew, velNew[0],velNew[1],velNew[2] );
1140                                                 df[l] = QCELL(lev, i,j,k,workSet, l);
1141                                         }
1142                                         const LbmFloat Qo = getLesNoneqTensorCoeff(df,feqOld);
1143                                         const LbmFloat oldOmega = getLesOmega(levOldOmega[lev], mLevel[lev].lcsmago,Qo);
1144                                         const LbmFloat newOmega = getLesOmega(mLevel[lev].omega,mLevel[lev].lcsmago,Qo);
1145                                         //newOmega = mLevel[lev].omega; // FIXME debug test
1146
1147                                         //LbmFloat dfScaleFac = (newSteptime/1.0)/(levOldStepsize[lev]/levOldOmega[lev]);
1148                                         const LbmFloat dfScale = (newSteptime/newOmega)/(levOldStepsize[lev]/oldOmega);
1149                                         //dfScale = dfScaleFac/newOmega;
1150                                         
1151                                         for(int l=0; l<cDfNum; l++) {
1152                                                 // org scaling
1153                                                 //df = eqOld + (df-eqOld)*dfScale; df *= (eqNew/eqOld); // non-eq. scaling, important
1154                                                 // new scaling
1155                                                 LbmFloat dfn = feqNew[l] + (df[l]-feqOld[l])*dfScale*feqNew[l]/feqOld[l]; // non-eq. scaling, important
1156                                                 //df = eqNew + (df-eqOld)*dfScale; // modified ig scaling, no real difference?
1157                                                 QCELL(lev, i,j,k,workSet, l) = dfn;
1158                                         }
1159
1160                                         if(RFLAG(lev,i,j,k, workSet) & CFInter) {
1161                                                 //if(workSet==mLevel[lev].setCurr) 
1162                                                 LbmFloat area = 1.0;
1163                                                 if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux);
1164                                                 massTOld += QCELL(lev, i,j,k,workSet, dMass) * area;
1165                                                 volTOld += QCELL(lev, i,j,k,workSet, dFfrac);
1166
1167                                                 // wrong... QCELL(i,j,k,workSet, dMass] = (QCELL(i,j,k,workSet, dFfrac]*rhoNew);
1168                                                 QCELL(lev, i,j,k,workSet, dMass) = (QCELL(lev, i,j,k,workSet, dMass)/rhoOld*rhoNew);
1169                                                 QCELL(lev, i,j,k,workSet, dFfrac) = (QCELL(lev, i,j,k,workSet, dMass)/rhoNew);
1170
1171                                                 //if(workSet==mLevel[lev].setCurr) 
1172                                                 massTNew += QCELL(lev, i,j,k,workSet, dMass);
1173                                                 volTNew += QCELL(lev, i,j,k,workSet, dFfrac);
1174                                         }
1175                                         if(RFLAG(lev,i,j,k, workSet) & CFFluid) { // DEBUG
1176                                                 if(RFLAG(lev,i,j,k, workSet) & (CFGrFromFine|CFGrFromCoarse)) { // DEBUG
1177                                                         // dont include 
1178                                                 } else {
1179                                                         LbmFloat area = 1.0;
1180                                                         if(lev!=mMaxRefine) area = QCELL(lev, i,j,k,workSet, dFlux) * mLevel[lev].lcellfactor;
1181                                                         //if(workSet==mLevel[lev].setCurr) 
1182                                                         massTOld += rhoOld*area;
1183                                                         //if(workSet==mLevel[lev].setCurr) 
1184                                                         massTNew += rhoNew*area;
1185                                                         volTOld += area;
1186                                                         volTNew += area;
1187                                                 }
1188                                         }
1189
1190                                 } // IJK
1191                         } // workSet
1192
1193                 } // lev
1194
1195                 if(!mSilent) {
1196                         debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE "<<mStepCnt<<
1197                                         " no"<<mTimeSwitchCounts<<" maxdt"<<mMaxTimestep<<
1198                                         " mindt"<<mMinTimestep<<" currdt"<<mLevel[mMaxRefine].timestep, 10);
1199                         debMsgStd("LbmFsgrSolver::step",DM_MSG,"REINIT DONE  masst:"<<massTNew<<","<<massTOld<<" org:"<<mCurrentMass<<"; "<<
1200                                         " volt:"<<volTNew<<","<<volTOld<<" org:"<<mCurrentVolume, 10);
1201                 } else {
1202                         debMsgStd("\nLbmOptSolver::step",DM_MSG,"Timestep changed by "<< (newdt/levOldStepsize[mMaxRefine]) <<" newDt:"<<newdt
1203                                         <<", oldDt:"<<levOldStepsize[mMaxRefine]<<" newOmega:"<<mOmega<<" gStar:"<<mpParam->getCurrentGStar()<<", step:"<<mStepCnt , 10);
1204                 }
1205         } // rescale?
1206         //NEWDEBCHECK("tt2");
1207         
1208         //errMsg("adaptTimestep","Warning - brute force rescale off!"); minCutoff = false; // DEBUG
1209         if(minCutoff) {
1210                 errMsg("adaptTimestep","Warning - performing Brute-Force rescale... (sim:"<<mName<<" step:"<<mStepCnt<<" newdt="<<desireddt<<" mindt="<<mpParam->getMinTimestep()<<") " );
1211                 //brute force resacle all the time?
1212
1213                 for(int lev=mMaxRefine; lev>=0 ; lev--) {
1214                 int rescs=0;
1215                 int wss = 0, wse = 1;
1216 #if COMPRESSGRIDS==1
1217                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
1218 #endif // COMPRESSGRIDS==1
1219                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
1220                 //for(int workSet = 0; workSet<=1; workSet++) {
1221                 FSGR_FORIJK1(lev) {
1222
1223                         //if( (RFLAG(lev, i,j,k, workSet) & CFFluid) || (RFLAG(lev, i,j,k, workSet) & CFInter) ) {
1224                         if( 
1225                                         (RFLAG(lev,i,j,k, workSet) & CFFluid) || 
1226                                         (RFLAG(lev,i,j,k, workSet) & CFInter) ||
1227                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 
1228                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 
1229                                         (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 
1230                                         ) {
1231                                 // these cells have to be scaled...
1232                         } else {
1233                                 continue;
1234                         }
1235
1236                         // collide on current set
1237                         LbmFloat rho, ux,uy,uz;
1238                         rho=0.0; ux =  uy = uz = 0.0;
1239                         for(int l=0; l<cDfNum; l++) {
1240                                 LbmFloat m = QCELL(lev, i, j, k, workSet, l); 
1241                                 rho += m;
1242                                 ux  += (dfDvecX[l]*m);
1243                                 uy  += (dfDvecY[l]*m); 
1244                                 uz  += (dfDvecZ[l]*m); 
1245                         } 
1246 #ifndef WIN32
1247                         if (!finite(rho)) {
1248                                 errMsg("adaptTimestep","Brute force non-finite rho at"<<PRINT_IJK);  // DEBUG!
1249                                 rho = 1.0;
1250                                 ux = uy = uz = 0.0;
1251                                 QCELL(lev, i, j, k, workSet, dMass) = 1.0;
1252                                 QCELL(lev, i, j, k, workSet, dFfrac) = 1.0;
1253                         }
1254 #endif // WIN32
1255
1256                         if( (ux*ux+uy*uy+uz*uz)> (allowMax*allowMax) ) {
1257                                 LbmFloat cfac = allowMax/sqrt(ux*ux+uy*uy+uz*uz);
1258                                 ux *= cfac;
1259                                 uy *= cfac;
1260                                 uz *= cfac;
1261                                 for(int l=0; l<cDfNum; l++) {
1262                                         QCELL(lev, i, j, k, workSet, l) = getCollideEq(l, rho, ux,uy,uz); }
1263                                 rescs++;
1264                                 debMsgDirect("B");
1265                         }
1266
1267                 } } 
1268                         //if(rescs>0) { errMsg("adaptTimestep","!!!!! Brute force rescaling was necessary !!!!!!!"); }
1269                         debMsgStd("adaptTimestep",DM_MSG,"Brute force rescale done. level:"<<lev<<" rescs:"<<rescs, 1);
1270                 //TTT mNumProblems += rescs; // add to problem display...
1271                 } // lev,set,ijk
1272
1273         } // try brute force rescale?
1274
1275         // time adap done...
1276 }
1277
1278
1279