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