- removed some unecessary files & code
[blender-staging.git] / intern / elbeem / intern / lbmfsgrsolver.h
1 /******************************************************************************
2  *
3  * El'Beem - the visual lattice boltzmann freesurface simulator
4  * All code distributed as part of El'Beem is covered by the version 2 of the 
5  * GNU General Public License. See the file COPYING for details.
6  * Copyright 2003-2005 Nils Thuerey
7  *
8  * Combined 2D/3D Lattice Boltzmann standard solver classes
9  *
10  *****************************************************************************/
11
12
13 #ifndef LBMFSGRSOLVER_H
14 #include "utilities.h"
15 #include "lbmdimensions.h"
16 #include "lbmfunctions.h"
17 #include "ntl_scene.h"
18 #include <stdio.h>
19
20 #if PARALLEL==1
21 #include <omp.h>
22 #endif // PARALLEL=1
23 #ifndef PARALLEL
24 #define PARALLEL 0
25 #endif // PARALLEL
26
27 // blender interface
28 #if ELBEEM_BLENDER==1
29 #include "SDL.h"
30 #include "SDL_thread.h"
31 #include "SDL_mutex.h"
32 extern "C" {
33         void simulateThreadIncreaseFrame(void);
34         extern SDL_mutex *globalBakeLock;
35         extern int        globalBakeState;
36         extern int        globalBakeFrame;
37 }
38 #endif // ELBEEM_BLENDER==1
39
40 #ifndef LBMMODEL_DEFINED
41 // force compiler error!
42 ERROR - define model first!
43 #endif // LBMMODEL_DEFINED
44
45
46 // general solver setting defines
47
48 //! debug coordinate accesses and the like? (much slower)
49 #define FSGR_STRICT_DEBUG 0
50
51 //! debug coordinate accesses and the like? (much slower)
52 #define FSGR_OMEGA_DEBUG 0
53
54 //! OPT3D quick LES on/off, only debug/benachmarking
55 #define USE_LES 1
56
57 //! order of interpolation (1/2)
58 #define INTORDER 2
59
60 //! order of interpolation (0=always current/1=interpolate/2=always other)
61 //#define TIMEINTORDER 0
62 // TODO remove interpol t param, also interTime
63
64 //! refinement border method (1 = small border / 2 = larger)
65 #define REFINEMENTBORDER 1
66
67 // use optimized 3D code?
68 #if LBMDIM==2
69 #define OPT3D false
70 #else
71 // determine with debugging...
72 #       if FSGR_STRICT_DEBUG==1
73 #               define OPT3D false
74 #       else // FSGR_STRICT_DEBUG==1
75 // usually switch optimizations for 3d on, when not debugging
76 #               define OPT3D true
77 // COMPRT
78 //#             define OPT3D false
79 #       endif // FSGR_STRICT_DEBUG==1
80 #endif
81
82 // enable/disable fine grid compression for finest level
83 #if LBMDIM==3
84 #define COMPRESSGRIDS 1
85 #else 
86 #define COMPRESSGRIDS 0
87 #endif 
88
89 //! threshold for level set fluid generation/isosurface
90 #define LS_FLUIDTHRESHOLD 0.5
91
92 //! invalid mass value for unused mass data
93 #define MASS_INVALID -1000.0
94
95 // empty/fill cells without fluid/empty NB's by inserting them into the full/empty lists?
96 #define FSGR_LISTTRICK          true
97 #define FSGR_LISTTTHRESHEMPTY   0.10
98 #define FSGR_LISTTTHRESHFULL    0.90
99 #define FSGR_MAGICNR            0.025
100 //0.04
101
102
103 // helper for comparing floats with epsilon
104 #define GFX_FLOATNEQ(x,y) ( ABS((x)-(y)) > (VECTOR_EPSILON) )
105 #define LBM_FLOATNEQ(x,y) ( ABS((x)-(y)) > (10.0*LBM_EPSILON) )
106
107
108 // macros for loops over all DFs
109 #define FORDF0 for(int l= 0; l< LBM_DFNUM; ++l)
110 #define FORDF1 for(int l= 1; l< LBM_DFNUM; ++l)
111 // and with different loop var to prevent shadowing
112 #define FORDF0M for(int m= 0; m< LBM_DFNUM; ++m)
113 #define FORDF1M for(int m= 1; m< LBM_DFNUM; ++m)
114
115 // aux. field indices (same for 2d)
116 #define dFfrac 19
117 #define dMass 20
118 #define dFlux 21
119 // max. no. of cell values for 3d
120 #define dTotalNum 22
121
122 // iso value defines
123 // border for marching cubes
124 #define ISOCORR 3
125
126
127 /*****************************************************************************/
128 /*! cell access classes */
129 template<typename D>
130 class UniformFsgrCellIdentifier : 
131         public CellIdentifierInterface 
132 {
133         public:
134                 //! which grid level?
135                 int level;
136                 //! location in grid
137                 int x,y,z;
138
139                 //! reset constructor
140                 UniformFsgrCellIdentifier() :
141                         x(0), y(0), z(0) { };
142
143                 // implement CellIdentifierInterface
144                 virtual string getAsString() {
145                         std::ostringstream ret;
146                         ret <<"{ i"<<x<<",j"<<y;
147                         if(D::cDimension>2) ret<<",k"<<z;
148                         ret <<" }";
149                         return ret.str();
150                 }
151
152                 virtual bool equal(CellIdentifierInterface* other) {
153                         //UniformFsgrCellIdentifier<D> *cid = dynamic_cast<UniformFsgrCellIdentifier<D> *>( other );
154                         UniformFsgrCellIdentifier<D> *cid = (UniformFsgrCellIdentifier<D> *)( other );
155                         if(!cid) return false;
156                         if( x==cid->x && y==cid->y && z==cid->z && level==cid->level ) return true;
157                         return false;
158                 }
159 };
160
161 //! information needed for each level in the simulation
162 class FsgrLevelData {
163 public:
164         int id; // level number
165
166         //! node size on this level (geometric, in world coordinates, not simulation units!) 
167         LbmFloat nodeSize;
168         //! node size on this level in simulation units 
169         LbmFloat simCellSize;
170         //! quadtree node relaxation parameter 
171         LbmFloat omega;
172         //! size this level was advanced to 
173         LbmFloat time;
174         //! size of a single lbm step in time units on this level 
175         LbmFloat stepsize;
176         //! step count
177         int lsteps;
178         //! gravity force for this level
179         LbmVec gravity;
180         //! level array 
181         LbmFloat *mprsCells[2];
182         CellFlagType *mprsFlags[2];
183
184         //! smago params and precalculated values
185         LbmFloat lcsmago;
186         LbmFloat lcsmago_sqr;
187         LbmFloat lcnu;
188
189         // LES statistics per level
190         double avgOmega;
191         double avgOmegaCnt;
192
193         //! current set of dist funcs 
194         int setCurr;
195         //! target/other set of dist funcs 
196         int setOther;
197
198         //! mass&volume for this level
199         LbmFloat lmass;
200         LbmFloat lvolume;
201         LbmFloat lcellfactor;
202
203         //! local storage of mSizes
204         int lSizex, lSizey, lSizez;
205         int lOffsx, lOffsy, lOffsz;
206
207 };
208
209
210
211 /*****************************************************************************/
212 /*! class for solving a LBM problem */
213 template<class D>
214 class LbmFsgrSolver : 
215         public /*? virtual */ D // this means, the solver is a lbmData object and implements the lbmInterface
216 {
217
218         public:
219                 //! Constructor 
220                 LbmFsgrSolver();
221                 //! Destructor 
222                 virtual ~LbmFsgrSolver();
223                 //! id string of solver
224                 virtual string getIdString() { return string("FsgrSolver[") + D::getIdString(); }
225
226                 //! initilize variables fom attribute list 
227                 virtual void parseAttrList();
228                 //! Initialize omegas and forces on all levels (for init/timestep change)
229                 void initLevelOmegas();
230                 //! finish the init with config file values (allocate arrays...) 
231                 virtual bool initialize( ntlTree* /*tree*/, vector<ntlGeometryObject*>* /*objects*/ );
232
233 #if LBM_USE_GUI==1
234                 //! show simulation info (implement SimulationObject pure virtual func)
235                 virtual void debugDisplay(fluidDispSettings *set);
236 #endif
237                 
238                 
239                 // implement CellIterator<UniformFsgrCellIdentifier> interface
240                 typedef UniformFsgrCellIdentifier<typename D::LbmCellContents> stdCellId;
241                 virtual CellIdentifierInterface* getFirstCell( );
242                 virtual void advanceCell( CellIdentifierInterface* );
243                 virtual bool noEndCell( CellIdentifierInterface* );
244                 virtual void deleteCellIterator( CellIdentifierInterface** );
245                 virtual CellIdentifierInterface* getCellAt( ntlVec3Gfx pos );
246                 virtual int        getCellSet      ( CellIdentifierInterface* );
247                 virtual ntlVec3Gfx getCellOrigin   ( CellIdentifierInterface* );
248                 virtual ntlVec3Gfx getCellSize     ( CellIdentifierInterface* );
249                 virtual int        getCellLevel    ( CellIdentifierInterface* );
250                 virtual LbmFloat   getCellDensity  ( CellIdentifierInterface* ,int set);
251                 virtual LbmVec     getCellVelocity ( CellIdentifierInterface* ,int set);
252                 virtual LbmFloat   getCellDf       ( CellIdentifierInterface* ,int set, int dir);
253                 virtual LbmFloat   getCellMass     ( CellIdentifierInterface* ,int set);
254                 virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set);
255                 virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set);
256                 virtual LbmFloat   getEquilDf      ( int );
257                 virtual int        getDfNum        ( );
258                 // convert pointers
259                 stdCellId* convertBaseCidToStdCid( CellIdentifierInterface* basecid);
260
261                 //! perform geometry init (if switched on) 
262                 bool initGeometryFlags();
263                 //! init part for all freesurface testcases 
264                 void initFreeSurfaces();
265                 //! init density gradient if enabled
266                 void initStandingFluidGradient();
267
268                 /*! init a given cell with flag, density, mass and equilibrium dist. funcs */
269                 inline void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass);
270                 inline void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
271                 inline void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
272
273                 /*! perform a single LBM step */
274                 virtual void step() { stepMain(); }
275                 void stepMain();
276                 void fineAdvance();
277                 void coarseAdvance(int lev);
278                 void coarseCalculateFluxareas(int lev);
279                 // coarsen a given level (returns true if sth. was changed)
280                 bool performRefinement(int lev);
281                 bool performCoarsening(int lev);
282                 //void oarseInterpolateToFineSpaceTime(int lev,LbmFloat t);
283                 void interpolateFineFromCoarse(int lev,LbmFloat t);
284                 void coarseRestrictFromFine(int lev);
285
286                 /*! init particle positions */
287                 virtual int initParticles(ParticleTracer *partt);
288                 /*! move all particles */
289                 virtual void advanceParticles(ParticleTracer *partt );
290
291
292                 /*! debug object display (used e.g. for preview surface) */
293                 virtual vector<ntlGeometryObject*> getDebugObjects();
294
295                 //! access the fillfrac field (for viz)
296                 inline float getFillFrac(int i, int j, int k);
297
298                 //! retrieve the fillfrac field ready to display
299                 void getIsofieldWeighted(float *iso);
300                 void getIsofield(float *iso){ return getIsofieldWeighted(iso); }
301                 //! for raytracing, preprocess
302                 void prepareVisualization( void );
303
304                 // rt interface
305                 void addDrop(bool active, float mx, float my);
306                 void initDrop(float mx, float my);
307                 void printCellStats();
308                 int checkGfxEndTime(); // {return 9;};
309                 //! get gfx geo setup id
310                 int getGfxGeoSetup() { return mGfxGeoSetup; }
311
312                 /*! type for cells */
313                 typedef typename D::LbmCell LbmCell;
314                 
315         protected:
316
317                 //! internal quick print function (for debugging) 
318                 void printLbmCell(int level, int i, int j, int k,int set);
319                 // debugging use CellIterator interface to mark cell
320                 void debugMarkCellCall(int level, int vi,int vj,int vk);
321
322                 void mainLoop(int lev);
323                 void adaptTimestep();
324                 //! init mObjectSpeeds for current parametrization
325                 void recalculateObjectSpeeds();
326                 //! flag reinit step - always works on finest grid!
327                 void reinitFlags( int workSet );
328                 //! mass dist weights
329                 LbmFloat getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l);
330                 //! add point to mListNewInter list
331                 inline void addToNewInterList( int ni, int nj, int nk );        
332                 //! cell is interpolated from coarse level (inited into set, source sets are determined by t)
333                 inline void interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet,bool markNbs);
334
335                 //! minimal and maximal z-coords (for 2D/3D loops)
336                 int getForZMinBnd() { return 0; }
337                 int getForZMin1()   { 
338                         if(D::cDimension==2) return 0;
339                         return 1; 
340                 }
341
342                 int getForZMaxBnd(int lev) { 
343                         if(D::cDimension==2) return 1;
344                         return mLevel[lev].lSizez -0;
345                 }
346                 int getForZMax1(int lev)   { 
347                         if(D::cDimension==2) return 1;
348                         return mLevel[lev].lSizez -1;
349                 }
350
351
352                 // member vars
353
354                 //! mass calculated during streaming step
355                 LbmFloat mCurrentMass;
356                 LbmFloat mCurrentVolume;
357                 LbmFloat mInitialMass;
358
359                 //! count problematic cases, that occured so far...
360                 int mNumProblems;
361
362                 // average mlsups, count how many so far...
363                 double mAvgMLSUPS;
364                 double mAvgMLSUPSCnt;
365
366                 //! Mcubes object for surface reconstruction 
367                 IsoSurface *mpPreviewSurface;
368                 float mSmoothSurface;
369                 float mSmoothNormals;
370                 
371                 //! use time adaptivity? 
372                 bool mTimeAdap;
373
374                 //! output surface preview? if >0 yes, and use as reduzed size 
375                 int mOutputSurfacePreview;
376                 LbmFloat mPreviewFactor;
377                 //! fluid vol height
378                 LbmFloat mFVHeight;
379                 LbmFloat mFVArea;
380                 bool mUpdateFVHeight;
381
382                 //! require some geo setup from the viz?
383                 int mGfxGeoSetup;
384                 //! force quit for gfx
385                 LbmFloat mGfxEndTime;
386                 //! smoother surface initialization?
387                 int mInitSurfaceSmoothing;
388
389                 int mTimestepReduceLock;
390                 int mTimeSwitchCounts;
391                 //! total simulation time so far 
392                 LbmFloat mSimulationTime;
393                 //! smallest and largest step size so far 
394                 LbmFloat mMinStepTime, mMaxStepTime;
395                 //! track max. velocity
396                 LbmFloat mMxvx, mMxvy, mMxvz, mMaxVlen;
397
398                 //! list of the cells to empty at the end of the step 
399                 vector<LbmPoint> mListEmpty;
400                 //! list of the cells to make fluid at the end of the step 
401                 vector<LbmPoint> mListFull;
402                 //! list of new interface cells to init
403         vector<LbmPoint> mListNewInter;
404                 //! class for handling redist weights in reinit flag function
405                 class lbmFloatSet {
406                         public:
407                                 LbmFloat val[dTotalNum];
408                                 LbmFloat numNbs;
409                 };
410                 //! normalized vectors for all neighboring cell directions (for e.g. massdweight calc)
411                 LbmVec mDvecNrm[27];
412                 
413                 
414                 //! debugging
415                 bool checkSymmetry(string idstring);
416                 //! symmetric init?
417                 //bool mStartSymm;
418                 //! kepp track of max/min no. of filled cells
419                 int mMaxNoCells, mMinNoCells;
420                 long long int mAvgNumUsedCells;
421
422                 //! for interactive - how to drop drops?
423                 int mDropMode;
424                 LbmFloat mDropSize;
425                 LbmVec mDropSpeed;
426                 //! dropping variables
427                 bool mDropping;
428                 LbmFloat mDropX, mDropY;
429                 LbmFloat mDropHeight;
430                 //! precalculated objects speeds for current parametrization
431                 vector<LbmVec> mObjectSpeeds;
432
433                 //! get isofield weights
434                 int mIsoWeightMethod;
435                 float mIsoWeight[27];
436
437                 // grid coarsening vars
438                 
439                 /*! vector for the data for each level */
440 #               define MAX_LEV 5
441                 FsgrLevelData mLevel[MAX_LEV];
442
443                 /*! minimal and maximal refinement levels */
444                 int mMaxRefine;
445
446                 /*! df scale factors for level up/down */
447                 LbmFloat mDfScaleUp, mDfScaleDown;
448
449                 /*! precomputed cell area values */
450                 LbmFloat mFsgrCellArea[27];
451
452                 /*! LES C_smago paramter for finest grid */
453                 float mInitialCsmago;
454                 /*! LES stats for non OPT3D */
455                 LbmFloat mDebugOmegaRet;
456
457                 //! fluid stats
458                 int mNumInterdCells;
459                 int mNumInvIfCells;
460                 int mNumInvIfTotal;
461                 int mNumFsgrChanges;
462
463                 //! debug function to disable standing f init
464                 int mDisableStandingFluidInit;
465                 //! debug function to force tadap syncing
466                 int mForceTadapRefine;
467
468
469                 // strict debug interface
470 #               if FSGR_STRICT_DEBUG==1
471                 int debLBMGI(int level, int ii,int ij,int ik, int is);
472                 CellFlagType& debRFLAG(int level, int xx,int yy,int zz,int set);
473                 CellFlagType& debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir);
474                 CellFlagType& debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir);
475                 int debLBMQI(int level, int ii,int ij,int ik, int is, int l);
476                 LbmFloat& debQCELL(int level, int xx,int yy,int zz,int set,int l);
477                 LbmFloat& debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l);
478                 LbmFloat& debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l);
479                 LbmFloat* debRACPNT(int level,  int ii,int ij,int ik, int is );
480                 LbmFloat& debRAC(LbmFloat* s,int l);
481 #               endif // FSGR_STRICT_DEBUG==1
482 };
483
484
485
486 /*****************************************************************************/
487 // relaxation_macros
488
489
490
491 // cell mark debugging
492 #if FSGR_STRICT_DEBUG==10
493 #define debugMarkCell(lev,x,y,z) \
494         errMsg("debugMarkCell",D::mName<<" step: "<<D::mStepCnt<<" lev:"<<(lev)<<" marking "<<PRINT_VEC((x),(y),(z))<<" line "<< __LINE__ ); \
495         debugMarkCellCall((lev),(x),(y),(z));
496 #else // FSGR_STRICT_DEBUG==1
497 #define debugMarkCell(lev,x,y,z) \
498         debugMarkCellCall((lev),(x),(y),(z));
499 #endif // FSGR_STRICT_DEBUG==1
500
501
502 // flag array defines -----------------------------------------------------------------------------------------------
503
504 // lbm testsolver get index define
505 #define _LBMGI(level, ii,ij,ik, is) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
506
507 //! flag array acces macro
508 #define _RFLAG(level,xx,yy,zz,set) mLevel[level].mprsFlags[set][ LBMGI((level),(xx),(yy),(zz),(set)) ]
509 #define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set) ]
510 #define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set) ]
511
512 // array data layouts
513 // standard array layout  -----------------------------------------------------------------------------------------------
514 #define ALSTRING "Standard Array Layout"
515 //#define _LBMQI(level, ii,ij,ik, is, lunused) ( ((is)*mLevel[level].lOffsz) + (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
516 #define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
517 #define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)])
518 #define _QCELL_NB(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[dir],(yy)+D::dfVecY[dir],(zz)+D::dfVecZ[dir],set, l)*dTotalNum +(l)])
519 #define _QCELL_NBINV(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+D::dfVecX[D::dfInv[dir]],(yy)+D::dfVecY[D::dfInv[dir]],(zz)+D::dfVecZ[D::dfInv[dir]],set, l)*dTotalNum +(l)])
520
521 #define QCELLSTEP dTotalNum
522 #define _RACPNT(level, ii,ij,ik, is )  &QCELL(level,ii,ij,ik,is,0)
523 #define _RAC(s,l) (s)[(l)]
524
525 // standard arrays
526 #define CSRC_C    RAC(ccel                                , dC )
527 #define CSRC_E    RAC(ccel + (-1)             *(dTotalNum), dE )
528 #define CSRC_W    RAC(ccel + (+1)             *(dTotalNum), dW )
529 #define CSRC_N    RAC(ccel + (-mLevel[lev].lOffsx)        *(dTotalNum), dN )
530 #define CSRC_S    RAC(ccel + (+mLevel[lev].lOffsx)        *(dTotalNum), dS )
531 #define CSRC_NE   RAC(ccel + (-mLevel[lev].lOffsx-1)      *(dTotalNum), dNE)
532 #define CSRC_NW   RAC(ccel + (-mLevel[lev].lOffsx+1)      *(dTotalNum), dNW)
533 #define CSRC_SE   RAC(ccel + (+mLevel[lev].lOffsx-1)      *(dTotalNum), dSE)
534 #define CSRC_SW   RAC(ccel + (+mLevel[lev].lOffsx+1)      *(dTotalNum), dSW)
535 #define CSRC_T    RAC(ccel + (-mLevel[lev].lOffsy)        *(dTotalNum), dT )
536 #define CSRC_B    RAC(ccel + (+mLevel[lev].lOffsy)        *(dTotalNum), dB )
537 #define CSRC_ET   RAC(ccel + (-mLevel[lev].lOffsy-1)      *(dTotalNum), dET)
538 #define CSRC_EB   RAC(ccel + (+mLevel[lev].lOffsy-1)      *(dTotalNum), dEB)
539 #define CSRC_WT   RAC(ccel + (-mLevel[lev].lOffsy+1)      *(dTotalNum), dWT)
540 #define CSRC_WB   RAC(ccel + (+mLevel[lev].lOffsy+1)      *(dTotalNum), dWB)
541 #define CSRC_NT   RAC(ccel + (-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNT)
542 #define CSRC_NB   RAC(ccel + (+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *(dTotalNum), dNB)
543 #define CSRC_ST   RAC(ccel + (-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dST)
544 #define CSRC_SB   RAC(ccel + (+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *(dTotalNum), dSB)
545
546 #define XSRC_C(x)    RAC(ccel + (x)                 *dTotalNum, dC )
547 #define XSRC_E(x)    RAC(ccel + ((x)-1)             *dTotalNum, dE )
548 #define XSRC_W(x)    RAC(ccel + ((x)+1)             *dTotalNum, dW )
549 #define XSRC_N(x)    RAC(ccel + ((x)-mLevel[lev].lOffsx)        *dTotalNum, dN )
550 #define XSRC_S(x)    RAC(ccel + ((x)+mLevel[lev].lOffsx)        *dTotalNum, dS )
551 #define XSRC_NE(x)   RAC(ccel + ((x)-mLevel[lev].lOffsx-1)      *dTotalNum, dNE)
552 #define XSRC_NW(x)   RAC(ccel + ((x)-mLevel[lev].lOffsx+1)      *dTotalNum, dNW)
553 #define XSRC_SE(x)   RAC(ccel + ((x)+mLevel[lev].lOffsx-1)      *dTotalNum, dSE)
554 #define XSRC_SW(x)   RAC(ccel + ((x)+mLevel[lev].lOffsx+1)      *dTotalNum, dSW)
555 #define XSRC_T(x)    RAC(ccel + ((x)-mLevel[lev].lOffsy)        *dTotalNum, dT )
556 #define XSRC_B(x)    RAC(ccel + ((x)+mLevel[lev].lOffsy)        *dTotalNum, dB )
557 #define XSRC_ET(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy-1)      *dTotalNum, dET)
558 #define XSRC_EB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy-1)      *dTotalNum, dEB)
559 #define XSRC_WT(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy+1)      *dTotalNum, dWT)
560 #define XSRC_WB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy+1)      *dTotalNum, dWB)
561 #define XSRC_NT(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNT)
562 #define XSRC_NB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy-mLevel[lev].lOffsx) *dTotalNum, dNB)
563 #define XSRC_ST(x)   RAC(ccel + ((x)-mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dST)
564 #define XSRC_SB(x)   RAC(ccel + ((x)+mLevel[lev].lOffsy+mLevel[lev].lOffsx) *dTotalNum, dSB)
565
566
567
568 #if FSGR_STRICT_DEBUG==1
569
570 #define LBMGI(level,ii,ij,ik, is)                 debLBMGI(level,ii,ij,ik, is)         
571 #define RFLAG(level,xx,yy,zz,set)                 debRFLAG(level,xx,yy,zz,set)            
572 #define RFLAG_NB(level,xx,yy,zz,set, dir)         debRFLAG_NB(level,xx,yy,zz,set, dir)    
573 #define RFLAG_NBINV(level,xx,yy,zz,set, dir)      debRFLAG_NBINV(level,xx,yy,zz,set, dir) 
574
575 #define LBMQI(level,ii,ij,ik, is, l)              debLBMQI(level,ii,ij,ik, is, l)         
576 #define QCELL(level,xx,yy,zz,set,l)               debQCELL(level,xx,yy,zz,set,l)         
577 #define QCELL_NB(level,xx,yy,zz,set, dir,l)       debQCELL_NB(level,xx,yy,zz,set, dir,l)
578 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    debQCELL_NBINV(level,xx,yy,zz,set, dir,l)
579 #define RACPNT(level, ii,ij,ik, is )              debRACPNT(level, ii,ij,ik, is )          
580 #define RAC(s,l)                                                debRAC(s,l)                  
581
582 #else // FSGR_STRICT_DEBUG==1
583
584 #define LBMGI(level,ii,ij,ik, is)                 _LBMGI(level,ii,ij,ik, is)         
585 #define RFLAG(level,xx,yy,zz,set)                 _RFLAG(level,xx,yy,zz,set)            
586 #define RFLAG_NB(level,xx,yy,zz,set, dir)         _RFLAG_NB(level,xx,yy,zz,set, dir)    
587 #define RFLAG_NBINV(level,xx,yy,zz,set, dir)      _RFLAG_NBINV(level,xx,yy,zz,set, dir) 
588
589 #define LBMQI(level,ii,ij,ik, is, l)              _LBMQI(level,ii,ij,ik, is, l)         
590 #define QCELL(level,xx,yy,zz,set,l)               _QCELL(level,xx,yy,zz,set,l)         
591 #define QCELL_NB(level,xx,yy,zz,set, dir,l)       _QCELL_NB(level,xx,yy,zz,set, dir, l)
592 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    _QCELL_NBINV(level,xx,yy,zz,set, dir,l)
593 #define RACPNT(level, ii,ij,ik, is )              _RACPNT(level, ii,ij,ik, is )          
594 #define RAC(s,l)                                  _RAC(s,l)                  
595
596 #endif // FSGR_STRICT_DEBUG==1
597
598 // general defines -----------------------------------------------------------------------------------------------
599
600 #define TESTFLAG(flag, compflag) ((flag & compflag)==compflag)
601
602 #if LBMDIM==2
603 #define dC 0
604 #define dN 1
605 #define dS 2
606 #define dE 3
607 #define dW 4
608 #define dNE 5
609 #define dNW 6
610 #define dSE 7
611 #define dSW 8
612 #define LBM_DFNUM 9
613 #else
614 // direction indices
615 #define dC 0
616 #define dN 1
617 #define dS 2
618 #define dE 3
619 #define dW 4
620 #define dT 5
621 #define dB 6
622 #define dNE 7
623 #define dNW 8
624 #define dSE 9
625 #define dSW 10
626 #define dNT 11
627 #define dNB 12
628 #define dST 13
629 #define dSB 14
630 #define dET 15
631 #define dEB 16
632 #define dWT 17
633 #define dWB 18
634 #define LBM_DFNUM 19
635 #endif
636 //? #define dWB 18
637
638 // default init for dFlux values
639 #define FLUX_INIT 0.5f * (float)(D::cDfNum)
640
641 // only for non DF dir handling!
642 #define dNET 19
643 #define dNWT 20
644 #define dSET 21
645 #define dSWT 22
646 #define dNEB 23
647 #define dNWB 24
648 #define dSEB 25
649 #define dSWB 26
650
651 //! fill value for boundary cells
652 #define BND_FILL 0.0
653
654 #define DFL1 (1.0/ 3.0)
655 #define DFL2 (1.0/18.0)
656 #define DFL3 (1.0/36.0)
657
658 #define OMEGA(l) mLevel[(l)].omega
659
660 #define EQC (  DFL1*(rho - usqr))
661 #define EQN (  DFL2*(rho + uy*(4.5*uy + 3.0) - usqr))
662 #define EQS (  DFL2*(rho + uy*(4.5*uy - 3.0) - usqr))
663 #define EQE (  DFL2*(rho + ux*(4.5*ux + 3.0) - usqr))
664 #define EQW (  DFL2*(rho + ux*(4.5*ux - 3.0) - usqr))
665 #define EQT (  DFL2*(rho + uz*(4.5*uz + 3.0) - usqr))
666 #define EQB (  DFL2*(rho + uz*(4.5*uz - 3.0) - usqr))
667                     
668 #define EQNE ( DFL3*(rho + (+ux+uy)*(4.5*(+ux+uy) + 3.0) - usqr))
669 #define EQNW ( DFL3*(rho + (-ux+uy)*(4.5*(-ux+uy) + 3.0) - usqr))
670 #define EQSE ( DFL3*(rho + (+ux-uy)*(4.5*(+ux-uy) + 3.0) - usqr))
671 #define EQSW ( DFL3*(rho + (-ux-uy)*(4.5*(-ux-uy) + 3.0) - usqr))
672 #define EQNT ( DFL3*(rho + (+uy+uz)*(4.5*(+uy+uz) + 3.0) - usqr))
673 #define EQNB ( DFL3*(rho + (+uy-uz)*(4.5*(+uy-uz) + 3.0) - usqr))
674 #define EQST ( DFL3*(rho + (-uy+uz)*(4.5*(-uy+uz) + 3.0) - usqr))
675 #define EQSB ( DFL3*(rho + (-uy-uz)*(4.5*(-uy-uz) + 3.0) - usqr))
676 #define EQET ( DFL3*(rho + (+ux+uz)*(4.5*(+ux+uz) + 3.0) - usqr))
677 #define EQEB ( DFL3*(rho + (+ux-uz)*(4.5*(+ux-uz) + 3.0) - usqr))
678 #define EQWT ( DFL3*(rho + (-ux+uz)*(4.5*(-ux+uz) + 3.0) - usqr))
679 #define EQWB ( DFL3*(rho + (-ux-uz)*(4.5*(-ux-uz) + 3.0) - usqr))
680
681
682 // this is a bit ugly, but necessary for the CSRC_ access...
683 #define MSRC_C    m[dC ]
684 #define MSRC_N    m[dN ]
685 #define MSRC_S    m[dS ]
686 #define MSRC_E    m[dE ]
687 #define MSRC_W    m[dW ]
688 #define MSRC_T    m[dT ]
689 #define MSRC_B    m[dB ]
690 #define MSRC_NE   m[dNE]
691 #define MSRC_NW   m[dNW]
692 #define MSRC_SE   m[dSE]
693 #define MSRC_SW   m[dSW]
694 #define MSRC_NT   m[dNT]
695 #define MSRC_NB   m[dNB]
696 #define MSRC_ST   m[dST]
697 #define MSRC_SB   m[dSB]
698 #define MSRC_ET   m[dET]
699 #define MSRC_EB   m[dEB]
700 #define MSRC_WT   m[dWT]
701 #define MSRC_WB   m[dWB]
702
703 // this is a bit ugly, but necessary for the ccel local access...
704 #define CCEL_C    RAC(ccel, dC )
705 #define CCEL_N    RAC(ccel, dN )
706 #define CCEL_S    RAC(ccel, dS )
707 #define CCEL_E    RAC(ccel, dE )
708 #define CCEL_W    RAC(ccel, dW )
709 #define CCEL_T    RAC(ccel, dT )
710 #define CCEL_B    RAC(ccel, dB )
711 #define CCEL_NE   RAC(ccel, dNE)
712 #define CCEL_NW   RAC(ccel, dNW)
713 #define CCEL_SE   RAC(ccel, dSE)
714 #define CCEL_SW   RAC(ccel, dSW)
715 #define CCEL_NT   RAC(ccel, dNT)
716 #define CCEL_NB   RAC(ccel, dNB)
717 #define CCEL_ST   RAC(ccel, dST)
718 #define CCEL_SB   RAC(ccel, dSB)
719 #define CCEL_ET   RAC(ccel, dET)
720 #define CCEL_EB   RAC(ccel, dEB)
721 #define CCEL_WT   RAC(ccel, dWT)
722 #define CCEL_WB   RAC(ccel, dWB)
723 // for coarse to fine interpol access
724 #define CCELG_C(f)    (RAC(ccel, dC )*mGaussw[(f)])
725 #define CCELG_N(f)    (RAC(ccel, dN )*mGaussw[(f)])
726 #define CCELG_S(f)    (RAC(ccel, dS )*mGaussw[(f)])
727 #define CCELG_E(f)    (RAC(ccel, dE )*mGaussw[(f)])
728 #define CCELG_W(f)    (RAC(ccel, dW )*mGaussw[(f)])
729 #define CCELG_T(f)    (RAC(ccel, dT )*mGaussw[(f)])
730 #define CCELG_B(f)    (RAC(ccel, dB )*mGaussw[(f)])
731 #define CCELG_NE(f)   (RAC(ccel, dNE)*mGaussw[(f)])
732 #define CCELG_NW(f)   (RAC(ccel, dNW)*mGaussw[(f)])
733 #define CCELG_SE(f)   (RAC(ccel, dSE)*mGaussw[(f)])
734 #define CCELG_SW(f)   (RAC(ccel, dSW)*mGaussw[(f)])
735 #define CCELG_NT(f)   (RAC(ccel, dNT)*mGaussw[(f)])
736 #define CCELG_NB(f)   (RAC(ccel, dNB)*mGaussw[(f)])
737 #define CCELG_ST(f)   (RAC(ccel, dST)*mGaussw[(f)])
738 #define CCELG_SB(f)   (RAC(ccel, dSB)*mGaussw[(f)])
739 #define CCELG_ET(f)   (RAC(ccel, dET)*mGaussw[(f)])
740 #define CCELG_EB(f)   (RAC(ccel, dEB)*mGaussw[(f)])
741 #define CCELG_WT(f)   (RAC(ccel, dWT)*mGaussw[(f)])
742 #define CCELG_WB(f)   (RAC(ccel, dWB)*mGaussw[(f)])
743
744
745 #if PARALLEL==1
746 #define CSMOMEGA_STATS(dlev, domega) 
747 #else // PARALLEL==1
748 #if FSGR_OMEGA_DEBUG==1
749 #define CSMOMEGA_STATS(dlev, domega) \
750         mLevel[dlev].avgOmega += domega; mLevel[dlev].avgOmegaCnt+=1.0; 
751 #else // FSGR_OMEGA_DEBUG==1
752 #define CSMOMEGA_STATS(dlev, domega) 
753 #endif // FSGR_OMEGA_DEBUG==1
754 #endif // PARALLEL==1
755
756
757 // used for main loops and grav init
758 // source set
759 #define SRCS(l) mLevel[(l)].setCurr
760 // target set
761 #define TSET(l) mLevel[(l)].setOther
762
763
764 // complete default stream&collide, 2d/3d
765 /* read distribution funtions of adjacent cells = sweep step */ 
766 #if OPT3D==false 
767
768 #if FSGR_STRICT_DEBUG==1
769 #define MARKCELLCHECK \
770         debugMarkCell(lev,i,j,k); D::mPanic=1;
771 #define STREAMCHECK(ni,nj,nk,nl) \
772         if((m[l] < -1.0) || (m[l]>1.0)) {\
773                 errMsg("STREAMCHECK","Invalid streamed DF l"<<l<<" value:"<<m[l]<<" at "<<PRINT_IJK<<" from "<<PRINT_VEC(ni,nj,nk)<<" nl"<<(nl)<<\
774                                 " nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther)  ); \
775                 MARKCELLCHECK; \
776         }
777 #define COLLCHECK \
778         if( (rho>2.0) || (rho<-1.0) || (ABS(ux)>1.0) || (ABS(uy)>1.0) |(ABS(uz)>1.0) ) {\
779                 errMsg("COLLCHECK","Invalid collision values r:"<<rho<<" u:"PRINT_VEC(ux,uy,uz)<<" at? "<<PRINT_IJK ); \
780                 MARKCELLCHECK; \
781         }
782 #else
783 #define STREAMCHECK(ni,nj,nk,nl) 
784 #define COLLCHECK
785 #endif
786
787 // careful ux,uy,uz need to be inited before!
788
789 #define DEFAULT_STREAM \
790                 m[dC] = RAC(ccel,dC); \
791                 FORDF1 { \
792                         if(NBFLAG( D::dfInv[l] )&CFBnd) { \
793                                 m[l] = RAC(ccel, D::dfInv[l] ); \
794                                 STREAMCHECK(i,j,k, D::dfInv[l]); \
795                         } else { \
796                                 m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
797                                 STREAMCHECK(i+D::dfVecX[D::dfInv[l]], j+D::dfVecY[D::dfInv[l]],k+D::dfVecZ[D::dfInv[l]], l); \
798                         } \
799                 }   
800
801 // careful ux,uy,uz need to be inited before!
802 #define DEFAULT_COLLIDE \
803                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago, &mDebugOmegaRet ); \
804                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
805                         FORDF0 { RAC(tcel,l) = m[l]; }   \
806                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
807                         COLLCHECK;
808 #define OPTIMIZED_STREAMCOLLIDE \
809                         m[0] = RAC(ccel,0); \
810                         FORDF1 { /* df0 is set later on... */ \
811                                 /* FIXME CHECK INV ? */\
812                                 if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); D::mPanic=1;  \
813                                 } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
814                                 STREAMCHECK(i+D::dfVecX[D::dfInv[l]], j+D::dfVecY[D::dfInv[l]],k+D::dfVecZ[D::dfInv[l]], l); \
815                         }   \
816                         rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
817                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
818                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago , &mDebugOmegaRet  ); \
819                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
820                         FORDF0 { RAC(tcel,l) = m[l]; } \
821                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
822                         COLLCHECK;
823
824 #else  // 3D, opt OPT3D==true
825
826 #define DEFAULT_STREAM \
827                 m[dC] = RAC(ccel,dC); \
828                 /* explicit streaming */ \
829                 if((!nbored & CFBnd)) { \
830                         /* no boundary near?, no real speed diff.? */ \
831                         m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
832                         m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
833                         m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
834                         m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
835                         m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
836                         m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
837                 } else { \
838                         /* explicit streaming */ \
839                         if(NBFLAG(dS )&CFBnd) { m[dN ] = RAC(ccel,dS ); } else { m[dN ] = CSRC_N ; } \
840                         if(NBFLAG(dN )&CFBnd) { m[dS ] = RAC(ccel,dN ); } else { m[dS ] = CSRC_S ; } \
841                         if(NBFLAG(dW )&CFBnd) { m[dE ] = RAC(ccel,dW ); } else { m[dE ] = CSRC_E ; } \
842                         if(NBFLAG(dE )&CFBnd) { m[dW ] = RAC(ccel,dE ); } else { m[dW ] = CSRC_W ; } \
843                         if(NBFLAG(dB )&CFBnd) { m[dT ] = RAC(ccel,dB ); } else { m[dT ] = CSRC_T ; } \
844                         if(NBFLAG(dT )&CFBnd) { m[dB ] = RAC(ccel,dT ); } else { m[dB ] = CSRC_B ; } \
845                         \
846                         if(NBFLAG(dSW)&CFBnd) { m[dNE] = RAC(ccel,dSW); } else { m[dNE] = CSRC_NE; } \
847                         if(NBFLAG(dSE)&CFBnd) { m[dNW] = RAC(ccel,dSE); } else { m[dNW] = CSRC_NW; } \
848                         if(NBFLAG(dNW)&CFBnd) { m[dSE] = RAC(ccel,dNW); } else { m[dSE] = CSRC_SE; } \
849                         if(NBFLAG(dNE)&CFBnd) { m[dSW] = RAC(ccel,dNE); } else { m[dSW] = CSRC_SW; } \
850                         if(NBFLAG(dSB)&CFBnd) { m[dNT] = RAC(ccel,dSB); } else { m[dNT] = CSRC_NT; } \
851                         if(NBFLAG(dST)&CFBnd) { m[dNB] = RAC(ccel,dST); } else { m[dNB] = CSRC_NB; } \
852                         if(NBFLAG(dNB)&CFBnd) { m[dST] = RAC(ccel,dNB); } else { m[dST] = CSRC_ST; } \
853                         if(NBFLAG(dNT)&CFBnd) { m[dSB] = RAC(ccel,dNT); } else { m[dSB] = CSRC_SB; } \
854                         if(NBFLAG(dWB)&CFBnd) { m[dET] = RAC(ccel,dWB); } else { m[dET] = CSRC_ET; } \
855                         if(NBFLAG(dWT)&CFBnd) { m[dEB] = RAC(ccel,dWT); } else { m[dEB] = CSRC_EB; } \
856                         if(NBFLAG(dEB)&CFBnd) { m[dWT] = RAC(ccel,dEB); } else { m[dWT] = CSRC_WT; } \
857                         if(NBFLAG(dET)&CFBnd) { m[dWB] = RAC(ccel,dET); } else { m[dWB] = CSRC_WB; } \
858                 } 
859
860
861
862 #define COLL_CALCULATE_DFEQ(dstarray) \
863                         dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \
864                         dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \
865                         dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \
866                         dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \
867                         dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \
868                         dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB; 
869 #define COLL_CALCULATE_NONEQTENSOR(csolev, srcArray ) \
870                         lcsmqadd  = (srcArray##NE - lcsmeq[ dNE ]); \
871                         lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \
872                         lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \
873                         lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
874                         lcsmqo = (lcsmqadd*    lcsmqadd); \
875                         lcsmqadd  = (srcArray##ET - lcsmeq[  dET ]); \
876                         lcsmqadd -= (srcArray##EB - lcsmeq[  dEB ]); \
877                         lcsmqadd -= (srcArray##WT - lcsmeq[  dWT ]); \
878                         lcsmqadd += (srcArray##WB - lcsmeq[  dWB ]); \
879                         lcsmqo += (lcsmqadd*    lcsmqadd); \
880                         lcsmqadd  = (srcArray##NT - lcsmeq[  dNT ]); \
881                         lcsmqadd -= (srcArray##NB - lcsmeq[  dNB ]); \
882                         lcsmqadd -= (srcArray##ST - lcsmeq[  dST ]); \
883                         lcsmqadd += (srcArray##SB - lcsmeq[  dSB ]); \
884                         lcsmqo += (lcsmqadd*    lcsmqadd); \
885                         lcsmqo *= 2.0; \
886                         lcsmqadd  = (srcArray##E  -  lcsmeq[ dE  ]); \
887                         lcsmqadd += (srcArray##W  -  lcsmeq[ dW  ]); \
888                         lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
889                         lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
890                         lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
891                         lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
892                         lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
893                         lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
894                         lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
895                         lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
896                         lcsmqo += (lcsmqadd*    lcsmqadd); \
897                         lcsmqadd  = (srcArray##N  -  lcsmeq[ dN  ]); \
898                         lcsmqadd += (srcArray##S  -  lcsmeq[ dS  ]); \
899                         lcsmqadd += (srcArray##NE -  lcsmeq[ dNE ]); \
900                         lcsmqadd += (srcArray##NW -  lcsmeq[ dNW ]); \
901                         lcsmqadd += (srcArray##SE -  lcsmeq[ dSE ]); \
902                         lcsmqadd += (srcArray##SW -  lcsmeq[ dSW ]); \
903                         lcsmqadd += (srcArray##NT  - lcsmeq[ dNT ]); \
904                         lcsmqadd += (srcArray##NB  - lcsmeq[ dNB ]); \
905                         lcsmqadd += (srcArray##ST  - lcsmeq[ dST ]); \
906                         lcsmqadd += (srcArray##SB  - lcsmeq[ dSB ]); \
907                         lcsmqo += (lcsmqadd*    lcsmqadd); \
908                         lcsmqadd  = (srcArray##T  -  lcsmeq[ dT  ]); \
909                         lcsmqadd += (srcArray##B  -  lcsmeq[ dB  ]); \
910                         lcsmqadd += (srcArray##NT -  lcsmeq[ dNT ]); \
911                         lcsmqadd += (srcArray##NB -  lcsmeq[ dNB ]); \
912                         lcsmqadd += (srcArray##ST -  lcsmeq[ dST ]); \
913                         lcsmqadd += (srcArray##SB -  lcsmeq[ dSB ]); \
914                         lcsmqadd += (srcArray##ET  - lcsmeq[ dET ]); \
915                         lcsmqadd += (srcArray##EB  - lcsmeq[ dEB ]); \
916                         lcsmqadd += (srcArray##WT  - lcsmeq[ dWT ]); \
917                         lcsmqadd += (srcArray##WB  - lcsmeq[ dWB ]); \
918                         lcsmqo += (lcsmqadd*    lcsmqadd); \
919                         lcsmqo = sqrt(lcsmqo); /* FIXME check effect of sqrt*/ \
920
921 //                      COLL_CALCULATE_CSMOMEGAVAL(csolev, lcsmomega); 
922
923 // careful - need lcsmqo 
924 #define COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega ) \
925                         dstomega =  1.0/\
926                                         ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*(\
927                                                         -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \
928                                                         / (6.0*mLevel[(csolev)].lcsmago_sqr)) \
929                                                 ) +0.5 ); 
930
931 #define DEFAULT_COLLIDE_LES \
932                         rho = + MSRC_C  + MSRC_N  \
933                                 + MSRC_S  + MSRC_E  \
934                                 + MSRC_W  + MSRC_T  \
935                                 + MSRC_B  + MSRC_NE \
936                                 + MSRC_NW + MSRC_SE \
937                                 + MSRC_SW + MSRC_NT \
938                                 + MSRC_NB + MSRC_ST \
939                                 + MSRC_SB + MSRC_ET \
940                                 + MSRC_EB + MSRC_WT \
941                                 + MSRC_WB; \
942                         \
943                         ux += MSRC_E - MSRC_W \
944                                 + MSRC_NE - MSRC_NW \
945                                 + MSRC_SE - MSRC_SW \
946                                 + MSRC_ET + MSRC_EB \
947                                 - MSRC_WT - MSRC_WB ;  \
948                         \
949                         uy += MSRC_N - MSRC_S \
950                                 + MSRC_NE + MSRC_NW \
951                                 - MSRC_SE - MSRC_SW \
952                                 + MSRC_NT + MSRC_NB \
953                                 - MSRC_ST - MSRC_SB ;  \
954                         \
955                         uz += MSRC_T - MSRC_B \
956                                 + MSRC_NT - MSRC_NB \
957                                 + MSRC_ST - MSRC_SB \
958                                 + MSRC_ET - MSRC_EB \
959                                 + MSRC_WT - MSRC_WB ;  \
960                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
961                         COLL_CALCULATE_DFEQ(lcsmeq); \
962                         COLL_CALCULATE_NONEQTENSOR(lev, MSRC_)\
963                         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
964                         CSMOMEGA_STATS(lev,lcsmomega); \
965                         \
966                         RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
967                         \
968                         RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ]; \
969                         RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ]; \
970                         RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
971                         RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ]; \
972                         RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ]; \
973                         RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
974                         \
975                         RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
976                         RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
977                         RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
978                         RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
979                         RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
980                         RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
981                         RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
982                         RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
983                         RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
984                         RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
985                         RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
986                         RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; 
987
988 #define DEFAULT_COLLIDE_NOLES \
989                         rho = + MSRC_C  + MSRC_N  \
990                                 + MSRC_S  + MSRC_E  \
991                                 + MSRC_W  + MSRC_T  \
992                                 + MSRC_B  + MSRC_NE \
993                                 + MSRC_NW + MSRC_SE \
994                                 + MSRC_SW + MSRC_NT \
995                                 + MSRC_NB + MSRC_ST \
996                                 + MSRC_SB + MSRC_ET \
997                                 + MSRC_EB + MSRC_WT \
998                                 + MSRC_WB; \
999                         \
1000                         ux += MSRC_E - MSRC_W \
1001                                 + MSRC_NE - MSRC_NW \
1002                                 + MSRC_SE - MSRC_SW \
1003                                 + MSRC_ET + MSRC_EB \
1004                                 - MSRC_WT - MSRC_WB ;  \
1005                         \
1006                         uy += MSRC_N - MSRC_S \
1007                                 + MSRC_NE + MSRC_NW \
1008                                 - MSRC_SE - MSRC_SW \
1009                                 + MSRC_NT + MSRC_NB \
1010                                 - MSRC_ST - MSRC_SB ;  \
1011                         \
1012                         uz += MSRC_T - MSRC_B \
1013                                 + MSRC_NT - MSRC_NB \
1014                                 + MSRC_ST - MSRC_SB \
1015                                 + MSRC_ET - MSRC_EB \
1016                                 + MSRC_WT - MSRC_WB ;  \
1017                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
1018                         \
1019                         RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C  + OMEGA(lev)*EQC ; \
1020                         \
1021                         RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N  + OMEGA(lev)*EQN ; \
1022                         RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S  + OMEGA(lev)*EQS ; \
1023                         RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E  + OMEGA(lev)*EQE ; \
1024                         RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W  + OMEGA(lev)*EQW ; \
1025                         RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T  + OMEGA(lev)*EQT ; \
1026                         RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B  + OMEGA(lev)*EQB ; \
1027                         \
1028                         RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \
1029                         RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \
1030                         RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \
1031                         RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \
1032                         RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \
1033                         RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \
1034                         RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \
1035                         RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \
1036                         RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \
1037                         RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \
1038                         RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \
1039                         RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB; 
1040
1041
1042
1043 #define OPTIMIZED_STREAMCOLLIDE_LES \
1044                         /* only surrounded by fluid cells...!, so safe streaming here... */ \
1045                         m[dC ] = CSRC_C ; \
1046                         m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
1047                         m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
1048                         m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
1049                         m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
1050                         m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
1051                         m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
1052                         \
1053                         rho = MSRC_C  + MSRC_N + MSRC_S  + MSRC_E + MSRC_W  + MSRC_T  \
1054                                 + MSRC_B  + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \
1055                                 + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \
1056                         ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \
1057                                 + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB + mLevel[lev].gravity[0];  \
1058                         uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \
1059                                 + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB + mLevel[lev].gravity[1];  \
1060                         uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \
1061                                 + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB + mLevel[lev].gravity[2];  \
1062                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
1063                         COLL_CALCULATE_DFEQ(lcsmeq); \
1064                         COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \
1065                         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
1066                         CSMOMEGA_STATS(lev,lcsmomega); \
1067                         \
1068                         RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C  + lcsmomega*EQC ; \
1069                         RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N  + lcsmomega*lcsmeq[ dN ];  \
1070                         RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S  + lcsmomega*lcsmeq[ dS ];  \
1071                         RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E  + lcsmomega*lcsmeq[ dE ]; \
1072                         RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W  + lcsmomega*lcsmeq[ dW ];  \
1073                         RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T  + lcsmomega*lcsmeq[ dT ];  \
1074                         RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B  + lcsmomega*lcsmeq[ dB ]; \
1075                         \
1076                         RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE];  \
1077                         RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW];  \
1078                         RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE];  \
1079                         RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
1080                         \
1081                         RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT];  \
1082                         RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB];  \
1083                         RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST];  \
1084                         RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
1085                         \
1086                         RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET];  \
1087                         RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
1088                         RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT];  \
1089                         RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB];  \
1090
1091 #define OPTIMIZED_STREAMCOLLIDE_UNUSED \
1092                         /* only surrounded by fluid cells...!, so safe streaming here... */ \
1093                         rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T  \
1094                                 + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
1095                                 + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
1096                         ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
1097                                 + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB + mLevel[lev].gravity[0];  \
1098                         uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
1099                                 + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB + mLevel[lev].gravity[1];  \
1100                         uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
1101                                 + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB + mLevel[lev].gravity[2];  \
1102                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
1103                         COLL_CALCULATE_DFEQ(lcsmeq); \
1104                         COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \
1105                         COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
1106                         \
1107                         RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C  + lcsmomega*EQC ; \
1108                         RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N  + lcsmomega*lcsmeq[ dN ];  \
1109                         RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S  + lcsmomega*lcsmeq[ dS ];  \
1110                         RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E  + lcsmomega*lcsmeq[ dE ]; \
1111                         RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W  + lcsmomega*lcsmeq[ dW ];  \
1112                         RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T  + lcsmomega*lcsmeq[ dT ];  \
1113                         RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B  + lcsmomega*lcsmeq[ dB ]; \
1114                         \
1115                         RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE];  \
1116                         RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW];  \
1117                         RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE];  \
1118                         RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \
1119                         \
1120                         RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT];  \
1121                         RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB];  \
1122                         RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST];  \
1123                         RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \
1124                         \
1125                         RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET];  \
1126                         RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \
1127                         RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT];  \
1128                         RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB];  \
1129
1130 #define OPTIMIZED_STREAMCOLLIDE_NOLES \
1131                         /* only surrounded by fluid cells...!, so safe streaming here... */ \
1132                         rho = CSRC_C  + CSRC_N + CSRC_S  + CSRC_E + CSRC_W  + CSRC_T  \
1133                                 + CSRC_B  + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
1134                                 + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
1135                         ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
1136                                 + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB + mLevel[lev].gravity[0];  \
1137                         uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
1138                                 + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB + mLevel[lev].gravity[1];  \
1139                         uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
1140                                 + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB + mLevel[lev].gravity[2];  \
1141                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
1142                         RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C  + OMEGA(lev)*EQC ; \
1143                         RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N  + OMEGA(lev)*EQN ;  \
1144                         RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S  + OMEGA(lev)*EQS ;  \
1145                         RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E  + OMEGA(lev)*EQE ; \
1146                         RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W  + OMEGA(lev)*EQW ;  \
1147                         RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T  + OMEGA(lev)*EQT ;  \
1148                         RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B  + OMEGA(lev)*EQB ; \
1149                          \
1150                         RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE;  \
1151                         RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW;  \
1152                         RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE;  \
1153                         RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \
1154                          \
1155                         RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT;  \
1156                         RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB;  \
1157                         RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST;  \
1158                         RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \
1159                          \
1160                         RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET;  \
1161                         RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \
1162                         RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT;  \
1163                         RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB;  \
1164
1165
1166 // debug version1
1167 #define STREAMCHECK(ni,nj,nk,nl) 
1168 #define COLLCHECK
1169 #define OPTIMIZED_STREAMCOLLIDE_DEBUG \
1170                         m[0] = RAC(ccel,0); \
1171                         FORDF1 { /* df0 is set later on... */ \
1172                                 if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); D::mPanic=1;  \
1173                                 } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
1174                                 STREAMCHECK(i+D::dfVecX[D::dfInv[l]], j+D::dfVecY[D::dfInv[l]],k+D::dfVecZ[D::dfInv[l]], l); \
1175                         }   \
1176                         rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
1177                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
1178                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago , &mDebugOmegaRet  ); \
1179                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
1180                         FORDF0 { RAC(tcel,l) = m[l]; } \
1181                         usqr = 1.5 * (ux*ux + uy*uy + uz*uz);  \
1182                         COLLCHECK;
1183
1184
1185
1186 // more debugging
1187 /*DEBUG \
1188                         m[0] = RAC(ccel,0); \
1189                         FORDF1 { \
1190                                 if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); D::mPanic=1;  \
1191                                 } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
1192                         }   \
1193 errMsg("T","QSDM at %d,%d,%d  lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lcsmomega ); \
1194                         rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
1195                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
1196                         D::collideArrays( m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago  , &mDebugOmegaRet ); \
1197                         CSMOMEGA_STATS(lev,mDebugOmegaRet); \
1198                         */
1199 #if USE_LES==1
1200 #define DEFAULT_COLLIDE DEFAULT_COLLIDE_LES
1201 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_LES
1202 #else 
1203 #define DEFAULT_COLLIDE DEFAULT_COLLIDE_NOLES
1204 #define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_NOLES
1205 #endif
1206
1207 #endif
1208
1209 #define USQRMAXCHECK(Cusqr,Cux,Cuy,Cuz,  CmMaxVlen,CmMxvx,CmMxvy,CmMxvz) \
1210                         if(Cusqr>CmMaxVlen) { \
1211                                 CmMxvx = Cux; CmMxvy = Cuy; CmMxvz = Cuz; CmMaxVlen = Cusqr; \
1212                         } /* stats */ 
1213
1214
1215
1216 // loops over _all_ cells (including boundary layer)
1217 #define FSGR_FORIJK_BOUNDS(leveli) \
1218         for(int k= getForZMinBnd(); k< getForZMaxBnd(leveli); ++k) \
1219    for(int j=0;j<mLevel[leveli].lSizey-0;++j) \
1220     for(int i=0;i<mLevel[leveli].lSizex-0;++i) \
1221         
1222 // loops over _only inner_ cells 
1223 #define FSGR_FORIJK1(leveli) \
1224         for(int k= getForZMin1(); k< getForZMax1(leveli); ++k) \
1225    for(int j=1;j<mLevel[leveli].lSizey-1;++j) \
1226     for(int i=1;i<mLevel[leveli].lSizex-1;++i) \
1227
1228 // relaxation_macros end
1229
1230
1231 /******************************************************************************
1232  * Lbm Constructor
1233  *****************************************************************************/
1234 template<class D>
1235 LbmFsgrSolver<D>::LbmFsgrSolver() :
1236         D(),
1237         mCurrentMass(0.0), mCurrentVolume(0.0),
1238         mNumProblems(0), 
1239         mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
1240         mpPreviewSurface(NULL), 
1241         mSmoothSurface(0.0), mSmoothNormals(0.0),
1242         mTimeAdap(false), 
1243         mOutputSurfacePreview(0), mPreviewFactor(0.25),
1244         mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
1245         mGfxGeoSetup(0), mGfxEndTime(-1.0), mInitSurfaceSmoothing(0),
1246         mTimestepReduceLock(0),
1247         mTimeSwitchCounts(0), 
1248         mSimulationTime(0.0),
1249         mMinStepTime(0.0), mMaxStepTime(0.0),
1250         mMaxNoCells(0), mMinNoCells(0), mAvgNumUsedCells(0),
1251         mDropMode(1), mDropSize(0.15), mDropSpeed(0.0),
1252         mDropping(false),
1253         mDropX(0.0), mDropY(0.0), mDropHeight(0.8),
1254         mIsoWeightMethod(2),
1255         mMaxRefine(1), 
1256         mDfScaleUp(-1.0), mDfScaleDown(-1.0),
1257         mInitialCsmago(0.04), mDebugOmegaRet(0.0),
1258         mNumInvIfTotal(0), mNumFsgrChanges(0),
1259         mDisableStandingFluidInit(0),
1260         mForceTadapRefine(-1)
1261 {
1262   // not much to do here... 
1263         D::mpIso = new IsoSurface( D::mIsoValue, false );
1264
1265   // init equilibrium dist. func 
1266   LbmFloat rho=1.0;
1267   FORDF0 {
1268                 D::dfEquil[l] = D::getCollideEq( l,rho,  0.0, 0.0, 0.0);
1269   }
1270
1271         // init LES
1272         int odm = 0;
1273         for(int m=0; m<D::cDimension; m++) { 
1274                 for(int l=0; l<D::cDfNum; l++) { 
1275                         D::lesCoeffDiag[m][l] = 
1276                         D::lesCoeffOffdiag[m][l] = 0.0;
1277                 }
1278         }
1279         for(int m=0; m<D::cDimension; m++) { 
1280                 for(int n=0; n<D::cDimension; n++) { 
1281                         for(int l=1; l<D::cDfNum; l++) { 
1282                                 LbmFloat em;
1283                                 switch(m) {
1284                                         case 0: em = D::dfDvecX[l]; break;
1285                                         case 1: em = D::dfDvecY[l]; break;
1286                                         case 2: em = D::dfDvecZ[l]; break;
1287                                         default: em = -1.0; errFatal("SMAGO1","err m="<<m, SIMWORLD_GENERICERROR);
1288                                 }
1289                                 LbmFloat en;
1290                                 switch(n) {
1291                                         case 0: en = D::dfDvecX[l]; break;
1292                                         case 1: en = D::dfDvecY[l]; break;
1293                                         case 2: en = D::dfDvecZ[l]; break;
1294                                         default: en = -1.0; errFatal("SMAGO2","err n="<<n, SIMWORLD_GENERICERROR);
1295                                 }
1296                                 const LbmFloat coeff = em*en;
1297                                 if(m==n) {
1298                                         D::lesCoeffDiag[m][l] = coeff;
1299                                 } else {
1300                                         if(m>n) {
1301                                                 D::lesCoeffOffdiag[odm][l] = coeff;
1302                                         }
1303                                 }
1304                         }
1305
1306                         if(m==n) {
1307                         } else {
1308                                 if(m>n) odm++;
1309                         }
1310                 }
1311         }
1312
1313         mDvecNrm[0] = LbmVec(0.0);
1314   FORDF1 {
1315                 mDvecNrm[l] = getNormalized( 
1316                         LbmVec(D::dfDvecX[D::dfInv[l]], D::dfDvecY[D::dfInv[l]], D::dfDvecZ[D::dfInv[l]] ) 
1317                         ) * -1.0; 
1318         }
1319
1320         addDrop(false,0,0);
1321 }
1322
1323 /*****************************************************************************/
1324 /* Destructor */
1325 /*****************************************************************************/
1326 template<class D>
1327 LbmFsgrSolver<D>::~LbmFsgrSolver()
1328 {
1329   if(!D::mInitDone){ debugOut("LbmFsgrSolver::LbmFsgrSolver : not inited...",0); return; }
1330
1331 #if COMPRESSGRIDS==1
1332         delete mLevel[mMaxRefine].mprsCells[1];
1333         mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
1334 #endif // COMPRESSGRIDS==1
1335
1336         for(int i=0; i<=mMaxRefine; i++) {
1337                 for(int s=0; s<2; s++) {
1338                         if(mLevel[i].mprsCells[s]) delete [] mLevel[i].mprsCells[s];
1339                         if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s];
1340                 }
1341         }
1342         delete D::mpIso;
1343         if(mpPreviewSurface) delete mpPreviewSurface;
1344
1345         // always output performance estimate
1346         debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
1347   if(!D::mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
1348 }
1349
1350
1351
1352
1353 /******************************************************************************
1354  * initilize variables fom attribute list 
1355  *****************************************************************************/
1356 template<class D>
1357 void 
1358 LbmFsgrSolver<D>::parseAttrList()
1359 {
1360         LbmSolverInterface::parseStdAttrList();
1361
1362         string matIso("default");
1363         matIso = D::mpAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
1364         D::mpIso->setMaterialName( matIso );
1365         mOutputSurfacePreview = D::mpAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false );
1366         mTimeAdap = D::mpAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
1367
1368         mIsoWeightMethod= D::mpAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
1369         mInitSurfaceSmoothing = D::mpAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
1370         mSmoothSurface = D::mpAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
1371         mSmoothNormals = D::mpAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
1372
1373         mInitialCsmago = D::mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
1374         // deprecated!
1375         float mInitialCsmagoCoarse = 0.0;
1376         mInitialCsmagoCoarse = D::mpAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
1377 #if USE_LES==1
1378 #else // USE_LES==1
1379         debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
1380         mInitialCsmago = 0.0;
1381 #endif // USE_LES==1
1382
1383         // refinement
1384         mMaxRefine  = D::mpAttrs->readInt("maxrefine",  mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", true);
1385         mDisableStandingFluidInit = D::mpAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
1386         mForceTadapRefine = D::mpAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
1387
1388         // demo mode settings
1389         mDropMode = D::mpAttrs->readInt("dropmode", mDropMode, "SimulationLbm","mDropMode", false );
1390         mDropSize = D::mpAttrs->readFloat("dropsize", mDropSize, "SimulationLbm","mDropSize", false );
1391         mDropHeight = D::mpAttrs->readFloat("dropheight", mDropHeight, "SimulationLbm","mDropHeight", false );
1392         mDropSpeed = vec2L( D::mpAttrs->readVec3d("dropspeed", ntlVec3d(0.0), "SimulationLbm","mDropSpeed", false ) );
1393         if( (mDropMode>2) || (mDropMode<-1) ) mDropMode=1;
1394         mGfxGeoSetup = D::mpAttrs->readInt("gfxgeosetup", mGfxGeoSetup, "SimulationLbm","mGfxGeoSetup", false );
1395         mGfxEndTime = D::mpAttrs->readFloat("gfxendtime", mGfxEndTime, "SimulationLbm","mGfxEndTime", false );
1396         mFVHeight = D::mpAttrs->readFloat("fvolheight", mFVHeight, "SimulationLbm","mFVHeight", false );
1397         mFVArea   = D::mpAttrs->readFloat("fvolarea", mFVArea, "SimulationLbm","mFArea", false );
1398
1399 }
1400
1401
1402 /******************************************************************************
1403  * Initialize omegas and forces on all levels (for init/timestep change)
1404  *****************************************************************************/
1405 template<class D>
1406 void 
1407 LbmFsgrSolver<D>::initLevelOmegas()
1408 {
1409         // no explicit settings
1410         D::mOmega = D::mpParam->calculateOmega();
1411         D::mGravity = vec2L( D::mpParam->calculateGravity() );
1412         D::mSurfaceTension = D::mpParam->calculateSurfaceTension(); // unused
1413
1414         if(mInitialCsmago<=0.0) {
1415                 if(OPT3D==1) {
1416                         errFatal("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version...",SIMWORLD_INITERROR); 
1417                         return;
1418                 }
1419         }
1420
1421         // use Tau instead of Omega for calculations
1422         { // init base level
1423                 int i = mMaxRefine;
1424                 mLevel[i].omega    = D::mOmega;
1425                 mLevel[i].stepsize = D::mpParam->getStepTime();
1426                 mLevel[i].lcsmago = mInitialCsmago; //CSMAGO_INITIAL;
1427                 mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
1428                 mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
1429         }
1430
1431         // init all sub levels
1432         for(int i=mMaxRefine-1; i>=0; i--) {
1433                 //mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5;
1434                 double nomega = 0.5 * (  (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5;
1435                 nomega                = 1.0/nomega;
1436                 mLevel[i].omega       = (LbmFloat)nomega;
1437                 mLevel[i].stepsize    = 2.0 * mLevel[i+1].stepsize;
1438                 //mLevel[i].lcsmago     = mLevel[i+1].lcsmago*mCsmagoRefineMultiplier;
1439                 //if(mLevel[i].lcsmago>1.0) mLevel[i].lcsmago = 1.0;
1440                 //if(strstr(D::getName().c_str(),"Debug")){ 
1441                 //mLevel[i].lcsmago = mLevel[mMaxRefine].lcsmago; // DEBUG
1442                 // if(strstr(D::getName().c_str(),"Debug")) mLevel[i].lcsmago = mLevel[mMaxRefine].lcsmago * (LbmFloat)(mMaxRefine-i)*0.5+1.0; 
1443                 //if(strstr(D::getName().c_str(),"Debug")) mLevel[i].lcsmago = mLevel[mMaxRefine].lcsmago * ((LbmFloat)(mMaxRefine-i)*1.0 + 1.0 ); 
1444                 //if(strstr(D::getName().c_str(),"Debug")) mLevel[i].lcsmago = 0.99;
1445                 mLevel[i].lcsmago = mInitialCsmago;
1446                 mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
1447                 mLevel[i].lcnu        = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
1448         }
1449         
1450         // for lbgk
1451         mLevel[ mMaxRefine ].gravity = D::mGravity / mLevel[ mMaxRefine ].omega;
1452         for(int i=mMaxRefine-1; i>=0; i--) {
1453                 // should be the same on all levels...
1454                 // for lbgk
1455                 mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega;
1456         }
1457
1458         // debug? invalidate old values...
1459         D::mGravity = -100.0;
1460         D::mOmega = -100.0;
1461
1462         for(int i=0; i<=mMaxRefine; i++) {
1463                 if(!D::mSilent) {
1464                         errMsg("LbmFsgrSolver", "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" offs:"<<mLevel[i].lOffsx<<","<<mLevel[i].lOffsy<<","<<mLevel[i].lOffsz 
1465                                         <<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
1466                                         <<" cmsagp:"<<mLevel[i].lcsmago<<", "
1467                                         << " ss"<<mLevel[i].stepsize<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
1468                 } else {
1469                         if(!D::mInitDone) {
1470                                 debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" "
1471                                                 <<"omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity , 5);
1472                         }
1473                 }
1474         }
1475         if(mMaxRefine>0) {
1476                 mDfScaleUp   = (mLevel[0  ].stepsize/mLevel[0+1].stepsize)* (1.0/mLevel[0  ].omega-1.0)/ (1.0/mLevel[0+1].omega-1.0); // yu
1477                 mDfScaleDown = (mLevel[0+1].stepsize/mLevel[0  ].stepsize)* (1.0/mLevel[0+1].omega-1.0)/ (1.0/mLevel[0  ].omega-1.0); // yu
1478         }
1479 }
1480
1481
1482 /******************************************************************************
1483  * Init Solver (values should be read from config file)
1484  *****************************************************************************/
1485 template<class D>
1486 bool 
1487 LbmFsgrSolver<D>::initialize( ntlTree* /*tree*/, vector<ntlGeometryObject*>* /*objects*/ )
1488 {
1489   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... (Layout:"<<ALSTRING<<") ",1);
1490
1491         // fix size inits to force cubic cells and mult4 level dimensions
1492         const int debugGridsizeInit = 1;
1493         mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)D::mSizex;
1494         int maxGridSize = D::mSizex; // get max size
1495         if(D::mSizey>maxGridSize) maxGridSize = D::mSizey;
1496         if(D::mSizez>maxGridSize) maxGridSize = D::mSizez;
1497         LbmFloat maxGeoSize = (D::mvGeoEnd[0]-D::mvGeoStart[0]); // get max size
1498         if((D::mvGeoEnd[1]-D::mvGeoStart[1])>maxGridSize) maxGeoSize = (D::mvGeoEnd[1]-D::mvGeoStart[1]);
1499         if((D::mvGeoEnd[2]-D::mvGeoStart[2])>maxGridSize) maxGeoSize = (D::mvGeoEnd[2]-D::mvGeoStart[2]);
1500         // FIXME better divide max geo size by corresponding resolution rather than max? no prob for rx==ry==rz though
1501         LbmFloat cellSize = (maxGeoSize / (LbmFloat)maxGridSize);
1502   if(debugGridsizeInit) debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Start:"<<D::mvGeoStart<<" End:"<<D::mvGeoEnd<<" maxS:"<<maxGeoSize<<" maxG:"<<maxGridSize<<" cs:"<<cellSize, 10);
1503         // force grid sizes according to geom. size, rounded
1504         D::mSizex = (int) ((D::mvGeoEnd[0]-D::mvGeoStart[0]) / cellSize +0.5);
1505         D::mSizey = (int) ((D::mvGeoEnd[1]-D::mvGeoStart[1]) / cellSize +0.5);
1506         D::mSizez = (int) ((D::mvGeoEnd[2]-D::mvGeoStart[2]) / cellSize +0.5);
1507         // match refinement sizes, round downwards to multiple of 4
1508         int sizeMask = 0;
1509         int maskBits = mMaxRefine;
1510         if(PARALLEL==1) maskBits+=2;
1511         for(int i=0; i<maskBits; i++) { sizeMask |= (1<<i); }
1512         sizeMask = ~sizeMask;
1513   if(debugGridsizeInit) debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Size X:"<<D::mSizex<<" Y:"<<D::mSizey<<" Z:"<<D::mSizez<<" m"<<convertCellFlagType2String(sizeMask) ,10);
1514         D::mSizex &= sizeMask;
1515         D::mSizey &= sizeMask;
1516         D::mSizez &= sizeMask;
1517         // force geom size to match rounded grid sizes
1518         D::mvGeoEnd[0] = D::mvGeoStart[0] + cellSize*(LbmFloat)D::mSizex;
1519         D::mvGeoEnd[1] = D::mvGeoStart[1] + cellSize*(LbmFloat)D::mSizey;
1520         D::mvGeoEnd[2] = D::mvGeoStart[2] + cellSize*(LbmFloat)D::mSizez;
1521
1522   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Final domain size X:"<<D::mSizex<<" Y:"<<D::mSizey<<" Z:"<<D::mSizez<<
1523                         ", Domain: "<<D::mvGeoStart<<":"<<D::mvGeoEnd<<", "<<(D::mvGeoEnd-D::mvGeoStart) ,2);
1524   //debMsgStd("LbmFsgrSolver::initialize",DM_MSG, ,2);
1525         D::mpParam->setSize(D::mSizex, D::mSizey, D::mSizez);
1526
1527   //debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Size X:"<<D::mSizex<<" Y:"<<D::mSizey<<" Z:"<<D::mSizez ,2);
1528
1529 #if ELBEEM_BLENDER!=1
1530   debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
1531                 <<"LBM_EPSILON="<<LBM_EPSILON       <<" "
1532                 <<"FSGR_STRICT_DEBUG="<<FSGR_STRICT_DEBUG       <<" "
1533                 <<"INTORDER="<<INTORDER        <<" "
1534                 <<"REFINEMENTBORDER="<<REFINEMENTBORDER        <<" "
1535                 <<"OPT3D="<<OPT3D        <<" "
1536                 <<"COMPRESSGRIDS="<<COMPRESSGRIDS<<" "
1537                 <<"LS_FLUIDTHRESHOLD="<<LS_FLUIDTHRESHOLD        <<" "
1538                 <<"MASS_INVALID="<<MASS_INVALID        <<" "
1539                 <<"FSGR_LISTTRICK="<<FSGR_LISTTRICK            <<" "
1540                 <<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY          <<" "
1541                 <<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL           <<" "
1542                 <<"FSGR_MAGICNR="<<FSGR_MAGICNR              <<" " 
1543                 <<"USE_LES="<<USE_LES              <<" " 
1544                 ,10);
1545 #endif // ELBEEM_BLENDER!=1
1546
1547         // perform 2D corrections...
1548         if(D::cDimension == 2) D::mSizez = 1;
1549
1550         D::mpParam->setSimulationMaxSpeed(0.0);
1551         if(mFVHeight>0.0) D::mpParam->setFluidVolumeHeight(mFVHeight);
1552         D::mpParam->setTadapLevels( mMaxRefine+1 );
1553
1554         if(mForceTadapRefine>mMaxRefine) {
1555                 D::mpParam->setTadapLevels( mForceTadapRefine+1 );
1556                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
1557         }
1558
1559         if(!D::mpParam->calculateAllMissingValues()) {
1560                 errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR);
1561                 return false;
1562         }
1563         // recalc objects speeds in geo init
1564
1565
1566
1567         // init vectors
1568         if(mMaxRefine >= MAX_LEV) {
1569                 errFatal("LbmFsgrSolver::initializeLbmGridref"," error: Too many levels!", SIMWORLD_INITERROR);
1570                 return false;
1571         }
1572         for(int i=0; i<=mMaxRefine; i++) {
1573                 mLevel[i].id = i;
1574                 mLevel[i].nodeSize = 0.0; 
1575                 mLevel[i].simCellSize = 0.0; 
1576                 mLevel[i].omega = 0.0; 
1577                 mLevel[i].time = 0.0; 
1578                 mLevel[i].stepsize = 1.0;
1579                 mLevel[i].gravity = LbmVec(0.0); 
1580                 mLevel[i].mprsCells[0] = NULL;
1581                 mLevel[i].mprsCells[1] = NULL;
1582                 mLevel[i].mprsFlags[0] = NULL;
1583                 mLevel[i].mprsFlags[1] = NULL;
1584
1585                 mLevel[i].avgOmega = 0.0; 
1586                 mLevel[i].avgOmegaCnt = 0.0;
1587         }
1588
1589         // init sizes
1590         mLevel[mMaxRefine].lSizex = D::mSizex;
1591         mLevel[mMaxRefine].lSizey = D::mSizey;
1592         mLevel[mMaxRefine].lSizez = D::mSizez;
1593         for(int i=mMaxRefine-1; i>=0; i--) {
1594                 mLevel[i].lSizex = mLevel[i+1].lSizex/2;
1595                 mLevel[i].lSizey = mLevel[i+1].lSizey/2;
1596                 mLevel[i].lSizez = mLevel[i+1].lSizez/2;
1597                 /*if( ((mLevel[i].lSizex % 4) != 0) || ((mLevel[i].lSizey % 4) != 0) || ((mLevel[i].lSizez % 4) != 0) ) {
1598                         errMsg("LbmFsgrSolver","Init: error invalid sizes on level "<<i<<" "<<PRINT_VEC(mLevel[i].lSizex,mLevel[i].lSizey,mLevel[i].lSizez) );
1599                         xit(1);
1600                 }// old QUAD handling */
1601         }
1602
1603         // estimate memory usage
1604         {
1605                 unsigned long int memCnt = 0;
1606                 unsigned long int rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
1607                 memCnt += sizeof(CellFlagType) * (rcellSize/dTotalNum +4) *2;
1608 #if COMPRESSGRIDS==0
1609                 memCnt += sizeof(LbmFloat) * (rcellSize +4) *2;
1610 #else // COMPRESSGRIDS==0
1611                 unsigned long int compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
1612                 memCnt += sizeof(LbmFloat) * (rcellSize+compressOffset +4);
1613 #endif // COMPRESSGRIDS==0
1614                 for(int i=mMaxRefine-1; i>=0; i--) {
1615                         rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
1616                         memCnt += sizeof(CellFlagType) * (rcellSize/dTotalNum +4) *2;
1617                         memCnt += sizeof(LbmFloat) * (rcellSize +4) *2;
1618                 }
1619                 double memd = memCnt;
1620                 char *sizeStr = "";
1621                 const double sfac = 1000.0;
1622                 if(memd>sfac){ memd /= sfac; sizeStr="KB"; }
1623                 if(memd>sfac){ memd /= sfac; sizeStr="MB"; }
1624                 if(memd>sfac){ memd /= sfac; sizeStr="GB"; }
1625                 if(memd>sfac){ memd /= sfac; sizeStr="TB"; }
1626                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Required Grid memory: "<< memd <<" "<< sizeStr<<" ",4);
1627         }
1628
1629         // safety check
1630         if(sizeof(CellFlagType) != CellFlagTypeSize) {
1631                 errFatal("LbmFsgrSolver::initialize","Fatal Error: CellFlagType has wrong size! Is:"<<sizeof(CellFlagType)<<", should be:"<<CellFlagTypeSize, SIMWORLD_GENERICERROR);
1632                 return false;
1633         }
1634
1635         mLevel[ mMaxRefine ].nodeSize = ((D::mvGeoEnd[0]-D::mvGeoStart[0]) / (LbmFloat)(D::mSizex));
1636         mLevel[ mMaxRefine ].simCellSize = D::mpParam->getCellSize();
1637         mLevel[ mMaxRefine ].lcellfactor = 1.0;
1638         unsigned long int rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
1639         // +4 for safety ?
1640         mLevel[ mMaxRefine ].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1641         mLevel[ mMaxRefine ].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1642
1643 #if COMPRESSGRIDS==0
1644         mLevel[ mMaxRefine ].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
1645         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
1646 #else // COMPRESSGRIDS==0
1647         unsigned long int compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
1648         mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +compressOffset +4 ];
1649         mLevel[ mMaxRefine ].mprsCells[0] = mLevel[ mMaxRefine ].mprsCells[1]+compressOffset;
1650         //errMsg("CGD","rcs:"<<rcellSize<<" cpff:"<<compressOffset<< " c0:"<<mLevel[ mMaxRefine ].mprsCells[0]<<" c1:"<<mLevel[ mMaxRefine ].mprsCells[1]<< " c0e:"<<(mLevel[ mMaxRefine ].mprsCells[0]+rcellSize)<<" c1:"<<(mLevel[ mMaxRefine ].mprsCells[1]+rcellSize)); // DEBUG
1651 #endif // COMPRESSGRIDS==0
1652
1653         LbmFloat lcfdimFac = 8.0;
1654         if(D::cDimension==2) lcfdimFac = 4.0;
1655         for(int i=mMaxRefine-1; i>=0; i--) {
1656                 mLevel[i].nodeSize = 2.0 * mLevel[i+1].nodeSize;
1657                 mLevel[i].simCellSize = 2.0 * mLevel[i+1].simCellSize;
1658                 mLevel[i].lcellfactor = mLevel[i+1].lcellfactor * lcfdimFac;
1659
1660                 if(D::cDimension==2){ mLevel[i].lSizez = 1; } // 2D
1661                 rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
1662                 mLevel[i].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1663                 mLevel[i].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
1664                 mLevel[i].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
1665                 mLevel[i].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
1666         }
1667
1668         // init sizes for _all_ levels
1669         for(int i=mMaxRefine; i>=0; i--) {
1670                 mLevel[i].lOffsx = mLevel[i].lSizex;
1671                 mLevel[i].lOffsy = mLevel[i].lOffsx*mLevel[i].lSizey;
1672                 mLevel[i].lOffsz = mLevel[i].lOffsy*mLevel[i].lSizez;
1673         mLevel[i].setCurr  = 0;
1674         mLevel[i].setOther = 1;
1675         mLevel[i].lsteps = 0;
1676         mLevel[i].lmass = 0.0;
1677         mLevel[i].lvolume = 0.0;
1678         }
1679
1680         // calc omega, force for all levels
1681         initLevelOmegas();
1682         mMinStepTime = D::mpParam->getStepTime();
1683         mMaxStepTime = D::mpParam->getStepTime();
1684
1685         // init isosurf
1686         D::mpIso->setIsolevel( D::mIsoValue );
1687         // approximate feature size with mesh resolution
1688         float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
1689         D::mpIso->setSmoothSurface( mSmoothSurface * featureSize );
1690         D::mpIso->setSmoothNormals( mSmoothNormals * featureSize );
1691
1692         // init iso weight values mIsoWeightMethod
1693         int wcnt = 0;
1694         float totw = 0.0;
1695         for(int ak=-1;ak<=1;ak++) 
1696                 for(int aj=-1;aj<=1;aj++) 
1697                         for(int ai=-1;ai<=1;ai++)  {
1698                                 switch(mIsoWeightMethod) {
1699                                 case 1: // light smoothing
1700                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
1701                                         break;
1702                                 case 2: // very light smoothing
1703                                         mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
1704                                         mIsoWeight[wcnt] *= mIsoWeight[wcnt];
1705                                         break;
1706                                 case 3: // no smoothing
1707                                         if(ai==0 && aj==0 && ak==0) mIsoWeight[wcnt] = 1.0;
1708                                         else mIsoWeight[wcnt] = 0.0;
1709                                         break;
1710                                 default: // strong smoothing (=0)
1711                                         mIsoWeight[wcnt] = 1.0;
1712                                         break;
1713                                 }
1714                                 totw += mIsoWeight[wcnt];
1715                                 wcnt++;
1716                         }
1717         for(int i=0; i<27; i++) mIsoWeight[i] /= totw;
1718
1719         LbmVec isostart = vec2L(D::mvGeoStart);
1720         LbmVec isoend   = vec2L(D::mvGeoEnd);
1721         int twodOff = 0; // 2d slices
1722         if(D::cDimension==2) {
1723                 LbmFloat sn,se;
1724                 sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(D::mSizex+1.0))*0.5;
1725                 se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(D::mSizex+1.0))*0.5;
1726                 isostart[2] = sn;
1727                 isoend[2] = se;
1728                 twodOff = 2;
1729         }
1730         //errMsg(" SETISO ", " "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(D::mSizex+1.0))*0.5)<<" "<<(LbmFloat)(D::mSizex+1.0)<<" " );
1731         D::mpIso->setStart( vec2G(isostart) );
1732         D::mpIso->setEnd(   vec2G(isoend) );
1733         LbmVec isodist = isoend-isostart;
1734         D::mpIso->initializeIsosurface( D::mSizex+2, D::mSizey+2, D::mSizez+2+twodOff, vec2G(isodist) );
1735         for(int ak=0;ak<D::mSizez+2+twodOff;ak++) 
1736                 for(int aj=0;aj<D::mSizey+2;aj++) 
1737                         for(int ai=0;ai<D::mSizex+2;ai++) { *D::mpIso->getData(ai,aj,ak) = 0.0; }
1738
1739   /* init array (set all invalid first) */
1740         for(int lev=0; lev<=mMaxRefine; lev++) {
1741                 FSGR_FORIJK_BOUNDS(lev) {
1742                         RFLAG(lev,i,j,k,0) = RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
1743                         initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); 
1744                 }
1745         }
1746
1747         // init defaults
1748         mAvgNumUsedCells = 0;
1749         D::mFixMass= 0.0;
1750
1751   /* init boundaries */
1752   debugOut("LbmFsgrSolver::initialize : Boundary init...",10);
1753
1754
1755         // use the density init?
1756         initGeometryFlags();
1757         D::initGenericTestCases();
1758         
1759         // new - init noslip 1 everywhere...
1760         // half fill boundary cells?
1761
1762   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1763     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1764                         initEmptyCell(mMaxRefine, i,0,k, CFBnd, 0.0, BND_FILL); 
1765                         initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, CFBnd, 0.0, BND_FILL); 
1766     }
1767
1768         if(D::cDimension == 3) {
1769                 // only for 3D
1770                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1771                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1772                                 initEmptyCell(mMaxRefine, i,j,0, CFBnd, 0.0, BND_FILL); 
1773                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, CFBnd, 0.0, BND_FILL); 
1774                         }
1775         }
1776
1777   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1778     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1779                         initEmptyCell(mMaxRefine, 0,j,k, CFBnd, 0.0, BND_FILL); 
1780                         initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, CFBnd, 0.0, BND_FILL); 
1781                         // DEBUG BORDER!
1782                         //initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, CFBnd, 0.0, BND_FILL); 
1783     }
1784
1785         // TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
1786   /*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1787     for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1788                         initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, CFBnd, 0.0, BND_FILL); 
1789     }
1790   for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1791     for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1792                         initEmptyCell(mMaxRefine, i,1,k, CFBnd, 0.0, BND_FILL); 
1793     }
1794         // */
1795
1796         /*for(int ii=0; ii<(int)pow(2.0,mMaxRefine)-1; ii++) {
1797                 errMsg("BNDTESTSYMM","set "<<mLevel[mMaxRefine].lSizex-2-ii );
1798                 for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
1799                         for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
1800                                 initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, CFBnd, 0.0, BND_FILL);  // SYMM!? 2D?
1801                         }
1802                 for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
1803                         for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {      
1804                                 initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, CFBnd, 0.0, BND_FILL);   // SYMM!? 3D?
1805                         }
1806         }
1807         // Symmetry tests */
1808
1809         // prepare interface cells
1810         initFreeSurfaces();
1811         initStandingFluidGradient();
1812
1813         // perform first step to init initial mass
1814         mInitialMass = 0.0;
1815         int inmCellCnt = 0;
1816         FSGR_FORIJK1(mMaxRefine) {
1817                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid) ) {
1818                         LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0); 
1819                         FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
1820                         mInitialMass += fluidRho;
1821                         inmCellCnt ++;
1822                 } else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
1823                         mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
1824                         inmCellCnt ++;
1825                 }
1826         }
1827         mCurrentVolume = mCurrentMass = mInitialMass;
1828
1829         ParamVec cspv = D::mpParam->calculateCellSize();
1830         if(D::cDimension==2) cspv[2] = 1.0;
1831         inmCellCnt = 1;
1832         double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
1833         debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Initial Mass:"<<mInitialMass<<" normalized:"<<nrmMass, 3);
1834         mInitialMass = 0.0; // reset, and use actual value after first step
1835
1836         //mStartSymm = false;
1837 #if ELBEEM_BLENDER!=1
1838         if((D::cDimension==2)&&(D::mSizex<200)) {
1839                 if(!checkSymmetry("init")) {
1840                         errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
1841                 } else {
1842                         errMsg("LbmFsgrSolver::initialize","Symmetric init!");
1843                 }
1844         }
1845 #endif // ELBEEM_BLENDER!=1
1846         
1847
1848         // ----------------------------------------------------------------------
1849         // coarsen region
1850         myTime_t fsgrtstart = getTime(); 
1851         for(int lev=mMaxRefine-1; lev>=0; lev--) {
1852                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
1853                 performRefinement(lev);
1854                 performCoarsening(lev);
1855                 coarseRestrictFromFine(lev);
1856                 performRefinement(lev);
1857                 performCoarsening(lev);
1858                 coarseRestrictFromFine(lev);
1859         }
1860         D::markedClearList();
1861         myTime_t fsgrtend = getTime(); 
1862         if(!D::mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< ((fsgrtend-fsgrtstart)/(double)1000.0)<<"s), changes:"<<mNumFsgrChanges , 10 ); }
1863         mNumFsgrChanges = 0;
1864
1865         for(int l=0; l<D::cDirNum; l++) { 
1866                 LbmFloat area = 0.5 * 0.5 *0.5;
1867                 if(D::cDimension==2) area = 0.5 * 0.5;
1868
1869                 if(D::dfVecX[l]!=0) area *= 0.5;
1870                 if(D::dfVecY[l]!=0) area *= 0.5;
1871                 if(D::dfVecZ[l]!=0) area *= 0.5;
1872                 mFsgrCellArea[l] = area;
1873         } // l
1874         /*for(int lev=0; lev<mMaxRefine; lev++) {
1875         FSGR_FORIJK_BOUNDS(lev) {
1876                 if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
1877                         if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
1878                                 LbmFloat totArea = mFsgrCellArea[0]; // for l=0
1879                                 for(int l=1; l<D::cDirNum; l++) { 
1880                                         int ni=(2*i)+D::dfVecX[l], nj=(2*j)+D::dfVecY[l], nk=(2*k)+D::dfVecZ[l];
1881                                         if(RFLAG(lev+1, ni,nj,nk, mLevel[lev+1].setCurr)&
1882                                                         (CFGrFromCoarse|CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
1883                                                         //(CFUnused|CFEmpty) //? (CFBnd|CFEmpty|CFGrFromCoarse|CFUnused)
1884                                                         ) { 
1885                                                 //LbmFloat area = 0.25; if(D::dfVecX[l]!=0) area *= 0.5; if(D::dfVecY[l]!=0) area *= 0.5; if(D::dfVecZ[l]!=0) area *= 0.5;
1886                                                 totArea += mFsgrCellArea[l];
1887                                         }
1888                                 } // l
1889                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = totArea;
1890                         } else if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFEmpty) {
1891                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 1.0;
1892                         } else {
1893                                 QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) = 0.0;
1894                         }
1895                         errMsg("DFINI"," at l"<<lev<<" "<<PRINT_IJK<<" v:"<<QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) );
1896                 }
1897         } } // */
1898
1899         // now really done...
1900   debugOut("LbmFsgrSolver::initialize : Init done ...",10);
1901         D::mInitDone = 1;
1902
1903         // make sure both sets are ok
1904         // copy from other to curr
1905         for(int lev=0; lev<=mMaxRefine; lev++) {
1906         FSGR_FORIJK_BOUNDS(lev) {
1907                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1908         } } // first copy flags */
1909
1910
1911         
1912         if(mOutputSurfacePreview) {
1913                 if(D::cDimension==2) {
1914                         errFatal("LbmFsgrSolver::init","No preview in 2D allowed!",SIMWORLD_INITERROR); return false;
1915                 }
1916
1917                 //int previewSize = mOutputSurfacePreview;
1918                 // same as normal one, but use reduced size
1919                 mpPreviewSurface = new IsoSurface( D::mIsoValue, false );
1920                 mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
1921                 mpPreviewSurface->setIsolevel( D::mIsoValue );
1922                 // usually dont display for rendering
1923                 mpPreviewSurface->setVisible( false );
1924
1925                 mpPreviewSurface->setStart( vec2G(isostart) );
1926                 mpPreviewSurface->setEnd(   vec2G(isoend) );
1927                 LbmVec pisodist = isoend-isostart;
1928                 mpPreviewSurface->initializeIsosurface( (int)(mPreviewFactor*D::mSizex)+2, (int)(mPreviewFactor*D::mSizey)+2, (int)(mPreviewFactor*D::mSizez)+2, vec2G(pisodist) );
1929                 //mpPreviewSurface->setName( D::getName() + "preview" );
1930                 mpPreviewSurface->setName( "preview" );
1931         
1932                 debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(mPreviewFactor*D::mSizex)<<","<<(mPreviewFactor*D::mSizey)<<","<<(mPreviewFactor*D::mSizez)<<" enabled",10);
1933         }
1934
1935 #if ELBEEM_BLENDER==1
1936         // make sure fill fracs are right for first surface generation
1937         stepMain();
1938 #endif // ELBEEM_BLENDER==1
1939
1940         // prepare once...
1941         prepareVisualization();
1942         // copy again for stats counting
1943         for(int lev=0; lev<=mMaxRefine; lev++) {
1944         FSGR_FORIJK_BOUNDS(lev) {
1945                 RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
1946         } } // first copy flags */
1947
1948         /*{ int lev=mMaxRefine;
1949                 FSGR_FORIJK_BOUNDS(lev) { // COMPRT deb out
1950                 debMsgDirect("\n x="<<PRINT_IJK);
1951                 for(int l=0; l<D::cDfNum; l++) {
1952                         debMsgDirect(" df="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, l) );
1953                 }
1954                 debMsgDirect(" m="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, dMass) );
1955                 debMsgDirect(" f="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, dFfrac) );
1956                 debMsgDirect(" x="<< QCELL(lev, i,j,k,mLevel[lev].setCurr, dFlux) );
1957         } } // COMPRT ON */
1958         return true;
1959 }
1960
1961
1962
1963 /*****************************************************************************/
1964 /*! perform geometry init (if switched on) */
1965 /*****************************************************************************/
1966 template<class D>
1967 bool 
1968 LbmFsgrSolver<D>::initGeometryFlags() {
1969         int level = mMaxRefine;
1970         myTime_t geotimestart = getTime(); 
1971         ntlGeometryObject *pObj;
1972         // getCellSize (due to forced cubes, use x values)
1973         ntlVec3Gfx dvec( (D::mvGeoEnd[0]-D::mvGeoStart[0])/ ((LbmFloat)D::mSizex*2.0));
1974         // real cell size from now on...
1975         dvec *= 2.0; 
1976         ntlVec3Gfx nodesize = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
1977         dvec = nodesize;
1978         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< D::mGeoInitId <<") v"<<dvec,3);
1979
1980         /* set interface cells */
1981         D::initGeoTree(D::mGeoInitId);
1982         ntlVec3Gfx maxIniVel = vec2G( D::mpParam->calculateLattVelocityFromRw( vec2P(D::getGeoMaxInitialVelocity()) ));
1983         D::mpParam->setSimulationMaxSpeed( norm(maxIniVel) + norm(mLevel[level].gravity) );
1984         LbmFloat allowMax = D::mpParam->getTadapMaxSpeed();  // maximum allowed velocity
1985         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxIniVel <<", allowed Max="<<allowMax ,5);
1986         if(D::mpParam->getSimulationMaxSpeed() > allowMax) {
1987                 // similar to adaptTimestep();
1988                 LbmFloat nextmax = D::mpParam->getSimulationMaxSpeed();
1989                 LbmFloat newdt = D::mpParam->getStepTime() * (allowMax / nextmax); // newtr
1990                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< D::mpParam->getStepTime() <<" ",5);
1991                 D::mpParam->setDesiredStepTime( newdt );
1992                 D::mpParam->calculateAllMissingValues( D::mSilent );
1993                 maxIniVel = vec2G( D::mpParam->calculateLattVelocityFromRw( vec2P(D::getGeoMaxInitialVelocity()) ));
1994                 debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxIniVel,5);
1995         }
1996         recalculateObjectSpeeds();
1997
1998         ntlVec3Gfx pos,iniPos; // position of current cell
1999         LbmFloat rhomass = 0.0;
2000         int savedNodes = 0;
2001         int OId = -1;
2002         gfxReal distance;
2003
2004         // 2d display as rectangles
2005         if(D::cDimension==2) {
2006                 dvec[2] = 0.0; 
2007                 iniPos =(D::mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (D::mvGeoEnd[2]-D::mvGeoStart[2])*0.5 ))-(dvec*0.0);
2008                 //iniPos =(D::mvGeoStart + ntlVec3Gfx( 0.0 ))+dvec;
2009         } else {
2010                 iniPos =(D::mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
2011                 iniPos[2] = D::mvGeoStart[2] + dvec[2]*getForZMin1();
2012         }
2013
2014
2015         // first init boundary conditions
2016 #define GETPOS(i,j,k) \
2017                                                 ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
2018                                                 iniPos[1]+ dvec[1]*(gfxReal)(j), \
2019                                                 iniPos[2]+ dvec[2]*(gfxReal)(k) )
2020         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
2021                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
2022                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
2023                                 CellFlagType ntype = CFInvalid;
2024                                 if(D::geoInitCheckPointInside( GETPOS(i,j,k) , FGI_ALLBOUNDS, OId, distance)) {
2025                                         pObj = (*D::mpGiObjects)[OId];
2026                                         switch( pObj->getGeoInitType() ){
2027                                         case FGI_MBNDINFLOW:  
2028                                                 rhomass = 1.0;
2029                                                 ntype = CFFluid|CFMbndInflow; 
2030                                                 break;
2031                                         case FGI_MBNDOUTFLOW: 
2032                                                 rhomass = 0.0;
2033                                                 ntype = CFEmpty|CFMbndOutflow; 
2034                                                 break;
2035                                         default:
2036                                                 rhomass = BND_FILL;
2037                                                 ntype = CFBnd; break;
2038                                         }
2039                                 }
2040                                 if(ntype != CFInvalid) {
2041                                         // initDefaultCell
2042                                         if((ntype == CFMbndInflow) || (ntype == CFMbndOutflow) ) {
2043                                                 ntype |= (OId<<24);
2044                                         }
2045
2046                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
2047                                 }
2048
2049                                 // walk along x until hit for following inits
2050                                 if(distance<=-1.0) { distance = 100.0; }
2051                                 if(distance>0.0) {
2052                                         gfxReal dcnt=dvec[0];
2053                                         while(( dcnt< distance )&&(i+1<mLevel[level].lSizex-1)) {
2054                                                 dcnt += dvec[0]; i++;
2055                                                 savedNodes++;
2056                                                 if(ntype != CFInvalid) {
2057                                                         // rhomass are still inited from above
2058                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
2059                                                 }
2060                                         }
2061                                 } 
2062                                 // */
2063
2064                         } 
2065                 } 
2066         } // zmax
2067
2068
2069         // now init fluid layer
2070         for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
2071                 for(int j=1;j<mLevel[level].lSizey-1;j++) {
2072                         for(int i=1;i<mLevel[level].lSizex-1;i++) {
2073                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
2074
2075                                 CellFlagType ntype = CFInvalid;
2076                                 int inits = 0;
2077                                 if(D::geoInitCheckPointInside( GETPOS(i,j,k) , FGI_FLUID, OId, distance)) {
2078                                         ntype = CFFluid;
2079                                 }
2080                                 if(ntype != CFInvalid) {
2081                                         // initDefaultCell
2082                                         rhomass = 1.0;
2083                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
2084                                         inits++;
2085                                 }
2086
2087                                 // walk along x until hit for following inits
2088                                 if(distance<=-1.0) { distance = 100.0; }
2089                                 if(distance>0.0) {
2090                                         gfxReal dcnt=dvec[0];
2091                                         while((dcnt< distance )&&(i+1<mLevel[level].lSizex-1)) {
2092                                                 dcnt += dvec[0]; i++;
2093                                                 savedNodes++;
2094                                                 if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
2095                                                 if(ntype != CFInvalid) {
2096                                                         // rhomass are still inited from above
2097                                                         initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
2098                                                         inits++;
2099                                                 }
2100                                         }
2101                                 } // distance>0
2102                                 
2103                         } 
2104                 } 
2105         } // zmax
2106
2107         D::freeGeoTree();
2108         myTime_t geotimeend = getTime(); 
2109         debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< ((geotimeend-geotimestart)/(double)1000.0)<<"s,"<<savedNodes<<") " , 10 ); 
2110         //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
2111         return true;
2112 }
2113
2114 /*****************************************************************************/
2115 /* init part for all freesurface testcases */
2116 template<class D>
2117 void 
2118 LbmFsgrSolver<D>::initFreeSurfaces() {
2119         double interfaceFill = 0.45;   // filling level of interface cells
2120
2121         // set interface cells 
2122         FSGR_FORIJK1(mMaxRefine) {
2123
2124                 /*if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFEmpty )) {
2125                         int initInter = 0;
2126                         // check for neighboring fluid cells 
2127                         FORDF1 {
2128                                 if( TESTFLAG( RFLAG_NBINV(mMaxRefine, i, j, k,  mLevel[mMaxRefine].setCurr,l), CFFluid ) ) {
2129                                         initInter = 1;
2130                                 }
2131                         }
2132
2133                         if(initInter) {
2134                                 initEmptyCell(mMaxRefine, i,j,k, CFInter, 1.0, interfaceFill);
2135                         }
2136
2137                 } // empty cells  OLD */
2138
2139                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid )) {
2140                         int initInter = 0; // check for neighboring empty cells 
2141                         FORDF1 {
2142                                 if( TESTFLAG( RFLAG_NBINV(mMaxRefine, i, j, k,  mLevel[mMaxRefine].setCurr,l), CFEmpty ) ) {
2143                                         initInter = 1;
2144                                 }
2145                         }
2146                         if(initInter) {
2147                                 QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr, dMass) = 
2148                                         //QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setOther, dMass) =  // COMPRT OFF
2149                                         interfaceFill;
2150                                 RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) = RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setOther) = CFInter;
2151                         }
2152                 }
2153         }
2154
2155         // remove invalid interface cells 
2156         FSGR_FORIJK1(mMaxRefine) {
2157                 // remove invalid interface cells 
2158                 if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
2159                         int delit = 0;
2160                         int NBs = 0; // neighbor flags or'ed 
2161                         int noEmptyNB = 1;
2162
2163                         FORDF1 {
2164                                 if( TESTFLAG( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l ), CFEmpty ) ) {
2165                                         noEmptyNB = 0;
2166                                 }
2167                                 NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
2168                         }
2169
2170                         // remove cells with no fluid or interface neighbors
2171                         if((NBs & CFFluid)==0) { delit = 1; }
2172                         if((NBs & CFInter)==0) { delit = 1; }
2173
2174                         // remove cells with no empty neighbors
2175                         if(noEmptyNB) { delit = 2; }
2176
2177                         // now we can remove the cell 
2178                         if(delit==1) {
2179                                 initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
2180                         }
2181                         if(delit==2) {
2182                                 initEmptyCell(mMaxRefine, i,j,k, CFFluid, 1.0, 1.0);
2183                         }
2184                 } // interface 
2185         }
2186
2187         // another brute force init, make sure the fill values are right...
2188         // and make sure both sets are equal
2189         for(int lev=0; lev<=mMaxRefine; lev++) {
2190         FSGR_FORIJK_BOUNDS(lev) {
2191                 if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) { 
2192                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL;
2193                         continue;
2194                 }
2195                 if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) { 
2196                         QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0;
2197                         continue;
2198                 }
2199         } }
2200
2201         // ----------------------------------------------------------------------
2202         // smoother surface...
2203         if(mInitSurfaceSmoothing>0) {
2204                 debMsgStd("Surface Smoothing init", DM_MSG, "Performing "<<(mInitSurfaceSmoothing)<<" smoothing steps ",10);
2205 #if COMPRESSGRIDS==1
2206                 errFatal("NYI","COMPRESSGRIDS mInitSurfaceSmoothing",SIMWORLD_INITERROR); return;
2207 #endif // COMPRESSGRIDS==0
2208         }
2209         for(int s=0; s<mInitSurfaceSmoothing; s++) {
2210                 FSGR_FORIJK1(mMaxRefine) {
2211                         if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFInter) ) {
2212                                 LbmFloat mass = 0.0;
2213                                 //LbmFloat nbdiv;
2214                                 FORDF0 {
2215                                         int ni=i+D::dfVecX[l], nj=j+D::dfVecY[l], nk=k+D::dfVecZ[l];
2216                                         if( RFLAG(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr) & CFFluid ){
2217                                                 mass += 1.0;
2218                                         }
2219                                         if( RFLAG(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr) & CFInter ){
2220                                                 mass += QCELL(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, dMass);
2221                                         }
2222                                         //nbdiv+=1.0;
2223                                 }
2224
2225                                 //errMsg(" I ", PRINT_IJK<<" m"<<mass );
2226                                 QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setOther, dMass) = (mass/19.0);
2227                                 QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setOther, dFfrac) = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setOther, dMass);
2228                         }
2229                 }
2230
2231                 mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
2232                 mLevel[mMaxRefine].setCurr ^= 1;
2233         }
2234         // copy back...?
2235
2236 }
2237
2238 /*****************************************************************************/
2239 /* init part for all freesurface testcases */
2240 template<class D>
2241 void 
2242 LbmFsgrSolver<D>::initStandingFluidGradient() {
2243         // ----------------------------------------------------------------------
2244         // standing fluid preinit
2245         const int debugStandingPreinit = 0;
2246         int haveStandingFluid = 0;
2247
2248 #define STANDFLAGCHECK(iindex) \
2249                                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || \
2250                                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  \
2251                                         if((iindex)>1) { haveStandingFluid=(iindex); } \
2252                                         j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=D::getForZMaxBnd(); \
2253                                         continue; \
2254                                 } 
2255         int gravIndex[3] = {0,0,0};
2256         int gravDir[3] = {1,1,1};
2257         int maxGravComp = 1; // by default y
2258         int gravComp1 = 0; // by default y
2259         int gravComp2 = 2; // by default y
2260         if( ABS(mLevel[mMaxRefine].gravity[0]) > ABS(mLevel[mMaxRefine].gravity[1]) ){ maxGravComp = 0; gravComp1=1; gravComp2=2; }
2261         if( ABS(mLevel[mMaxRefine].gravity[2]) > ABS(mLevel[mMaxRefine].gravity[0]) ){ maxGravComp = 2; gravComp1=0; gravComp2=1; }
2262
2263         int gravIMin[3] = { 0 , 0 , 0 };
2264         int gravIMax[3] = {
2265                 mLevel[mMaxRefine].lSizex + 0,
2266                 mLevel[mMaxRefine].lSizey + 0,
2267                 mLevel[mMaxRefine].lSizez + 0 };
2268         if(LBMDIM==2) gravIMax[2] = 1;
2269
2270         //int gravDir = 1;
2271         if( mLevel[mMaxRefine].gravity[maxGravComp] > 0.0 ) {
2272                 // swap directions
2273                 int i=maxGravComp;
2274                 int tmp = gravIMin[i];
2275                 gravIMin[i] = gravIMax[i] - 1;
2276                 gravIMax[i] = tmp - 1;
2277                 gravDir[i] = -1;
2278         }
2279 #define PRINTGDIRS \
2280         errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \
2281         errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \
2282         errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] ); 
2283         // _PRINTGDIRS;
2284
2285         bool gravAbort = false;
2286 #define GRAVLOOP \
2287         gravAbort=false; \
2288         for(gravIndex[2]= gravIMin[2];     (gravIndex[2]!=gravIMax[2])&&(!gravAbort);  gravIndex[2] += gravDir[2]) \
2289                 for(gravIndex[1]= gravIMin[1];   (gravIndex[1]!=gravIMax[1])&&(!gravAbort);  gravIndex[1] += gravDir[1]) \
2290                         for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort);  gravIndex[0] += gravDir[0]) 
2291
2292         GRAVLOOP {
2293                 int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2294                 //if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) {debMsgStd("Standing fluid preinit", DM_MSG, "fluidheightinit check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 ); }
2295                 //STANDFLAGCHECK(gravIndex[maxGravComp]);
2296                 if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || 
2297                                 ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){  
2298                         int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
2299                         if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
2300                         //if(gravIndex[maxGravComp]>1)  
2301                         if(fluidHeight>1) 
2302                         {
2303                                 haveStandingFluid = fluidHeight; //gravIndex[maxGravComp]; 
2304                                 gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
2305                         }
2306                         gravAbort = true; continue; 
2307                 } 
2308         } // GRAVLOOP
2309         // _PRINTGDIRS;
2310
2311         LbmFloat fluidHeight;
2312         //if(gravDir>0) { fluidHeight = (LbmFloat)haveStandingFluid;
2313         //} else { fluidHeight = (LbmFloat)haveStandingFluid; }
2314         fluidHeight = (LbmFloat)(ABS(gravIMax[maxGravComp]-gravIMin[maxGravComp]));
2315         if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "fheight="<<fluidHeight<<" min="<<PRINT_VEC(gravIMin[0],gravIMin[1],       gravIMin[2])<<" max="<<PRINT_VEC(gravIMax[0], gravIMax[1],gravIMax[2])<<
2316                         " mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
2317                                 
2318         if(mDisableStandingFluidInit) {
2319                 debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
2320                 haveStandingFluid=0;
2321         }
2322
2323         // copy flags and init , as no flags will be changed during grav init
2324         // also important for corasening later on
2325         const int lev = mMaxRefine;
2326         CellFlagType nbflag[LBM_DFNUM], nbored; 
2327         for(int k=D::getForZMinBnd();k<D::getForZMaxBnd();++k) {
2328                 for(int j=0;j<mLevel[lev].lSizey-0;++j) {
2329                         for(int i=0;i<mLevel[lev].lSizex-0;++i) {
2330                                 if( (RFLAG(lev, i,j,k,SRCS(lev)) & (CFFluid)) ) {
2331                                         nbored = 0;
2332                                         FORDF1 {
2333                                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
2334                                                 nbored |= nbflag[l];
2335                                         } 
2336                                         if(nbored&CFBnd) {
2337                                                 RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
2338                                         } else {
2339                                                 RFLAG(lev, i,j,k,SRCS(lev)) |= CFNoBndFluid;
2340                                         }
2341                                 }
2342                                 RFLAG(lev, i,j,k,TSET(lev)) = RFLAG(lev, i,j,k,SRCS(lev));
2343         } } }
2344
2345         if(haveStandingFluid) {
2346                 int rhoworkSet = mLevel[lev].setCurr;
2347                 myTime_t timestart = getTime(); // FIXME use user time here?
2348 #if OPT3D==true 
2349                 LbmFloat lcsmqadd, lcsmqo, lcsmeq[LBM_DFNUM], lcsmomega;
2350 #endif // OPT3D==true 
2351
2352                 GRAVLOOP {
2353                         int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
2354                         //debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
2355                         if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
2356                                         ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){ 
2357                                 //gravAbort = true; 
2358                                 continue;
2359                         }
2360
2361                         LbmFloat rho = 1.0;
2362                         // 1/6 velocity from denisty gradient, 1/2 for delta of two cells
2363                         rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) * 
2364                                 (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega); 
2365                         if(debugStandingPreinit) 
2366                                 if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) { 
2367                                         errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK); 
2368                                 }
2369
2370                         if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
2371                                         (RFLAG(lev, i,j,k, rhoworkSet) & CFInter) ) {
2372                                 FORDF0 { QCELL(lev, i,j,k, rhoworkSet, l) *= rho; }
2373                                 QCELL(lev, i,j,k, rhoworkSet, dMass) *= rho;
2374                         }
2375
2376                 } // GRAVLOOP
2377                 debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited", 8);
2378                 
2379                 int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
2380                 preinitSteps = (haveStandingFluid>>2); // not much use...?
2381                 //preinitSteps = 4; // DEBUG!!!!
2382                 //D::mInitDone = 1; // GRAVTEST
2383                 //preinitSteps = 0;
2384                 debMsgNnl("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
2385                 for(int s=0; s<preinitSteps; s++) {
2386                         int workSet = SRCS(lev); //mLevel[lev].setCurr;
2387                         int otherSet = TSET(lev); //mLevel[lev].setOther;
2388                         debMsgDirect(".");
2389                         if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "s="<<s<<" curset="<<workSet<<" srcs"<<SRCS(lev), 10);
2390                         LbmFloat *ccel;
2391                         LbmFloat *tcel;
2392                         LbmFloat m[LBM_DFNUM];
2393
2394                 // grav loop not necessary here
2395 #define NBFLAG(l) (nbflag[(l)])
2396                 LbmFloat rho, ux,uy,uz, usqr; 
2397                 int kstart=D::getForZMinBnd(), kend=D::getForZMaxBnd();
2398 #if COMPRESSGRIDS==0
2399                 for(int k=kstart;k<kend;++k) {
2400 #else // COMPRESSGRIDS==0
2401                 int kdir = 1; // COMPRT ON
2402                 if(mLevel[lev].setCurr==1) {
2403                         kdir = -1;
2404                         int temp = kend;
2405                         kend = kstart-1;
2406                         kstart = temp-1;
2407                 } // COMPRT
2408                 for(int k=kstart;k!=kend;k+=kdir) {
2409
2410                 //errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir ); // debug
2411 #endif // COMPRESSGRIDS==0
2412
2413                 for(int j=0;j<mLevel[lev].lSizey-0;++j) {
2414                 for(int i=0;i<mLevel[lev].lSizex-0;++i) {
2415                                 const CellFlagType currFlag = RFLAG(lev, i,j,k,workSet);
2416                                 if( (currFlag & (CFEmpty|CFBnd)) ) continue;
2417                                 ccel = RACPNT(lev, i,j,k,workSet); 
2418                                 tcel = RACPNT(lev, i,j,k,otherSet);
2419
2420                                 if( (currFlag & (CFInter)) ) {
2421                                         // copy all values
2422                                         for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
2423                                         continue;
2424                                 }
2425
2426                                 if( (currFlag & CFNoBndFluid)) {
2427                                         OPTIMIZED_STREAMCOLLIDE;
2428                                 } else {
2429                                         FORDF1 {
2430                                                 nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
2431                                         } 
2432                                         DEFAULT_STREAM;
2433                                         ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; 
2434                                         DEFAULT_COLLIDE;
2435                                 }
2436                                 for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
2437                         } } } // GRAVLOOP
2438
2439                         mLevel[lev].setOther = mLevel[lev].setCurr;
2440                         mLevel[lev].setCurr ^= 1;
2441                 }
2442                 //D::mInitDone = 0; // GRAVTEST
2443                 // */
2444
2445                 myTime_t timeend = getTime();
2446                 debMsgDirect(" done, "<<((timeend-timestart)/(double)1000.0)<<"s \n");
2447 #undef  NBFLAG
2448         }
2449 }
2450
2451
2452
2453 /*****************************************************************************/
2454 /* init a given cell with flag, density, mass and equilibrium dist. funcs */
2455
2456 template<class D>
2457 void LbmFsgrSolver<D>::changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
2458         CellFlagType pers = RFLAG(level,xx,yy,zz,set) & CFPersistMask;
2459         RFLAG(level,xx,yy,zz,set) = newflag | pers;
2460 }
2461
2462 template<class D>
2463 void 
2464 LbmFsgrSolver<D>::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass) {
2465   /* init eq. dist funcs */
2466         LbmFloat *ecel;
2467         int workSet = mLevel[level].setCurr;
2468
2469         ecel = RACPNT(level, i,j,k, workSet);
2470         FORDF0 { RAC(ecel, l) = D::dfEquil[l] * rho; }
2471         RAC(ecel, dMass) = mass;
2472         RAC(ecel, dFfrac) = mass/rho;
2473         RAC(ecel, dFlux) = FLUX_INIT;
2474         //RFLAG(level, i,j,k, workSet)= flag;
2475         changeFlag(level, i,j,k, workSet, flag);
2476
2477   workSet ^= 1;
2478         changeFlag(level, i,j,k, workSet, flag);
2479         return;
2480 }
2481
2482 template<class D>
2483 void 
2484 LbmFsgrSolver<D>::initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel) {
2485         LbmFloat *ecel;
2486         int workSet = mLevel[level].setCurr;
2487
2488         ecel = RACPNT(level, i,j,k, workSet);
2489         FORDF0 { RAC(ecel, l) = D::getCollideEq(l, rho,vel[0],vel[1],vel[2]); }
2490         RAC(ecel, dMass) = mass;
2491         RAC(ecel, dFfrac) = mass/rho;
2492         RAC(ecel, dFlux) = FLUX_INIT;
2493         //RFLAG(level, i,j,k, workSet) = flag;
2494         changeFlag(level, i,j,k, workSet, flag);
2495
2496   workSet ^= 1;
2497         changeFlag(level, i,j,k, workSet, flag);
2498         return;
2499 }
2500
2501 template<class D>
2502 bool 
2503 LbmFsgrSolver<D>::checkSymmetry(string idstring)
2504 {
2505         bool erro = false;
2506         bool symm = true;
2507         int msgs = 0;
2508         const int maxMsgs = 10;
2509         const bool markCells = false;
2510
2511         //for(int lev=0; lev<=mMaxRefine; lev++) {
2512         { int lev = mMaxRefine;
2513
2514         // no point if not symm.
2515         if( (mLevel[lev].lSizex==mLevel[lev].lSizey) && (mLevel[lev].lSizex==mLevel[lev].lSizez)) {
2516                 // ok
2517         } else {
2518                 return false;
2519         }
2520
2521         for(int s=0; s<2; s++) {
2522         FSGR_FORIJK1(lev) {
2523                 if(i<(mLevel[lev].lSizex/2)) {
2524                         int inb = (mLevel[lev].lSizey-1-i); 
2525
2526                         if(lev==mMaxRefine) inb -= 1;           // FSGR_SYMM_T
2527
2528                         if( RFLAG(lev, i,j,k,s) != RFLAG(lev, inb,j,k,s) ) { erro = true;
2529                                 if(D::cDimension==2) {
2530                                         if(msgs<maxMsgs) { msgs++;
2531                                                 errMsg("EFLAG", PRINT_IJK<<"s"<<s<<" flag "<<RFLAG(lev, i,j,k,s)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" flag "<<RFLAG(lev, inb,j,k,s) );
2532                                         }
2533                                 }
2534                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2535                                 symm = false;
2536                         }
2537                         if( LBM_FLOATNEQ(QCELL(lev, i,j,k,s, dMass), QCELL(lev, inb,j,k,s, dMass)) ) { erro = true;
2538                                 if(D::cDimension==2) {
2539                                         if(msgs<maxMsgs) { msgs++;
2540                                                 //debMsgDirect(" mass1 "<<QCELL(lev, i,j,k,s, dMass)<<" mass2 "<<QCELL(lev, inb,j,k,s, dMass) <<std::endl);
2541                                                 errMsg("EMASS", PRINT_IJK<<"s"<<s<<" mass "<<QCELL(lev, i,j,k,s, dMass)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" mass "<<QCELL(lev, inb,j,k,s, dMass) );
2542                                         }
2543                                 }
2544                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2545                                 symm = false;
2546                         }
2547
2548                         LbmFloat nbrho = QCELL(lev, i,j,k, s, dC);
2549                         FORDF1 { nbrho += QCELL(lev, i,j,k, s, l); }
2550                         LbmFloat otrho = QCELL(lev, inb,j,k, s, dC);
2551                         FORDF1 { otrho += QCELL(lev, inb,j,k, s, l); }
2552                         if( LBM_FLOATNEQ(nbrho, otrho) ) { erro = true;
2553                                 if(D::cDimension==2) {
2554                                         if(msgs<maxMsgs) { msgs++;
2555                                                 //debMsgDirect(" rho 1 "<<nbrho <<" rho 2 "<<otrho  <<std::endl);
2556                                                 errMsg("ERHO ", PRINT_IJK<<"s"<<s<<" rho  "<<nbrho <<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" rho  "<<otrho  );
2557                                         }
2558                                 }
2559                                 if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
2560                                 symm = false;
2561                         }
2562                 }
2563         } }
2564         } // lev
2565         LbmFloat maxdiv =0.0;
2566         if(erro) {
2567                 errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
2568                 //if(D::cDimension==2) D::mPanic = true; 
2569                 //return false;
2570         } else {
2571                 errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
2572         }
2573         // all ok...
2574         return symm;
2575 }// */
2576
2577
2578 /*****************************************************************************/
2579 /*! debug object display */
2580 /*****************************************************************************/
2581 template<class D>
2582 vector<ntlGeometryObject*> LbmFsgrSolver<D>::getDebugObjects() { 
2583         vector<ntlGeometryObject*> debo; 
2584         if(mOutputSurfacePreview) {
2585                 debo.push_back( mpPreviewSurface );
2586         }
2587         return debo; 
2588 }
2589
2590 /*****************************************************************************/
2591 /*! perform a single LBM step */
2592 /*****************************************************************************/
2593
2594 template<class D>
2595 void 
2596 LbmFsgrSolver<D>::stepMain()
2597 {
2598 #if ELBEEM_BLENDER==1
2599                 // update gui display
2600                 SDL_mutexP(globalBakeLock);
2601                 if(globalBakeState<0) {
2602                         // this means abort... cause panic
2603                         D::mPanic = 1;
2604                         errMsg("LbmFsgrSolver::step","Got abort signal from GUI, causing panic, aborting...");
2605                 }
2606                 SDL_mutexV(globalBakeLock);
2607 #endif // ELBEEM_BLENDER==1
2608         D::markedClearList(); // DMC clearMarkedCellsList
2609         if(mDropping) {
2610                 initDrop(mDropX, mDropY);
2611         }
2612         if(mGfxGeoSetup==6) {
2613                 // xobs init hack
2614                 if(mSimulationTime<0.400) {
2615                         if((mSimulationTime>0.25) && (mSimulationTime<0.325)) {
2616                                 // stop shortly...
2617                                 mDropping = false;
2618                         } else {
2619                                 initDrop(0.0, 1.0);
2620                         }
2621                 } else {
2622                         mDropping=false;
2623                 }
2624         }
2625
2626         // safety check, counter reset
2627         D::mNumUsedCells = 0;
2628         mNumInterdCells = 0;
2629         mNumInvIfCells = 0;
2630
2631   //debugOutNnl("LbmFsgrSolver::step : "<<D::mStepCnt, 10);
2632   if(!D::mSilent){ debMsgNnl("LbmFsgrSolver::step", DM_MSG, D::mName<<" cnt:"<<D::mStepCnt<<"  ", 10); }
2633         //debMsgDirect(  "LbmFsgrSolver::step : "<<D::mStepCnt<<" ");
2634         myTime_t timestart = getTime();
2635         //myTime_t timestart = 0;
2636         //if(mStartSymm) { checkSymmetry("step1"); } // DEBUG 
2637
2638         // important - keep for tadap
2639         mCurrentMass = D::mFixMass; // reset here for next step
2640         mCurrentVolume = 0.0;
2641         
2642         //stats
2643         mMaxVlen = mMxvz = mMxvy = mMxvx = 0.0;
2644
2645         //change to single step advance!
2646         int levsteps = 0;
2647         int dsbits = D::mStepCnt ^ (D::mStepCnt-1);
2648         //errMsg("S"," step:"<<D::mStepCnt<<" s-1:"<<(D::mStepCnt-1)<<" xf:"<<convertCellFlagType2String(dsbits));
2649         for(int lev=0; lev<=mMaxRefine; lev++) {
2650                 //if(! (D::mStepCnt&(1<<lev)) ) {
2651                 if( dsbits & (1<<(mMaxRefine-lev)) ) {
2652                         //errMsg("S"," l"<<lev);
2653
2654                         if(lev==mMaxRefine) {
2655                                 // always advance fine level...
2656                                 fineAdvance();
2657                         } else {
2658                                 performRefinement(lev);
2659                                 performCoarsening(lev);
2660                                 coarseRestrictFromFine(lev);
2661                                 coarseAdvance(lev);
2662                         }
2663 #if FSGR_OMEGA_DEBUG==1
2664                         errMsg("LbmFsgrSolver::step","LES stats l="<<lev<<" omega="<<mLevel[lev].omega<<" avgOmega="<< (mLevel[lev].avgOmega/mLevel[lev].avgOmegaCnt) );
2665                         mLevel[lev].avgOmega = 0.0; mLevel[lev].avgOmegaCnt = 0.0;
2666 #endif // FSGR_OMEGA_DEBUG==1
2667                         levsteps++;
2668                 }
2669                 mCurrentMass   += mLevel[lev].lmass;
2670                 mCurrentVolume += mLevel[lev].lvolume;
2671         }
2672
2673   // prepare next step
2674         D::mStepCnt++;
2675
2676
2677         // some dbugging output follows
2678         // calculate MLSUPS
2679         myTime_t timeend = getTime();
2680
2681         D::mNumUsedCells += mNumInterdCells; // count both types for MLSUPS
2682         mAvgNumUsedCells += D::mNumUsedCells;
2683         D::mMLSUPS = (D::mNumUsedCells / ((timeend-timestart)/(double)1000.0) ) / (1000000);
2684         if(D::mMLSUPS>10000){ D::mMLSUPS = -1; }
2685         else { mAvgMLSUPS += D::mMLSUPS; mAvgMLSUPSCnt += 1.0; } // track average mlsups
2686         
2687         LbmFloat totMLSUPS = ( ((mLevel[mMaxRefine].lSizex-2)*(mLevel[mMaxRefine].lSizey-2)*(getForZMax1(mMaxRefine)-getForZMin1())) / ((timeend-timestart)/(double)1000.0) ) / (1000000);
2688         if(totMLSUPS>10000) totMLSUPS = -1;
2689         mNumInvIfTotal += mNumInvIfCells; // debug
2690
2691   // do some formatting 
2692   if(!D::mSilent){ 
2693                 string sepStr(""); // DEBUG
2694                 debMsgDirect( 
2695                         "mlsups(curr:"<<D::mMLSUPS<<
2696                         " avg:"<<(mAvgMLSUPS/mAvgMLSUPSCnt)<<"), "<< sepStr<<
2697                         " totcls:"<<(D::mNumUsedCells)<< sepStr<<
2698                         " avgcls:"<< (int)(mAvgNumUsedCells/(long long int)D::mStepCnt)<< sepStr<<
2699                         " intd:"<<mNumInterdCells<< sepStr<<
2700                         " invif:"<<mNumInvIfCells<< sepStr<<
2701                         " invift:"<<mNumInvIfTotal<< sepStr<<
2702                         " fsgrcs:"<<mNumFsgrChanges<< sepStr<<
2703                         " filled:"<<D::mNumFilledCells<<", emptied:"<<D::mNumEmptiedCells<< sepStr<<
2704                         " mMxv:"<<mMxvx<<","<<mMxvy<<","<<mMxvz<<", tscnts:"<<mTimeSwitchCounts<< sepStr<<
2705                         /*" rhoMax:"<<mRhoMax<<", rhoMin:"<<mRhoMin<<", vlenMax:"<<mMaxVlen<<", "*/
2706                         " probs:"<<mNumProblems<< sepStr<<
2707                         " simt:"<<mSimulationTime<< sepStr<<
2708                         " for '"<<D::mName<<"' " );
2709
2710                 //wrong?
2711                 //debMsgDirect(", dccd:"<< mCurrentMass<<"/"<<mCurrentVolume<<"(fix:"<<D::mFixMass<<",ini:"<<mInitialMass<<") ");
2712                 debMsgDirect(std::endl);
2713                 debMsgDirect(D::mStepCnt<<": dccd="<< mCurrentMass<<"/"<<mCurrentVolume<<"(fix="<<D::mFixMass<<",ini="<<mInitialMass<<") ");
2714                 debMsgDirect(std::endl);
2715
2716                 // nicer output
2717                 debMsgDirect(std::endl); // 
2718                 //debMsgStd(" ",DM_MSG," ",10);
2719         } else {
2720                 debMsgDirect(".");
2721                 //if((mStepCnt%10)==9) debMsgDirect("\n");
2722         }
2723
2724         if(D::mStepCnt==1) {
2725                 mMinNoCells = mMaxNoCells = D::mNumUsedCells;
2726         } else {
2727                 if(D::mNumUsedCells>mMaxNoCells) mMaxNoCells = D::mNumUsedCells;
2728                 if(D::mNumUsedCells<mMinNoCells) mMinNoCells = D::mNumUsedCells;
2729         }
2730         
2731         // mass scale test
2732         if((mMaxRefine>0)&&(mInitialMass>0.0)) {
2733                 LbmFloat mscale = mInitialMass/mCurrentMass;
2734
2735                 mscale = 1.0;
2736                 const LbmFloat dchh = 0.001;
2737                 if(mCurrentMass<mInitialMass) mscale = 1.0+dchh;
2738                 if(mCurrentMass>mInitialMass) mscale = 1.0-dchh;
2739
2740                 // use mass rescaling?
2741                 // with float precision this seems to be nonsense...
2742                 const bool MREnable = false;
2743
2744                 const int MSInter = 2;
2745                 static int mscount = 0;
2746                 if( (MREnable) && ((mLevel[0].lsteps%MSInter)== (MSInter-1)) && ( ABS( (mInitialMass/mCurrentMass)-1.0 ) > 0.01) && ( dsbits & (1<<(mMaxRefine-0)) ) ){
2747                         // example: FORCE RESCALE MASS! ini:1843.5, cur:1817.6, f=1.01425 step:22153 levstep:5539 msc:37
2748                         // mass rescale MASS RESCALE check
2749                         errMsg("MDTDD","\n\n");
2750                         errMsg("MDTDD","FORCE RESCALE MASS! "
2751                                         <<"ini:"<<mInitialMass<<", cur:"<<mCurrentMass<<", f="<<ABS(mInitialMass/mCurrentMass)
2752                                         <<" step:"<<D::mStepCnt<<" levstep:"<<mLevel[0].lsteps<<" msc:"<<mscount<<" "
2753                                         );
2754                         errMsg("MDTDD","\n\n");
2755
2756                         mscount++;
2757                         for(int lev=mMaxRefine; lev>=0 ; lev--) {
2758                                 //for(int workSet = 0; workSet<=1; workSet++) {
2759                                 int wss = 0;
2760                                 int wse = 1;
2761 #if COMPRESSGRIDS==1
2762                                 if(lev== mMaxRefine) wss = wse = mLevel[lev].setCurr;
2763 #endif // COMPRESSGRIDS==1
2764                                 for(int workSet = wss; workSet<=wse; workSet++) { // COMPRT
2765
2766                                         FSGR_FORIJK1(lev) {
2767                                                 if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) 
2768                                                         ) {
2769
2770                                                         FORDF0 { QCELL(lev, i,j,k,workSet, l) *= mscale; }
2771                                                         QCELL(lev, i,j,k,workSet, dMass) *= mscale;
2772                                                         QCELL(lev, i,j,k,workSet, dFfrac) *= mscale;
2773
2774                                                 } else {
2775                                                         continue;
2776                                                 }
2777                                         }
2778                                 }
2779                                 mLevel[lev].lmass *= mscale;
2780                         }
2781                 } 
2782
2783                 mCurrentMass *= mscale;
2784         }// if mass scale test */
2785         else {
2786                 // use current mass after full step for initial setting
2787                 if((mMaxRefine>0)&&(mInitialMass<=0.0) && (levsteps == (mMaxRefine+1))) {
2788                         mInitialMass = mCurrentMass;
2789                         debMsgStd("MDTDD",DM_NOTIFY,"Second Initial Mass Init: "<<mInitialMass, 2);
2790                 }
2791         }
2792
2793         // one of the last things to do - adapt timestep
2794         // was in fineAdvance before... 
2795         if(mTimeAdap) {
2796                 adaptTimestep();
2797         } // time adaptivity
2798
2799         // debug - raw dump of ffrac values
2800         /*if((mStepCnt%200)==1){
2801                 std::ostringstream name;
2802                 name <<"fill_" << mStepCnt <<".dump";
2803                 FILE *file = fopen(name.str().c_str(),"w");
2804                 for(int k= getForZMinBnd(mMaxRefine); k< getForZMaxBnd(mMaxRefine); ++k) 
2805                  for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++) 
2806                         for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
2807                                 float val = QCELL(mMaxRefine,i,j,k, mLevel[mMaxRefine].setCurr,dFfrac);
2808                                 fwrite( &val, sizeof(val), 1, file);
2809                                 //errMsg("W", PRINT_IJK<<" val:"<<val);
2810                         }
2811                 fclose(file);
2812         } // */
2813
2814         /*
2815         if(1) { // DEBUG
2816                 const int lev = mMaxRefine;
2817                 int workSet = mLevel[lev].setCurr;
2818                 FSGR_FORIJK1(lev) {
2819                         if( 
2820                                         (RFLAG(lev,i,j,k, workSet) & CFFluid) || 
2821                                         (RFLAG(lev,i,j,k, workSet) & CFInter) ||
2822                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromCoarse) || 
2823                                         (RFLAG(lev,i,j,k, workSet) & CFGrFromFine) || 
2824                                         (RFLAG(lev,i,j,k, workSet) & CFGrNorm) 
2825                                         ) {
2826                                 // these cells have to be scaled...
2827                         } else {
2828                                 continue;
2829                         }
2830
2831                         // collide on current set
2832                         LbmFloat rho, ux,uy,uz;
2833                         rho=0.0; ux =  uy = uz = 0.0;
2834                         for(int l=0; l<D::cDfNum; l++) {
2835                                 LbmFloat m = QCELL(lev, i, j, k, workSet, l); 
2836                                 rho += m;
2837                                 ux  += (D::dfDvecX[l]*m);
2838                                 uy  += (D::dfDvecY[l]*m); 
2839                                 uz  += (D::dfDvecZ[l]*m); 
2840                         } 
2841                         //errMsg("DEBUG"," "<<PRINT_IJK <<" rho="<<rho<<" vel="<<PRINT_VEC(ux,uy,uz) );
2842                         errMsg(stdout,"D %d,%d rho=%+'.5f vel=[%+'.5f,%+'.5f] \n", i,j, rho, ux,uy );
2843                 }
2844         }       // DEBUG */
2845
2846 }
2847
2848 template<class D>
2849 void 
2850 LbmFsgrSolver<D>::fineAdvance()
2851 {
2852         // do the real thing...
2853         mainLoop( mMaxRefine );
2854         if(mUpdateFVHeight) {
2855                 // warning assume -Y gravity...
2856                 mFVHeight = mCurrentMass*mFVArea/((LbmFloat)(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizez));
2857                 if(mFVHeight<1.0) mFVHeight = 1.0;
2858                 D::mpParam->setFluidVolumeHeight(mFVHeight);
2859         }
2860
2861         // advance time before timestep change
2862         mSimulationTime += D::mpParam->getStepTime();
2863         // time adaptivity
2864         D::mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
2865         //if(mStartSymm) { checkSymmetry("step2"); } // DEBUG 
2866         if(!D::mSilent){ errMsg("fineAdvance"," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps) ); }
2867
2868         // update other set
2869   mLevel[mMaxRefine].setOther   = mLevel[mMaxRefine].setCurr;
2870   mLevel[mMaxRefine].setCurr   ^= 1;
2871   mLevel[mMaxRefine].lsteps++;
2872
2873         // flag init... (work on current set, to simplify flag checks)
2874         reinitFlags( mLevel[mMaxRefine].setCurr );
2875         //if((D::cDimension==2)&&(mStartSymm)) { checkSymmetry("step3"); } // DEBUG 
2876         if(!D::mSilent){ errMsg("fineAdvance"," flags reinit on set "<< mLevel[mMaxRefine].setCurr ); }
2877 }
2878
2879
2880 /*****************************************************************************/
2881 //! coarse/fine step functions
2882 /*****************************************************************************/
2883
2884 // access to own dfs during step (may be changed to local array)
2885 #define MYDF(l) RAC(ccel, l)
2886
2887 template<class D>
2888 void 
2889 LbmFsgrSolver<D>::mainLoop(int lev)
2890 {
2891         // loops over _only inner_ cells  -----------------------------------------------------------------------------------
2892         LbmFloat calcCurrentMass = 0.0; //mCurrentMass;
2893         LbmFloat calcCurrentVolume = 0.0; //mCurrentVolume;
2894         int      calcCellsFilled = D::mNumFilledCells;
2895         int      calcCellsEmptied = D::mNumEmptiedCells;
2896         int      calcNumUsedCells = D::mNumUsedCells;
2897
2898 #if PARALLEL==1
2899 #include "paraloop.h"
2900 #else // PARALLEL==1
2901   { // main loop region
2902         int kstart=D::getForZMin1(), kend=D::getForZMax1();
2903 #define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
2904 #endif // PARALLEL==1
2905
2906
2907         // local to loop
2908         CellFlagType nbflag[LBM_DFNUM]; 
2909 #define NBFLAG(l) nbflag[(l)]
2910         // */
2911         
2912         LbmFloat *ccel = NULL;
2913         LbmFloat *tcel = NULL;
2914         int oldFlag;
2915         int newFlag;
2916         int nbored;
2917         LbmFloat m[LBM_DFNUM];
2918         LbmFloat rho, ux, uy, uz, tmp, usqr;
2919         LbmFloat mass, change;
2920         usqr = tmp = 0.0;