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