intern/elbeem/intern/solver_class.h

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