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