|
Blender
V2.59
|
00001 00004 /****************************************************************************** 00005 * 00006 * El'Beem - the visual lattice boltzmann freesurface simulator 00007 * All code distributed as part of El'Beem is covered by the version 2 of the 00008 * GNU General Public License. See the file COPYING for details. 00009 * Copyright 2003-2006 Nils Thuerey 00010 * 00011 * Combined 2D/3D Lattice Boltzmann standard solver classes 00012 * 00013 *****************************************************************************/ 00014 00015 00016 #ifndef LBM_SOLVERCLASS_H 00017 #define LBM_SOLVERCLASS_H 00018 00019 #include "utilities.h" 00020 #include "solver_interface.h" 00021 #include "ntl_ray.h" 00022 #include <stdio.h> 00023 00024 #if PARALLEL==1 00025 #include <omp.h> 00026 #endif // PARALLEL=1 00027 #ifndef PARALLEL 00028 #define PARALLEL 0 00029 #endif // PARALLEL 00030 00031 00032 // general solver setting defines 00033 00035 // might be enabled by compilation 00036 #ifndef FSGR_STRICT_DEBUG 00037 #define FSGR_STRICT_DEBUG 0 00038 #endif // FSGR_STRICT_DEBUG 00039 00041 #define FSGR_OMEGA_DEBUG 0 00042 00044 #define USE_LES 1 00045 00047 //#define TIMEINTORDER 0 00048 // TODO remove interpol t param, also interTime 00049 00050 // use optimized 3D code? 00051 #if LBMDIM==2 00052 #define OPT3D 0 00053 #else 00054 // determine with debugging... 00055 # if FSGR_STRICT_DEBUG==1 00056 # define OPT3D 0 00057 # else // FSGR_STRICT_DEBUG==1 00058 // usually switch optimizations for 3d on, when not debugging 00059 # define OPT3D 1 00060 // COMPRT 00061 //# define OPT3D 0 00062 # endif // FSGR_STRICT_DEBUG==1 00063 #endif 00064 00066 #define MASS_INVALID -1000.0 00067 00068 // empty/fill cells without fluid/empty NB's by inserting them into the full/empty lists? 00069 #define FSGR_LISTTRICK 1 00070 #define FSGR_LISTTTHRESHEMPTY 0.10 00071 #define FSGR_LISTTTHRESHFULL 0.90 00072 #define FSGR_MAGICNR 0.025 00073 //0.04 00074 00076 #define FSGR_MAXNOOFLEVELS 5 00077 00078 // enable/disable fine grid compression for finest level 00079 // make sure this is same as useGridComp in calculateMemreqEstimate 00080 #if LBMDIM==3 00081 #define COMPRESSGRIDS 1 00082 #else 00083 #define COMPRESSGRIDS 0 00084 #endif 00085 00086 // helper for comparing floats with epsilon 00087 #define GFX_FLOATNEQ(x,y) ( ABS((x)-(y)) > (VECTOR_EPSILON) ) 00088 #define LBM_FLOATNEQ(x,y) ( ABS((x)-(y)) > (10.0*LBM_EPSILON) ) 00089 00090 00091 // macros for loops over all DFs 00092 #define FORDF0 for(int l= 0; l< LBM_DFNUM; ++l) 00093 #define FORDF1 for(int l= 1; l< LBM_DFNUM; ++l) 00094 // and with different loop var to prevent shadowing 00095 #define FORDF0M for(int m= 0; m< LBM_DFNUM; ++m) 00096 #define FORDF1M for(int m= 1; m< LBM_DFNUM; ++m) 00097 00098 // iso value defines 00099 // border for marching cubes 00100 #define ISOCORR 3 00101 00102 #define LBM_INLINED inline 00103 00104 // sirdude fix for solaris 00105 #if !defined(linux) && defined(sun) 00106 #include "ieeefp.h" 00107 #ifndef expf 00108 #define expf(x) exp((double)(x)) 00109 #endif 00110 #endif 00111 00112 #include "solver_control.h" 00113 00114 #if LBM_INCLUDE_TESTSOLVERS==1 00115 #include "solver_test.h" 00116 #endif // LBM_INCLUDE_TESTSOLVERS==1 00117 00118 /*****************************************************************************/ 00120 class UniformFsgrCellIdentifier : 00121 public CellIdentifierInterface , public LbmCellContents 00122 { 00123 public: 00125 int level; 00127 int x,y,z; 00128 00130 UniformFsgrCellIdentifier() : 00131 x(0), y(0), z(0) { }; 00132 00133 // implement CellIdentifierInterface 00134 virtual string getAsString() { 00135 std::ostringstream ret; 00136 ret <<"{ i"<<x<<",j"<<y; 00137 if(LBMDIM>2) ret<<",k"<<z; 00138 ret <<" }"; 00139 return ret.str(); 00140 } 00141 00142 virtual bool equal(CellIdentifierInterface* other) { 00143 UniformFsgrCellIdentifier *cid = (UniformFsgrCellIdentifier *)( other ); 00144 if(!cid) return false; 00145 if( x==cid->x && y==cid->y && z==cid->z && level==cid->level ) return true; 00146 return false; 00147 } 00148 }; 00149 00151 class FsgrLevelData { 00152 public: 00153 int id; // level number 00154 00156 LbmFloat nodeSize; 00158 LbmFloat simCellSize; 00160 LbmFloat omega; 00162 LbmFloat time; 00164 LbmFloat timestep; 00166 int lsteps; 00168 LbmVec gravity; 00170 LbmFloat *mprsCells[2]; 00171 CellFlagType *mprsFlags[2]; 00172 00174 LbmFloat lcsmago; 00175 LbmFloat lcsmago_sqr; 00176 LbmFloat lcnu; 00177 00178 // LES statistics per level 00179 double avgOmega; 00180 double avgOmegaCnt; 00181 00183 int setCurr; 00185 int setOther; 00186 00188 LbmFloat lmass; 00189 LbmFloat lvolume; 00190 LbmFloat lcellfactor; 00191 00193 int lSizex, lSizey, lSizez; 00194 int lOffsx, lOffsy, lOffsz; 00195 00196 }; 00197 00198 00199 00200 /*****************************************************************************/ 00202 class LbmFsgrSolver : 00203 public LbmSolverInterface // this means, the solver is a lbmData object and implements the lbmInterface 00204 { 00205 00206 public: 00208 LbmFsgrSolver(); 00210 virtual ~LbmFsgrSolver(); 00211 00213 virtual void parseAttrList(); 00215 void initLevelOmegas(); 00216 00217 // multi step solver init 00219 virtual bool initializeSolverMemory(); 00221 virtual bool initializeSolverGrids(); 00223 virtual bool initializeSolverPostinit(); 00224 00226 virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename); 00227 00228 # if LBM_USE_GUI==1 00229 00230 virtual void debugDisplay(int set); 00231 # endif 00232 00233 // implement CellIterator<UniformFsgrCellIdentifier> interface 00234 typedef UniformFsgrCellIdentifier stdCellId; 00235 virtual CellIdentifierInterface* getFirstCell( ); 00236 virtual void advanceCell( CellIdentifierInterface* ); 00237 virtual bool noEndCell( CellIdentifierInterface* ); 00238 virtual void deleteCellIterator( CellIdentifierInterface** ); 00239 virtual CellIdentifierInterface* getCellAt( ntlVec3Gfx pos ); 00240 virtual int getCellSet ( CellIdentifierInterface* ); 00241 virtual ntlVec3Gfx getCellOrigin ( CellIdentifierInterface* ); 00242 virtual ntlVec3Gfx getCellSize ( CellIdentifierInterface* ); 00243 virtual int getCellLevel ( CellIdentifierInterface* ); 00244 virtual LbmFloat getCellDensity ( CellIdentifierInterface* ,int set); 00245 virtual LbmVec getCellVelocity ( CellIdentifierInterface* ,int set); 00246 virtual LbmFloat getCellDf ( CellIdentifierInterface* ,int set, int dir); 00247 virtual LbmFloat getCellMass ( CellIdentifierInterface* ,int set); 00248 virtual LbmFloat getCellFill ( CellIdentifierInterface* ,int set); 00249 virtual CellFlagType getCellFlag ( CellIdentifierInterface* ,int set); 00250 virtual LbmFloat getEquilDf ( int ); 00251 virtual ntlVec3Gfx getVelocityAt (float x, float y, float z); 00252 // convert pointers 00253 stdCellId* convertBaseCidToStdCid( CellIdentifierInterface* basecid); 00254 00256 bool initGeometryFlags(); 00258 void initFreeSurfaces(); 00260 void initStandingFluidGradient(); 00261 00263 LBM_INLINED void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass); 00264 LBM_INLINED void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel); 00265 LBM_INLINED void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag); 00266 LBM_INLINED void forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag); 00267 LBM_INLINED void initInterfaceVars(int level, int i,int j,int k,int workSet, bool initMass); 00269 void interpolateCellValues(int level,int ei,int ej,int ek,int workSet, LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz); 00270 00272 void stepMain(); 00274 void fineAdvance(); 00276 void coarseAdvance(int lev); 00278 void coarseCalculateFluxareas(int lev); 00279 // adaptively coarsen grid 00280 bool adaptGrid(int lev); 00281 // restrict fine grid DFs to coarse grid 00282 void coarseRestrictFromFine(int lev); 00283 00284 /* simulation object interface, just calls stepMain */ 00285 virtual void step(); 00287 virtual int initParticles(); 00289 virtual void advanceParticles(); 00291 void handleObstacleParticle(ParticleObject *p); 00298 virtual vector<ntlGeometryObject*> getDebugObjects(); 00299 00300 // gui/output debugging functions 00301 # if LBM_USE_GUI==1 00302 virtual void debugDisplayNode(int dispset, CellIdentifierInterface* cell ); 00303 virtual void lbmDebugDisplay(int dispset); 00304 virtual void lbmMarkedCellDisplay(); 00305 # endif // LBM_USE_GUI==1 00306 virtual void debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet=-1); 00307 00309 void prepareVisualization( void ); 00310 00311 /* surface generation settings */ 00312 virtual void setSurfGenSettings(short value); 00313 00314 protected: 00315 00317 void printLbmCell(int level, int i, int j, int k,int set); 00318 // debugging use CellIterator interface to mark cell 00319 void debugMarkCellCall(int level, int vi,int vj,int vk); 00320 00321 // loop over grid, stream&collide update 00322 void mainLoop(int lev); 00323 // change time step size 00324 void adaptTimestep(); 00326 void recalculateObjectSpeeds(); 00328 void initMovingObstacles(bool staticInit); 00330 void reinitFlags( int workSet ); 00332 LbmFloat getMassdWeight(bool dirForw, int i,int j,int k,int workSet, int l); 00334 void computeFluidSurfaceNormal(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret); 00335 void computeFluidSurfaceNormalAcc(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret); 00336 void computeObstacleSurfaceNormal(LbmFloat *cell, CellFlagType *cellflag, LbmFloat *snret); 00337 void computeObstacleSurfaceNormalAcc(int i,int j,int k, LbmFloat *snret); 00339 LBM_INLINED void addToNewInterList( int ni, int nj, int nk ); 00341 void interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet,bool markNbs); 00342 void coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet); 00343 00345 LBM_INLINED int getForZMinBnd(); 00346 LBM_INLINED int getForZMin1(); 00347 LBM_INLINED int getForZMaxBnd(int lev); 00348 LBM_INLINED int getForZMax1(int lev); 00349 LBM_INLINED bool checkDomainBounds(int lev,int i,int j,int k); 00350 LBM_INLINED bool checkDomainBoundsPos(int lev,LbmVec pos); 00351 00352 // touch grid and flags once 00353 void preinitGrids(); 00354 // one relaxation step for standing fluid 00355 void standingFluidPreinit(); 00356 00357 00358 // member vars 00359 00361 LbmFloat mCurrentMass; 00362 LbmFloat mCurrentVolume; 00363 LbmFloat mInitialMass; 00364 00366 int mNumProblems; 00367 00368 // average mlsups, count how many so far... 00369 double mAvgMLSUPS; 00370 double mAvgMLSUPSCnt; 00371 00373 IsoSurface *mpPreviewSurface; 00374 00376 bool mTimeAdap; 00378 bool mForceTimeStepReduce; 00379 00381 LbmFloat mFVHeight; 00382 LbmFloat mFVArea; 00383 bool mUpdateFVHeight; 00384 00386 LbmFloat mGfxEndTime; 00388 int mInitSurfaceSmoothing; 00390 // each flag switches side on off, fssgNoObs is for obstacle sides 00391 // -1 equals all on 00392 typedef enum { 00393 fssgNormal = 0, 00394 fssgNoNorth = 1, 00395 fssgNoSouth = 2, 00396 fssgNoEast = 4, 00397 fssgNoWest = 8, 00398 fssgNoTop = 16, 00399 fssgNoBottom = 32, 00400 fssgNoObs = 64 00401 } fsSurfaceGen; 00402 int mFsSurfGenSetting; 00403 00405 int mTimestepReduceLock; 00407 int mTimeSwitchCounts; 00408 // only switch of maxvel is higher for several steps... 00409 int mTimeMaxvelStepCnt; 00410 00412 LbmFloat mSimulationTime, mLastSimTime; 00414 LbmFloat mMinTimestep, mMaxTimestep; 00416 LbmFloat mMxvx, mMxvy, mMxvz, mMaxVlen; 00417 00419 vector<LbmPoint> mListEmpty; 00421 vector<LbmPoint> mListFull; 00423 vector<LbmPoint> mListNewInter; 00425 class lbmFloatSet { 00426 public: 00427 LbmFloat val[dTotalNum]; 00428 LbmFloat numNbs; 00429 }; 00431 LbmVec mDvecNrm[27]; 00432 00433 00435 bool checkSymmetry(string idstring); 00437 int mMaxNoCells, mMinNoCells; 00438 LONGINT mAvgNumUsedCells; 00439 00441 vector<LbmVec> mObjectSpeeds; 00443 vector<LbmFloat> mObjectPartslips; 00445 vector<LbmFloat> mObjectMassMovnd; 00446 00448 vector<ntlVec3Gfx> mMOIVertices; 00449 vector<ntlVec3Gfx> mMOIVerticesOld; 00450 vector<ntlVec3Gfx> mMOINormals; 00451 00453 int mIsoWeightMethod; 00454 float mIsoWeight[27]; 00455 00456 // grid coarsening vars 00457 00459 FsgrLevelData mLevel[FSGR_MAXNOOFLEVELS]; 00460 00462 int mMaxRefine; 00463 00465 LbmFloat mDfScaleUp, mDfScaleDown; 00466 00468 LbmFloat mFsgrCellArea[27]; 00470 LbmFloat mGaussw[27]; 00471 00473 float mInitialCsmago; 00475 LbmFloat mDebugOmegaRet; 00477 LbmFloat mLastOmega; 00478 LbmVec mLastGravity; 00479 00481 int mNumInterdCells; 00482 int mNumInvIfCells; 00483 int mNumInvIfTotal; 00484 int mNumFsgrChanges; 00485 00487 int mDisableStandingFluidInit; 00489 bool mInit2dYZ; 00491 int mForceTadapRefine; 00493 int mCutoff; 00494 00495 // strict debug interface 00496 # if FSGR_STRICT_DEBUG==1 00497 int debLBMGI(int level, int ii,int ij,int ik, int is); 00498 CellFlagType& debRFLAG(int level, int xx,int yy,int zz,int set); 00499 CellFlagType& debRFLAG_NB(int level, int xx,int yy,int zz,int set, int dir); 00500 CellFlagType& debRFLAG_NBINV(int level, int xx,int yy,int zz,int set, int dir); 00501 int debLBMQI(int level, int ii,int ij,int ik, int is, int l); 00502 LbmFloat& debQCELL(int level, int xx,int yy,int zz,int set,int l); 00503 LbmFloat& debQCELL_NB(int level, int xx,int yy,int zz,int set, int dir,int l); 00504 LbmFloat& debQCELL_NBINV(int level, int xx,int yy,int zz,int set, int dir,int l); 00505 LbmFloat* debRACPNT(int level, int ii,int ij,int ik, int is ); 00506 LbmFloat& debRAC(LbmFloat* s,int l); 00507 # endif // FSGR_STRICT_DEBUG==1 00508 00509 LbmControlData *mpControl; 00510 00511 void initCpdata(); 00512 void handleCpdata(); 00513 void cpDebugDisplay(int dispset); 00514 00515 bool mUseTestdata; 00516 # if LBM_INCLUDE_TESTSOLVERS==1 00517 // test functions 00518 LbmTestdata *mpTest; 00519 void initTestdata(); 00520 void destroyTestdata(); 00521 void handleTestdata(); 00522 void set3dHeight(int ,int ); 00523 00524 int mMpNum,mMpIndex; 00525 int mOrgSizeX; 00526 LbmFloat mOrgStartX; 00527 LbmFloat mOrgEndX; 00528 void mrSetup(); 00529 void mrExchange(); 00530 void mrIsoExchange(); 00531 LbmFloat mrInitTadap(LbmFloat max); 00532 void gcFillBuffer( LbmGridConnector *gc, int *retSizeCnt, const int *bdfs); 00533 void gcUnpackBuffer(LbmGridConnector *gc, int *retSizeCnt, const int *bdfs); 00534 public: 00535 // needed for testdata 00536 void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz); 00537 // mptest 00538 int getMpIndex() { return mMpIndex; }; 00539 # endif // LBM_INCLUDE_TESTSOLVERS==1 00540 00541 // former LbmModelLBGK functions 00542 // relaxation funtions - implemented together with relax macros 00543 static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz); 00544 static inline LbmFloat getCollideEq(int l, LbmFloat rho, LbmFloat ux, LbmFloat uy, LbmFloat uz); 00545 inline LbmFloat getLesNoneqTensorCoeff( LbmFloat df[], LbmFloat feq[] ); 00546 inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo); 00547 inline void collideArrays( int lev, int i, int j, int k, // position - more for debugging 00548 LbmFloat df[], LbmFloat &outrho, // out only! 00549 // velocity modifiers (returns actual velocity!) 00550 LbmFloat &mux, LbmFloat &muy, LbmFloat &muz, 00551 LbmFloat omega, LbmVec gravity, LbmFloat csmago, 00552 LbmFloat *newOmegaRet, LbmFloat *newQoRet); 00553 00554 00555 // former LBM models 00557 # define STCON static const 00558 00559 # if LBMDIM==3 00560 00562 virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D3Q19]"); } 00563 00565 STCON int cDimension; 00566 00567 // Wi factors for collide step 00568 STCON LbmFloat cCollenZero; 00569 STCON LbmFloat cCollenOne; 00570 STCON LbmFloat cCollenSqrtTwo; 00571 00573 STCON LbmFloat cMagicNr2; 00574 STCON LbmFloat cMagicNr2Neg; 00575 STCON LbmFloat cMagicNr; 00576 STCON LbmFloat cMagicNrNeg; 00577 00579 STCON int cDfNum; 00581 STCON int cDirNum; 00582 00584 typedef enum { 00585 cDirInv= -1, 00586 cDirC = 0, 00587 cDirN = 1, 00588 cDirS = 2, 00589 cDirE = 3, 00590 cDirW = 4, 00591 cDirT = 5, 00592 cDirB = 6, 00593 cDirNE = 7, 00594 cDirNW = 8, 00595 cDirSE = 9, 00596 cDirSW = 10, 00597 cDirNT = 11, 00598 cDirNB = 12, 00599 cDirST = 13, 00600 cDirSB = 14, 00601 cDirET = 15, 00602 cDirEB = 16, 00603 cDirWT = 17, 00604 cDirWB = 18 00605 } dfDir; 00606 00607 /* Vector Order 3D: 00608 * 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 00609 * 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 00610 * 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 00611 * 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 00612 */ 00613 00616 STCON char* dfString[ 19 ]; 00617 00619 STCON int dfNorm[ 19 ]; 00620 00622 STCON int dfInv[ 19 ]; 00623 00625 STCON int dfRefX[ 19 ]; 00627 STCON int dfRefY[ 19 ]; 00629 STCON int dfRefZ[ 19 ]; 00630 00632 STCON int dfVecX[ 27 ]; 00633 STCON int dfVecY[ 27 ]; 00634 STCON int dfVecZ[ 27 ]; 00635 00637 STCON LbmFloat dfDvecX[ 27 ]; 00638 STCON LbmFloat dfDvecY[ 27 ]; 00639 STCON LbmFloat dfDvecZ[ 27 ]; 00640 00642 STCON int princDirX[ 2*3 ]; 00643 STCON int princDirY[ 2*3 ]; 00644 STCON int princDirZ[ 2*3 ]; 00645 00647 STCON LbmFloat dfLength[ 19 ]; 00648 00650 static LbmFloat dfEquil[ dTotalNum ]; 00651 00653 static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ]; 00654 static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ]; 00655 00656 # else // end LBMDIM==3 , LBMDIM==2 00657 00659 virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D2Q9]"); } 00660 00662 STCON int cDimension; 00663 00665 STCON LbmFloat cCollenZero; 00666 STCON LbmFloat cCollenOne; 00667 STCON LbmFloat cCollenSqrtTwo; 00668 00670 STCON LbmFloat cMagicNr2; 00671 STCON LbmFloat cMagicNr2Neg; 00672 STCON LbmFloat cMagicNr; 00673 STCON LbmFloat cMagicNrNeg; 00674 00676 STCON int cDfNum; 00677 STCON int cDirNum; 00678 00680 typedef enum { 00681 cDirInv= -1, 00682 cDirC = 0, 00683 cDirN = 1, 00684 cDirS = 2, 00685 cDirE = 3, 00686 cDirW = 4, 00687 cDirNE = 5, 00688 cDirNW = 6, 00689 cDirSE = 7, 00690 cDirSW = 8 00691 } dfDir; 00692 00693 /* Vector Order 2D: 00694 * 0 1 2 3 4 5 6 7 8 00695 * 0, 0,0, 1,-1, 1,-1,1,-1 00696 * 0, 1,-1, 0,0, 1,1,-1,-1 */ 00697 00698 /* name of the dist. function 00699 only for nicer output */ 00700 STCON char* dfString[ 9 ]; 00701 00702 /* index of normal dist func, not used so far?... */ 00703 STCON int dfNorm[ 9 ]; 00704 00705 /* index of inverse dist func, not fast, but useful... */ 00706 STCON int dfInv[ 9 ]; 00707 00708 /* index of x reflected dist func for free slip, not valid for all DFs... */ 00709 STCON int dfRefX[ 9 ]; 00710 /* index of x reflected dist func for free slip, not valid for all DFs... */ 00711 STCON int dfRefY[ 9 ]; 00712 /* index of x reflected dist func for free slip, not valid for all DFs... */ 00713 STCON int dfRefZ[ 9 ]; 00714 00715 /* dist func vectors */ 00716 STCON int dfVecX[ 9 ]; 00717 STCON int dfVecY[ 9 ]; 00718 /* Z, 2D values are all 0! */ 00719 STCON int dfVecZ[ 9 ]; 00720 00721 /* arrays as before with doubles */ 00722 STCON LbmFloat dfDvecX[ 9 ]; 00723 STCON LbmFloat dfDvecY[ 9 ]; 00724 /* Z, 2D values are all 0! */ 00725 STCON LbmFloat dfDvecZ[ 9 ]; 00726 00728 STCON int princDirX[ 2*2 ]; 00729 STCON int princDirY[ 2*2 ]; 00730 STCON int princDirZ[ 2*2 ]; 00731 00732 /* vector lengths */ 00733 STCON LbmFloat dfLength[ 9 ]; 00734 00735 /* equilibrium distribution functions, precalculated = getCollideEq(i, 0,0,0,0) */ 00736 static LbmFloat dfEquil[ dTotalNum ]; 00737 00739 static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ]; 00740 static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ]; 00741 00742 # endif // LBMDIM==2 00743 }; 00744 00745 #undef STCON 00746 00747 00748 /*****************************************************************************/ 00749 // relaxation_macros 00750 00751 00752 00753 // cell mark debugging 00754 #if FSGR_STRICT_DEBUG==10 00755 #define debugMarkCell(lev,x,y,z) \ 00756 errMsg("debugMarkCell",this->mName<<" step: "<<this->mStepCnt<<" lev:"<<(lev)<<" marking "<<PRINT_VEC((x),(y),(z))<<" line "<< __LINE__ ); \ 00757 debugMarkCellCall((lev),(x),(y),(z)); 00758 #else // FSGR_STRICT_DEBUG==1 00759 #define debugMarkCell(lev,x,y,z) \ 00760 debugMarkCellCall((lev),(x),(y),(z)); 00761 #endif // FSGR_STRICT_DEBUG==1 00762 00763 00764 // flag array defines ----------------------------------------------------------------------------------------------- 00765 00766 // lbm testsolver get index define, note - ignores is (set) as flag 00767 // array is only a single entry 00768 #define _LBMGI(level, ii,ij,ik, is) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) ) 00769 00771 #define _RFLAG(level,xx,yy,zz,set) mLevel[level].mprsFlags[set][ LBMGI((level),(xx),(yy),(zz),(set)) ] 00772 #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) ] 00773 #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) ] 00774 00775 // array handling ----------------------------------------------------------------------------------------------- 00776 00777 #define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) ) 00778 #define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)]) 00779 #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)]) 00780 #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)]) 00781 00782 #define QCELLSTEP dTotalNum 00783 #define _RACPNT(level, ii,ij,ik, is ) &QCELL(level,ii,ij,ik,is,0) 00784 #define _RAC(s,l) (s)[(l)] 00785 00786 00787 #if FSGR_STRICT_DEBUG==1 00788 00789 #define LBMGI(level,ii,ij,ik, is) debLBMGI(level,ii,ij,ik, is) 00790 #define RFLAG(level,xx,yy,zz,set) debRFLAG(level,xx,yy,zz,set) 00791 #define RFLAG_NB(level,xx,yy,zz,set, dir) debRFLAG_NB(level,xx,yy,zz,set, dir) 00792 #define RFLAG_NBINV(level,xx,yy,zz,set, dir) debRFLAG_NBINV(level,xx,yy,zz,set, dir) 00793 00794 #define LBMQI(level,ii,ij,ik, is, l) debLBMQI(level,ii,ij,ik, is, l) 00795 #define QCELL(level,xx,yy,zz,set,l) debQCELL(level,xx,yy,zz,set,l) 00796 #define QCELL_NB(level,xx,yy,zz,set, dir,l) debQCELL_NB(level,xx,yy,zz,set, dir,l) 00797 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l) debQCELL_NBINV(level,xx,yy,zz,set, dir,l) 00798 #define RACPNT(level, ii,ij,ik, is ) debRACPNT(level, ii,ij,ik, is ) 00799 #define RAC(s,l) debRAC(s,l) 00800 00801 #else // FSGR_STRICT_DEBUG==1 00802 00803 #define LBMGI(level,ii,ij,ik, is) _LBMGI(level,ii,ij,ik, is) 00804 #define RFLAG(level,xx,yy,zz,set) _RFLAG(level,xx,yy,zz,set) 00805 #define RFLAG_NB(level,xx,yy,zz,set, dir) _RFLAG_NB(level,xx,yy,zz,set, dir) 00806 #define RFLAG_NBINV(level,xx,yy,zz,set, dir) _RFLAG_NBINV(level,xx,yy,zz,set, dir) 00807 00808 #define LBMQI(level,ii,ij,ik, is, l) _LBMQI(level,ii,ij,ik, is, l) 00809 #define QCELL(level,xx,yy,zz,set,l) _QCELL(level,xx,yy,zz,set,l) 00810 #define QCELL_NB(level,xx,yy,zz,set, dir,l) _QCELL_NB(level,xx,yy,zz,set, dir, l) 00811 #define QCELL_NBINV(level,xx,yy,zz,set, dir,l) _QCELL_NBINV(level,xx,yy,zz,set, dir,l) 00812 #define RACPNT(level, ii,ij,ik, is ) _RACPNT(level, ii,ij,ik, is ) 00813 #define RAC(s,l) _RAC(s,l) 00814 00815 #endif // FSGR_STRICT_DEBUG==1 00816 00817 // general defines ----------------------------------------------------------------------------------------------- 00818 00819 // replace TESTFLAG 00820 #define FLAGISEXACT(flag, compflag) ((flag & compflag)==compflag) 00821 00822 #if LBMDIM==2 00823 #define dC 0 00824 #define dN 1 00825 #define dS 2 00826 #define dE 3 00827 #define dW 4 00828 #define dNE 5 00829 #define dNW 6 00830 #define dSE 7 00831 #define dSW 8 00832 #else 00833 // direction indices 00834 #define dC 0 00835 #define dN 1 00836 #define dS 2 00837 #define dE 3 00838 #define dW 4 00839 #define dT 5 00840 #define dB 6 00841 #define dNE 7 00842 #define dNW 8 00843 #define dSE 9 00844 #define dSW 10 00845 #define dNT 11 00846 #define dNB 12 00847 #define dST 13 00848 #define dSB 14 00849 #define dET 15 00850 #define dEB 16 00851 #define dWT 17 00852 #define dWB 18 00853 #endif 00854 //? #define dWB 18 00855 00856 // default init for dFlux values 00857 #define FLUX_INIT 0.5f * (float)(this->cDfNum) 00858 00859 // only for non DF dir handling! 00860 #define dNET 19 00861 #define dNWT 20 00862 #define dSET 21 00863 #define dSWT 22 00864 #define dNEB 23 00865 #define dNWB 24 00866 #define dSEB 25 00867 #define dSWB 26 00868 00870 #define BND_FILL 0.0 00871 00872 #define DFL1 (1.0/ 3.0) 00873 #define DFL2 (1.0/18.0) 00874 #define DFL3 (1.0/36.0) 00875 00876 // loops over _all_ cells (including boundary layer) 00877 #define FSGR_FORIJK_BOUNDS(leveli) \ 00878 for(int k= getForZMinBnd(); k< getForZMaxBnd(leveli); ++k) \ 00879 for(int j=0;j<mLevel[leveli].lSizey-0;++j) \ 00880 for(int i=0;i<mLevel[leveli].lSizex-0;++i) \ 00881 00882 // loops over _only inner_ cells 00883 #define FSGR_FORIJK1(leveli) \ 00884 for(int k= getForZMin1(); k< getForZMax1(leveli); ++k) \ 00885 for(int j=1;j<mLevel[leveli].lSizey-1;++j) \ 00886 for(int i=1;i<mLevel[leveli].lSizex-1;++i) \ 00887 00888 // relaxation_macros end 00889 00890 00891 00892 /******************************************************************************/ 00894 /******************************************************************************/ 00895 00897 inline LbmFloat LbmFsgrSolver::getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz) { 00898 return ((ux)*dfDvecX[l]+(uy)*dfDvecY[l]+(uz)*dfDvecZ[l]); 00899 }; 00900 00902 inline LbmFloat LbmFsgrSolver::getCollideEq(int l, LbmFloat rho, LbmFloat ux, LbmFloat uy, LbmFloat uz) { 00903 #if FSGR_STRICT_DEBUG==1 00904 if((l<0)||(l>LBM_DFNUM)) { errFatal("LbmFsgrSolver::getCollideEq","Invalid DFEQ call "<<l, SIMWORLD_PANIC ); /* no access to mPanic here */ } 00905 #endif // FSGR_STRICT_DEBUG==1 00906 LbmFloat tmp = getVelVecLen(l,ux,uy,uz); 00907 return( dfLength[l] *( 00908 + rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz)) 00909 + 3.0 *tmp 00910 + 9.0/2.0 *(tmp*tmp) ) 00911 ); // */ 00912 }; 00913 00914 /*****************************************************************************/ 00915 /* init a given cell with flag, density, mass and equilibrium dist. funcs */ 00916 00917 void LbmFsgrSolver::forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) { 00918 // also overwrite persisting flags 00919 // function is useful for tracking accesses... 00920 RFLAG(level,xx,yy,zz,set) = newflag; 00921 } 00922 void LbmFsgrSolver::changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) { 00923 CellFlagType pers = RFLAG(level,xx,yy,zz,set) & CFPersistMask; 00924 RFLAG(level,xx,yy,zz,set) = newflag | pers; 00925 } 00926 00927 void 00928 LbmFsgrSolver::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass) { 00929 /* init eq. dist funcs */ 00930 LbmFloat *ecel; 00931 int workSet = mLevel[level].setCurr; 00932 00933 ecel = RACPNT(level, i,j,k, workSet); 00934 FORDF0 { RAC(ecel, l) = this->dfEquil[l] * rho; } 00935 RAC(ecel, dMass) = mass; 00936 RAC(ecel, dFfrac) = mass/rho; 00937 RAC(ecel, dFlux) = FLUX_INIT; 00938 changeFlag(level, i,j,k, workSet, flag); 00939 00940 workSet ^= 1; 00941 changeFlag(level, i,j,k, workSet, flag); 00942 return; 00943 } 00944 00945 void 00946 LbmFsgrSolver::initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel) { 00947 LbmFloat *ecel; 00948 int workSet = mLevel[level].setCurr; 00949 00950 ecel = RACPNT(level, i,j,k, workSet); 00951 FORDF0 { RAC(ecel, l) = getCollideEq(l, rho,vel[0],vel[1],vel[2]); } 00952 RAC(ecel, dMass) = mass; 00953 RAC(ecel, dFfrac) = mass/rho; 00954 RAC(ecel, dFlux) = FLUX_INIT; 00955 changeFlag(level, i,j,k, workSet, flag); 00956 00957 workSet ^= 1; 00958 changeFlag(level, i,j,k, workSet, flag); 00959 return; 00960 } 00961 00962 int LbmFsgrSolver::getForZMinBnd() { 00963 return 0; 00964 } 00965 int LbmFsgrSolver::getForZMin1() { 00966 if(LBMDIM==2) return 0; 00967 return 1; 00968 } 00969 00970 int LbmFsgrSolver::getForZMaxBnd(int lev) { 00971 if(LBMDIM==2) return 1; 00972 return mLevel[lev].lSizez -0; 00973 } 00974 int LbmFsgrSolver::getForZMax1(int lev) { 00975 if(LBMDIM==2) return 1; 00976 return mLevel[lev].lSizez -1; 00977 } 00978 00979 bool LbmFsgrSolver::checkDomainBounds(int lev,int i,int j,int k) { 00980 if(i<0) return false; 00981 if(j<0) return false; 00982 if(k<0) return false; 00983 if(i>mLevel[lev].lSizex-1) return false; 00984 if(j>mLevel[lev].lSizey-1) return false; 00985 if(k>mLevel[lev].lSizez-1) return false; 00986 return true; 00987 } 00988 bool LbmFsgrSolver::checkDomainBoundsPos(int lev,LbmVec pos) { 00989 const int i= (int)pos[0]; 00990 if(i<0) return false; 00991 if(i>mLevel[lev].lSizex-1) return false; 00992 const int j= (int)pos[1]; 00993 if(j<0) return false; 00994 if(j>mLevel[lev].lSizey-1) return false; 00995 const int k= (int)pos[2]; 00996 if(k<0) return false; 00997 if(k>mLevel[lev].lSizez-1) return false; 00998 return true; 00999 } 01000 01001 void LbmFsgrSolver::initInterfaceVars(int level, int i,int j,int k,int workSet, bool initMass) { 01002 LbmFloat *ccel = &QCELL(level ,i,j,k, workSet,0); 01003 LbmFloat nrho = 0.0; 01004 FORDF0 { nrho += RAC(ccel,l); } 01005 if(initMass) { 01006 RAC(ccel,dMass) = nrho; 01007 RAC(ccel, dFfrac) = 1.; 01008 } else { 01009 // preinited, e.g. from reinitFlags 01010 RAC(ccel, dFfrac) = RAC(ccel, dMass)/nrho; 01011 RAC(ccel, dFlux) = FLUX_INIT; 01012 } 01013 } 01014 01015 01016 #endif 01017 01018