← Back to branch summary

~siretart/ubuntu/utopic/blender/libav10

« back to all changes in this revision

Viewing changes to intern/elbeem/intern/solver_interface.h

Tags: upstream-2.40
Import upstream version 2.40

Show diffs side-by-side

added added

removed removed

Lines of Context:
intern/elbeem/intern/solver_interface.h
 
1
/******************************************************************************
 
2
 *
 
3
 * El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
 
4
 * All code distributed as part of El'Beem is covered by the version 2 of the 
 
5
 * GNU General Public License. See the file COPYING for details.
 
6
 * Copyright 2003-2005 Nils Thuerey
 
7
 *
 
8
 * Header for Combined 2D/3D Lattice Boltzmann Interface Class
 
9
 * 
 
10
 *****************************************************************************/
 
11
#ifndef LBMINTERFACE_H
 
12
#define LBMINTERFACE_H
 
13
 
 
14
//! include gui support?
 
15
#ifndef NOGUI
 
16
#define LBM_USE_GUI      1
 
17
#else
 
18
#define LBM_USE_GUI      0
 
19
#endif
 
20
 
 
21
#if LBM_USE_GUI==1
 
22
#define USE_GLUTILITIES
 
23
// for debug display
 
24
#include <GL/gl.h>
 
25
#include "../gui/guifuncs.h"
 
26
#endif
 
27
 
 
28
#include <sstream>
 
29
#include "utilities.h"
 
30
#include "ntl_bsptree.h"
 
31
#include "ntl_geometryobject.h"
 
32
#include "ntl_rndstream.h"
 
33
#include "parametrizer.h"
 
34
#include "attributes.h"
 
35
#include "particletracer.h"
 
36
#include "isosurface.h"
 
37
 
 
38
// use which fp-precision for LBM? 1=float, 2=double
 
39
#ifdef PRECISION_LBM_SINGLE
 
40
#define LBM_PRECISION 1
 
41
#else
 
42
#ifdef PRECISION_LBM_DOUBLE
 
43
#define LBM_PRECISION 2
 
44
#else
 
45
// default to floats
 
46
#define LBM_PRECISION 1
 
47
#endif
 
48
#endif
 
49
 
 
50
#if LBM_PRECISION==1
 
51
/* low precision for LBM solver */
 
52
typedef float LbmFloat;
 
53
typedef ntlVec3f LbmVec;
 
54
#define LBM_EPSILON (1e-5)
 
55
#else
 
56
/* standard precision for LBM solver */
 
57
typedef double LbmFloat;
 
58
typedef ntlVec3d LbmVec;
 
59
#define LBM_EPSILON (1e-10)
 
60
#endif
 
61
 
 
62
// long integer, needed e.g. for memory calculations
 
63
#ifndef USE_MSVC6FIXES
 
64
#define LONGINT long long int
 
65
#else
 
66
#define LONGINT _int64
 
67
#endif
 
68
 
 
69
 
 
70
// conversions (lbm and parametrizer)
 
71
template<class T> inline LbmVec     vec2L(T v) { return LbmVec(v[0],v[1],v[2]); }
 
72
template<class T> inline ParamVec   vec2P(T v) { return ParamVec(v[0],v[1],v[2]); }
 
73
 
 
74
 
 
75
// bubble id type
 
76
typedef int BubbleId;
 
77
 
 
78
// for both short int/char
 
79
#define CFUnused              (1<< 0)
 
80
#define CFEmpty               (1<< 1)
 
81
#define CFBnd                 (1<< 2)
 
82
#define CFBndNoslip           (1<< 3)
 
83
#define CFBndFreeslip         (1<< 4)
 
84
#define CFBndPartslip         (1<< 5)
 
85
// force symmetry for flag reinit
 
86
#define CFNoInterpolSrc       (1<< 6) 
 
87
#define CFFluid               (1<< 7)
 
88
#define CFInter               (1<< 8)
 
89
#define CFNoNbFluid           (1<< 9)
 
90
#define CFNoNbEmpty           (1<<10)
 
91
#define CFNoDelete            (1<<11)
 
92
#define CFNoBndFluid          (1<<12)
 
93
        
 
94
//! refinement tags
 
95
// cell treated normally on coarser grids
 
96
#define CFGrNorm              (1<<13)
 
97
// border cells to be interpolated from finer grid
 
98
#define CFGrFromFine          (1<<14)
 
99
#define CFGrFromCoarse        (1<<15)
 
100
#define CFGrCoarseInited      (1<<16)
 
101
// 32k aux border marker 
 
102
#define CFGrToFine            (1<<17)
 
103
#define CFMbndInflow          (1<<18)
 
104
#define CFMbndOutflow         (1<<19)
 
105
 
 
106
// debug/helper type
 
107
#define CFIgnore              (1<<20)
 
108
 
 
109
// above 24 is used to encode in/outflow object type
 
110
#define CFPersistMask (0xFF000000 | CFMbndInflow | CFMbndOutflow)
 
111
 
 
112
// nk
 
113
#define CFInvalid             (CellFlagType)(1<<31)
 
114
 
 
115
// use 32bit flag types
 
116
#ifdef __x86_64__
 
117
 typedef int cfINT32;
 
118
#else
 
119
 typedef long cfINT32;
 
120
#endif // defined (_IA64)  
 
121
#define CellFlagType cfINT32
 
122
#define CellFlagTypeSize 4
 
123
 
 
124
 
 
125
// aux. field indices (same for 2d)
 
126
#define dFfrac 19
 
127
#define dMass 20
 
128
#define dFlux 21
 
129
// max. no. of cell values for 3d
 
130
#define dTotalNum 22
 
131
 
 
132
 
 
133
/*****************************************************************************/
 
134
/*! a single lbm cell */
 
135
/*  the template is only needed for 
 
136
 *  dimension dependend constants e.g. 
 
137
 *  number of df's in model */
 
138
template<typename D>
 
139
class LbmCellTemplate {
 
140
        public:
 
141
                LbmFloat     df[ 27 ]; // be on the safe side here...
 
142
        LbmFloat     rho;
 
143
                LbmVec       vel;
 
144
        LbmFloat     mass;
 
145
                CellFlagType flag;
 
146
                BubbleId     bubble;
 
147
        LbmFloat     ffrac;
 
148
 
 
149
                //! test if a flag is set 
 
150
                inline bool test(CellFlagType t) {
 
151
                        return ((flag & t)==t);
 
152
                }
 
153
                //! test if any of the given flags is set 
 
154
                inline bool testAny(CellFlagType t) {
 
155
                        return ((flag & t)!=0);
 
156
                }
 
157
                //! test if the cell is empty 
 
158
                inline bool isEmpty() {
 
159
                        return (flag == CFEmpty);
 
160
                }
 
161
 
 
162
                //! init default values for a certain flag type
 
163
                inline void initDefaults(CellFlagType type) {
 
164
                        flag = type;
 
165
                        vel = LbmVec(0.0);
 
166
                        for(int l=0; l<D::cDfNum;l++) df[l] = D::dfEquil[l];
 
167
                                
 
168
                        if(type & CFFluid) {
 
169
                                rho = mass = ffrac = 1.0;
 
170
                                bubble = -1;
 
171
                        }
 
172
                        else if(type & CFInter) {
 
173
                                rho = mass = ffrac = 0.0;
 
174
                                bubble = 0;
 
175
                        }
 
176
                        else if(type & CFBnd) {
 
177
                                rho = mass = ffrac = 0.0;
 
178
                                bubble = -1;
 
179
                        }
 
180
                        else if(type & CFEmpty) {
 
181
                                rho = mass = ffrac = 0.0;
 
182
                                bubble = 0;
 
183
                        } else {
 
184
                                // ?
 
185
                                rho = mass = ffrac = 0.0;
 
186
                                bubble = -1;
 
187
                        }
 
188
                }
 
189
 
 
190
                //TODO add init method?
 
191
};
 
192
 
 
193
 
 
194
/* struct for the coordinates of a cell in the grid */
 
195
typedef struct {
 
196
  int x,y,z;
 
197
} LbmPoint;
 
198
 
 
199
 
 
200
/* struct for the coordinates of a cell in the grid */
 
201
typedef struct {
 
202
        char active;            // bubble in use, oder may be overwritten?
 
203
  LbmFloat volume;          // volume of this bubble (0 vor atmosphere)
 
204
        LbmFloat mass;            // "mass" of bubble 
 
205
        int i,j,k;              // index of a cell in the bubble
 
206
} LbmBubble;
 
207
 
 
208
 
 
209
 
 
210
 
 
211
//! choose which data to display
 
212
#define FLUIDDISPINVALID    0
 
213
#define FLUIDDISPNothing    1
 
214
#define FLUIDDISPCelltypes  2
 
215
#define FLUIDDISPVelocities 3
 
216
#define FLUIDDISPCellfills  4
 
217
#define FLUIDDISPDensity    5
 
218
#define FLUIDDISPGrid       6
 
219
#define FLUIDDISPSurface    7
 
220
 
 
221
//! settings for a debug display
 
222
typedef struct fluidDispSettings_T {
 
223
        int            type;  // what to display
 
224
        bool           on;    // display enabled?
 
225
        float          scale; // additional scale param
 
226
} fluidDispSettings;
 
227
 
 
228
 
 
229
 
 
230
/*****************************************************************************/
 
231
//! cell identifier interface
 
232
class CellIdentifierInterface {
 
233
        public:
 
234
                //! reset constructor
 
235
                CellIdentifierInterface():mEnd(false) { };
 
236
                //! virtual destructor
 
237
                virtual ~CellIdentifierInterface() {};
 
238
 
 
239
                //! return node as string (with some basic info)
 
240
                virtual string getAsString() = 0;
 
241
 
 
242
                //! compare cids
 
243
                virtual bool equal(CellIdentifierInterface* other) = 0;
 
244
 
 
245
                //! set/get end flag for grid traversal (not needed for marked cells)
 
246
                inline void setEnd(bool set){ mEnd = set;  }
 
247
                inline bool getEnd( )       { return mEnd; }
 
248
 
 
249
                //! has the grid been traversed?
 
250
                bool mEnd;
 
251
 
 
252
};
 
253
 
 
254
 
 
255
 
 
256
/*****************************************************************************/
 
257
/*! class defining abstract function interface */
 
258
/*  has to provide iterating functionality */
 
259
class LbmSolverInterface 
 
260
{
 
261
        public:
 
262
                //! Constructor 
 
263
                LbmSolverInterface();
 
264
                //! Destructor 
 
265
                virtual ~LbmSolverInterface() { };
 
266
                //! id string of solver
 
267
                virtual string getIdString() = 0;
 
268
                //! dimension of solver
 
269
                virtual int getDimension() = 0;
 
270
 
 
271
                /*! finish the init with config file values (allocate arrays...) */
 
272
                virtual bool initializeSolver() =0; //( ntlTree *tree, vector<ntlGeometryObject*> *objects ) = 0;
 
273
 
 
274
                /*! parse a boundary flag string */
 
275
                CellFlagType readBoundaryFlagInt(string name, int defaultValue, string source,string target, bool needed);
 
276
                /*! parse standard attributes */
 
277
                void parseStdAttrList();
 
278
                /*! initilize variables fom attribute list (should at least call std parse) */
 
279
                virtual void parseAttrList() = 0;
 
280
 
 
281
                virtual void step() = 0;
 
282
                virtual void prepareVisualization() { /* by default off */ };
 
283
 
 
284
                /*! particle handling */
 
285
                virtual int initParticles(ParticleTracer *partt) = 0;
 
286
                virtual void advanceParticles(ParticleTracer *partt ) = 0;
 
287
                /*! get surface object (NULL if no surface) */
 
288
                ntlGeometryObject* getSurfaceGeoObj() { return mpIso; }
 
289
 
 
290
                /*! debug object display */
 
291
                virtual vector<ntlGeometryObject*> getDebugObjects() { vector<ntlGeometryObject*> empty(0); return empty; }
 
292
 
 
293
#if LBM_USE_GUI==1
 
294
                /*! show simulation info */
 
295
                virtual void debugDisplay(fluidDispSettings *) = 0;
 
296
#endif
 
297
 
 
298
                /*! init tree for certain geometry init */
 
299
                void initGeoTree(int id);
 
300
                /*! destroy tree etc. when geometry init done */
 
301
                void freeGeoTree();
 
302
                /*! check for a certain flag type at position org (needed for e.g. quadtree refinement) */
 
303
                bool geoInitCheckPointInside(ntlVec3Gfx org, int flags, int &OId, gfxReal &distance);
 
304
                bool geoInitCheckPointInside(ntlVec3Gfx org, ntlVec3Gfx dir, int flags, int &OId, gfxReal &distance, 
 
305
                                const gfxReal halfCellsize, bool &thinHit, bool recurse);
 
306
                /*! set render globals, for scene/tree access */
 
307
                void setRenderGlobals(ntlRenderGlobals *glob) { mpGlob = glob; };
 
308
                /*! get max. velocity of all objects to initialize as fluid regions */
 
309
                ntlVec3Gfx getGeoMaxInitialVelocity();
 
310
 
 
311
                /* rt interface functions */
 
312
                unsigned int getIsoVertexCount()  { return mpIso->getIsoVertexCount(); }
 
313
                unsigned int getIsoIndexCount()   { return mpIso->getIsoIndexCount(); }
 
314
                char* getIsoVertexArray()         { return mpIso->getIsoVertexArray(); }
 
315
                unsigned int *getIsoIndexArray()  { return mpIso->getIsoIndexArray(); }
 
316
                void triangulateSurface()         { mpIso->triangulate(); }
 
317
 
 
318
                /* access functions */
 
319
 
 
320
                /*! return grid sizes */
 
321
                int getSizeX( void ) { return mSizex; }
 
322
                int getSizeY( void ) { return mSizey; }
 
323
                int getSizeZ( void ) { return mSizez; }
 
324
                /*! return grid sizes */
 
325
                void setSizeX( int ns ) { mSizex = ns; }
 
326
                void setSizeY( int ns ) { mSizey = ns; }
 
327
                void setSizeZ( int ns ) { mSizez = ns; }
 
328
                /*! access fluid only simulation flag */
 
329
                void setAllfluid(bool set) { mAllfluid=set; }
 
330
                bool getAllfluid()         { return mAllfluid; }
 
331
 
 
332
                /*! set attr list pointer */
 
333
                void setAttrList(AttributeList *set) { mpAttrs = set; }
 
334
                /*! Returns the attribute list pointer */
 
335
                inline AttributeList *getAttributeList() { return mpAttrs; }
 
336
 
 
337
                /*! set parametrizer pointer */
 
338
                inline void setParametrizer(Parametrizer *set) { mpParam = set; }
 
339
                /*! get parametrizer pointer */
 
340
                inline Parametrizer *getParametrizer() { return mpParam; }
 
341
 
 
342
                /*! set density gradient init from e.g. init test cases */
 
343
                inline void setInitDensityGradient(bool set) { mInitDensityGradient = set; }
 
344
 
 
345
                /*! access geometry start vector */
 
346
                inline void setGeoStart(ntlVec3Gfx set) { mvGeoStart = set; }
 
347
                inline ntlVec3Gfx getGeoStart() const   { return mvGeoStart; }
 
348
 
 
349
                /*! access geometry end vector */
 
350
                inline void setGeoEnd(ntlVec3Gfx set)   { mvGeoEnd = set; }
 
351
                inline ntlVec3Gfx getGeoEnd() const     { return mvGeoEnd; }
 
352
 
 
353
                /*! access geo init vars */
 
354
                inline void setGeoInitId(int set)       { mGeoInitId = set; }
 
355
                inline int getGeoInitId() const { return mGeoInitId; }
 
356
 
 
357
                /*! access name string */
 
358
                inline void setName(string set) { mName = set; }
 
359
                inline string getName() const   { return mName; }
 
360
 
 
361
                /*! access string for node info debugging output */
 
362
                inline string getNodeInfoString() const { return mNodeInfoString; }
 
363
 
 
364
                /*! get panic flag */
 
365
                inline bool getPanic() { return mPanic; }
 
366
 
 
367
                //! set silent mode?
 
368
                inline void setSilent(bool set){ mSilent = set; }
 
369
 
 
370
                //! set amount of surface/normal smoothing
 
371
                inline void setSmoothing(float setss,float setns){ mSmoothSurface=setss; mSmoothNormals=setns; }
 
372
                //! set desired refinement
 
373
                inline void setPreviewSize(int set){ mOutputSurfacePreview = set; }
 
374
                //! set desired refinement
 
375
                inline void setRefinementDesired(int set){ mRefinementDesired = set; }
 
376
 
 
377
 
 
378
                // cell iterator interface
 
379
                
 
380
                // cell id type
 
381
                typedef CellIdentifierInterface* CellIdentifier;
 
382
 
 
383
                //! cell iteration methods
 
384
                virtual CellIdentifierInterface* getFirstCell( ) = 0;
 
385
                virtual void advanceCell( CellIdentifierInterface* ) = 0;
 
386
                virtual bool noEndCell( CellIdentifierInterface* ) = 0;
 
387
                //! clean up iteration, this should be called, when the iteration is not completely finished
 
388
                virtual void deleteCellIterator( CellIdentifierInterface** ) = 0;
 
389
 
 
390
                //! find cell at a given position (returns NULL if not in domain)
 
391
                virtual CellIdentifierInterface* getCellAt( ntlVec3Gfx pos ) = 0;
 
392
 
 
393
                //! return node information
 
394
                virtual int        getCellSet      ( CellIdentifierInterface* ) = 0;
 
395
                virtual ntlVec3Gfx getCellOrigin   ( CellIdentifierInterface* ) = 0;
 
396
                virtual ntlVec3Gfx getCellSize     ( CellIdentifierInterface* ) = 0;
 
397
                virtual int        getCellLevel    ( CellIdentifierInterface* ) = 0;
 
398
                virtual LbmFloat   getCellDensity  ( CellIdentifierInterface*,int ) = 0;
 
399
                virtual LbmVec     getCellVelocity ( CellIdentifierInterface*,int ) = 0;
 
400
                /*! get equilibrium distribution functions */
 
401
                virtual LbmFloat   getEquilDf      ( int ) = 0;
 
402
                /*! get number of distribution functions */
 
403
                virtual int        getDfNum        ( ) = 0;
 
404
                /*! redundant cell functions */
 
405
                virtual LbmFloat   getCellDf       ( CellIdentifierInterface* ,int set, int dir) = 0;
 
406
                virtual LbmFloat   getCellMass     ( CellIdentifierInterface* ,int set) = 0;
 
407
                virtual LbmFloat   getCellFill     ( CellIdentifierInterface* ,int set) = 0;
 
408
                virtual CellFlagType getCellFlag   ( CellIdentifierInterface* ,int set) = 0;
 
409
 
 
410
                // gui/output debugging functions
 
411
#if LBM_USE_GUI==1
 
412
                virtual void debugDisplayNode(fluidDispSettings *dispset, CellIdentifier cell ) = 0;
 
413
                virtual void lbmDebugDisplay(fluidDispSettings *dispset) = 0;
 
414
                virtual void lbmMarkedCellDisplay() = 0;
 
415
#endif // LBM_USE_GUI==1
 
416
                virtual void debugPrintNodeInfo(CellIdentifier cell, int forceSet=-1) = 0;
 
417
 
 
418
                // debugging cell marker functions
 
419
 
 
420
                //! add cell to mMarkedCells list
 
421
                void addCellToMarkedList( CellIdentifierInterface *cid );
 
422
                //! marked cell iteration methods
 
423
                CellIdentifierInterface* markedGetFirstCell( );
 
424
                CellIdentifierInterface* markedAdvanceCell();
 
425
                void markedClearList();
 
426
 
 
427
 
 
428
        protected:
 
429
 
 
430
                /*! abort simulation on error... */
 
431
                bool mPanic;
 
432
 
 
433
 
 
434
                /*! Size of the array in x,y,z direction */
 
435
                int mSizex, mSizey, mSizez;
 
436
                /*! only fluid in sim? */
 
437
                bool mAllfluid;
 
438
 
 
439
 
 
440
                /*! step counter */
 
441
                int mStepCnt;
 
442
 
 
443
                /*! mass change from one step to the next, for extreme cases fix globally */
 
444
                LbmFloat mFixMass;
 
445
 
 
446
                // deprecated param vars
 
447
                /*! omega for lbm */
 
448
                LbmFloat mOmega;
 
449
                /*! gravity strength in neg. z direction */
 
450
                LbmVec mGravity;
 
451
                /*! Surface tension of the fluid */
 
452
                LbmFloat mSurfaceTension;
 
453
 
 
454
 
 
455
                /* boundary inits */
 
456
                CellFlagType mBoundaryEast, mBoundaryWest, 
 
457
                  mBoundaryNorth, mBoundarySouth, 
 
458
                  mBoundaryTop, mBoundaryBottom;
 
459
 
 
460
                /*! initialization from config file done? */
 
461
                int mInitDone;
 
462
 
 
463
                /*! init density gradient? */
 
464
                bool mInitDensityGradient;
 
465
 
 
466
                /*! pointer to the attribute list */
 
467
                AttributeList *mpAttrs;
 
468
 
 
469
                /*! get parameters from this parametrize in finishInit */
 
470
                Parametrizer *mpParam;
 
471
 
 
472
                /*! number of particles lost so far */
 
473
                int mNumParticlesLost;
 
474
                /*! number of particles lost so far */
 
475
                int mNumInvalidDfs;
 
476
                /*! no of filled/emptied cells per time step */
 
477
                int mNumFilledCells, mNumEmptiedCells;
 
478
                /*! counter number of used cells for performance */
 
479
                int mNumUsedCells;
 
480
                /*! MLSUPS counter */
 
481
                LbmFloat mMLSUPS;
 
482
                /*! debug - velocity output scaling factor */
 
483
                LbmFloat mDebugVelScale;
 
484
                /*! string for node info debugging output */
 
485
                string mNodeInfoString;
 
486
 
 
487
 
 
488
                /*! an own random stream */
 
489
                ntlRandomStream mRandom;
 
490
 
 
491
 
 
492
                // geo init vars
 
493
                // TODO deprecate SimulationObject vars
 
494
 
 
495
                /*! for display - start and end vectors for geometry */
 
496
                ntlVec3Gfx mvGeoStart, mvGeoEnd;
 
497
 
 
498
                /*! perform accurate geometry init? */
 
499
                bool mAccurateGeoinit;
 
500
 
 
501
                /*! name of this lbm object (for debug output) */
 
502
                string mName;
 
503
 
 
504
                //! Mcubes object for surface reconstruction 
 
505
                IsoSurface *mpIso;
 
506
                /*! isolevel value for marching cubes surface reconstruction */
 
507
                LbmFloat mIsoValue;
 
508
 
 
509
                //! debug output?
 
510
                bool mSilent;
 
511
 
 
512
                /*! geometry init id */
 
513
                int mGeoInitId;
 
514
                /*! tree object for geomerty initialization */
 
515
                ntlTree *mpGiTree;
 
516
                /*! object vector for geo init */
 
517
                vector<ntlGeometryObject*> *mpGiObjects;
 
518
                /*! inside which objects? */
 
519
                vector<int> mGiObjInside;
 
520
                /*! inside which objects? */
 
521
                vector<gfxReal> mGiObjDistance;
 
522
                vector<gfxReal> mGiObjSecondDist;
 
523
                /*! remember globals */
 
524
                ntlRenderGlobals *mpGlob;
 
525
                
 
526
                //! use refinement/coarsening?
 
527
                int mRefinementDesired;
 
528
 
 
529
                //! output surface preview? if >0 yes, and use as reduzed size 
 
530
                int mOutputSurfacePreview;
 
531
                LbmFloat mPreviewFactor;
 
532
 
 
533
                /* enable surface and normals smoothing? */
 
534
                float mSmoothSurface;
 
535
                float mSmoothNormals;
 
536
 
 
537
                // list for marked cells
 
538
                vector<CellIdentifierInterface *> mMarkedCells;
 
539
                int mMarkedCellIndex;
 
540
};
 
541
 
 
542
 
 
543
//! shorten static const definitions
 
544
#define STCON static const
 
545
 
 
546
 
 
547
/*****************************************************************************/
 
548
/*! class for solver templating - 3D implementation */
 
549
class LbmD3Q19 {
 
550
 
 
551
        public:
 
552
 
 
553
                // constructor, init interface
 
554
                LbmD3Q19() {};
 
555
                // virtual destructor 
 
556
                virtual ~LbmD3Q19() {};
 
557
                //! id string of solver
 
558
                string getIdString() { return string("3D"); }
 
559
 
 
560
                //! how many dimensions?
 
561
                STCON int cDimension;
 
562
 
 
563
                // Wi factors for collide step 
 
564
                STCON LbmFloat cCollenZero;
 
565
                STCON LbmFloat cCollenOne;
 
566
                STCON LbmFloat cCollenSqrtTwo;
 
567
 
 
568
                //! threshold value for filled/emptied cells 
 
569
                STCON LbmFloat cMagicNr2;
 
570
                STCON LbmFloat cMagicNr2Neg;
 
571
                STCON LbmFloat cMagicNr;
 
572
                STCON LbmFloat cMagicNrNeg;
 
573
 
 
574
                //! size of a single set of distribution functions 
 
575
                STCON int    cDfNum;
 
576
                //! direction vector contain vecs for all spatial dirs, even if not used for LBM model
 
577
                STCON int    cDirNum;
 
578
 
 
579
                //! distribution functions directions 
 
580
                typedef enum {
 
581
                         cDirInv=  -1,
 
582
                         cDirC  =  0,
 
583
                         cDirN  =  1,
 
584
                         cDirS  =  2,
 
585
                         cDirE  =  3,
 
586
                         cDirW  =  4,
 
587
                         cDirT  =  5,
 
588
                         cDirB  =  6,
 
589
                         cDirNE =  7,
 
590
                         cDirNW =  8,
 
591
                         cDirSE =  9,
 
592
                         cDirSW = 10,
 
593
                         cDirNT = 11,
 
594
                         cDirNB = 12,
 
595
                         cDirST = 13,
 
596
                         cDirSB = 14,
 
597
                         cDirET = 15,
 
598
                         cDirEB = 16,
 
599
                         cDirWT = 17,
 
600
                         cDirWB = 18
 
601
                } dfDir;
 
602
 
 
603
                /* Vector Order 3D:
 
604
                 *  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
 
605
                 *  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
 
606
                 *  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
 
607
                 *  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
 
608
                 */
 
609
 
 
610
                /*! name of the dist. function 
 
611
                         only for nicer output */
 
612
                STCON char* dfString[ 19 ];
 
613
 
 
614
                /*! index of normal dist func, not used so far?... */
 
615
                STCON int dfNorm[ 19 ];
 
616
 
 
617
                /*! index of inverse dist func, not fast, but useful... */
 
618
                STCON int dfInv[ 19 ];
 
619
 
 
620
                /*! index of x reflected dist func for free slip, not valid for all DFs... */
 
621
                STCON int dfRefX[ 19 ];
 
622
                /*! index of x reflected dist func for free slip, not valid for all DFs... */
 
623
                STCON int dfRefY[ 19 ];
 
624
                /*! index of x reflected dist func for free slip, not valid for all DFs... */
 
625
                STCON int dfRefZ[ 19 ];
 
626
 
 
627
                /*! dist func vectors */
 
628
                STCON int dfVecX[ 27 ];
 
629
                STCON int dfVecY[ 27 ];
 
630
                STCON int dfVecZ[ 27 ];
 
631
 
 
632
                /*! arrays as before with doubles */
 
633
                STCON LbmFloat dfDvecX[ 27 ];
 
634
                STCON LbmFloat dfDvecY[ 27 ];
 
635
                STCON LbmFloat dfDvecZ[ 27 ];
 
636
 
 
637
                /*! principal directions */
 
638
                STCON int princDirX[ 2*3 ];
 
639
                STCON int princDirY[ 2*3 ];
 
640
                STCON int princDirZ[ 2*3 ];
 
641
 
 
642
                /*! vector lengths */
 
643
                STCON LbmFloat dfLength[ 19 ];
 
644
 
 
645
                /*! equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
 
646
                static LbmFloat dfEquil[ 19 ];
 
647
 
 
648
                /*! arrays for les model coefficients */
 
649
                static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ];
 
650
                static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ];
 
651
 
 
652
}; // LbmData3D
 
653
 
 
654
 
 
655
 
 
656
/*****************************************************************************/
 
657
//! class for solver templating - 2D implementation 
 
658
class LbmD2Q9 {
 
659
        
 
660
        public:
 
661
 
 
662
                // constructor, init interface
 
663
                LbmD2Q9() {};
 
664
                // virtual destructor 
 
665
                virtual ~LbmD2Q9() {};
 
666
                //! id string of solver
 
667
                string getIdString() { return string("2D"); }
 
668
 
 
669
                //! how many dimensions?
 
670
                STCON int cDimension;
 
671
 
 
672
                //! Wi factors for collide step 
 
673
                STCON LbmFloat cCollenZero;
 
674
                STCON LbmFloat cCollenOne;
 
675
                STCON LbmFloat cCollenSqrtTwo;
 
676
 
 
677
                //! threshold value for filled/emptied cells 
 
678
                STCON LbmFloat cMagicNr2;
 
679
                STCON LbmFloat cMagicNr2Neg;
 
680
                STCON LbmFloat cMagicNr;
 
681
                STCON LbmFloat cMagicNrNeg;
 
682
 
 
683
                //! size of a single set of distribution functions 
 
684
                STCON int    cDfNum;
 
685
                STCON int    cDirNum;
 
686
 
 
687
                //! distribution functions directions 
 
688
                typedef enum {
 
689
                         cDirInv=  -1,
 
690
                         cDirC  =  0,
 
691
                         cDirN  =  1,
 
692
                         cDirS  =  2,
 
693
                         cDirE  =  3,
 
694
                         cDirW  =  4,
 
695
                         cDirNE =  5,
 
696
                         cDirNW =  6,
 
697
                         cDirSE =  7,
 
698
                         cDirSW =  8
 
699
                } dfDir;
 
700
 
 
701
                /* Vector Order 2D:
 
702
                 * 0  1 2  3  4  5  6 7  8
 
703
                 * 0, 0,0, 1,-1, 1,-1,1,-1 
 
704
                 * 0, 1,-1, 0,0, 1,1,-1,-1  */
 
705
 
 
706
                /* name of the dist. function 
 
707
                         only for nicer output */
 
708
                STCON char* dfString[ 9 ];
 
709
 
 
710
                /* index of normal dist func, not used so far?... */
 
711
                STCON int dfNorm[ 9 ];
 
712
 
 
713
                /* index of inverse dist func, not fast, but useful... */
 
714
                STCON int dfInv[ 9 ];
 
715
 
 
716
                /* index of x reflected dist func for free slip, not valid for all DFs... */
 
717
                STCON int dfRefX[ 9 ];
 
718
                /* index of x reflected dist func for free slip, not valid for all DFs... */
 
719
                STCON int dfRefY[ 9 ];
 
720
                /* index of x reflected dist func for free slip, not valid for all DFs... */
 
721
                STCON int dfRefZ[ 9 ];
 
722
 
 
723
                /* dist func vectors */
 
724
                STCON int dfVecX[ 9 ];
 
725
                STCON int dfVecY[ 9 ];
 
726
                /* Z, 2D values are all 0! */
 
727
                STCON int dfVecZ[ 9 ];
 
728
 
 
729
                /* arrays as before with doubles */
 
730
                STCON LbmFloat dfDvecX[ 9 ];
 
731
                STCON LbmFloat dfDvecY[ 9 ];
 
732
                /* Z, 2D values are all 0! */
 
733
                STCON LbmFloat dfDvecZ[ 9 ];
 
734
 
 
735
                /*! principal directions */
 
736
                STCON int princDirX[ 2*2 ];
 
737
                STCON int princDirY[ 2*2 ];
 
738
                STCON int princDirZ[ 2*2 ];
 
739
 
 
740
                /* vector lengths */
 
741
                STCON LbmFloat dfLength[ 9 ];
 
742
 
 
743
                /* equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */
 
744
                static LbmFloat dfEquil[ 9 ];
 
745
 
 
746
                /*! arrays for les model coefficients */
 
747
                static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ];
 
748
                static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ];
 
749
 
 
750
}; // LbmData3D
 
751
 
 
752
 
 
753
 
 
754
// lbmdimensions
 
755
 
 
756
// not needed hereafter
 
757
#undef STCON
 
758
 
 
759
 
 
760
 
 
761
/*****************************************************************************/
 
762
//! class for solver templating - lbgk (srt) model implementation 
 
763
template<class DQ>
 
764
class LbmModelLBGK : public DQ , public LbmSolverInterface {
 
765
        public:
 
766
 
 
767
                /*! type for cells contents, needed for cell id interface */
 
768
                typedef DQ LbmCellContents;
 
769
                /*! type for cells */
 
770
                typedef LbmCellTemplate< LbmCellContents > LbmCell;
 
771
 
 
772
                // constructor
 
773
                LbmModelLBGK() : DQ(), LbmSolverInterface() {};
 
774
                // virtual destructor 
 
775
                virtual ~LbmModelLBGK() {};
 
776
                //! id string of solver
 
777
                string getIdString() { return DQ::getIdString() + string("lbgk]"); }
 
778
 
 
779
                /*! calculate length of velocity vector */
 
780
                static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz) {
 
781
                        return ((ux)*DQ::dfDvecX[l]+(uy)*DQ::dfDvecY[l]+(uz)*DQ::dfDvecZ[l]);
 
782
                };
 
783
 
 
784
                /*! calculate equilibrium DF for given values */
 
785
                static inline LbmFloat getCollideEq(int l, LbmFloat rho,  LbmFloat ux, LbmFloat uy, LbmFloat uz) {
 
786
                        LbmFloat tmp = getVelVecLen(l,ux,uy,uz); 
 
787
                        return( DQ::dfLength[l] *( 
 
788
                                                + rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) 
 
789
                                                + 3.0 *tmp 
 
790
                                                + 9.0/2.0 *(tmp*tmp) )
 
791
                                        );
 
792
                };
 
793
 
 
794
                
 
795
                /*! relaxation LES functions */
 
796
                inline LbmFloat getLesNoneqTensorCoeff(
 
797
                                LbmFloat df[],                          
 
798
                                LbmFloat feq[] ) {
 
799
                        LbmFloat Qo = 0.0;
 
800
                        for(int m=0; m< ((DQ::cDimension*DQ::cDimension)-DQ::cDimension)/2 ; m++) { 
 
801
                                LbmFloat qadd = 0.0;
 
802
                                for(int l=1; l<DQ::cDfNum; l++) { 
 
803
                                        if(DQ::lesCoeffOffdiag[m][l]==0.0) continue;
 
804
                                        qadd += DQ::lesCoeffOffdiag[m][l]*(df[l]-feq[l]);
 
805
                                }
 
806
                                Qo += (qadd*qadd);
 
807
                        }
 
808
                        Qo *= 2.0; // off diag twice
 
809
                        for(int m=0; m<DQ::cDimension; m++) { 
 
810
                                LbmFloat qadd = 0.0;
 
811
                                for(int l=1; l<DQ::cDfNum; l++) { 
 
812
                                        if(DQ::lesCoeffDiag[m][l]==0.0) continue;
 
813
                                        qadd += DQ::lesCoeffDiag[m][l]*(df[l]-feq[l]);
 
814
                                }
 
815
                                Qo += (qadd*qadd);
 
816
                        }
 
817
                        Qo = sqrt(Qo);
 
818
                        return Qo;
 
819
                }
 
820
                inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo) {
 
821
                        const LbmFloat tau = 1.0/omega;
 
822
                        const LbmFloat nu = (2.0*tau-1.0) * (1.0/6.0);
 
823
                        const LbmFloat C = csmago;
 
824
                        const LbmFloat Csqr = C*C;
 
825
                        LbmFloat S = -nu + sqrt( nu*nu + 18.0*Csqr*Qo ) / (6.0*Csqr);
 
826
                        return( 1.0/( 3.0*( nu+Csqr*S ) +0.5 ) );
 
827
                }
 
828
 
 
829
                // "normal" collision
 
830
                inline void collideArrays(LbmFloat df[],                                
 
831
                                LbmFloat &outrho, // out only!
 
832
                                // velocity modifiers (returns actual velocity!)
 
833
                                LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 
 
834
                                LbmFloat omega, LbmFloat csmago, 
 
835
                                LbmFloat *newOmegaRet, LbmFloat *newQoRet
 
836
                        ) {
 
837
                        LbmFloat rho=df[0]; 
 
838
                        LbmFloat ux = mux;
 
839
                        LbmFloat uy = muy;
 
840
                        LbmFloat uz = muz; 
 
841
                        for(int l=1; l<DQ::cDfNum; l++) { 
 
842
                                rho += df[l]; 
 
843
                                ux  += (DQ::dfDvecX[l]*df[l]); 
 
844
                                uy  += (DQ::dfDvecY[l]*df[l]);  
 
845
                                uz  += (DQ::dfDvecZ[l]*df[l]);  
 
846
                        }  
 
847
                        LbmFloat feq[19];
 
848
                        for(int l=0; l<DQ::cDfNum; l++) { 
 
849
                                feq[l] = getCollideEq(l,rho,ux,uy,uz); 
 
850
                        }
 
851
 
 
852
                        LbmFloat omegaNew;
 
853
                        LbmFloat Qo = 0.0;
 
854
                        if(csmago>0.0) {
 
855
                                Qo = getLesNoneqTensorCoeff(df,feq);
 
856
                                omegaNew = getLesOmega(omega,csmago,Qo);
 
857
                        } else {
 
858
                                omegaNew = omega; // smago off...
 
859
                        }
 
860
                        if(newOmegaRet) *newOmegaRet = omegaNew; // return value for stats
 
861
                        if(newQoRet)    *newQoRet = Qo; // return value of non-eq. stress tensor
 
862
 
 
863
                        for(int l=0; l<DQ::cDfNum; l++) { 
 
864
                                df[l] = (1.0-omegaNew ) * df[l] + omegaNew * feq[l]; 
 
865
                        }  
 
866
 
 
867
                        mux = ux;
 
868
                        muy = uy;
 
869
                        muz = uz;
 
870
                        outrho = rho;
 
871
                };
 
872
 
 
873
}; // LBGK
 
874
 
 
875
#ifdef LBMMODEL_DEFINED
 
876
// force compiler error!
 
877
ERROR - Dont include several LBM models at once...
 
878
#endif
 
879
#define LBMMODEL_DEFINED 1
 
880
 
 
881
 
 
882
typedef LbmModelLBGK<  LbmD2Q9 > LbmBGK2D;
 
883
typedef LbmModelLBGK< LbmD3Q19 > LbmBGK3D;
 
884
 
 
885
 
 
886
// helper function to create consistent grid resolutions
 
887
void initGridSizes(int &mSizex, int &mSizey, int &mSizez,
 
888
                ntlVec3Gfx &mvGeoStart, ntlVec3Gfx &mvGeoEnd, 
 
889
                int mMaxRefine, bool parallel);
 
890
void calculateMemreqEstimate(int resx,int resy,int resz, int refine,
 
891
                double *reqret, string *reqstr);
 
892
 
 
893
//! helper function to convert flag to string (for debuggin)
 
894
string convertCellFlagType2String( CellFlagType flag );
 
895
string convertSingleFlag2String(CellFlagType cflag);
 
896
 
 
897
#endif // LBMINTERFACE_H