elbeem update:
[blender.git] / intern / elbeem / intern / solver_main.cpp
1 /******************************************************************************
2  *
3  * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
4  * Copyright 2003,2004,2005 Nils Thuerey
5  *
6  * Standard LBM Factory implementation
7  *
8  *****************************************************************************/
9
10 #include "solver_class.h"
11 #include "solver_relax.h"
12 #include "particletracer.h"
13 #include "loop_tools.h"
14
15 /*****************************************************************************/
16 /*! perform a single LBM step */
17 /*****************************************************************************/
18
19 double globdfcnt;
20 double globdfavg[19];
21 double globdfmax[19];
22 double globdfmin[19];
23
24 void LbmFsgrSolver::step() { 
25         initLevelOmegas();
26         stepMain(); 
27
28         // intlbm test
29         if(0) {
30                 if(this->mStepCnt<5) {
31                         // init
32                         globdfcnt=0.;
33                         FORDF0{ 
34                                 globdfavg[l] = 0.;
35                                 globdfmax[l] = -1000.; //this->dfEquil[l];
36                                 globdfmin[l] = 1000.; // this->dfEquil[l];
37                         }
38                 } else {
39
40                         int workSet = mLevel[mMaxRefine].setCurr;
41                         for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k)  {
42                                 for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++)  {
43                                         for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
44                                                 //float val = 0.;
45                                                 if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) {
46                                                         //val = QCELL(mMaxRefine,i,j,k, workSet,dFfrac);
47                                                         FORDF0{ 
48                                                                 const double df = (double)QCELL(mMaxRefine,i,j,k, workSet,l);
49                                                                 globdfavg[l] += df;
50                                                                 if(df>globdfmax[l]) globdfmax[l] = df;
51                                                                 //if(df<=0.01) { errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l); }
52                                                                 if(df<globdfmin[l]) globdfmin[l] = df;
53                                                                 //errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l<<" "<<df<<" "<<globdfmin[l]); 
54                                                         }
55                                                         globdfcnt+=1.;
56                                                 }
57                                         }
58                                 }
59                         }
60                         if(this->mStepCnt%10==0) {
61                                 FORDF0{ errMsg("GLOBDF","l="<<l<<" avg="<<(globdfavg[l]/globdfcnt)<<"   max="<<globdfmax[l]<<"   min="<<globdfmin[l]<<"   "<<globdfcnt); }
62                         }
63
64                 }
65         }
66         // intlbm test */
67 }
68
69 void LbmFsgrSolver::stepMain()
70 {
71         this->markedClearList(); // DMC clearMarkedCellsList
72
73         // safety check, counter reset
74         this->mNumUsedCells = 0;
75         mNumInterdCells = 0;
76         mNumInvIfCells = 0;
77
78   //debugOutNnl("LbmFsgrSolver::step : "<<this->mStepCnt, 10);
79   if(!this->mSilent){ debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<<this->mStepCnt<<" t:"<<mSimulationTime, 10); }
80         //debMsgDirect(  "LbmFsgrSolver::step : "<<this->mStepCnt<<" ");
81         myTime_t timestart = getTime();
82         //myTime_t timestart = 0;
83         //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 
84
85         // time adapt
86         mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;
87
88         // init moving bc's, can change mMaxVlen
89         initMovingObstacles(false);
90 #if LBM_INCLUDE_TESTSOLVERS==1
91         handleCpdata();
92 #endif
93
94         // important - keep for tadap
95         LbmFloat lastMass = mCurrentMass;
96         mCurrentMass = this->mFixMass; // reset here for next step
97         mCurrentVolume = 0.0;
98         
99         //change to single step advance!
100         int levsteps = 0;
101         int dsbits = this->mStepCnt ^ (this->mStepCnt-1);
102         //errMsg("S"," step:"<<this->mStepCnt<<" s-1:"<<(this->mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
103         for(int lev=0; lev<=mMaxRefine; lev++) {
104                 //if(! (this->mStepCnt&(1<<lev)) ) {
105                 if( dsbits & (1<<(mMaxRefine-lev)) ) {
106                         //errMsg("S"," l"<<lev);
107
108                         if(lev==mMaxRefine) {
109                                 // always advance fine level...
110                                 fineAdvance();
111                         } else {
112                                 adaptGrid(lev);
113                                 coarseRestrictFromFine(lev);
114                                 coarseAdvance(lev);
115                         }
116 #if FSGR_OMEGA_DEBUG==1
117                         errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
118                         mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
119 #endif // FSGR_OMEGA_DEBUG==1
120                         levsteps++;
121                 }
122                 mCurrentMass   += mLevel[lev].lmass;
123                 mCurrentVolume += mLevel[lev].lvolume;
124         }
125
126   // prepare next step
127         this->mStepCnt++;
128
129
130         // some dbugging output follows
131         // calculate MLSUPS
132         myTime_t timeend = getTime();
133
134         this->mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
135         mAvgNumUsedCells += this->mNumUsedCells;
136         this->mMLSUPS = (this->mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000);
137         if(this->mMLSUPS>10000){ this->mMLSUPS = -1; }
138         else { mAvgMLSUPS += this->mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
139         
140         LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
141         if(totMLSUPS>10000) totMLSUPS = -1;
142         mNumInvIfTotal += mNumInvIfCells; // debug
143
144   // do some formatting 
145   if(!this->mSilent){ 
146                 string sepStr(""); // DEBUG
147                 int avgcls = (int)(mAvgNumUsedCells/(LONGINT)this->mStepCnt);
148         debMsgStd("LbmFsgrSolver::step", DM_MSG, this->mName<<" cnt:"<<this->mStepCnt<<" t:"<<mSimulationTime<<
149                 //debMsgDirect( 
150                         " mlsups(curr:"<<this->mMLSUPS<<
151                         " avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<<
152                         " totcls:"<<(this->mNumUsedCells)<< sepStr<<
153                         " avgcls:"<< avgcls<< sepStr<<
154                         " intd:"<<mNumInterdCells<< sepStr<<
155                         " invif:"<<mNumInvIfCells<< sepStr<<
156                         " invift:"<<mNumInvIfTotal<< sepStr<<
157                         " fsgrcs:"<<mNumFsgrChanges<< sepStr<<
158                         " filled:"<<this->mNumFilledCells<<", emptied:"<<this->mNumEmptiedCells<< sepStr<<
159                         " mMxv:"<<mMxvx<<","<<mMxvy<<","<<mMxvz<<", tscnts:"<<mTimeSwitchCounts<< sepStr<<
160                         " RWmxv:"<<ntlVec3Gfx(mMxvx,mMxvy,mMxvz)*(mLevel[mMaxRefine].simCellSize / mLevel[mMaxRefine].timestep)<<" "<< /* realworld vel output */
161                         " probs:"<<mNumProblems<< sepStr<<
162                         " simt:"<<mSimulationTime<< sepStr<<
163                         " for '"<<this->mName<<"' " , 10);
164         } else { debMsgDirect("."); }
165
166         if(this->mStepCnt==1) {
167                 mMinNoCells = mMaxNoCells = this->mNumUsedCells;
168         } else {
169                 if(this->mNumUsedCells>mMaxNoCells) mMaxNoCells = this->mNumUsedCells;
170                 if(this->mNumUsedCells<mMinNoCells) mMinNoCells = this->mNumUsedCells;
171         }
172         
173         // mass scale test
174         if((mMaxRefine>0)&&(mInitialMass>0.0)) {
175                 LbmFloat mscale = mInitialMass/mCurrentMass;
176
177                 mscale = 1.0;
178                 const LbmFloat dchh = 0.001;
179                 if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
180                 if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;
181
182                 // use mass rescaling?
183                 // with float precision this seems to be nonsense...
184                 const bool MREnable = false;
185
186                 const int MSInter = 2;
187                 static int mscount = 0;
188                 if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
189                         // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
190                         // mass rescale MASS RESCALE check
191                         errMsg("MDTDD","\n\n");
192                         errMsg("MDTDD","FORCE RESCALE MASS! "
193                                         <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
194                                         <<" step:"<<this->mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
195                                         );
196                         errMsg("MDTDD","\n\n");
197
198                         mscount++;
199                         for(int lev=mMaxRefine; lev>=0 ; lev--) {
200                                 //for(int workSet = 0; workSet<=1; workSet++) {
201                                 int wss = 0;
202                                 int wse = 1;
203 #if COMPRESSGRIDS==1
204                                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
205 #endif // COMPRESSGRIDS==1
206                                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
207
208                                         FSGR_FORIJK1(lev) {
209                                                 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
210                                                         ) {
211
212                                                         FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
213                                                         QCELL(lev, i,j,k,workSet, dMass) *= mscale;
214                                                         QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;
215
216                                                 } else {
217                                                         continue;
218                                                 }
219                                         }
220                                 }
221                                 mLevel[lev].lmass *= mscale;
222                         }
223                 } 
224
225                 mCurrentMass *= mscale;
226         }// if mass scale test */
227         else {
228                 // use current mass after full step for initial setting
229                 if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) {
230                         mInitialMass = mCurrentMass;
231                         debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2);
232                 }
233         }
234
235 #if LBM_INCLUDE_TESTSOLVERS==1
236         if((mUseTestdata)&&(this->mInitDone)) { handleTestdata(); }
237         testXchng();
238 #endif
239
240         // advance positions with current grid
241         advanceParticles();
242         if(mpParticles) {
243                 mpParticles->checkDumpTextPositions(mSimulationTime);
244                 mpParticles->checkTrails(mSimulationTime);
245         }
246
247         // one of the last things to do - adapt timestep
248         // was in fineAdvance before... 
249         if(mTimeAdap) {
250                 adaptTimestep();
251         } // time adaptivity
252
253
254 #ifndef WIN32
255         // good indicator for instabilities
256         if( (!finite(mMxvx)) || (!finite(mMxvy)) || (!finite(mMxvz)) ) { CAUSE_PANIC; }
257         if( (!finite(mCurrentMass)) || (!finite(mCurrentVolume)) ) { CAUSE_PANIC; }
258 #endif // WIN32
259
260         // output total step time
261         timeend = getTime();
262         double mdelta = (lastMass-mCurrentMass);
263         if(ABS(mdelta)<1e-12) mdelta=0.;
264         debMsgStd("LbmFsgrSolver::stepMain",DM_MSG,"step:"<<this->mStepCnt
265                         <<": dccd="<< mCurrentMass
266                         <<",d"<<mdelta
267                         <<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
268                         <<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
269                         <<" totst:"<<getTimeString(timeend-timestart), 3);
270         // nicer output
271         debMsgDirect(std::endl); 
272
273  //#endif // ELBEEM_PLUGIN!=1
274 }
275
276 #define NEWDEBCHECK(str) \
277         if(!this->mPanic){ FSGR_FORIJK_BOUNDS(mMaxRefine) { \
278                 if(RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr)&(CFFluid|CFInter)) { \
279                 for(int l=0;l<dTotalNum;l++) { \
280                         if(!finite(QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr,l))) { errMsg("NNOFIN"," "<<str<<" at "<<PRINT_IJK<<" l"<<l<<" "); }\
281                 }/*for*/ \
282                 }/*if*/ \
283         } }
284
285 void LbmFsgrSolver::fineAdvance()
286 {
287         // do the real thing...
288         //NEWDEBCHECK("t1");
289         mainLoop( mMaxRefine );
290         if(mUpdateFVHeight) {
291                 // warning assume -Y gravity...
292                 mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
293                 if(mFVHeight<1.0) mFVHeight = 1.0;
294                 this->mpParam->setFluidVolumeHeight(mFVHeight);
295         }
296
297         // advance time before timestep change
298         mSimulationTime += this->mpParam->getTimestep();
299         // time adaptivity
300         this->mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
301         //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
302         if(!this->mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); }
303
304         // update other set
305   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
306   mLevel[mMaxRefine].setCurr   ^= 1;
307   mLevel[mMaxRefine].lsteps++;
308
309         // flag init... (work on current set, to simplify flag checks)
310         reinitFlags( mLevel[mMaxRefine].setCurr );
311         if(!this->mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); }
312
313         // DEBUG VEL CHECK
314         if(0) {
315                 int lev = mMaxRefine;
316                 int workSet = mLevel[lev].setCurr;
317                 int mi=0,mj=0,mk=0;
318                 LbmFloat compRho=0.;
319                 LbmFloat compUx=0.;
320                 LbmFloat compUy=0.;
321                 LbmFloat compUz=0.;
322                 LbmFloat maxUlen=0.;
323                 LbmVec maxU(0.);
324                 LbmFloat maxRho=-100.;
325                 int ri=0,rj=0,rk=0;
326
327                 FSGR_FORIJK1(lev) {
328                         if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) ) {
329                                 compRho=QCELL(lev, i,j,k,workSet, dC);
330                                 compUx = compUy = compUz = 0.0;
331                                 for(int l=1; l<this->cDfNum; l++) {
332                                         LbmFloat df = QCELL(lev, i,j,k,workSet, l);
333                                         compRho += df;
334                                         compUx  += (this->dfDvecX[l]*df);
335                                         compUy  += (this->dfDvecY[l]*df); 
336                                         compUz  += (this->dfDvecZ[l]*df); 
337                                 } 
338                                 LbmVec u(compUx,compUy,compUz);
339                                 LbmFloat nu = norm(u);
340                                 if(nu>maxUlen) {
341                                         maxU = u;
342                                         maxUlen = nu;
343                                         mi=i; mj=j; mk=k;
344                                 }
345                                 if(compRho>maxRho) {
346                                         maxRho=compRho;
347                                         ri=i; rj=j; rk=k;
348                                 }
349                         } else {
350                                 continue;
351                         }
352                 }
353
354                 errMsg("MAXVELCHECK"," at "<<PRINT_VEC(mi,mj,mk)<<" norm:"<<maxUlen<<" u:"<<maxU);
355                 errMsg("MAXRHOCHECK"," at "<<PRINT_VEC(ri,rj,rk)<<" rho:"<<maxRho);
356                 printLbmCell(lev, 30,36,23, -1);
357         } // DEBUG VEL CHECK
358
359 }
360
361
362
363 // fine step defines
364
365 // access to own dfs during step (may be changed to local array)
366 #define MYDF(l) RAC(ccel, l)
367
368 // drop model definitions
369 #define RWVEL_THRESH 1.5
370 #define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
371
372 #if LBMDIM==3
373                         // normal
374 #define SLOWDOWNREGION (this->mSizez/4)
375 #else // LBMDIM==2
376                         // off
377 #define SLOWDOWNREGION 10 
378 #endif // LBMDIM==2
379 #define P_LCSMQO 0.01
380
381 /*****************************************************************************/
382 //! fine step function
383 /*****************************************************************************/
384 void 
385 LbmFsgrSolver::mainLoop(int lev)
386 {
387         // loops over _only inner_ cells  -----------------------------------------------------------------------------------
388         
389         // slow surf regions smooth (if below)
390         const LbmFloat smoothStrength = 0.0; //0.01;
391         const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
392         const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
393
394         const int cutMin  = 1;
395         const int cutConst = mCutoff+2;
396
397
398 #       if LBM_INCLUDE_TESTSOLVERS==1
399         // 3d region off... quit
400         if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { return; }
401 #       endif // ELBEEM_PLUGIN!=1
402         
403   // main loop region
404         const bool doReduce = true;
405         const int gridLoopBound=1;
406         GRID_REGION_INIT();
407 #if PARALLEL==1
408 #include "paraloopstart.h"
409         GRID_REGION_START();
410 #else // PARALLEL==1
411         GRID_REGION_START();
412 #endif // PARALLEL==1
413
414         // local to main
415         CellFlagType nbflag[LBM_DFNUM]; 
416         int oldFlag, newFlag, nbored;
417         LbmFloat m[LBM_DFNUM];
418         LbmFloat rho, ux, uy, uz, tmp, usqr;
419
420         // smago vars
421         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
422         
423         // ifempty cell conversion flags
424         bool iffilled, ifemptied;
425         LbmFloat nbfracs[LBM_DFNUM]; // ffracs of neighbors
426         int recons[LBM_DFNUM];   // reconstruct this DF?
427         int numRecons;           // how many are reconstructed?
428
429         LbmFloat mass=0., change=0., lcsmqo=0.;
430         rho= ux= uy= uz= usqr= tmp= 0.; 
431         lcsmqadd = lcsmomega = 0.;
432         FORDF0{ lcsmeq[l] = 0.; }
433
434         // ---
435         // now stream etc.
436         // use //template functions for 2D/3D
437
438         GRID_LOOP_START();
439                 // loop start
440                 // stream from current set to other, then collide and store
441                 //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
442
443 #               if FSGR_STRICT_DEBUG==1
444                 // safety pointer check
445                 rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
446                 if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) || 
447                     (&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
448                         errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<<
449                                         RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<","<<RFLAG(lev, i,j,k,mLevel[lev].setOther)<<" but is "<<
450                                         (*pFlagSrc)<<","<<(*pFlagDst) <<",  pointers "<<
451           (int)(&RFLAG(lev, i,j,k,mLevel[lev].setCurr))<<","<<(int)(&RFLAG(lev, i,j,k,mLevel[lev].setOther))<<" but is "<<
452           (int)(pFlagSrc)<<","<<(int)(pFlagDst)<<" "
453                                         ); 
454                         CAUSE_PANIC;
455                 }       
456                 if( (&QCELL(lev, i,j,k,mLevel[lev].setCurr,0) != ccel) || 
457                     (&QCELL(lev, i,j,k,mLevel[lev].setOther,0) != tcel) ) {
458                         errMsg("LbmFsgrSolver::mainLoop","Err cellp "<<PRINT_IJK<<"="<<
459           (int)(&QCELL(lev, i,j,k,mLevel[lev].setCurr,0))<<","<<(int)(&QCELL(lev, i,j,k,mLevel[lev].setOther,0))<<" but is "<<
460           (int)(ccel)<<","<<(int)(tcel)<<" "
461                                         ); 
462                         CAUSE_PANIC;
463                 }       
464 #               endif
465                 oldFlag = *pFlagSrc;
466                 
467                 // old INTCFCOARSETEST==1
468                 if( (oldFlag & (CFGrFromCoarse)) ) { 
469                         if(( this->mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
470                                 FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
471                         } else {
472                                 interpolateCellFromCoarse( lev, i,j,k, TSET(lev), 0.0, CFFluid|CFGrFromCoarse, false);
473                                 calcNumUsedCells++;
474                         }
475                         continue; // interpolateFineFromCoarse test!
476                 } // interpolateFineFromCoarse test! 
477         
478                 if(oldFlag & (CFMbndInflow)) {
479                         // fluid & if are ok, fill if later on
480                         int isValid = oldFlag & (CFFluid|CFInter);
481                         const LbmFloat iniRho = 1.0;
482                         const int OId = oldFlag>>24;
483                         if(!isValid) {
484                                 // make new if cell
485                                 const LbmVec vel(mObjectSpeeds[OId]);
486                                 // TODO add OPT3D treatment
487                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
488                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
489                                 RAC(tcel, dFlux) = FLUX_INIT;
490                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
491                                 calcCurrentMass += iniRho; 
492                                 calcCurrentVolume += 1.0; 
493                                 calcNumUsedCells++;
494                                 mInitialMass += iniRho;
495                                 // dont treat cell until next step
496                                 continue;
497                         } 
498                 } 
499                 else  // these are exclusive
500                 if(oldFlag & (CFMbndOutflow)) {
501                         int isnotValid = oldFlag & (CFFluid);
502                         if(isnotValid) {
503                                 // remove fluid cells, shouldnt be here anyway
504                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
505                                 mInitialMass -= fluidRho;
506                                 const LbmFloat iniRho = 0.0;
507                                 RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
508                                 RAC(tcel, dFlux) = FLUX_INIT;
509                                 changeFlag(lev, i,j,k, TSET(lev), CFInter);
510
511                                 // same as ifemptied for if below
512                                 LbmPoint oemptyp; oemptyp.flag = 0;
513                                 oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
514                                 LIST_EMPTY(oemptyp);
515                                 calcCellsEmptied++;
516                                 continue;
517                         }
518                 }
519
520                 if(oldFlag & (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)) { 
521                         *pFlagDst = oldFlag;
522                         continue;
523                 }
524                 /*if( oldFlag & CFNoBndFluid ) {  // TEST ME FASTER?
525                         OPTIMIZED_STREAMCOLLIDE; PERFORM_USQRMAXCHECK;
526                         RAC(tcel,dFfrac) = 1.0; 
527                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
528                         calcCurrentMass += rho; calcCurrentVolume += 1.0;
529                         calcNumUsedCells++;
530                         continue;
531                 }// TEST ME FASTER? */
532
533                 // only neighbor flags! not own flag
534                 nbored = 0;
535                 
536 #if OPT3D==0
537                 FORDF1 {
538                         nbflag[l] = RFLAG_NB(lev, i,j,k,SRCS(lev),l);
539                         nbored |= nbflag[l];
540                 } 
541 #else
542                 nbflag[dSB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dSB];
543                 nbflag[dWB] = *(pFlagSrc + (-mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWB];
544                 nbflag[ dB] = *(pFlagSrc + (-mLevel[lev].lOffsy)); nbored |= nbflag[dB];
545                 nbflag[dEB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dEB];
546                 nbflag[dNB] = *(pFlagSrc + (-mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNB];
547
548                 nbflag[dSW] = *(pFlagSrc + (-mLevel[lev].lOffsx+-1)); nbored |= nbflag[dSW];
549                 nbflag[ dS] = *(pFlagSrc + (-mLevel[lev].lOffsx)); nbored |= nbflag[dS];
550                 nbflag[dSE] = *(pFlagSrc + (-mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dSE];
551
552                 nbflag[ dW] = *(pFlagSrc + (-1)); nbored |= nbflag[dW];
553                 nbflag[ dE] = *(pFlagSrc + ( 1)); nbored |= nbflag[dE];
554
555                 nbflag[dNW] = *(pFlagSrc + ( mLevel[lev].lOffsx+-1)); nbored |= nbflag[dNW];
556           nbflag[ dN] = *(pFlagSrc + ( mLevel[lev].lOffsx)); nbored |= nbflag[dN];
557                 nbflag[dNE] = *(pFlagSrc + ( mLevel[lev].lOffsx+ 1)); nbored |= nbflag[dNE];
558
559                 nbflag[dST] = *(pFlagSrc + ( mLevel[lev].lOffsy+-mLevel[lev].lOffsx)); nbored |= nbflag[dST];
560                 nbflag[dWT] = *(pFlagSrc + ( mLevel[lev].lOffsy+-1)); nbored |= nbflag[dWT];
561                 nbflag[ dT] = *(pFlagSrc + ( mLevel[lev].lOffsy)); nbored |= nbflag[dT];
562                 nbflag[dET] = *(pFlagSrc + ( mLevel[lev].lOffsy+ 1)); nbored |= nbflag[dET];
563                 nbflag[dNT] = *(pFlagSrc + ( mLevel[lev].lOffsy+ mLevel[lev].lOffsx)); nbored |= nbflag[dNT];
564                 // */
565 #endif
566
567                 // pointer to destination cell
568                 calcNumUsedCells++;
569
570                 // FLUID cells 
571                 if( oldFlag & CFFluid ) { 
572                         // only standard fluid cells (with nothing except fluid as nbs
573
574                         if(oldFlag&CFMbndInflow) {
575                                 // force velocity for inflow, necessary to have constant direction of flow
576                                 // FIXME , test also set interface cells?
577                                 const int OId = oldFlag>>24;
578                                 //? DEFAULT_STREAM;
579                                 //const LbmFloat fluidRho = 1.0;
580                                 // for submerged inflows, streaming would have to be performed...
581                                 LbmFloat fluidRho = m[0]; FORDF1 { fluidRho += m[l]; }
582                                 const LbmVec vel(mObjectSpeeds[OId]);
583                                 ux=vel[0], uy=vel[1], uz=vel[2]; 
584                                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
585                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); }
586                                 rho = fluidRho;
587                                 //errMsg("INFLOW_DEBUG","std at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
588                         } else {
589                                 if(nbored&CFBnd) {
590                                         DEFAULT_STREAM;
591                                         //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
592                                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
593                                         oldFlag &= (~CFNoBndFluid);
594                                 } else {
595                                         // do standard stream/collide
596                                         OPTIMIZED_STREAMCOLLIDE;
597                                         oldFlag |= CFNoBndFluid;
598                                 } 
599                         }
600
601                         PERFORM_USQRMAXCHECK;
602                         // "normal" fluid cells
603                         RAC(tcel,dFfrac) = 1.0; 
604                         *pFlagDst = (CellFlagType)oldFlag; // newFlag;
605                         calcCurrentMass += rho; 
606                         calcCurrentVolume += 1.0;
607                         continue;
608                 }
609                 
610                 newFlag  = oldFlag;
611                 // make sure here: always check which flags to really unset...!
612                 newFlag = newFlag & (~( 
613                                         CFNoNbFluid|CFNoNbEmpty| CFNoDelete 
614                                         | CFNoInterpolSrc
615                                         | CFNoBndFluid
616                                         ));
617                 if(!(nbored&CFBndNoslip)) { //NEWSURFT NEWSURFTNOS
618                         newFlag |= CFNoBndFluid;
619                 }
620                 /*if(!(nbored&CFBnd)) { //NEWSURFT NEWSURFTNOS
621                         // explicitly check for noslip neighbors
622                         bool hasnoslipnb = false;
623                         FORDF1 { if((nbflag[l]&CFBnd)&&(nbflag[l]&CFBndNoslip)) hasnoslipnb=true; }
624                         if(!hasnoslipnb) newFlag |= CFNoBndFluid; 
625                 } // */
626
627                 // store own dfs and mass
628                 mass = RAC(ccel,dMass);
629
630                 // WARNING - only interface cells arrive here!
631                 // read distribution funtions of adjacent cells = stream step
632                 DEFAULT_STREAM;
633
634                 if((nbored & CFFluid)==0) { newFlag |= CFNoNbFluid; mNumInvIfCells++; }
635                 if((nbored & CFEmpty)==0) { newFlag |= CFNoNbEmpty; mNumInvIfCells++; }
636
637                 // calculate mass exchange for interface cells 
638                 LbmFloat myfrac = RAC(ccel,dFfrac);
639                 if(myfrac<0.) myfrac=0.; // NEWSURFT
640 #               define nbdf(l) m[ this->dfInv[(l)] ]
641
642                 // update mass 
643                 // only do boundaries for fluid cells, and interface cells without
644                 // any fluid neighbors (assume that interface cells _with_ fluid
645                 // neighbors are affected enough by these) 
646                 // which Df's have to be reconstructed? 
647                 // for fluid cells - just the f_i difference from streaming to empty cells  ----
648                 numRecons = 0;
649                 bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
650                 //onlyBndnb = false; // DEBUG test off
651
652                 FORDF1 { // dfl loop
653                         recons[l] = 0;
654                         nbfracs[l] = 0.0;
655                         // finally, "normal" interface cells ----
656                         if( nbflag[l]&(CFFluid|CFBnd) ) { // NEWTEST! FIXME check!!!!!!!!!!!!!!!!!!
657                                 change = nbdf(l) - MYDF(l);
658                         }
659                         // interface cells - distuingish cells that shouldn't fill/empty 
660                         else if( nbflag[l] & CFInter ) {
661                                 
662                                 LbmFloat mynbfac,nbnbfac;
663                                 // NEW TEST t1
664                                 // t2 -> off
665                                 if((oldFlag&CFNoBndFluid)&&(nbflag[l]&CFNoBndFluid)) {
666                                         mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
667                                         nbnbfac = 1.0/mynbfac;
668                                         onlyBndnb = false;
669                                 } else {
670                                         mynbfac = nbnbfac = 1.0; // switch calc flux off
671                                         goto changeDefault;  // NEWSURFT
672                                         //change = 0.; goto changeDone;  // NEWSURFT
673                                 }
674                                 //mynbfac = nbnbfac = 1.0; // switch calc flux off t3
675
676                                 // perform interface case handling
677                                 if ((oldFlag|nbflag[l])&(CFNoNbFluid|CFNoNbEmpty)) {
678                                 switch (oldFlag&(CFNoNbFluid|CFNoNbEmpty)) {
679                                         case 0: 
680                                                 // we are a normal cell so... 
681                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
682                                                         case CFNoNbFluid: 
683                                                                 // just fill current cell = empty neighbor 
684                                                                 change = nbnbfac*nbdf(l) ; goto changeDone; 
685                                                         case CFNoNbEmpty: 
686                                                                 // just empty current cell = fill neighbor 
687                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
688                                                 }
689                                                 break;
690
691                                         case CFNoNbFluid: 
692                                                 // we dont have fluid nb's so...
693                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
694                                                         case 0: 
695                                                         case CFNoNbEmpty: 
696                                                                 // we have no fluid nb's -> just empty
697                                                                 change = - mynbfac*MYDF(l) ; goto changeDone; 
698                                                 }
699                                                 break;
700
701                                         case CFNoNbEmpty: 
702                                                 // we dont have empty nb's so...
703                                                 switch (nbflag[l]&(CFNoNbFluid|CFNoNbEmpty)) {
704                                                         case 0: 
705                                                         case CFNoNbFluid: 
706                                                                 // we have no empty nb's -> just fill
707                                                                 change = nbnbfac*nbdf(l); goto changeDone; 
708                                                 }
709                                                 break;
710                                 }} // inter-inter exchange
711
712                         changeDefault: ;
713                                 // just do normal mass exchange...
714                                 change = ( nbnbfac*nbdf(l) - mynbfac*MYDF(l) ) ;
715                         changeDone: ;
716                                 nbfracs[l] = QCELL_NB(lev, i,j,k, SRCS(lev),l, dFfrac);
717                                 if(nbfracs[l]<0.) nbfracs[l] = 0.; // NEWSURFT
718                                 change *=  (myfrac + nbfracs[l]) * 0.5;
719                         } // the other cell is interface
720
721                         // last alternative - reconstruction in this direction
722                         else {
723                                 // empty + bnd case
724                                 recons[l] = 1; 
725                                 numRecons++;
726                                 change = 0.0; 
727                         }
728
729                         // modify mass at SRCS
730                         mass += change;
731                 } // l
732                 // normal interface, no if empty/fluid
733
734                 LbmFloat nv1,nv2;
735                 LbmFloat nx,ny,nz;
736
737                 if(nbflag[dE] &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0;
738                 if(nbflag[dW] &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0;
739                 nx = 0.5* (nv2-nv1);
740                 if(nbflag[dN] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
741                 if(nbflag[dS] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
742                 ny = 0.5* (nv2-nv1);
743 #               if LBMDIM==3
744                 if(nbflag[dT] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
745                 if(nbflag[dB] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
746                 nz = 0.5* (nv2-nv1);
747 #               else // LBMDIM==3
748                 nz = 0.0;
749 #               endif // LBMDIM==3
750
751                 if( (ABS(nx)+ABS(ny)+ABS(nz)) > LBM_EPSILON) {
752                         // normal ok and usable...
753                         FORDF1 {
754                                 if( (this->dfDvecX[l]*nx + this->dfDvecY[l]*ny + this->dfDvecZ[l]*nz)  // dot Dvec,norml
755                                                 > LBM_EPSILON) {
756                                         recons[l] = 2; 
757                                         numRecons++;
758                                 } 
759                         }
760                 }
761
762                 // calculate macroscopic cell values
763                 LbmFloat oldUx, oldUy, oldUz;
764                 LbmFloat oldRho; // OLD rho = ccel->rho;
765 #               define REFERENCE_PRESSURE 1.0 // always atmosphere...
766 #               if OPT3D==0
767                 oldRho=RAC(ccel,0);
768                 oldUx = oldUy = oldUz = 0.0;
769                 for(int l=1; l<this->cDfNum; l++) {
770                         oldRho += RAC(ccel,l);
771                         oldUx  += (this->dfDvecX[l]*RAC(ccel,l));
772                         oldUy  += (this->dfDvecY[l]*RAC(ccel,l)); 
773                         oldUz  += (this->dfDvecZ[l]*RAC(ccel,l)); 
774                 } 
775                 // reconstruct dist funcs from empty cells
776                 FORDF1 {
777                         if(recons[ l ]) {
778                                 m[ this->dfInv[l] ] = 
779                                         this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) + 
780                                         this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz) 
781                                         - MYDF( l );
782                         }
783                 }
784                 usqr = 1.5 * (oldUx*oldUx + oldUy*oldUy + oldUz*oldUz); // needed later on
785 #               else // OPT3D==0
786                 oldRho = + RAC(ccel,dC)  + RAC(ccel,dN )
787                                 + RAC(ccel,dS ) + RAC(ccel,dE )
788                                 + RAC(ccel,dW ) + RAC(ccel,dT )
789                                 + RAC(ccel,dB ) + RAC(ccel,dNE)
790                                 + RAC(ccel,dNW) + RAC(ccel,dSE)
791                                 + RAC(ccel,dSW) + RAC(ccel,dNT)
792                                 + RAC(ccel,dNB) + RAC(ccel,dST)
793                                 + RAC(ccel,dSB) + RAC(ccel,dET)
794                                 + RAC(ccel,dEB) + RAC(ccel,dWT)
795                                 + RAC(ccel,dWB);
796
797                 oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
798                                 + RAC(ccel,dNE) - RAC(ccel,dNW)
799                                 + RAC(ccel,dSE) - RAC(ccel,dSW)
800                                 + RAC(ccel,dET) + RAC(ccel,dEB)
801                                 - RAC(ccel,dWT) - RAC(ccel,dWB);
802
803                 oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
804                                 + RAC(ccel,dNE) + RAC(ccel,dNW)
805                                 - RAC(ccel,dSE) - RAC(ccel,dSW)
806                                 + RAC(ccel,dNT) + RAC(ccel,dNB)
807                                 - RAC(ccel,dST) - RAC(ccel,dSB);
808
809                 oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
810                                 + RAC(ccel,dNT) - RAC(ccel,dNB)
811                                 + RAC(ccel,dST) - RAC(ccel,dSB)
812                                 + RAC(ccel,dET) - RAC(ccel,dEB)
813                                 + RAC(ccel,dWT) - RAC(ccel,dWB);
814
815                 // now reconstruction
816                 ux=oldUx, uy=oldUy, uz=oldUz;  // no local vars, only for usqr
817                 rho = REFERENCE_PRESSURE;
818                 usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
819                 if(recons[dN ]) { m[dS ] = EQN  + EQS  - MYDF(dN ); }
820                 if(recons[dS ]) { m[dN ] = EQS  + EQN  - MYDF(dS ); }
821                 if(recons[dE ]) { m[dW ] = EQE  + EQW  - MYDF(dE ); }
822                 if(recons[dW ]) { m[dE ] = EQW  + EQE  - MYDF(dW ); }
823                 if(recons[dT ]) { m[dB ] = EQT  + EQB  - MYDF(dT ); }
824                 if(recons[dB ]) { m[dT ] = EQB  + EQT  - MYDF(dB ); }
825                 if(recons[dNE]) { m[dSW] = EQNE + EQSW - MYDF(dNE); }
826                 if(recons[dNW]) { m[dSE] = EQNW + EQSE - MYDF(dNW); }
827                 if(recons[dSE]) { m[dNW] = EQSE + EQNW - MYDF(dSE); }
828                 if(recons[dSW]) { m[dNE] = EQSW + EQNE - MYDF(dSW); }
829                 if(recons[dNT]) { m[dSB] = EQNT + EQSB - MYDF(dNT); }
830                 if(recons[dNB]) { m[dST] = EQNB + EQST - MYDF(dNB); }
831                 if(recons[dST]) { m[dNB] = EQST + EQNB - MYDF(dST); }
832                 if(recons[dSB]) { m[dNT] = EQSB + EQNT - MYDF(dSB); }
833                 if(recons[dET]) { m[dWB] = EQET + EQWB - MYDF(dET); }
834                 if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
835                 if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
836                 if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
837 #               endif           
838
839
840                 // inflow bc handling
841                 if(oldFlag & (CFMbndInflow)) {
842                         // fill if cells in inflow region
843                         if(myfrac<0.5) { 
844                                 mass += 0.25; 
845                                 mInitialMass += 0.25;
846                         }
847                         const int OId = oldFlag>>24;
848                         const LbmVec vel(mObjectSpeeds[OId]);
849                         ux=vel[0], uy=vel[1], uz=vel[2]; 
850                         //? usqr = 1.5 * (ux*ux + uy*uy + uz*uz);
851                         //FORDF0 { RAC(tcel, l) = this->getCollideEq(l, fluidRho,ux,uy,uz); } rho = fluidRho;
852                         rho = REFERENCE_PRESSURE;
853                         FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
854                         //errMsg("INFLOW_DEBUG","if at "<<PRINT_IJK<<" v="<<vel<<" rho="<<rho);
855                 }  else { 
856                         // NEWSURFT, todo optimize!
857                         if(onlyBndnb) { //if(usqr<0.001*0.001) {
858                                 rho = ux = uy = uz = 0.;
859                                 FORDF0{ 
860                                         rho += m[l]; 
861                                         ux  += (this->dfDvecX[l]*m[l]); 
862                                         uy  += (this->dfDvecY[l]*m[l]);  
863                                         uz  += (this->dfDvecZ[l]*m[l]);  
864                                 }
865                                 FORDF0 { RAC(tcel, l) = this->getCollideEq(l, rho,ux,uy,uz); }
866                         } else {// NEWSURFT */
867                                 // normal collide
868                                 // mass streaming done... do normal collide
869                                 LbmVec grav = mLevel[lev].gravity*mass;
870                                 DEFAULT_COLLIDEG(grav);
871                                 PERFORM_USQRMAXCHECK;
872                                 // rho init from default collide necessary for fill/empty check below
873                         } // test
874                 }
875
876                 // testing..., particle generation
877                 // also check oldFlag for CFNoNbFluid, test
878                 // for inflow no pargen test // NOBUBBB!
879                 if((this->mInitDone) //&&(mUseTestdata) 
880                                 //&& (!((oldFlag|newFlag)&CFNoNbEmpty)) 
881                                 && (!((oldFlag|newFlag)&CFNoDelete)) 
882                                 && (this->mPartGenProb>0.0)) {
883                         bool doAdd = true;
884                         bool bndOk=true;
885                         if( (i<cutMin)||(i>this->mSizex-cutMin)||
886                                         (j<cutMin)||(j>this->mSizey-cutMin)||
887                                         (k<cutMin)||(k>this->mSizez-cutMin) ) { bndOk=false; }
888                         if(!bndOk) doAdd=false;
889                         
890                         LbmFloat realWorldFac = (mLevel[lev].simCellSize / mLevel[lev].timestep);
891                         LbmFloat rux = (ux * realWorldFac);
892                         LbmFloat ruy = (uy * realWorldFac);
893                         LbmFloat ruz = (uz * realWorldFac);
894                         LbmFloat rl = norm(ntlVec3Gfx(rux,ruy,ruz));
895                         // WHMOD
896
897                         LbmFloat prob = (rand()/(RAND_MAX+1.0));
898                         LbmFloat basethresh = this->mPartGenProb*lcsmqo*rl;
899
900                         // reduce probability in outer region?
901                         const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst;
902                         const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst;
903                         LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord);
904                         LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord);
905                         if(pifac<0.) pifac=0.;
906                         if(pjfac<0.) pjfac=0.;
907
908                         //if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
909                         if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
910                                 // add
911                         } else {
912                                 doAdd = false; // dont...
913                         }
914
915                         // "wind" disturbance
916                         // use realworld relative velocity here instead?
917                         if( (doAdd && 
918                                         ((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
919                                         ||(k>this->mSizez-SLOWDOWNREGION)   ) {
920                                 LbmFloat nuz = uz;
921                                 if(k>this->mSizez-SLOWDOWNREGION) {
922                                         // special case
923                                         LbmFloat zfac = (LbmFloat)( k-(this->mSizez-SLOWDOWNREGION) );
924                                         zfac /= (LbmFloat)(SLOWDOWNREGION);
925                                         nuz += (1.0) * zfac; // check max speed? OFF?
926                                         //errMsg("TOPT"," at "<<PRINT_IJK<<" zfac"<<zfac);
927                                 } else {
928                                         // normal probability
929                                         //? LbmFloat fac = P_LCSMQO-lcsmqo;
930                                         //? jdf *= fac;
931                                 }
932                                 FORDF1 {
933                                         const LbmFloat jdf = 0.05 * (rand()/(RAND_MAX+1.0));
934                                         // TODO  use wind velocity?
935                                         if(jdf>0.025) {
936                                         const LbmFloat add =  this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf;
937                                         RAC(tcel,l) += add; }
938                                 }
939                                 //errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
940                         } // wind disturbance
941
942                         if(usqr<0.0001) doAdd=false;   // TODO check!?
943                         // if outside, and 20% above sea level, delete, TODO really check level?
944                         //if((!bndOk)&&((LbmFloat)k>pTest->mFluidHeight*1.5)) { doAdd=true; } // FORCEDISSOLVE
945                         //if(this->mStepCnt>700) errMsg("DFJITT"," at "<<PRINT_IJK<<"rwl:"<<rl<<"  usqr:"<<usqr <<" qo:"<<lcsmqo<<" add="<<doAdd );
946
947                         if( (doAdd)  ) { // ADD DROP
948                                 LbmFloat len = norm(LbmVec(ux,uy,uz));
949                                 // WHMOD
950                                 //for(int s=0; s<10; s++) { // multiple parts
951                                 for(int s=0; s<1; s++) { // one part!
952                                 //LbmFloat prob = this->mPartGenProb * 0.02* (rand()/(RAND_MAX+1.0));
953                                 const LbmFloat posjitter = 1.0;
954                                 const LbmFloat posjitteroffs = posjitter*-0.5;
955                                 LbmFloat jpx = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
956                                 LbmFloat jpy = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
957                                 LbmFloat jpz = posjitteroffs+ posjitter* (rand()/(RAND_MAX+1.0));
958
959                                 const LbmFloat jitterstr = 0.1;
960                                 const LbmFloat jitteroffs = jitterstr*-0.5;
961                                 LbmFloat jx = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
962                                 LbmFloat jy = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
963                                 LbmFloat jz = jitteroffs+ jitterstr* (rand()/(RAND_MAX+1.0));
964
965                                 // average normal & velocity 
966                                 // -> mostly along velocity dir, many into surface
967                                 LbmVec pv = (LbmVec(nx+jx,ny+jy,nz+jz)*0.75 + getNormalized(LbmVec(ux,uy,uz)) )*0.35; 
968                                 normalize(pv);
969
970                                 LbmFloat srci = i+0.5+jpx; // TEST? + (pv[0]*1.41);
971                                 LbmFloat srcj = j+0.5+jpy; // TEST? + (pv[1]*1.41);
972                                 LbmFloat srck = k+0.5+jpz; // TEST? + (pv[2]*1.41);
973                                 int type=0;
974                                 //if((s%3)!=2) {} else { type=PART_FLOAT; }
975                                 //type = PART_DROP;
976                                 type = PART_INTER;
977                                 // drop
978                                 /*srci += (pv[0]*1.41);
979                                 srcj += (pv[1]*1.41);
980                                 srck += (pv[2]*1.41);
981                                 if(!(RFLAG(lev, (int)(srci),(int)(srcj),(int)(srck),SRCS(lev)) &CFEmpty)) continue; // only add in good direction */
982
983                                 pv *= len;
984                                 LbmFloat size = 1.0+ 9.0* (rand()/(RAND_MAX+1.0));
985
986                                 mpParticles->addParticle(srci, srcj, srck); //i+0.5+jpx,j+0.5+jpy,k+0.5+jpz);
987                                 mpParticles->getLast()->setVel(pv[0],pv[1],pv[2]);
988                                 //? mpParticles->getLast()->advanceVel(); // advance a bit outwards
989                                 mpParticles->getLast()->setStatus(PART_IN);
990                                 mpParticles->getLast()->setType(type);
991                                 //if((s%3)==2) mpParticles->getLast()->setType(PART_FLOAT);
992                                 mpParticles->getLast()->setSize(size);
993                                 //errMsg("NEWPART"," at "<<PRINT_IJK<<"   u="<<norm(LbmVec(ux,uy,uz)) <<" RWu="<<norm(LbmVec(rux,ruy,ruz))<<" add"<<doAdd<<" pvel="<<norm(pv) );
994                                 //mass -= size*0.001; // NTEST!
995                                 mass -= size*0.0020; // NTEST!
996 #                               if LBMDIM==2
997                                 mpParticles->getLast()->setVel(pv[0],pv[1],0.0);
998                                 mpParticles->getLast()->setPos(ntlVec3Gfx(srci,srcj,0.5));
999 #                               endif // LBMDIM==2
1000                 //errMsg("PIT","NEW pit"<<mpParticles->getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<<convertFlags2String(mpParticles->getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel()  );
1001                                 } // multiple parts
1002                         } // doAdd
1003                 } // */
1004
1005
1006                 // interface cell filled or emptied?
1007                 iffilled = ifemptied = 0;
1008                 // interface cells empty/full?, WARNING: to mark these cells, better do it at the end of reinitCellFlags
1009                 // interface cell if full?
1010                 if( (mass) >= (rho * (1.0+FSGR_MAGICNR)) ) { iffilled = 1; }
1011                 // interface cell if empty?
1012                 if( (mass) <= (rho * (   -FSGR_MAGICNR)) ) { ifemptied = 1; }
1013
1014                 if(oldFlag & (CFMbndOutflow)) {
1015                         mInitialMass -= mass;
1016                         mass = myfrac = 0.0;
1017                         iffilled = 0; ifemptied = 1;
1018                 }
1019
1020                 // looks much nicer... LISTTRICK
1021 #               if FSGR_LISTTRICK==1
1022                 //if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
1023                 if(newFlag&CFNoBndFluid) { // test NEW TEST
1024                         if(!iffilled) {
1025                                 // remove cells independent from amount of change...
1026                                 if( (oldFlag & CFNoNbEmpty)&&(newFlag & CFNoNbEmpty)&&
1027                                                 ( (mass>(rho*FSGR_LISTTTHRESHFULL))  || ((nbored&CFInter)==0)  )) { 
1028                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","filled "<<PRINT_IJK); };
1029                                         iffilled = 1; 
1030                                 } 
1031                         }
1032                         if(!ifemptied) {
1033                                 if( (oldFlag & CFNoNbFluid)&&(newFlag & CFNoNbFluid)&&
1034                                                 ( (mass<(rho*FSGR_LISTTTHRESHEMPTY)) || ((nbored&CFInter)==0)  )) { 
1035                                         //if((nbored&CFInter)==0){ errMsg("NBORED!CFINTER","emptied "<<PRINT_IJK); };
1036                                         ifemptied = 1; 
1037                                 } 
1038                         }
1039                 } // nobndfluid only */
1040 #               endif
1041                 //iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
1042                 
1043
1044                 // now that all dfs are known, handle last changes
1045                 if(iffilled) {
1046                         LbmPoint filledp; filledp.flag=0;
1047                         if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1;  // NEWSURFT
1048                         filledp.x = i; filledp.y = j; filledp.z = k;
1049                         LIST_FULL(filledp);
1050                         //this->mNumFilledCells++; // DEBUG
1051                         calcCellsFilled++;
1052                 }
1053                 else if(ifemptied) {
1054                         LbmPoint emptyp; emptyp.flag=0;
1055                         if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; //  NEWSURFT
1056                         emptyp.x = i; emptyp.y = j; emptyp.z = k;
1057                         LIST_EMPTY(emptyp);
1058                         //this->mNumEmptiedCells++; // DEBUG
1059                         calcCellsEmptied++;
1060                 } 
1061                 // dont cutoff values -> better cell conversions
1062                 RAC(tcel,dFfrac)   = (mass/rho);
1063
1064                 // init new flux value
1065                 float flux = FLUX_INIT; // dxqn on
1066                 if(newFlag&CFNoBndFluid) {
1067                         //flux = 50.0; // extreme on
1068                         for(int nn=1; nn<this->cDfNum; nn++) { 
1069                                 if(nbflag[nn] & (CFFluid|CFInter|CFBnd)) { flux += this->dfLength[nn]; }
1070                         }
1071                         // optical hack - smooth slow moving
1072                         // surface regions
1073                         if(usqr< sssUsqrLimit) {
1074                         for(int nn=1; nn<this->cDfNum; nn++) { 
1075                                 if(nbfracs[nn]!=0.0) {
1076                                         LbmFloat curSmooth = (sssUsqrLimit-usqr)*sssUsqrLimitInv;
1077                                         if(curSmooth>1.0) curSmooth=1.0;
1078                                         flux *= (1.0+ smoothStrength*curSmooth * (nbfracs[nn]-myfrac)) ;
1079                                 }
1080                         } }
1081                         // NEW TEST */
1082                 }
1083                 // flux = FLUX_INIT; // calc flux off
1084                 QCELL(lev, i,j,k,TSET(lev), dFlux) = flux; // */
1085
1086                 // perform mass exchange with streamed values
1087                 QCELL(lev, i,j,k,TSET(lev), dMass) = mass; // MASST
1088                 // set new flag 
1089                 *pFlagDst = (CellFlagType)newFlag;
1090                 calcCurrentMass += mass; 
1091                 calcCurrentVolume += RAC(tcel,dFfrac);
1092
1093                 // interface cell handling done...
1094
1095 //#if PARALLEL==1
1096 //#include "paraloopend.h"
1097         //GRID_REGION_END();
1098 //#else // PARALLEL==1
1099         //GRID_LOOPREG_END();
1100         //GRID_REGION_END();
1101 //#endif // PARALLEL==1
1102 #if PARALLEL!=1
1103         GRID_LOOPREG_END();
1104 #else // PARALLEL==1
1105 #include "paraloopend.h" // = GRID_LOOPREG_END();
1106 #endif // PARALLEL==1
1107
1108         // write vars from computations to class
1109         mLevel[lev].lmass    = calcCurrentMass;
1110         mLevel[lev].lvolume  = calcCurrentVolume;
1111         this->mNumFilledCells  = calcCellsFilled;
1112         this->mNumEmptiedCells = calcCellsEmptied;
1113         this->mNumUsedCells = calcNumUsedCells;
1114 }
1115
1116
1117
1118 void 
1119 LbmFsgrSolver::preinitGrids()
1120 {
1121         const int lev = mMaxRefine;
1122         const bool doReduce = false;
1123         const int gridLoopBound=0;
1124
1125         // touch both grids
1126         for(int s=0; s<2; s++) {
1127         
1128         GRID_REGION_INIT();
1129 #if PARALLEL==1
1130 #include "paraloopstart.h"
1131 #endif // PARALLEL==1
1132         GRID_REGION_START();
1133         GRID_LOOP_START();
1134                 FORDF0{ RAC(ccel,l) = 0.; }
1135                 *pFlagSrc =0;
1136                 *pFlagDst =0;
1137                 //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
1138 #if PARALLEL!=1
1139         GRID_LOOPREG_END();
1140 #else // PARALLEL==1
1141 #include "paraloopend.h" // = GRID_LOOPREG_END();
1142 #endif // PARALLEL==1
1143         //GRID_REGION_END();
1144         // TEST! */
1145
1146         /* dummy remove warnings */ 
1147         calcCurrentMass = calcCurrentVolume = 0.;
1148         calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
1149         
1150         // change grid
1151   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
1152   mLevel[mMaxRefine].setCurr   ^= 1;
1153         }
1154 }
1155
1156 void 
1157 LbmFsgrSolver::standingFluidPreinit()
1158 {
1159         const int lev = mMaxRefine;
1160         const bool doReduce = false;
1161         const int gridLoopBound=1;
1162
1163         GRID_REGION_INIT();
1164 #if PARALLEL==1
1165 #include "paraloopstart.h"
1166 #endif // PARALLEL==1
1167         GRID_REGION_START();
1168
1169         LbmFloat rho, ux,uy,uz, usqr; 
1170         CellFlagType nbflag[LBM_DFNUM];
1171         LbmFloat m[LBM_DFNUM];
1172         LbmFloat lcsmqo;
1173 #       if OPT3D==1 
1174         LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
1175         CellFlagType nbored=0;
1176 #       endif // OPT3D==true 
1177
1178         GRID_LOOP_START();
1179                 //FORDF0{ RAC(ccel,l) = 0.; }
1180                 //*pFlagSrc =0;
1181                 //*pFlagDst =0;
1182                 //errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
1183                 const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
1184                 if( (currFlag & (CFEmpty|CFBnd)) ) continue;
1185                 //ccel = RACPNT(lev, i,j,k,workSet); 
1186                 //tcel = RACPNT(lev, i,j,k,otherSet);
1187
1188                 if( (currFlag & (CFInter)) ) {
1189                         // copy all values
1190                         for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
1191                         continue;
1192                 }
1193
1194                 if( (currFlag & CFNoBndFluid)) {
1195                         OPTIMIZED_STREAMCOLLIDE;
1196                 } else {
1197                         FORDF1 {
1198                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
1199                         } 
1200                         DEFAULT_STREAM;
1201                         //ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
1202                         DEFAULT_COLLIDEG(mLevel[lev].gravity);
1203                 }
1204                 for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
1205 #if PARALLEL!=1
1206         GRID_LOOPREG_END();
1207 #else // PARALLEL==1
1208 #include "paraloopend.h" // = GRID_LOOPREG_END();
1209 #endif // PARALLEL==1
1210         //GRID_REGION_END();
1211         // TEST! */
1212
1213         /* dummy remove warnings */ 
1214         calcCurrentMass = calcCurrentVolume = 0.;
1215         calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
1216         
1217         // change grid
1218   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
1219   mLevel[mMaxRefine].setCurr   ^= 1;
1220 }
1221
1222
1223 /******************************************************************************
1224  * work on lists from updateCellMass to reinit cell flags
1225  *****************************************************************************/
1226
1227 LbmFloat LbmFsgrSolver::getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l) {
1228         //return 0.0; // test
1229         int level = mMaxRefine;
1230         LbmFloat *ccel = RACPNT(level, i,j,k, workSet);
1231
1232         LbmFloat nx,ny,nz, nv1,nv2;
1233         if(RFLAG_NB(level,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = RAC((ccel+QCELLSTEP ),dFfrac); } else nv1 = 0.0;
1234         if(RFLAG_NB(level,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = RAC((ccel-QCELLSTEP ),dFfrac); } else nv2 = 0.0;
1235         nx = 0.5* (nv2-nv1);
1236         if(RFLAG_NB(level,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
1237         if(RFLAG_NB(level,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
1238         ny = 0.5* (nv2-nv1);
1239 #if LBMDIM==3
1240         if(RFLAG_NB(level,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
1241         if(RFLAG_NB(level,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[level].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
1242         nz = 0.5* (nv2-nv1);
1243 #else //LBMDIM==3
1244         nz = 0.0;
1245 #endif //LBMDIM==3
1246         LbmFloat scal = mDvecNrm[l][0]*nx + mDvecNrm[l][1]*ny + mDvecNrm[l][2]*nz;
1247
1248         LbmFloat ret = 1.0;
1249         // forward direction, add mass (for filling cells):
1250         if(dirForw) {
1251                 if(scal<LBM_EPSILON) ret = 0.0;
1252                 else ret = scal;
1253         } else {
1254                 // backward for emptying
1255                 if(scal>-LBM_EPSILON) ret = 0.0;
1256                 else ret = scal * -1.0;
1257         }
1258         //errMsg("massd", PRINT_IJK<<" nv"<<nvel<<" : ret="<<ret ); //xit(1); //VECDEB
1259         return ret;
1260 }
1261
1262 void LbmFsgrSolver::addToNewInterList( int ni, int nj, int nk ) {
1263 #if FSGR_STRICT_DEBUG==10
1264         // dangerous, this can change the simulation...
1265   /*for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1266        iter != mListNewInter.end(); iter++ ) {
1267     if(ni!=iter->x) continue;
1268     if(nj!=iter->y) continue;
1269     if(nk!=iter->z) continue;
1270                 // all 3 values match... skip point
1271                 return;
1272         } */
1273 #endif // FSGR_STRICT_DEBUG==1
1274         // store point
1275         LbmPoint newinter; newinter.flag = 0;
1276         newinter.x = ni; newinter.y = nj; newinter.z = nk;
1277         mListNewInter.push_back(newinter);
1278 }
1279
1280 void LbmFsgrSolver::reinitFlags( int workSet ) { 
1281         // reinitCellFlags OLD mods:
1282         // add all to intel list?
1283         // check ffrac for new cells
1284         // new if cell inits (last loop)
1285         // vweights handling
1286
1287         const int debugFlagreinit = 0;
1288         
1289         // some things need to be read/modified on the other set
1290         int otherSet = (workSet^1);
1291         // fixed level on which to perform 
1292         int workLev = mMaxRefine;
1293
1294   /* modify interface cells from lists */
1295   /* mark filled interface cells as fluid, or emptied as empty */
1296         /* count neighbors and distribute excess mass to interface neighbor cells
1297    * problems arise when there are no interface neighbors anymore
1298          * then just distribute to any fluid neighbors...
1299          */
1300
1301         // for symmetry, first init all neighbor cells */
1302         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1303        iter != mListFull.end(); iter++ ) {
1304     int i=iter->x, j=iter->y, k=iter->z;
1305                 if(debugFlagreinit) errMsg("FULL", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
1306     FORDF1 {
1307                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1308                         //if((LBMDIM>2)&&( (ni<=0) || (nj<=0) || (nk<=0) || (ni>=mLevel[workLev].lSizex-1) || (nj>=mLevel[workLev].lSizey-1) || (nk>=mLevel[workLev].lSizez-1) )) {
1309                         if( (ni<=0) || (nj<=0) || 
1310                                   (ni>=mLevel[workLev].lSizex-1) ||
1311                                   (nj>=mLevel[workLev].lSizey-1) 
1312 #                                       if LBMDIM==3
1313                                   || (nk<=0) || (nk>=mLevel[workLev].lSizez-1) 
1314 #                                       endif // LBMDIM==1
1315                                  ) {
1316                                 continue; } // new bc, dont treat cells on boundary NEWBC
1317       if( RFLAG(workLev, ni,nj,nk, workSet) & CFEmpty ){
1318                                 
1319                                 // preinit speed, get from average surrounding cells
1320                                 // interpolate from non-workset to workset, sets are handled in function
1321
1322                                 // new and empty interface cell, dont change old flag here!
1323                                 addToNewInterList(ni,nj,nk);
1324
1325                                 LbmFloat avgrho = 0.0;
1326                                 LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
1327                                 interpolateCellValues(workLev,ni,nj,nk,workSet, avgrho,avgux,avguy,avguz);
1328
1329                                 // careful with l's...
1330                                 FORDF0M { 
1331                                         QCELL(workLev,ni,nj,nk, workSet, m) = this->getCollideEq( m,avgrho,  avgux, avguy, avguz ); 
1332                                         //QCELL(workLev,ni,nj,nk, workSet, l) = avgnbdf[l]; // CHECK FIXME test?
1333                                 }
1334                                 //errMsg("FNEW", PRINT_VEC(ni,nj,nk)<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<<avgrho<<" vel"<<PRINT_VEC(avgux,avguy,avguz) ); // DEBUG SYMM
1335                                 QCELL(workLev,ni,nj,nk, workSet, dMass) = 0.0; //?? new
1336                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) = 0.0; //?? new
1337                                 //RFLAG(workLev,ni,nj,nk,workSet) = (CellFlagType)(CFInter|CFNoInterpolSrc);
1338                                 changeFlag(workLev,ni,nj,nk,workSet, (CFInter|CFNoInterpolSrc));
1339                                 if(debugFlagreinit) errMsg("NEWE", PRINT_IJK<<" newif "<<PRINT_VEC(ni,nj,nk)<<" rho"<<avgrho<<" vel("<<avgux<<","<<avguy<<","<<avguz<<") " );
1340       }
1341                         /* prevent surrounding interface cells from getting removed as empty cells 
1342                          * (also cells that are not newly inited) */
1343       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1344                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete);
1345                                 changeFlag(workLev,ni,nj,nk, workSet, (RFLAG(workLev,ni,nj,nk, workSet) | CFNoDelete));
1346                                 // also add to list...
1347                                 addToNewInterList(ni,nj,nk);
1348                         } // NEW?
1349     }
1350
1351                 // NEW? no extra loop...
1352                 //RFLAG(workLev,i,j,k, workSet) = CFFluid;
1353                 changeFlag(workLev,i,j,k, workSet,CFFluid);
1354         }
1355
1356         /* remove empty interface cells that are not allowed to be removed anyway
1357          * this is important, otherwise the dreaded cell-type-flickering can occur! */
1358   for(int n=0; n<(int)mListEmpty.size(); n++) {
1359     int i=mListEmpty[n].x, j=mListEmpty[n].y, k=mListEmpty[n].z;
1360                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)) {
1361                         // treat as "new inter"
1362                         addToNewInterList(i,j,k);
1363                         // remove entry
1364                         if(debugFlagreinit) errMsg("EMPT REMOVED!!!", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) ); // DEBUG SYMM
1365                         if(n<(int)mListEmpty.size()-1) mListEmpty[n] = mListEmpty[mListEmpty.size()-1]; 
1366                         mListEmpty.pop_back();
1367                         n--; 
1368                 }
1369         } 
1370
1371
1372         /* problems arise when adjacent cells empty&fill ->
1373                  let fill cells+surrounding interface cells have the higher importance */
1374   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1375        iter != mListEmpty.end(); iter++ ) {
1376     int i=iter->x, j=iter->y, k=iter->z;
1377                 if((RFLAG(workLev,i,j,k, workSet)&(CFInter|CFNoDelete)) == (CFInter|CFNoDelete)){ errMsg("A"," ARGHARGRAG "); } // DEBUG
1378                 if(debugFlagreinit) errMsg("EMPT", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" rho"<< QCELL(workLev, i,j,k, workSet, 0) );
1379
1380                 /* set surrounding fluid cells to interface cells */
1381     FORDF1 {
1382                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1383       if( RFLAG(workLev,ni,nj,nk, workSet) & CFFluid){
1384                                 // init fluid->interface 
1385                                 //RFLAG(workLev,ni,nj,nk, workSet) = (CellFlagType)(CFInter); 
1386                                 changeFlag(workLev,ni,nj,nk, workSet, CFInter); 
1387                                 /* new mass = current density */
1388                                 LbmFloat nbrho = QCELL(workLev,ni,nj,nk, workSet, dC);
1389                 for(int rl=1; rl< this->cDfNum ; ++rl) { nbrho += QCELL(workLev,ni,nj,nk, workSet, rl); }
1390                                 QCELL(workLev,ni,nj,nk, workSet, dMass) =  nbrho; 
1391                                 QCELL(workLev,ni,nj,nk, workSet, dFfrac) =  1.0; 
1392
1393                                 // store point
1394                                 addToNewInterList(ni,nj,nk);
1395       }
1396       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter){
1397                                 // test, also add to list...
1398                                 addToNewInterList(ni,nj,nk);
1399                         } // NEW?
1400     }
1401
1402                 /* for symmetry, set our flag right now */
1403                 changeFlag(workLev,i,j,k, workSet, CFEmpty);
1404                 // mark cell not be changed mass... - not necessary, not in list anymore anyway!
1405         } // emptylist
1406
1407
1408         
1409         // precompute weights to get rid of order dependancies
1410         vector<lbmFloatSet> vWeights;
1411         vWeights.resize( mListFull.size() + mListEmpty.size() );
1412         int weightIndex = 0;
1413   int nbCount = 0;
1414         LbmFloat nbWeights[LBM_DFNUM];
1415         LbmFloat nbTotWeights = 0.0;
1416         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1417        iter != mListFull.end(); iter++ ) {
1418     int i=iter->x, j=iter->y, k=iter->z;
1419     nbCount = 0; nbTotWeights = 0.0;
1420     FORDF1 {
1421                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1422       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1423                                 nbCount++;
1424                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1425                                 else nbWeights[l] = getMassdWeight(1,i,j,k,workSet,l); // NEWSURFT
1426                                 nbTotWeights += nbWeights[l];
1427       } else {
1428                                 nbWeights[l] = -100.0; // DEBUG;
1429                         }
1430     }
1431                 if(nbCount>0) { 
1432                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1433         vWeights[weightIndex].val[0] = nbTotWeights;
1434         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1435         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1436                 } else { 
1437         vWeights[weightIndex].numNbs = 0.0;
1438                 }
1439                 weightIndex++;
1440         }
1441   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1442        iter != mListEmpty.end(); iter++ ) {
1443     int i=iter->x, j=iter->y, k=iter->z;
1444     nbCount = 0; nbTotWeights = 0.0;
1445     FORDF1 {
1446                         int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1447       if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1448                                 nbCount++;
1449                                 if(iter->flag&1) nbWeights[l] = 1.; // NEWSURFT
1450                                 else nbWeights[l] = getMassdWeight(0,i,j,k,workSet,l); // NEWSURFT
1451                                 nbTotWeights += nbWeights[l];
1452       } else {
1453                                 nbWeights[l] = -100.0; // DEBUG;
1454                         }
1455     }
1456                 if(nbCount>0) { 
1457                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeights);
1458         vWeights[weightIndex].val[0] = nbTotWeights;
1459         FORDF1 { vWeights[weightIndex].val[l] = nbWeights[l]; }
1460         vWeights[weightIndex].numNbs = (LbmFloat)nbCount;
1461                 } else { 
1462         vWeights[weightIndex].numNbs = 0.0;
1463                 }
1464                 weightIndex++;
1465         } 
1466         weightIndex = 0;
1467         
1468
1469         /* process full list entries, filled cells are done after this loop */
1470         for( vector<LbmPoint>::iterator iter=mListFull.begin();
1471        iter != mListFull.end(); iter++ ) {
1472     int i=iter->x, j=iter->y, k=iter->z;
1473
1474                 LbmFloat myrho = QCELL(workLev,i,j,k, workSet, dC);
1475     FORDF1 { myrho += QCELL(workLev,i,j,k, workSet, l); } // QCELL.rho
1476
1477     LbmFloat massChange = QCELL(workLev,i,j,k, workSet, dMass) - myrho;
1478                 if(vWeights[weightIndex].numNbs>0.0) {
1479                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1480                         //errMsg("FF  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1481                         FORDF1 {
1482                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1483         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1484                                         LbmFloat change = -1.0;
1485                                         if(nbTotWeightsp>0.0) {
1486                                                 //change = massChange * ( nbWeights[l]/nbTotWeightsp );
1487                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1488                                         } else {
1489                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1490                                         }
1491                                         QCELL(workLev,ni,nj,nk, workSet, dMass) += change;
1492                                 }
1493                         }
1494                         massChange = 0.0;
1495                 } else {
1496                         // Problem! no interface neighbors
1497                         this->mFixMass += massChange;
1498                         //TTT mNumProblems++;
1499                         //errMsg(" FULL PROBLEM ", PRINT_IJK<<" "<<this->mFixMass);
1500                 }
1501                 weightIndex++;
1502
1503     // already done? RFLAG(workLev,i,j,k, workSet) = CFFluid;
1504     QCELL(workLev,i,j,k, workSet, dMass) = myrho; // should be rho... but unused?
1505     QCELL(workLev,i,j,k, workSet, dFfrac) = 1.0; // should be rho... but unused?
1506   } // fulllist
1507
1508
1509         /* now, finally handle the empty cells - order is important, has to be after
1510          * full cell handling */
1511   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1512        iter != mListEmpty.end(); iter++ ) {
1513     int i=iter->x, j=iter->y, k=iter->z;
1514
1515     LbmFloat massChange = QCELL(workLev, i,j,k, workSet, dMass);
1516                 if(vWeights[weightIndex].numNbs>0.0) {
1517                         const LbmFloat nbTotWeightsp = vWeights[weightIndex].val[0];
1518                         //errMsg("EE  I", PRINT_IJK<<" "<<weightIndex<<" "<<nbTotWeightsp);
1519                         FORDF1 {
1520                                 int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
1521         if( RFLAG(workLev,ni,nj,nk, workSet) & CFInter) {
1522                                         LbmFloat change = -1.0;
1523                                         if(nbTotWeightsp>0.0) {
1524                                                 change = massChange * ( vWeights[weightIndex].val[l]/nbTotWeightsp );
1525                                         } else {
1526                                                 change = (LbmFloat)(massChange/vWeights[weightIndex].numNbs);
1527                                         }
1528                                         QCELL(workLev, ni,nj,nk, workSet, dMass) += change;
1529                                 }
1530                         }
1531                         massChange = 0.0;
1532                 } else {
1533                         // Problem! no interface neighbors
1534                         this->mFixMass += massChange;
1535                         //TTT mNumProblems++;
1536                         //errMsg(" EMPT PROBLEM ", PRINT_IJK<<" "<<this->mFixMass);
1537                 }
1538                 weightIndex++;
1539                 
1540                 // finally... make it empty 
1541     // already done? RFLAG(workLev,i,j,k, workSet) = CFEmpty;
1542     QCELL(workLev,i,j,k, workSet, dMass) = 0.0;
1543     QCELL(workLev,i,j,k, workSet, dFfrac) = 0.0;
1544         }
1545   for( vector<LbmPoint>::iterator iter=mListEmpty.begin();
1546        iter != mListEmpty.end(); iter++ ) {
1547     int i=iter->x, j=iter->y, k=iter->z;
1548     changeFlag(workLev,i,j,k, otherSet, CFEmpty);
1549         } 
1550
1551
1552         // check if some of the new interface cells can be removed again 
1553         // never happens !!! not necessary
1554         // calculate ffrac for new IF cells NEW
1555
1556         // how many are really new interface cells?
1557         int numNewIf = 0;
1558   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1559        iter != mListNewInter.end(); iter++ ) {
1560     int i=iter->x, j=iter->y, k=iter->z;
1561                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1562                         continue; 
1563                         // FIXME remove from list?
1564                 }
1565                 numNewIf++;
1566         }
1567
1568         // redistribute mass, reinit flags
1569         if(debugFlagreinit) errMsg("NEWIF", "total:"<<mListNewInter.size());
1570         float newIfFac = 1.0/(LbmFloat)numNewIf;
1571   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1572        iter != mListNewInter.end(); iter++ ) {
1573     int i=iter->x, j=iter->y, k=iter->z;
1574                 if((i<=0) || (j<=0) || 
1575                          (i>=mLevel[workLev].lSizex-1) ||
1576                          (j>=mLevel[workLev].lSizey-1) ||
1577                          ((LBMDIM==3) && ((k<=0) || (k>=mLevel[workLev].lSizez-1) ) )
1578                          ) {
1579                         continue; } // new bc, dont treat cells on boundary NEWBC
1580                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { 
1581                         //errMsg("???"," "<<PRINT_IJK);
1582                         continue; 
1583                 } // */
1584
1585     QCELL(workLev,i,j,k, workSet, dMass) += (this->mFixMass * newIfFac);
1586                 
1587                 int nbored = 0;
1588                 FORDF1 { nbored |= RFLAG_NB(workLev, i,j,k, workSet,l); }
1589                 if(!(nbored & CFBndNoslip)) { RFLAG(workLev,i,j,k, workSet) |= CFNoBndFluid; }
1590                 if(!(nbored & CFFluid))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbFluid; }
1591                 if(!(nbored & CFEmpty))     { RFLAG(workLev,i,j,k, workSet) |= CFNoNbEmpty; }
1592
1593                 if(!(RFLAG(workLev,i,j,k, otherSet)&CFInter)) {
1594                         RFLAG(workLev,i,j,k, workSet) = (CellFlagType)(RFLAG(workLev,i,j,k, workSet) | CFNoDelete);
1595                 }
1596                 if(debugFlagreinit) errMsg("NEWIF", PRINT_IJK<<" mss"<<QCELL(workLev, i,j,k, workSet, dMass) <<" f"<< convertCellFlagType2String(RFLAG(workLev,i,j,k, workSet))<<" wl"<<workLev );
1597         }
1598
1599         // reinit fill fraction
1600   for( vector<LbmPoint>::iterator iter=mListNewInter.begin();
1601        iter != mListNewInter.end(); iter++ ) {
1602     int i=iter->x, j=iter->y, k=iter->z;
1603                 if(!(RFLAG(workLev,i,j,k, workSet)&CFInter)) { continue; }
1604
1605                 LbmFloat nrho = 0.0;
1606                 FORDF0 { nrho += QCELL(workLev, i,j,k, workSet, l); }
1607     QCELL(workLev,i,j,k, workSet, dFfrac) = QCELL(workLev,i,j,k, workSet, dMass)/nrho;
1608     QCELL(workLev,i,j,k, workSet, dFlux) = FLUX_INIT;
1609         }
1610
1611         if(mListNewInter.size()>0){ 
1612                 //errMsg("FixMassDisted"," fm:"<<this->mFixMass<<" nif:"<<mListNewInter.size() );
1613                 this->mFixMass = 0.0; 
1614         }
1615
1616         // empty lists for next step
1617         mListFull.clear();
1618         mListEmpty.clear();
1619         mListNewInter.clear();
1620 } // reinitFlags
1621
1622
1623