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