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