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