merge from trunk #37722
[blender-staging.git] / intern / elbeem / intern / solver_class.h
1 /** \file elbeem/intern/solver_class.h
2  *  \ingroup elbeem
3  */
4 /******************************************************************************
5  *
6  * El'Beem - the visual lattice boltzmann freesurface simulator
7  * All code distributed as part of El'Beem is covered by the version 2 of the 
8  * GNU General Public License. See the file COPYING for details.
9  * Copyright 2003-2006 Nils Thuerey
10  *
11  * Combined 2D/3D Lattice Boltzmann standard solver classes
12  *
13  *****************************************************************************/
14
15
16 #ifndef LBM_SOLVERCLASS_H
17 #define LBM_SOLVERCLASS_H
18
19 #include "utilities.h"
20 #include "solver_interface.h"
21 #include "ntl_ray.h"
22 #include <stdio.h>
23
24 #if PARALLEL==1
25 #include <omp.h>
26 #endif // PARALLEL=1
27 #ifndef PARALLEL
28 #define PARALLEL 0
29 #endif // PARALLEL
30
31
32 // general solver setting defines
33
34 //! debug coordinate accesses and the like? (much slower)
35 //  might be enabled by compilation
36 #ifndef FSGR_STRICT_DEBUG
37 #define FSGR_STRICT_DEBUG 0
38 #endif // FSGR_STRICT_DEBUG
39
40 //! debug coordinate accesses and the like? (much slower)
41 #define FSGR_OMEGA_DEBUG 0
42
43 //! OPT3D quick LES on/off, only debug/benchmarking
44 #define USE_LES 1
45
46 //! order of interpolation (0=always current/1=interpolate/2=always other)
47 //#define TIMEINTORDER 0
48 // TODO remove interpol t param, also interTime
49
50 // use optimized 3D code?
51 #if LBMDIM==2
52 #define OPT3D 0
53 #else
54 // determine with debugging...
55 #       if FSGR_STRICT_DEBUG==1
56 #               define OPT3D 0
57 #       else // FSGR_STRICT_DEBUG==1
58 // usually switch optimizations for 3d on, when not debugging
59 #               define OPT3D 1
60 // COMPRT
61 //#             define OPT3D 0
62 #       endif // FSGR_STRICT_DEBUG==1
63 #endif
64
65 //! invalid mass value for unused mass data
66 #define MASS_INVALID -1000.0
67
68 // empty/fill cells without fluid/empty NB's by inserting them into the full/empty lists?
69 #define FSGR_LISTTRICK          1
70 #define FSGR_LISTTTHRESHEMPTY   0.10
71 #define FSGR_LISTTTHRESHFULL    0.90
72 #define FSGR_MAGICNR            0.025
73 //0.04
74
75 //! maxmimum no. of grid levels
76 #define FSGR_MAXNOOFLEVELS 5
77
78 // enable/disable fine grid compression for finest level
79 // make sure this is same as useGridComp in calculateMemreqEstimate
80 #if LBMDIM==3
81 #define COMPRESSGRIDS 1
82 #else 
83 #define COMPRESSGRIDS 0
84 #endif 
85
86 // helper for comparing floats with epsilon
87 #define GFX_FLOATNEQ(x,y) ( ABS((x)-(y)) > (VECTOR_EPSILON) )
88 #define LBM_FLOATNEQ(x,y) ( ABS((x)-(y)) > (10.0*LBM_EPSILON) )
89
90
91 // macros for loops over all DFs
92 #define FORDF0 for(int l= 0; l< LBM_DFNUM; ++l)
93 #define FORDF1 for(int l= 1; l< LBM_DFNUM; ++l)
94 // and with different loop var to prevent shadowing
95 #define FORDF0M for(int m= 0; m< LBM_DFNUM; ++m)
96 #define FORDF1M for(int m= 1; m< LBM_DFNUM; ++m)
97
98 // iso value defines
99 // border for marching cubes
100 #define ISOCORR 3
101
102 #define LBM_INLINED  inline
103
104 // sirdude fix for solaris
105 #if !defined(linux) && defined(sun)
106 #include "ieeefp.h"
107 #ifndef expf
108 #define expf(x) exp((double)(x))
109 #endif
110 #endif
111
112 #include "solver_control.h"
113
114 #if LBM_INCLUDE_TESTSOLVERS==1
115 #include "solver_test.h"
116 #endif // LBM_INCLUDE_TESTSOLVERS==1
117
118 /*****************************************************************************/
119 /*! cell access classes */
120 class UniformFsgrCellIdentifier : 
121         public CellIdentifierInterface , public LbmCellContents
122 {
123         public:
124                 //! which grid level?
125                 int level;
126                 //! location in grid
127                 int x,y,z;
128
129                 //! reset constructor
130                 UniformFsgrCellIdentifier() :
131                         x(0), y(0), z(0) { };
132
133                 // implement CellIdentifierInterface
134                 virtual string getAsString() {
135                         std::ostringstream ret;
136                         ret <<"{ i"<<x<<",j"<<y;
137                         if(LBMDIM>2) ret<<",k"<<z;
138                         ret <<" }";
139                         return ret.str();
140                 }
141
142                 virtual bool equal(CellIdentifierInterface* other) {
143                         UniformFsgrCellIdentifier *cid = (UniformFsgrCellIdentifier *)( other );
144                         if(!cid) return false;
145                         if( x==cid->x && y==cid->y && z==cid->z && level==cid->level ) return true;
146                         return false;
147                 }
148 };
149
150 //! information needed for each level in the simulation
151 class FsgrLevelData {
152 public:
153         int id; // level number
154
155         //! node size on this level (geometric, in world coordinates, not simulation units!) 
156         LbmFloat nodeSize;
157         //! node size on this level in simulation units 
158         LbmFloat simCellSize;
159         //! quadtree node relaxation parameter 
160         LbmFloat omega;
161         //! size this level was advanced to 
162         LbmFloat time;
163         //! size of a single lbm step in time units on this level 
164         LbmFloat timestep;
165         //! step count
166         int lsteps;
167         //! gravity force for this level
168         LbmVec gravity;
169         //! level array 
170         LbmFloat *mprsCells[2];
171         CellFlagType *mprsFlags[2];
172
173         //! smago params and precalculated values
174         LbmFloat lcsmago;
175         LbmFloat lcsmago_sqr;
176         LbmFloat lcnu;
177
178         // LES statistics per level
179         double avgOmega;
180         double avgOmegaCnt;
181
182         //! current set of dist funcs 
183         int setCurr;
184         //! target/other set of dist funcs 
185         int setOther;
186
187         //! mass&volume for this level
188         LbmFloat lmass;
189         LbmFloat lvolume;
190         LbmFloat lcellfactor;
191
192         //! local storage of mSizes
193         int lSizex, lSizey, lSizez;
194         int lOffsx, lOffsy, lOffsz;
195
196 };
197
198
199
200 /*****************************************************************************/
201 /*! class for solving a LBM problem */
202 class LbmFsgrSolver : 
203         public LbmSolverInterface // this means, the solver is a lbmData object and implements the lbmInterface
204 {
205
206         public:
207                 //! Constructor 
208                 LbmFsgrSolver();
209                 //! Destructor 
210                 virtual ~LbmFsgrSolver();
211
212                 //! initilize variables fom attribute list 
213                 virtual void parseAttrList();
214                 //! Initialize omegas and forces on all levels (for init/timestep change)
215                 void initLevelOmegas();
216
217                 // multi step solver init
218                 /*! finish the init with config file values (allocate arrays...) */
219                 virtual bool initializeSolverMemory();
220                 /*! init solver arrays */
221                 virtual bool initializeSolverGrids();
222                 /*! prepare actual simulation start, setup viz etc */
223                 virtual bool initializeSolverPostinit();
224
225                 //! notify object that dump is in progress (e.g. for field dump) 
226                 virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename);
227
228 #               if LBM_USE_GUI==1
229                 //! show simulation info (implement LbmSolverInterface pure virtual func)
230                 virtual void debugDisplay(int set);
231 #               endif
232                 
233                 // implement CellIterator<UniformFsgrCellIdentifier> interface
234                 typedef UniformFsgrCellIdentifier stdCellId;
235                 virtual CellIdentifierInterface* getFirstCell( );
236                 virtual void advanceCell( CellIdentifierInterface* );
237                 virtual bool noEndCell( CellIdentifierInterface* );
238                 virtual void deleteCellIterator( CellIdentifierInterface** );
239                 virtual CellIdentifierInterface* getCellAt( ntlVec3Gfx pos );
240                 virtual int        getCellSet      ( CellIdentifierInterface* );
241                 virtual ntlVec3Gfx getCellOrigin   ( CellIdentifierInterface* );
242                 virtual ntlVec3Gfx getCellSize     ( CellIdentifierInterface* );
243                 virtual int        getCellLevel    ( CellIdentifierInterface* );
244                 virtual LbmFloat   getCellDensity  ( CellIdentifierInterface* ,int set);
245                 virtual LbmVec     getCellVelocity ( CellIdentifierInterface* ,int set);
246                 virtual LbmFloat   getCellDf       ( CellIdentifierInterface* ,int set, int dir);
247                 virtual LbmFloat   getCellMass     ( CellIdentifierInterface* ,int set);
248                 virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set);
249                 virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set);
250                 virtual LbmFloat   getEquilDf      ( int );
251                 virtual ntlVec3Gfx getVelocityAt   (float x, float y, float z);
252                 // convert pointers
253                 stdCellId* convertBaseCidToStdCid( CellIdentifierInterface* basecid);
254
255                 //! perform geometry init (if switched on) 
256                 bool initGeometryFlags();
257                 //! init part for all freesurface testcases 
258                 void initFreeSurfaces();
259                 //! init density gradient if enabled
260                 void initStandingFluidGradient();
261
262                 /*! init a given cell with flag, density, mass and equilibrium dist. funcs */
263                 LBM_INLINED void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass);
264                 LBM_INLINED void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
265                 LBM_INLINED void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
266                 LBM_INLINED void forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
267                 LBM_INLINED void initInterfaceVars(int level, int i,int j,int k,int workSet, bool initMass);
268                 //! interpolate velocity and density at a given position
269                 void interpolateCellValues(int level,int ei,int ej,int ek,int workSet, LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz);
270
271                 /*! perform a single LBM step */
272                 void stepMain();
273                 //! advance fine grid
274                 void fineAdvance();
275                 //! advance coarse grid
276                 void coarseAdvance(int lev);
277                 //! update flux area values on coarse grids
278                 void coarseCalculateFluxareas(int lev);
279                 // adaptively coarsen grid
280                 bool adaptGrid(int lev);
281                 // restrict fine grid DFs to coarse grid
282                 void coarseRestrictFromFine(int lev);
283
284                 /* simulation object interface, just calls stepMain */
285                 virtual void step();
286                 /*! init particle positions */
287                 virtual int initParticles();
288                 /*! move all particles */
289                 virtual void advanceParticles();
290                 /*! move a particle at a boundary */
291                 void handleObstacleParticle(ParticleObject *p);
292                 /*! check whether to add particle 
293                 bool checkAddParticle();
294                 void performAddParticle();*/
295
296
297                 /*! debug object display (used e.g. for preview surface) */
298                 virtual vector<ntlGeometryObject*> getDebugObjects();
299         
300                 // gui/output debugging functions
301 #               if LBM_USE_GUI==1
302                 virtual void debugDisplayNode(int dispset, CellIdentifierInterface* cell );
303                 virtual void lbmDebugDisplay(int dispset);
304                 virtual void lbmMarkedCellDisplay();
305 #               endif // LBM_USE_GUI==1
306                 virtual void debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet=-1);
307
308                 //! for raytracing, preprocess
309                 void prepareVisualization( void );
310
311                 /* surface generation settings */
312                 virtual void setSurfGenSettings(short value);
313
314         protected:
315
316                 //! internal quick print function (for debugging) 
317                 void printLbmCell(int level, int i, int j, int k,int set);
318                 // debugging use CellIterator interface to mark cell
319                 void debugMarkCellCall(int level, int vi,int vj,int vk);
320                 
321                 // loop over grid, stream&collide update
322                 void mainLoop(int lev);
323                 // change time step size
324                 void adaptTimestep();
325                 //! init mObjectSpeeds for current parametrization
326                 void recalculateObjectSpeeds();
327                 //! init moving obstacles for next sim step sim 
328                 void initMovingObstacles(bool staticInit);
329                 //! flag reinit step - always works on finest grid!
330                 void reinitFlags( int workSet );
331                 //! mass dist weights
332                 LbmFloat getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l);
333                 //! compute surface normals: fluid, fluid more accurate, and for obstacles
334                 void computeFluidSurfaceNormal(LbmFloat *cell, CellFlagType *cellflag,    LbmFloat *snret);
335                 void computeFluidSurfaceNormalAcc(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret);
336                 void computeObstacleSurfaceNormal(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret);
337                 void computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret);
338                 //! add point to mListNewInter list
339                 LBM_INLINED void addToNewInterList( int ni, int nj, int nk );   
340                 //! cell is interpolated from coarse level (inited into set, source sets are determined by t)
341                 void interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet,bool markNbs);
342                 void coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet);
343
344                 //! minimal and maximal z-coords (for 2D/3D loops)
345                 LBM_INLINED int getForZMinBnd();
346                 LBM_INLINED int getForZMin1();
347                 LBM_INLINED int getForZMaxBnd(int lev);
348                 LBM_INLINED int getForZMax1(int lev);
349                 LBM_INLINED bool checkDomainBounds(int lev,int i,int j,int k);
350                 LBM_INLINED bool checkDomainBoundsPos(int lev,LbmVec pos);
351
352                 // touch grid and flags once
353                 void preinitGrids();
354                 // one relaxation step for standing fluid
355                 void standingFluidPreinit();
356
357
358                 // member vars
359
360                 //! mass calculated during streaming step
361                 LbmFloat mCurrentMass;
362                 LbmFloat mCurrentVolume;
363                 LbmFloat mInitialMass;
364
365                 //! count problematic cases, that occured so far...
366                 int mNumProblems;
367
368                 // average mlsups, count how many so far...
369                 double mAvgMLSUPS;
370                 double mAvgMLSUPSCnt;
371
372                 //! Mcubes object for surface reconstruction 
373                 IsoSurface *mpPreviewSurface;
374                 
375                 //! use time adaptivity? 
376                 bool mTimeAdap;
377                 //! force smaller timestep for next LBM step? (eg for mov obj)
378                 bool mForceTimeStepReduce;
379
380                 //! fluid vol height
381                 LbmFloat mFVHeight;
382                 LbmFloat mFVArea;
383                 bool mUpdateFVHeight;
384
385                 //! force quit for gfx
386                 LbmFloat mGfxEndTime;
387                 //! smoother surface initialization?
388                 int mInitSurfaceSmoothing;
389                 //! surface generation settings, default is all off (=0)
390                 //  each flag switches side on off,  fssgNoObs is for obstacle sides
391                 //  -1 equals all on
392                 typedef enum {
393                          fssgNormal   =  0,
394                          fssgNoNorth  =  1,
395                          fssgNoSouth  =  2,
396                          fssgNoEast   =  4,
397                          fssgNoWest   =  8,
398                          fssgNoTop    = 16,
399                          fssgNoBottom = 32,
400                          fssgNoObs    = 64
401                 } fsSurfaceGen;
402                 int mFsSurfGenSetting;
403
404                 //! lock time step down switching
405                 int mTimestepReduceLock;
406                 //! count no. of switches
407                 int mTimeSwitchCounts;
408                 // only switch of maxvel is higher for several steps...
409                 int mTimeMaxvelStepCnt;
410
411                 //! total simulation time so far 
412                 LbmFloat mSimulationTime, mLastSimTime;
413                 //! smallest and largest step size so far 
414                 LbmFloat mMinTimestep, mMaxTimestep;
415                 //! track max. velocity
416                 LbmFloat mMxvx, mMxvy, mMxvz, mMaxVlen;
417
418                 //! list of the cells to empty at the end of the step 
419                 vector<LbmPoint> mListEmpty;
420                 //! list of the cells to make fluid at the end of the step 
421                 vector<LbmPoint> mListFull;
422                 //! list of new interface cells to init
423         vector<LbmPoint> mListNewInter;
424                 //! class for handling redist weights in reinit flag function
425                 class lbmFloatSet {
426                         public:
427                                 LbmFloat val[dTotalNum];
428                                 LbmFloat numNbs;
429                 };
430                 //! normalized vectors for all neighboring cell directions (for e.g. massdweight calc)
431                 LbmVec mDvecNrm[27];
432                 
433                 
434                 //! debugging
435                 bool checkSymmetry(string idstring);
436                 //! kepp track of max/min no. of filled cells
437                 int mMaxNoCells, mMinNoCells;
438                 LONGINT mAvgNumUsedCells;
439
440                 //! precalculated objects speeds for current parametrization
441                 vector<LbmVec> mObjectSpeeds;
442                 //! partslip bc. values for obstacle boundary conditions
443                 vector<LbmFloat> mObjectPartslips;
444                 //! moving object mass boundary condition values
445                 vector<LbmFloat> mObjectMassMovnd;
446
447                 //! permanent movobj vert storage
448           vector<ntlVec3Gfx>  mMOIVertices;
449         vector<ntlVec3Gfx>  mMOIVerticesOld;
450           vector<ntlVec3Gfx>  mMOINormals;
451
452                 //! get isofield weights
453                 int mIsoWeightMethod;
454                 float mIsoWeight[27];
455
456                 // grid coarsening vars
457                 
458                 /*! vector for the data for each level */
459                 FsgrLevelData mLevel[FSGR_MAXNOOFLEVELS];
460
461                 /*! minimal and maximal refinement levels */
462                 int mMaxRefine;
463
464                 /*! df scale factors for level up/down */
465                 LbmFloat mDfScaleUp, mDfScaleDown;
466
467                 /*! precomputed cell area values */
468                 LbmFloat mFsgrCellArea[27];
469                 /*! restriction interpolation weights */
470                 LbmFloat mGaussw[27];
471
472                 /*! LES C_smago paramter for finest grid */
473                 float mInitialCsmago;
474                 /*! LES stats for non OPT3D */
475                 LbmFloat mDebugOmegaRet;
476                 /*! remember last init for animated params */
477                 LbmFloat mLastOmega;
478                 LbmVec   mLastGravity;
479
480                 //! fluid stats
481                 int mNumInterdCells;
482                 int mNumInvIfCells;
483                 int mNumInvIfTotal;
484                 int mNumFsgrChanges;
485
486                 //! debug function to disable standing f init
487                 int mDisableStandingFluidInit;
488                 //! init 2d with skipped Y/Z coords
489                 bool mInit2dYZ;
490                 //! debug function to force tadap syncing
491                 int mForceTadapRefine;
492                 //! border cutoff value
493                 int mCutoff;
494
495                 // strict debug interface
496 #               if FSGR_STRICT_DEBUG==1
497                 int debLBMGI(int level, int ii,int ij,int ik, int is);
498                 CellFlagType& debRFLAG(int level, int xx,int yy,int zz,int set);
499                 CellFlagType& debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir);
500                 CellFlagType& debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir);
501                 int debLBMQI(int level, int ii,int ij,int ik, int is, int l);
502                 LbmFloat& debQCELL(int level, int xx,int yy,int zz,int set,int l);
503                 LbmFloat& debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l);
504                 LbmFloat& debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l);
505                 LbmFloat* debRACPNT(int level,  int ii,int ij,int ik, int is );
506                 LbmFloat& debRAC(LbmFloat* s,int l);
507 #               endif // FSGR_STRICT_DEBUG==1
508
509                 LbmControlData *mpControl;
510
511                 void initCpdata();
512                 void handleCpdata();
513                 void cpDebugDisplay(int dispset); 
514
515                 bool mUseTestdata;
516 #               if LBM_INCLUDE_TESTSOLVERS==1
517                 // test functions
518                 LbmTestdata *mpTest;
519                 void initTestdata();
520                 void destroyTestdata();
521                 void handleTestdata();
522                 void set3dHeight(int ,int );
523
524                 int mMpNum,mMpIndex;
525                 int mOrgSizeX;
526                 LbmFloat mOrgStartX;
527                 LbmFloat mOrgEndX;
528                 void mrSetup();
529                 void mrExchange(); 
530                 void mrIsoExchange(); 
531                 LbmFloat mrInitTadap(LbmFloat max); 
532                 void gcFillBuffer(  LbmGridConnector *gc, int *retSizeCnt, const int *bdfs);
533                 void gcUnpackBuffer(LbmGridConnector *gc, int *retSizeCnt, const int *bdfs);
534         public:
535                 // needed for testdata
536                 void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz);
537                 // mptest
538                 int getMpIndex() { return mMpIndex; };
539 #               endif // LBM_INCLUDE_TESTSOLVERS==1
540
541                 // former LbmModelLBGK  functions
542                 // relaxation funtions - implemented together with relax macros
543                 static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz);
544                 static inline LbmFloat getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz);
545                 inline LbmFloat getLesNoneqTensorCoeff( LbmFloat df[],                          LbmFloat feq[] );
546                 inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo);
547                 inline void collideArrays( int lev, int i, int j, int k, // position - more for debugging
548                                 LbmFloat df[], LbmFloat &outrho, // out only!
549                                 // velocity modifiers (returns actual velocity!)
550                                 LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
551                                 LbmFloat omega, LbmVec gravity, LbmFloat csmago, 
552                                 LbmFloat *newOmegaRet, LbmFloat *newQoRet);
553
554
555                 // former LBM models
556                 //! shorten static const definitions
557 #               define STCON static const
558
559 #               if LBMDIM==3
560                 
561                 //! id string of solver
562                 virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D3Q19]"); }
563
564                 //! how many dimensions? UNUSED? replace by LBMDIM?
565                 STCON int cDimension;
566
567                 // Wi factors for collide step 
568                 STCON LbmFloat cCollenZero;
569                 STCON LbmFloat cCollenOne;
570                 STCON LbmFloat cCollenSqrtTwo;
571
572                 //! threshold value for filled/emptied cells 
573                 STCON LbmFloat cMagicNr2;
574                 STCON LbmFloat cMagicNr2Neg;
575                 STCON LbmFloat cMagicNr;
576                 STCON LbmFloat cMagicNrNeg;
577
578                 //! size of a single set of distribution functions 
579                 STCON int    cDfNum;
580                 //! direction vector contain vecs for all spatial dirs, even if not used for LBM model
581                 STCON int    cDirNum;
582
583                 //! distribution functions directions 
584                 typedef enum {
585                          cDirInv=  -1,
586                          cDirC  =  0,
587                          cDirN  =  1,
588                          cDirS  =  2,
589                          cDirE  =  3,
590                          cDirW  =  4,
591                          cDirT  =  5,
592                          cDirB  =  6,
593                          cDirNE =  7,
594                          cDirNW =  8,
595                          cDirSE =  9,
596                          cDirSW = 10,
597                          cDirNT = 11,
598                          cDirNB = 12,
599                          cDirST = 13,
600                          cDirSB = 14,
601                          cDirET = 15,
602                          cDirEB = 16,
603                          cDirWT = 17,
604                          cDirWB = 18
605                 } dfDir;
606
607                 /* Vector Order 3D:
608                  *  0   1  2   3  4   5  6       7  8  9 10  11 12 13 14  15 16 17 18     19 20 21 22  23 24 25 26
609                  *  0,  0, 0,  1,-1,  0, 0,      1,-1, 1,-1,  0, 0, 0, 0,  1, 1,-1,-1,     1,-1, 1,-1,  1,-1, 1,-1
610                  *  0,  1,-1,  0, 0,  0, 0,      1, 1,-1,-1,  1, 1,-1,-1,  0, 0, 0, 0,     1, 1,-1,-1,  1, 1,-1,-1
611                  *  0,  0, 0,  0, 0,  1,-1,      0, 0, 0, 0,  1,-1, 1,-1,  1,-1, 1,-1,     1, 1, 1, 1, -1,-1,-1,-1
612                  */
613
614                 /*! name of the dist. function 
615                          only for nicer output */
616                 STCON char* dfString[ 19 ];
617
618                 /*! index of normal dist func, not used so far?... */
619                 STCON int dfNorm[ 19 ];
620
621                 /*! index of inverse dist func, not fast, but useful... */
622                 STCON int dfInv[ 19 ];
623
624                 /*! index of x reflected dist func for free slip, not valid for all DFs... */
625                 STCON int dfRefX[ 19 ];
626                 /*! index of x reflected dist func for free slip, not valid for all DFs... */
627                 STCON int dfRefY[ 19 ];
628                 /*! index of x reflected dist func for free slip, not valid for all DFs... */
629                 STCON int dfRefZ[ 19 ];
630
631                 /*! dist func vectors */
632                 STCON int dfVecX[ 27 ];
633                 STCON int dfVecY[ 27 ];
634                 STCON int dfVecZ[ 27 ];
635
636                 /*! arrays as before with doubles */
637                 STCON LbmFloat dfDvecX[ 27 ];
638                 STCON LbmFloat dfDvecY[ 27 ];
639                 STCON LbmFloat dfDvecZ[ 27 ];
640
641                 /*! principal directions */
642                 STCON int princDirX[ 2*3 ];
643                 STCON int princDirY[ 2*3 ];
644                 STCON int princDirZ[ 2*3 ];
645
646                 /*! vector lengths */
647                 STCON LbmFloat dfLength[ 19 ];
648
649                 /*! equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
650                 static LbmFloat dfEquil[ dTotalNum ];
651
652                 /*! arrays for les model coefficients */
653                 static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ];
654                 static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ];
655
656 #               else // end LBMDIM==3 , LBMDIM==2
657                 
658                 //! id string of solver
659                 virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D2Q9]"); }
660
661                 //! how many dimensions?
662                 STCON int cDimension;
663
664                 //! Wi factors for collide step 
665                 STCON LbmFloat cCollenZero;
666                 STCON LbmFloat cCollenOne;
667                 STCON LbmFloat cCollenSqrtTwo;
668
669                 //! threshold value for filled/emptied cells 
670                 STCON LbmFloat cMagicNr2;
671                 STCON LbmFloat cMagicNr2Neg;
672                 STCON LbmFloat cMagicNr;
673                 STCON LbmFloat cMagicNrNeg;
674
675                 //! size of a single set of distribution functions 
676                 STCON int    cDfNum;
677                 STCON int    cDirNum;
678
679                 //! distribution functions directions 
680                 typedef enum {
681                          cDirInv=  -1,
682                          cDirC  =  0,
683                          cDirN  =  1,
684                          cDirS  =  2,
685                          cDirE  =  3,
686                          cDirW  =  4,
687                          cDirNE =  5,
688                          cDirNW =  6,
689                          cDirSE =  7,
690                          cDirSW =  8
691                 } dfDir;
692
693                 /* Vector Order 2D:
694                  * 0  1 2  3  4  5  6 7  8
695                  * 0, 0,0, 1,-1, 1,-1,1,-1 
696                  * 0, 1,-1, 0,0, 1,1,-1,-1  */
697
698                 /* name of the dist. function 
699                          only for nicer output */
700                 STCON char* dfString[ 9 ];
701
702                 /* index of normal dist func, not used so far?... */
703                 STCON int dfNorm[ 9 ];
704
705                 /* index of inverse dist func, not fast, but useful... */
706                 STCON int dfInv[ 9 ];
707
708                 /* index of x reflected dist func for free slip, not valid for all DFs... */
709                 STCON int dfRefX[ 9 ];
710                 /* index of x reflected dist func for free slip, not valid for all DFs... */
711                 STCON int dfRefY[ 9 ];
712                 /* index of x reflected dist func for free slip, not valid for all DFs... */
713                 STCON int dfRefZ[ 9 ];
714
715                 /* dist func vectors */
716                 STCON int dfVecX[ 9 ];
717                 STCON int dfVecY[ 9 ];
718                 /* Z, 2D values are all 0! */
719                 STCON int dfVecZ[ 9 ];
720
721                 /* arrays as before with doubles */
722                 STCON LbmFloat dfDvecX[ 9 ];
723                 STCON LbmFloat dfDvecY[ 9 ];
724                 /* Z, 2D values are all 0! */
725                 STCON LbmFloat dfDvecZ[ 9 ];
726
727                 /*! principal directions */
728                 STCON int princDirX[ 2*2 ];
729                 STCON int princDirY[ 2*2 ];
730                 STCON int princDirZ[ 2*2 ];
731
732                 /* vector lengths */
733                 STCON LbmFloat dfLength[ 9 ];
734
735                 /* equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
736                 static LbmFloat dfEquil[ dTotalNum ];
737
738                 /*! arrays for les model coefficients */
739                 static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ];
740                 static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ];
741
742 #               endif  // LBMDIM==2
743 };
744
745 #undef STCON
746
747
748 /*****************************************************************************/
749 // relaxation_macros
750
751
752
753 // cell mark debugging
754 #if FSGR_STRICT_DEBUG==10
755 #define debugMarkCell(lev,x,y,z) \
756         errMsg("debugMarkCell",this->mName<<" step: "<<this->mStepCnt<<" lev:"<<(lev)<<" marking "<<PRINT_VEC((x),(y),(z))<<" line "<< __LINE__ ); \
757         debugMarkCellCall((lev),(x),(y),(z));
758 #else // FSGR_STRICT_DEBUG==1
759 #define debugMarkCell(lev,x,y,z) \
760         debugMarkCellCall((lev),(x),(y),(z));
761 #endif // FSGR_STRICT_DEBUG==1
762
763
764 // flag array defines -----------------------------------------------------------------------------------------------
765
766 // lbm testsolver get index define, note - ignores is (set) as flag
767 // array is only a single entry
768 #define _LBMGI(level, ii,ij,ik, is) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
769
770 //! flag array acces macro
771 #define _RFLAG(level,xx,yy,zz,set) mLevel[level].mprsFlags[set][ LBMGI((level),(xx),(yy),(zz),(set)) ]
772 #define _RFLAG_NB(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set) ]
773 #define _RFLAG_NBINV(level,xx,yy,zz,set, dir) mLevel[level].mprsFlags[set][ LBMGI((level),(xx)+this->dfVecX[this->dfInv[dir]],(yy)+this->dfVecY[this->dfInv[dir]],(zz)+this->dfVecZ[this->dfInv[dir]],set) ]
774
775 // array handling  -----------------------------------------------------------------------------------------------
776
777 #define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
778 #define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)])
779 #define _QCELL_NB(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set, l)*dTotalNum +(l)])
780 #define _QCELL_NBINV(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+this->dfVecX[this->dfInv[dir]],(yy)+this->dfVecY[this->dfInv[dir]],(zz)+this->dfVecZ[this->dfInv[dir]],set, l)*dTotalNum +(l)])
781
782 #define QCELLSTEP dTotalNum
783 #define _RACPNT(level, ii,ij,ik, is )  &QCELL(level,ii,ij,ik,is,0)
784 #define _RAC(s,l) (s)[(l)]
785
786
787 #if FSGR_STRICT_DEBUG==1
788
789 #define LBMGI(level,ii,ij,ik, is)                 debLBMGI(level,ii,ij,ik, is)         
790 #define RFLAG(level,xx,yy,zz,set)                 debRFLAG(level,xx,yy,zz,set)            
791 #define RFLAG_NB(level,xx,yy,zz,set, dir)         debRFLAG_NB(level,xx,yy,zz,set, dir)    
792 #define RFLAG_NBINV(level,xx,yy,zz,set, dir)      debRFLAG_NBINV(level,xx,yy,zz,set, dir) 
793
794 #define LBMQI(level,ii,ij,ik, is, l)              debLBMQI(level,ii,ij,ik, is, l)         
795 #define QCELL(level,xx,yy,zz,set,l)               debQCELL(level,xx,yy,zz,set,l)         
796 #define QCELL_NB(level,xx,yy,zz,set, dir,l)       debQCELL_NB(level,xx,yy,zz,set, dir,l)
797 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    debQCELL_NBINV(level,xx,yy,zz,set, dir,l)
798 #define RACPNT(level, ii,ij,ik, is )              debRACPNT(level, ii,ij,ik, is )          
799 #define RAC(s,l)                                                debRAC(s,l)                  
800
801 #else // FSGR_STRICT_DEBUG==1
802
803 #define LBMGI(level,ii,ij,ik, is)                 _LBMGI(level,ii,ij,ik, is)         
804 #define RFLAG(level,xx,yy,zz,set)                 _RFLAG(level,xx,yy,zz,set)            
805 #define RFLAG_NB(level,xx,yy,zz,set, dir)         _RFLAG_NB(level,xx,yy,zz,set, dir)    
806 #define RFLAG_NBINV(level,xx,yy,zz,set, dir)      _RFLAG_NBINV(level,xx,yy,zz,set, dir) 
807
808 #define LBMQI(level,ii,ij,ik, is, l)              _LBMQI(level,ii,ij,ik, is, l)         
809 #define QCELL(level,xx,yy,zz,set,l)               _QCELL(level,xx,yy,zz,set,l)         
810 #define QCELL_NB(level,xx,yy,zz,set, dir,l)       _QCELL_NB(level,xx,yy,zz,set, dir, l)
811 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l)    _QCELL_NBINV(level,xx,yy,zz,set, dir,l)
812 #define RACPNT(level, ii,ij,ik, is )              _RACPNT(level, ii,ij,ik, is )          
813 #define RAC(s,l)                                  _RAC(s,l)                  
814
815 #endif // FSGR_STRICT_DEBUG==1
816
817 // general defines -----------------------------------------------------------------------------------------------
818
819 // replace TESTFLAG
820 #define FLAGISEXACT(flag, compflag)  ((flag & compflag)==compflag)
821
822 #if LBMDIM==2
823 #define dC 0
824 #define dN 1
825 #define dS 2
826 #define dE 3
827 #define dW 4
828 #define dNE 5
829 #define dNW 6
830 #define dSE 7
831 #define dSW 8
832 #else
833 // direction indices
834 #define dC 0
835 #define dN 1
836 #define dS 2
837 #define dE 3
838 #define dW 4
839 #define dT 5
840 #define dB 6
841 #define dNE 7
842 #define dNW 8
843 #define dSE 9
844 #define dSW 10
845 #define dNT 11
846 #define dNB 12
847 #define dST 13
848 #define dSB 14
849 #define dET 15
850 #define dEB 16
851 #define dWT 17
852 #define dWB 18
853 #endif
854 //? #define dWB 18
855
856 // default init for dFlux values
857 #define FLUX_INIT 0.5f * (float)(this->cDfNum)
858
859 // only for non DF dir handling!
860 #define dNET 19
861 #define dNWT 20
862 #define dSET 21
863 #define dSWT 22
864 #define dNEB 23
865 #define dNWB 24
866 #define dSEB 25
867 #define dSWB 26
868
869 //! fill value for boundary cells
870 #define BND_FILL 0.0
871
872 #define DFL1 (1.0/ 3.0)
873 #define DFL2 (1.0/18.0)
874 #define DFL3 (1.0/36.0)
875
876 // loops over _all_ cells (including boundary layer)
877 #define FSGR_FORIJK_BOUNDS(leveli) \
878         for(int k= getForZMinBnd(); k< getForZMaxBnd(leveli); ++k) \
879    for(int j=0;j<mLevel[leveli].lSizey-0;++j) \
880     for(int i=0;i<mLevel[leveli].lSizex-0;++i) \
881         
882 // loops over _only inner_ cells 
883 #define FSGR_FORIJK1(leveli) \
884         for(int k= getForZMin1(); k< getForZMax1(leveli); ++k) \
885    for(int j=1;j<mLevel[leveli].lSizey-1;++j) \
886     for(int i=1;i<mLevel[leveli].lSizex-1;++i) \
887
888 // relaxation_macros end
889
890
891
892 /******************************************************************************/
893 /*! equilibrium functions */
894 /******************************************************************************/
895
896 /*! calculate length of velocity vector */
897 inline LbmFloat LbmFsgrSolver::getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz) {
898         return ((ux)*dfDvecX[l]+(uy)*dfDvecY[l]+(uz)*dfDvecZ[l]);
899 };
900
901 /*! calculate equilibrium DF for given values */
902 inline LbmFloat LbmFsgrSolver::getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz) {
903 #if FSGR_STRICT_DEBUG==1
904         if((l<0)||(l>LBM_DFNUM)) { errFatal("LbmFsgrSolver::getCollideEq","Invalid DFEQ call "<<l, SIMWORLD_PANIC ); /* no access to mPanic here */     }
905 #endif // FSGR_STRICT_DEBUG==1
906         LbmFloat tmp = getVelVecLen(l,ux,uy,uz); 
907         return( dfLength[l] *( 
908                                 + rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) 
909                                 + 3.0 *tmp 
910                                 + 9.0/2.0 *(tmp*tmp) )
911                         ); // */
912 };
913
914 /*****************************************************************************/
915 /* init a given cell with flag, density, mass and equilibrium dist. funcs */
916
917 void LbmFsgrSolver::forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
918         // also overwrite persisting flags
919         // function is useful for tracking accesses...
920         RFLAG(level,xx,yy,zz,set) = newflag;
921 }
922 void LbmFsgrSolver::changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
923         CellFlagType pers = RFLAG(level,xx,yy,zz,set) & CFPersistMask;
924         RFLAG(level,xx,yy,zz,set) = newflag | pers;
925 }
926
927 void 
928 LbmFsgrSolver::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass) {
929   /* init eq. dist funcs */
930         LbmFloat *ecel;
931         int workSet = mLevel[level].setCurr;
932
933         ecel = RACPNT(level, i,j,k, workSet);
934         FORDF0 { RAC(ecel, l) = this->dfEquil[l] * rho; }
935         RAC(ecel, dMass) = mass;
936         RAC(ecel, dFfrac) = mass/rho;
937         RAC(ecel, dFlux) = FLUX_INIT;
938         changeFlag(level, i,j,k, workSet, flag);
939
940   workSet ^= 1;
941         changeFlag(level, i,j,k, workSet, flag);
942         return;
943 }
944
945 void 
946 LbmFsgrSolver::initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel) {
947         LbmFloat *ecel;
948         int workSet = mLevel[level].setCurr;
949
950         ecel = RACPNT(level, i,j,k, workSet);
951         FORDF0 { RAC(ecel, l) = getCollideEq(l, rho,vel[0],vel[1],vel[2]); }
952         RAC(ecel, dMass) = mass;
953         RAC(ecel, dFfrac) = mass/rho;
954         RAC(ecel, dFlux) = FLUX_INIT;
955         changeFlag(level, i,j,k, workSet, flag);
956
957   workSet ^= 1;
958         changeFlag(level, i,j,k, workSet, flag);
959         return;
960 }
961
962 int LbmFsgrSolver::getForZMinBnd() { 
963         return 0; 
964 }
965 int LbmFsgrSolver::getForZMin1()   { 
966         if(LBMDIM==2) return 0;
967         return 1; 
968 }
969
970 int LbmFsgrSolver::getForZMaxBnd(int lev) { 
971         if(LBMDIM==2) return 1;
972         return mLevel[lev].lSizez -0;
973 }
974 int LbmFsgrSolver::getForZMax1(int lev)   { 
975         if(LBMDIM==2) return 1;
976         return mLevel[lev].lSizez -1;
977 }
978
979 bool LbmFsgrSolver::checkDomainBounds(int lev,int i,int j,int k) { 
980         if(i<0) return false;
981         if(j<0) return false;
982         if(k<0) return false;
983         if(i>mLevel[lev].lSizex-1) return false;
984         if(j>mLevel[lev].lSizey-1) return false;
985         if(k>mLevel[lev].lSizez-1) return false;
986         return true;
987 }
988 bool LbmFsgrSolver::checkDomainBoundsPos(int lev,LbmVec pos) { 
989         const int i= (int)pos[0]; 
990         if(i<0) return false;
991         if(i>mLevel[lev].lSizex-1) return false;
992         const int j= (int)pos[1]; 
993         if(j<0) return false;
994         if(j>mLevel[lev].lSizey-1) return false;
995         const int k= (int)pos[2];
996         if(k<0) return false;
997         if(k>mLevel[lev].lSizez-1) return false;
998         return true;
999 }
1000
1001 void LbmFsgrSolver::initInterfaceVars(int level, int i,int j,int k,int workSet, bool initMass) {
1002         LbmFloat *ccel = &QCELL(level ,i,j,k, workSet,0);
1003         LbmFloat nrho = 0.0;
1004         FORDF0 { nrho += RAC(ccel,l); }
1005         if(initMass) {
1006                 RAC(ccel,dMass) = nrho;
1007         RAC(ccel, dFfrac) = 1.;
1008         } else {
1009                 // preinited, e.g. from reinitFlags
1010                 RAC(ccel, dFfrac) = RAC(ccel, dMass)/nrho;
1011                 RAC(ccel, dFlux) = FLUX_INIT;
1012         }
1013 }
1014
1015
1016 #endif
1017
1018