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