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