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