Page MenuHome
Differential D3744 Diff 11987 intern/elbeem/intern/solver_init.cpp

Changeset View

Standalone View

View Options

intern/elbeem/intern/solver_init.cpp

Context not available.
/*****************************************************************************/ /*****************************************************************************/
//! common variables //! common variables
/*****************************************************************************/ /*****************************************************************************/
/*! 3D implementation D3Q19 */ /*! 3D implementation D3Q19 */
Context not available.
//! how many dimensions? //! how many dimensions?
const int LbmFsgrSolver::cDimension = 3; const int LbmFsgrSolver::cDimension = 3;
// Wi factors for collide step // Wi factors for collide step
const LbmFloat LbmFsgrSolver::cCollenZero = (1.0/3.0); const LbmFloat LbmFsgrSolver::cCollenZero = (1.0/3.0);
const LbmFloat LbmFsgrSolver::cCollenOne = (1.0/18.0); const LbmFloat LbmFsgrSolver::cCollenOne = (1.0/18.0);
const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0); const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
//! threshold value for filled/emptied cells //! threshold value for filled/emptied cells
const LbmFloat LbmFsgrSolver::cMagicNr2 = 1.0005; const LbmFloat LbmFsgrSolver::cMagicNr2 = 1.0005;
const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005; const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
const LbmFloat LbmFsgrSolver::cMagicNr = 1.010001; const LbmFloat LbmFsgrSolver::cMagicNr = 1.010001;
const LbmFloat LbmFsgrSolver::cMagicNrNeg = -0.010001; const LbmFloat LbmFsgrSolver::cMagicNrNeg = -0.010001;
//! size of a single set of distribution functions //! size of a single set of distribution functions
const int LbmFsgrSolver::cDfNum = 19; const int LbmFsgrSolver::cDfNum = 19;
//! direction vector contain vecs for all spatial dirs, even if not used for LBM model //! direction vector contain vecs for all spatial dirs, even if not used for LBM model
const int LbmFsgrSolver::cDirNum = 27; const int LbmFsgrSolver::cDirNum = 27;
//const string LbmFsgrSolver::dfString[ cDfNum ] = { //const string LbmFsgrSolver::dfString[ cDfNum ] = {
const char* LbmFsgrSolver::dfString[ cDfNum ] = { const char* LbmFsgrSolver::dfString[ cDfNum ] = {
" C", " N"," S"," E"," W"," T"," B", " C", " N"," S"," E"," W"," T"," B",
"NE","NW","SE","SW", "NE","NW","SE","SW",
"NT","NB","ST","SB", "NT","NB","ST","SB",
"ET","EB","WT","WB" "ET","EB","WT","WB"
}; };
const int LbmFsgrSolver::dfNorm[ cDfNum ] = { const int LbmFsgrSolver::dfNorm[ cDfNum ] = {
cDirC, cDirN, cDirS, cDirE, cDirW, cDirT, cDirB, cDirC, cDirN, cDirS, cDirE, cDirW, cDirT, cDirB,
cDirNE, cDirNW, cDirSE, cDirSW, cDirNE, cDirNW, cDirSE, cDirSW,
cDirNT, cDirNB, cDirST, cDirSB, cDirNT, cDirNB, cDirST, cDirSB,
cDirET, cDirEB, cDirWT, cDirWB cDirET, cDirEB, cDirWT, cDirWB
}; };
const int LbmFsgrSolver::dfInv[ cDfNum ] = { const int LbmFsgrSolver::dfInv[ cDfNum ] = {
cDirC, cDirS, cDirN, cDirW, cDirE, cDirB, cDirT, cDirC, cDirS, cDirN, cDirW, cDirE, cDirB, cDirT,
cDirSW, cDirSE, cDirNW, cDirNE, cDirSW, cDirSE, cDirNW, cDirNE,
cDirSB, cDirST, cDirNB, cDirNT, cDirSB, cDirST, cDirNB, cDirNT,
cDirWB, cDirWT, cDirEB, cDirET cDirWB, cDirWT, cDirEB, cDirET
}; };
const int LbmFsgrSolver::dfRefX[ cDfNum ] = { const int LbmFsgrSolver::dfRefX[ cDfNum ] = {
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
cDirSE, cDirSW, cDirNE, cDirNW, cDirSE, cDirSW, cDirNE, cDirNW,
0, 0, 0, 0, 0, 0, 0, 0,
cDirEB, cDirET, cDirWB, cDirWT cDirEB, cDirET, cDirWB, cDirWT
}; };
const int LbmFsgrSolver::dfRefY[ cDfNum ] = { const int LbmFsgrSolver::dfRefY[ cDfNum ] = {
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
cDirNW, cDirNE, cDirSW, cDirSE, cDirNW, cDirNE, cDirSW, cDirSE,
cDirNB, cDirNT, cDirSB, cDirST, cDirNB, cDirNT, cDirSB, cDirST,
0, 0, 0, 0 0, 0, 0, 0
}; };
const int LbmFsgrSolver::dfRefZ[ cDfNum ] = { const int LbmFsgrSolver::dfRefZ[ cDfNum ] = {
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
cDirST, cDirSB, cDirNT, cDirNB, cDirST, cDirSB, cDirNT, cDirNB,
cDirWT, cDirWB, cDirET, cDirEB cDirWT, cDirWB, cDirET, cDirEB
}; };
Context not available.
// 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 // 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
// 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 // 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
const int LbmFsgrSolver::dfVecX[ cDirNum ] = { const int LbmFsgrSolver::dfVecX[ cDirNum ] = {
0, 0,0, 1,-1, 0,0, 0, 0,0, 1,-1, 0,0,
1,-1,1,-1, 1,-1,1,-1,
0,0,0,0, 0,0,0,0,
Context not available.
1,-1, 1,-1, 1,-1, 1,-1,
1,-1, 1,-1, 1,-1, 1,-1,
}; };
const int LbmFsgrSolver::dfVecY[ cDirNum ] = { const int LbmFsgrSolver::dfVecY[ cDirNum ] = {
0, 1,-1, 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,
Context not available.
1, 1,-1,-1, 1, 1,-1,-1,
1, 1,-1,-1 1, 1,-1,-1
}; };
const int LbmFsgrSolver::dfVecZ[ cDirNum ] = { const int LbmFsgrSolver::dfVecZ[ cDirNum ] = {
0, 0,0,0,0,1,-1, 0, 0,0,0,0,1,-1,
0,0,0,0, 0,0,0,0,
1,-1,1,-1, 1,-1,1,-1,
Context not available.
}; };
/* principal directions */ /* principal directions */
const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = { const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = {
1,-1, 0,0, 0,0 1,-1, 0,0, 0,0
}; };
const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = { const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = {
0,0, 1,-1, 0,0 0,0, 1,-1, 0,0
}; };
const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = { const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = {
0,0, 0,0, 1,-1 0,0, 0,0, 1,-1
}; };
Context not available.
LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ]; LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= { const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= {
cCollenZero, cCollenZero,
cCollenOne, cCollenOne, cCollenOne,
cCollenOne, cCollenOne, cCollenOne, cCollenOne, cCollenOne, cCollenOne,
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenOne, cCollenOne, cCollenOne,
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo,
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo,
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
}; };
Context not available.
//! how many dimensions? //! how many dimensions?
const int LbmFsgrSolver::cDimension = 2; const int LbmFsgrSolver::cDimension = 2;
//! Wi factors for collide step //! Wi factors for collide step
const LbmFloat LbmFsgrSolver::cCollenZero = (4.0/9.0); const LbmFloat LbmFsgrSolver::cCollenZero = (4.0/9.0);
const LbmFloat LbmFsgrSolver::cCollenOne = (1.0/9.0); const LbmFloat LbmFsgrSolver::cCollenOne = (1.0/9.0);
const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0); const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
//! threshold value for filled/emptied cells //! threshold value for filled/emptied cells
const LbmFloat LbmFsgrSolver::cMagicNr2 = 1.0005; const LbmFloat LbmFsgrSolver::cMagicNr2 = 1.0005;
const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005; const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
const LbmFloat LbmFsgrSolver::cMagicNr = 1.010001; const LbmFloat LbmFsgrSolver::cMagicNr = 1.010001;
const LbmFloat LbmFsgrSolver::cMagicNrNeg = -0.010001; const LbmFloat LbmFsgrSolver::cMagicNrNeg = -0.010001;
//! size of a single set of distribution functions //! size of a single set of distribution functions
const int LbmFsgrSolver::cDfNum = 9; const int LbmFsgrSolver::cDfNum = 9;
const int LbmFsgrSolver::cDirNum = 9; const int LbmFsgrSolver::cDirNum = 9;
//const string LbmFsgrSolver::dfString[ cDfNum ] = { //const string LbmFsgrSolver::dfString[ cDfNum ] = {
const char* LbmFsgrSolver::dfString[ cDfNum ] = { const char* LbmFsgrSolver::dfString[ cDfNum ] = {
" C", " C",
" N", " S", " E", " W", " N", " S", " E", " W",
"NE", "NW", "SE","SW" "NE", "NW", "SE","SW"
}; };
const int LbmFsgrSolver::dfNorm[ cDfNum ] = { const int LbmFsgrSolver::dfNorm[ cDfNum ] = {
cDirC, cDirC,
cDirN, cDirS, cDirE, cDirW, cDirN, cDirS, cDirE, cDirW,
cDirNE, cDirNW, cDirSE, cDirSW cDirNE, cDirNW, cDirSE, cDirSW
}; };
const int LbmFsgrSolver::dfInv[ cDfNum ] = { const int LbmFsgrSolver::dfInv[ cDfNum ] = {
cDirC, cDirC,
cDirS, cDirN, cDirW, cDirE, cDirS, cDirN, cDirW, cDirE,
cDirSW, cDirSE, cDirNW, cDirNE cDirSW, cDirSE, cDirNW, cDirNE
}; };
const int LbmFsgrSolver::dfRefX[ cDfNum ] = { const int LbmFsgrSolver::dfRefX[ cDfNum ] = {
0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
cDirSE, cDirSW, cDirNE, cDirNW cDirSE, cDirSW, cDirNE, cDirNW
}; };
const int LbmFsgrSolver::dfRefY[ cDfNum ] = { const int LbmFsgrSolver::dfRefY[ cDfNum ] = {
0, 0,
0, 0, 0, 0, 0, 0, 0, 0,
cDirNW, cDirNE, cDirSW, cDirSE cDirNW, cDirNE, cDirSW, cDirSE
}; };
const int LbmFsgrSolver::dfRefZ[ cDfNum ] = { const int LbmFsgrSolver::dfRefZ[ cDfNum ] = {
0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0 0, 0, 0, 0
}; };
// Vector Order 2D: // Vector Order 2D:
// 0 1 2 3 4 5 6 7 8 // 0 1 2 3 4 5 6 7 8
// 0, 0,0, 1,-1, 1,-1,1,-1 // 0, 0,0, 1,-1, 1,-1,1,-1
// 0, 1,-1, 0,0, 1,1,-1,-1 // 0, 1,-1, 0,0, 1,1,-1,-1
const int LbmFsgrSolver::dfVecX[ cDirNum ] = { const int LbmFsgrSolver::dfVecX[ cDirNum ] = {
0, 0,
0,0, 1,-1, 0,0, 1,-1,
1,-1,1,-1 1,-1,1,-1
}; };
const int LbmFsgrSolver::dfVecY[ cDirNum ] = { const int LbmFsgrSolver::dfVecY[ cDirNum ] = {
0, 0,
1,-1, 0,0, 1,-1, 0,0,
1,1,-1,-1 1,1,-1,-1
}; };
const int LbmFsgrSolver::dfVecZ[ cDirNum ] = { const int LbmFsgrSolver::dfVecZ[ cDirNum ] = {
0, 0,0,0,0, 0,0,0,0 0, 0,0,0,0, 0,0,0,0
}; };
const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = { const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
0, 0,
0,0, 1,-1, 0,0, 1,-1,
1,-1,1,-1 1,-1,1,-1
}; };
const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = { const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
0, 0,
1,-1, 0,0, 1,-1, 0,0,
1,1,-1,-1 1,1,-1,-1
}; };
const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = { const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
0, 0,0,0,0, 0,0,0,0 0, 0,0,0,0, 0,0,0,0
}; };
const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = { const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = {
1,-1, 0,0 1,-1, 0,0
}; };
const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = { const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = {
0,0, 1,-1 0,0, 1,-1
}; };
const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = { const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = {
0,0, 0,0 0,0, 0,0
}; };
Context not available.
LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ]; LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= { const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= {
cCollenZero, cCollenZero,
cCollenOne, cCollenOne, cCollenOne, cCollenOne, cCollenOne, cCollenOne, cCollenOne, cCollenOne,
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
}; };
Context not available.
LbmFsgrSolver::LbmFsgrSolver() : LbmFsgrSolver::LbmFsgrSolver() :
//D(), //D(),
mCurrentMass(0.0), mCurrentVolume(0.0), mCurrentMass(0.0), mCurrentVolume(0.0),
mNumProblems(0), mNumProblems(0),
mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0), mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
mpPreviewSurface(NULL), mpPreviewSurface(NULL),
mTimeAdap(true), mForceTimeStepReduce(false), mTimeAdap(true), mForceTimeStepReduce(false),
mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false), mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
mInitSurfaceSmoothing(0), mFsSurfGenSetting(0), mInitSurfaceSmoothing(0), mFsSurfGenSetting(0),
Context not available.
mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(), mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
mMOIVertices(), mMOIVerticesOld(), mMOINormals(), mMOIVertices(), mMOIVerticesOld(), mMOINormals(),
mIsoWeightMethod(1), mIsoWeightMethod(1),
mMaxRefine(1), mMaxRefine(1),
mDfScaleUp(-1.0), mDfScaleDown(-1.0), mDfScaleUp(-1.0), mDfScaleDown(-1.0),
mInitialCsmago(0.02), // set to 0.02 for mMaxRefine==0 below and default for fine level, coarser ones are 0.03 mInitialCsmago(0.02), // set to 0.02 for mMaxRefine==0 below and default for fine level, coarser ones are 0.03
mDebugOmegaRet(0.0), mDebugOmegaRet(0.0),
Context not available.
#endif // LBM_INCLUDE_TESTSOLVERS!=1 #endif // LBM_INCLUDE_TESTSOLVERS!=1
mpIso = new IsoSurface( mIsoValue ); mpIso = new IsoSurface( mIsoValue );
// init equilibrium dist. func // init equilibrium dist. func
LbmFloat rho=1.0; LbmFloat rho=1.0;
FORDF0 { FORDF0 {
dfEquil[l] = this->getCollideEq( l,rho, 0.0, 0.0, 0.0); dfEquil[l] = this->getCollideEq( l,rho, 0.0, 0.0, 0.0);
Context not available.
// init LES // init LES
int odm = 0; int odm = 0;
for(int m=0; m<LBMDIM; m++) { for(int m=0; m<LBMDIM; m++) {
for(int l=0; l<cDfNum; l++) { for(int l=0; l<cDfNum; l++) {
this->lesCoeffDiag[m][l] = this->lesCoeffDiag[m][l] =
this->lesCoeffOffdiag[m][l] = 0.0; this->lesCoeffOffdiag[m][l] = 0.0;
} }
} }
for(int m=0; m<LBMDIM; m++) { for(int m=0; m<LBMDIM; m++) {
for(int n=0; n<LBMDIM; n++) { for(int n=0; n<LBMDIM; n++) {
for(int l=1; l<cDfNum; l++) { for(int l=1; l<cDfNum; l++) {
LbmFloat em; LbmFloat em;
switch(m) { switch(m) {
case 0: em = dfDvecX[l]; break; case 0: em = dfDvecX[l]; break;
Context not available.
mDvecNrm[0] = LbmVec(0.0); mDvecNrm[0] = LbmVec(0.0);
FORDF1 { FORDF1 {
mDvecNrm[l] = getNormalized( mDvecNrm[l] = getNormalized(
LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] ) LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] )
) * -1.0; ) * -1.0;
} }
// calculate gauss weights for restriction // calculate gauss weights for restriction
Context not available.
errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!"); errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!");
#endif #endif
for(int n=0;(n<cDirNum); n++) { mGaussw[n] = 0.0; } for(int n=0;(n<cDirNum); n++) { mGaussw[n] = 0.0; }
//for(int n=0;(n<cDirNum); n++) { //for(int n=0;(n<cDirNum); n++) {
for(int n=0;(n<cDfNum); n++) { for(int n=0;(n<cDfNum); n++) {
const LbmFloat d = norm(LbmVec(dfVecX[n], dfVecY[n], dfVecZ[n])); const LbmFloat d = norm(LbmVec(dfVecX[n], dfVecY[n], dfVecZ[n]));
LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw ); LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw );
mGaussw[n] = w; mGaussw[n] = w;
totGaussw += w; totGaussw += w;
} }
for(int n=0;(n<cDirNum); n++) { for(int n=0;(n<cDirNum); n++) {
mGaussw[n] = mGaussw[n]/totGaussw; mGaussw[n] = mGaussw[n]/totGaussw;
} }
Context not available.
delete mpIso; delete mpIso;
if(mpPreviewSurface) delete mpPreviewSurface; if(mpPreviewSurface) delete mpPreviewSurface;
// cleanup done during scene deletion... // cleanup done during scene deletion...
if(mpControl) delete mpControl; if(mpControl) delete mpControl;
// always output performance estimate // always output performance estimate
Context not available.
/****************************************************************************** /******************************************************************************
* initilize variables fom attribute list * initilize variables fom attribute list
*****************************************************************************/ *****************************************************************************/
void LbmFsgrSolver::parseAttrList() void LbmFsgrSolver::parseAttrList()
{ {
Context not available.
if(mInitialCsmago<=0.0) { if(mInitialCsmago<=0.0) {
if(OPT3D==1) { if(OPT3D==1) {
errFatal("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version...",SIMWORLD_INITERROR); errFatal("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version...",SIMWORLD_INITERROR);
return; return;
} }
} }
LbmFloat fineCsmago = mInitialCsmago; LbmFloat fineCsmago = mInitialCsmago;
LbmFloat coarseCsmago = mInitialCsmago; LbmFloat coarseCsmago = mInitialCsmago;
LbmFloat maxFineCsmago1 = 0.026; LbmFloat maxFineCsmago1 = 0.026;
LbmFloat maxCoarseCsmago1 = 0.029; // try stabilizing LbmFloat maxCoarseCsmago1 = 0.029; // try stabilizing
LbmFloat maxFineCsmago2 = 0.028; LbmFloat maxFineCsmago2 = 0.028;
LbmFloat maxCoarseCsmago2 = 0.032; // try stabilizing some more LbmFloat maxCoarseCsmago2 = 0.032; // try stabilizing some more
if((mMaxRefine==1)&&(mInitialCsmago<maxFineCsmago1)) { if((mMaxRefine==1)&&(mInitialCsmago<maxFineCsmago1)) {
fineCsmago = maxFineCsmago1; fineCsmago = maxFineCsmago1;
Context not available.
fineCsmago = maxFineCsmago2; fineCsmago = maxFineCsmago2;
coarseCsmago = maxCoarseCsmago2; coarseCsmago = maxCoarseCsmago2;
} }
// use Tau instead of Omega for calculations // use Tau instead of Omega for calculations
{ // init base level { // init base level
Context not available.
mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago; mLevel[i].lcsmago_sqr = mLevel[i].lcsmago*mLevel[i].lcsmago;
mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0); mLevel[i].lcnu = (2.0* (1.0/mLevel[i].omega)-1.0) * (1.0/6.0);
} }
// for lbgk // for lbgk
mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega; mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega;
for(int i=mMaxRefine-1; i>=0; i--) { for(int i=mMaxRefine-1; i>=0; i--) {
Context not available.
for(int i=0; i<=mMaxRefine; i++) { for(int i=0; i<=mMaxRefine; i++) {
if(!mSilent) { if(!mSilent) {
errMsg("LbmFsgrSolver", "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" offs:"<<mLevel[i].lOffsx<<","<<mLevel[i].lOffsy<<","<<mLevel[i].lOffsz errMsg("LbmFsgrSolver", "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" offs:"<<mLevel[i].lOffsx<<","<<mLevel[i].lOffsy<<","<<mLevel[i].lOffsz
<<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", " <<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
<<" cmsagp:"<<mLevel[i].lcsmago<<", " <<" cmsagp:"<<mLevel[i].lcsmago<<", "
<< " ss"<<mLevel[i].timestep<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize ); << " ss"<<mLevel[i].timestep<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
Context not available.
} }
if(mFsSurfGenSetting==-1) { if(mFsSurfGenSetting==-1) {
// all on // all on
mFsSurfGenSetting = mFsSurfGenSetting =
fssgNormal | fssgNoNorth | fssgNoSouth | fssgNoEast | fssgNormal | fssgNoNorth | fssgNoSouth | fssgNoEast |
fssgNoWest | fssgNoTop | fssgNoBottom | fssgNoObs ; fssgNoWest | fssgNoTop | fssgNoBottom | fssgNoObs ;
} }
Context not available.
double sizeReduction = 1.0; double sizeReduction = 1.0;
double memEstFromFunc = -1.0; double memEstFromFunc = -1.0;
double memEstFine = -1.0; double memEstFine = -1.0;
string memreqStr(""); string memreqStr("");
bool firstMInit = true; bool firstMInit = true;
int minitTries=0; int minitTries=0;
while(!memOk) { while(!memOk) {
Context not available.
#endif // LBM_INCLUDE_TESTSOLVERS==1 #endif // LBM_INCLUDE_TESTSOLVERS==1
firstMInit=false; firstMInit=false;
calculateMemreqEstimate( mSizex, mSizey, mSizez, calculateMemreqEstimate( mSizex, mSizey, mSizez,
mMaxRefine, mFarFieldSize, &memEstFromFunc, &memEstFine, &memreqStr ); mMaxRefine, mFarFieldSize, &memEstFromFunc, &memEstFine, &memreqStr );
bool noLimit = false; bool noLimit = false;
double memLimit = 0.; double memLimit = 0.;
string memLimStr("-"); string memLimStr("-");
Context not available.
noLimit = true; noLimit = true;
} }
// restrict max. chunk of 1 mem block to 1GB for windos // restrict max. chunk of 1 mem block to 1GB for windows
bool memBlockAllocProblem = false; bool memBlockAllocProblem = false;
double maxDefaultMemChunk = 2.*1024.*1024.*1024.; double maxDefaultMemChunk = 2.*1024.*1024.*1024.;
//std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG //std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG
Context not available.
, 3 ); , 3 );
} else { } else {
memOk = true; memOk = true;
} }
} }
mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex; mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex;
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<< debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<<
", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<< ", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<<
Context not available.
<<"FSGR_LISTTRICK="<<FSGR_LISTTRICK <<" " <<"FSGR_LISTTRICK="<<FSGR_LISTTRICK <<" "
<<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY <<" " <<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY <<" "
<<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL <<" " <<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL <<" "
<<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" " <<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" "
<<"USE_LES="<<USE_LES <<" " <<"USE_LES="<<USE_LES <<" "
,10); ,10);
// perform 2D corrections... // perform 2D corrections...
Context not available.
// init vectors // init vectors
for(int i=0; i<=mMaxRefine; i++) { for(int i=0; i<=mMaxRefine; i++) {
mLevel[i].id = i; mLevel[i].id = i;
mLevel[i].nodeSize = 0.0; mLevel[i].nodeSize = 0.0;
mLevel[i].simCellSize = 0.0; mLevel[i].simCellSize = 0.0;
mLevel[i].omega = 0.0; mLevel[i].omega = 0.0;
mLevel[i].time = 0.0; mLevel[i].time = 0.0;
mLevel[i].timestep = 1.0; mLevel[i].timestep = 1.0;
mLevel[i].gravity = LbmVec(0.0); mLevel[i].gravity = LbmVec(0.0);
mLevel[i].mprsCells[0] = NULL; mLevel[i].mprsCells[0] = NULL;
mLevel[i].mprsCells[1] = NULL; mLevel[i].mprsCells[1] = NULL;
mLevel[i].mprsFlags[0] = NULL; mLevel[i].mprsFlags[0] = NULL;
mLevel[i].mprsFlags[1] = NULL; mLevel[i].mprsFlags[1] = NULL;
mLevel[i].avgOmega = 0.0; mLevel[i].avgOmega = 0.0;
mLevel[i].avgOmegaCnt = 0.0; mLevel[i].avgOmegaCnt = 0.0;
} }
Context not available.
errMsg("LbmFsgrSolver::initialize","Sanity Error - memory estimate is off! real:"<<ownMemCheck<<" vs. estimate:"<<memEstFromFunc ); errMsg("LbmFsgrSolver::initialize","Sanity Error - memory estimate is off! real:"<<ownMemCheck<<" vs. estimate:"<<memEstFromFunc );
} }
#endif // ELBEEM_PLUGIN!=1 #endif // ELBEEM_PLUGIN!=1
// init sizes for _all_ levels // init sizes for _all_ levels
for(int i=mMaxRefine; i>=0; i--) { for(int i=mMaxRefine; i>=0; i--) {
mLevel[i].lOffsx = mLevel[i].lSizex; mLevel[i].lOffsx = mLevel[i].lSizex;
Context not available.
// init iso weight values mIsoWeightMethod // init iso weight values mIsoWeightMethod
int wcnt = 0; int wcnt = 0;
float totw = 0.0; float totw = 0.0;
for(int ak=-1;ak<=1;ak++) for(int ak=-1;ak<=1;ak++)
for(int aj=-1;aj<=1;aj++) for(int aj=-1;aj<=1;aj++)
for(int ai=-1;ai<=1;ai++) { for(int ai=-1;ai<=1;ai++) {
switch(mIsoWeightMethod) { switch(mIsoWeightMethod) {
case 1: // light smoothing case 1: // light smoothing
Context not available.
mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) ); mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) );
// reset iso field // reset iso field
for(int ak=0;ak<isosz;ak++) for(int ak=0;ak<isosz;ak++)
for(int aj=0;aj<isosy;aj++) for(int aj=0;aj<isosy;aj++)
for(int ai=0;ai<isosx;ai++) { *mpIso->getData(ai,aj,ak) = 0.0; } for(int ai=0;ai<isosx;ai++) { *mpIso->getData(ai,aj,ak) = 0.0; }
Context not available.
FSGR_FORIJK_BOUNDS(lev) { FSGR_FORIJK_BOUNDS(lev) {
RFLAG(lev,i,j,k,0) = 0, RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage RFLAG(lev,i,j,k,0) = 0, RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
if(!mAllfluid) { if(!mAllfluid) {
initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0); initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0);
} else { } else {
initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0); initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0);
} }
} }
} }
Context not available.
mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) ); mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) );
//mpPreviewSurface->setName( getName() + "preview" ); //mpPreviewSurface->setName( getName() + "preview" );
mpPreviewSurface->setName( "preview" ); mpPreviewSurface->setName( "preview" );
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10); debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10);
} }
Context not available.
bool LbmFsgrSolver::initializeSolverGrids() { bool LbmFsgrSolver::initializeSolverGrids() {
/* init boundaries */ /* init boundaries */
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10); debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10);
// init obstacles, and reinit time step size // init obstacles, and reinit time step size
initGeometryFlags(); initGeometryFlags();
mLastSimTime = -1.0; mLastSimTime = -1.0;
// TODO check for invalid cells? nitGenericTestCases(); // TODO check for invalid cells? nitGenericTestCases();
// new - init noslip 1 everywhere... // new - init noslip 1 everywhere...
// half fill boundary cells? // half fill boundary cells?
Context not available.
} else if(mDomainBound.find(string("part")) != string::npos) { } else if(mDomainBound.find(string("part")) != string::npos) {
domainBoundType = CFBnd | CFBndPartslip; // part slip type domainBoundType = CFBnd | CFBndPartslip; // part slip type
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10); debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10);
} else { } else {
domainBoundType = CFBnd | CFBndNoslip; domainBoundType = CFBnd | CFBndNoslip;
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10); debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10);
} }
Context not available.
//for(int i=0; i<(int)(domainobj+0); i++) { //for(int i=0; i<(int)(domainobj+0); i++) {
//errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName()); //errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName());
//if((*mpGiObjects)[i] == mpIso) { //check... //if((*mpGiObjects)[i] == mpIso) { //check...
//} //}
//} //}
//errMsg("GEOIN"," dm "<<(domainBoundType>>24)); //errMsg("GEOIN"," dm "<<(domainBoundType>>24));
for(int k=0;k<mLevel[mMaxRefine].lSizez;k++) for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) { for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
initEmptyCell(mMaxRefine, i,0,k, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, i,0,k, domainBoundType, 0.0, BND_FILL);
initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, i,mLevel[mMaxRefine].lSizey-1,k, domainBoundType, 0.0, BND_FILL);
} }
for(int k=0;k<mLevel[mMaxRefine].lSizez;k++) for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) { for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
initEmptyCell(mMaxRefine, 0,j,k, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, 0,j,k, domainBoundType, 0.0, BND_FILL);
initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, domainBoundType, 0.0, BND_FILL);
// DEBUG BORDER! // DEBUG BORDER!
//initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); //initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL);
} }
if(LBMDIM == 3) { if(LBMDIM == 3) {
// only for 3D // only for 3D
for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) { for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
initEmptyCell(mMaxRefine, i,j,0, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, i,j,0, domainBoundType, 0.0, BND_FILL);
initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, domainBoundType, 0.0, BND_FILL);
} }
} }
// TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11 // TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
/*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++) /*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) { for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL);
} }
for(int k=0;k<mLevel[mMaxRefine].lSizez;k++) for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) { for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
initEmptyCell(mMaxRefine, i,1,k, domainBoundType, 0.0, BND_FILL); initEmptyCell(mMaxRefine, i,1,k, domainBoundType, 0.0, BND_FILL);
} }
// */ // */
Context not available.
initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, domainBoundType, 0.0, BND_FILL); // SYMM!? 2D? initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2-ii,j,k, domainBoundType, 0.0, BND_FILL); // SYMM!? 2D?
} }
for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) { for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, domainBoundType, 0.0, BND_FILL); // SYMM!? 3D? initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, domainBoundType, 0.0, BND_FILL); // SYMM!? 3D?
} }
} }
Context not available.
int cy1=cyo-(int)(vdist*sy); int cy1=cyo-(int)(vdist*sy);
int cy2=cyo+(int)(vdist*sy); int cy2=cyo+(int)(vdist*sy);
//for(int j=cy-sy;j<cy+sy;j++) for(int i=cx-sx;i<cx+sx;i++) { //for(int j=cy-sy;j<cy+sy;j++) for(int i=cx-sx;i<cx+sx;i++) {
for(int j=1;j<mLevel[level].lSizey-1;j++) for(int j=1;j<mLevel[level].lSizey-1;j++)
for(int i=1;i<mLevel[level].lSizex-1;i++) { for(int i=1;i<mLevel[level].lSizex-1;i++) {
LbmFloat d1 = norm(LbmVec(cx,cy1,0.)-LbmVec(i,j,0)); LbmFloat d1 = norm(LbmVec(cx,cy1,0.)-LbmVec(i,j,0));
Context not available.
int inmCellCnt = 0; int inmCellCnt = 0;
FSGR_FORIJK1(mMaxRefine) { FSGR_FORIJK1(mMaxRefine) {
if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) { if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) {
LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0); LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0);
FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); } FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
mInitialMass += fluidRho; mInitialMass += fluidRho;
inmCellCnt ++; inmCellCnt ++;
Context not available.
/*! prepare actual simulation start, setup viz etc */ /*! prepare actual simulation start, setup viz etc */
bool LbmFsgrSolver::initializeSolverPostinit() { bool LbmFsgrSolver::initializeSolverPostinit() {
// coarsen region // coarsen region
myTime_t fsgrtstart = getTime(); myTime_t fsgrtstart = getTime();
for(int lev=mMaxRefine-1; lev>=0; lev--) { for(int lev=mMaxRefine-1; lev>=0; lev--) {
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8); debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
adaptGrid(lev); adaptGrid(lev);
Context not available.
coarseRestrictFromFine(lev); coarseRestrictFromFine(lev);
} }
markedClearList(); markedClearList();
myTime_t fsgrtend = getTime(); myTime_t fsgrtend = getTime();
if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); } if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
mNumFsgrChanges = 0; mNumFsgrChanges = 0;
for(int l=0; l<cDirNum; l++) { for(int l=0; l<cDirNum; l++) {
LbmFloat area = 0.5 * 0.5 *0.5; LbmFloat area = 0.5 * 0.5 *0.5;
if(LBMDIM==2) area = 0.5 * 0.5; if(LBMDIM==2) area = 0.5 * 0.5;
Context not available.
if(k<=0) continue; \ if(k<=0) continue; \
if(k>=mLevel[level].lSizez-1) continue; \ if(k>=mLevel[level].lSizez-1) continue; \
#endif // LBMDIM #endif // LBMDIM
// calculate object velocity from vert arrays in objvel vec // calculate object velocity from vert arrays in objvel vec
#define OBJVEL_CALC \ #define OBJVEL_CALC \
Context not available.
/*****************************************************************************/ /*****************************************************************************/
//! init moving obstacles for next sim step sim //! init moving obstacles for next sim step sim
/*****************************************************************************/ /*****************************************************************************/
void LbmFsgrSolver::initMovingObstacles(bool staticInit) { void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
myTime_t monstart = getTime(); myTime_t monstart = getTime();
Context not available.
// 2d display as rectangles // 2d display as rectangles
ntlVec3Gfx iniPos(0.0); ntlVec3Gfx iniPos(0.0);
if(LBMDIM==2) { if(LBMDIM==2) {
dvec[2] = 1.0; dvec[2] = 1.0;
iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0); iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0);
} else { } else {
iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0); iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
} }
if( (int)mObjectMassMovnd.size() < numobjs) { if( (int)mObjectMassMovnd.size() < numobjs) {
for(int i=mObjectMassMovnd.size(); i<numobjs; i++) { for(int i=mObjectMassMovnd.size(); i<numobjs; i++) {
mObjectMassMovnd.push_back(0.); mObjectMassMovnd.push_back(0.);
} }
} }
// stats // stats
int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0; int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
int nbored; int nbored;
Context not available.
if(skip) continue; if(skip) continue;
debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<mLbmInitId, 10); debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" skip:"<<skip<<", static:"<<staticInit<<" anim:"<<obj->getIsAnimated()<<" gid:"<<obj->getGeoInitId()<<" simgid:"<<mLbmInitId, 10);
if( (obj->getGeoInitType()&FGI_ALLBOUNDS) || if( (obj->getGeoInitType()&FGI_ALLBOUNDS) ||
((obj->getGeoInitType()&FGI_FLUID) && staticInit) ) { ((obj->getGeoInitType()&FGI_FLUID) && staticInit) ) {
otype = ntype = CFInvalid; otype = ntype = CFInvalid;
switch(obj->getGeoInitType()) { switch(obj->getGeoInitType()) {
/* case FGI_BNDPART: // old, use noslip for moving part/free objs /* case FGI_BNDPART: // old, use noslip for moving part/free objs
case FGI_BNDFREE: case FGI_BNDFREE:
if(!staticInit) { if(!staticInit) {
errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() ); errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
otype = ntype = CFBnd|CFBndNoslip; otype = ntype = CFBnd|CFBndNoslip;
Context not available.
if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip; if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip;
if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip; if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
} }
break; break;
// off */ // off */
case FGI_BNDPART: rhomass = BND_FILL; case FGI_BNDPART: rhomass = BND_FILL;
otype = ntype = CFBnd|CFBndPartslip|(OId<<24); otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
Context not available.
case FGI_BNDNO: rhomass = BND_FILL; case FGI_BNDNO: rhomass = BND_FILL;
otype = ntype = CFBnd|CFBndNoslip|(OId<<24); otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
break; break;
case FGI_FLUID: case FGI_FLUID:
otype = ntype = CFFluid; otype = ntype = CFFluid;
break; break;
case FGI_MBNDINFLOW: case FGI_MBNDINFLOW:
otype = ntype = CFMbndInflow; otype = ntype = CFMbndInflow;
break; break;
case FGI_MBNDOUTFLOW: case FGI_MBNDOUTFLOW:
otype = ntype = CFMbndOutflow; otype = ntype = CFMbndOutflow;
break; break;
} }
int wasActive = ((obj->getGeoActive(sourceTime)>0.)? 1:0); int wasActive = ((obj->getGeoActive(sourceTime)>0.)? 1:0);
Context not available.
if(ntype&(CFBndFreeslip|CFBndPartslip)) { pNormals = &mMOINormals; } if(ntype&(CFBndFreeslip|CFBndPartslip)) { pNormals = &mMOINormals; }
mMOIVertices.clear(); mMOIVertices.clear();
if(obj->getMeshAnimated()) { if(obj->getMeshAnimated()) {
// do two full update // do two full update
// TODO tNormals handling!? // TODO tNormals handling!?
mMOIVerticesOld.clear(); mMOIVerticesOld.clear();
Context not available.
if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; } if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
CellFlagType rflagnb[27]; CellFlagType rflagnb[27];
LbmFloat massCheck = 0.; LbmFloat massCheck = 0.;
int massReinits=0; int massReinits=0;
bool fillCells = (mObjectMassMovnd[OId]<=-1.); bool fillCells = (mObjectMassMovnd[OId]<=-1.);
LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime); LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime);
// first pass - set new obs. cells // first pass - set new obs. cells
if(active) { if(active) {
Context not available.
//{ errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<", t="<<targetTime); } //{ errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<", t="<<targetTime); }
if(QCELL(level, i,j,k, workSet, dFlux)==targetTime) continue; if(QCELL(level, i,j,k, workSet, dFlux)==targetTime) continue;
monPoints++; monPoints++;
// check mass // check mass
if(RFLAG(level, i,j,k, workSet)&(CFFluid)) { if(RFLAG(level, i,j,k, workSet)&(CFFluid)) {
FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); } FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); }
Context not available.
rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l); rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
if(rflagnb[l]&(CFFluid|CFInter)) { if(rflagnb[l]&(CFFluid|CFInter)) {
rflagnb[l] &= (~CFNoBndFluid); // remove CFNoBndFluid flag rflagnb[l] &= (~CFNoBndFluid); // remove CFNoBndFluid flag
RFLAG_NB(level, i,j,k,workSet,l) &= rflagnb[l]; RFLAG_NB(level, i,j,k,workSet,l) &= rflagnb[l];
} }
} }
LbmFloat *dstCell = RACPNT(level, i,j,k,workSet); LbmFloat *dstCell = RACPNT(level, i,j,k,workSet);
RAC(dstCell,0) = 0.0; RAC(dstCell,0) = 0.0;
if(ntype&CFBndMoving) { if(ntype&CFBndMoving) {
OBJVEL_CALC; OBJVEL_CALC;
// compute fluid acceleration // compute fluid acceleration
FORDF1 { FORDF1 {
if(rflagnb[l]&(CFFluid|CFInter)) { if(rflagnb[l]&(CFFluid|CFInter)) {
const LbmFloat ux = dfDvecX[l]*objvel[0]; const LbmFloat ux = dfDvecX[l]*objvel[0];
const LbmFloat uy = dfDvecY[l]*objvel[1]; const LbmFloat uy = dfDvecY[l]*objvel[1];
const LbmFloat uz = dfDvecZ[l]*objvel[2]; const LbmFloat uz = dfDvecZ[l]*objvel[2];
LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); // LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); //
if(ntype&(CFBndFreeslip|CFBndPartslip)) { if(ntype&(CFBndFreeslip|CFBndPartslip)) {
// missing, diag mass checks... // missing, diag mass checks...
//if(l<=LBMDIM*2) factor *= 1.4142; //if(l<=LBMDIM*2) factor *= 1.4142;
Context not available.
//rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l); //rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
rflagnb[l] = RFLAG_NB(level, i,j,k,otherSet,l); // test use other set to not have loop dependance rflagnb[l] = RFLAG_NB(level, i,j,k,otherSet,l); // test use other set to not have loop dependance
nbored |= rflagnb[l]; nbored |= rflagnb[l];
} }
CellFlagType settype = CFInvalid; CellFlagType settype = CFInvalid;
if(nbored&CFFluid) { if(nbored&CFFluid) {
settype = CFInter|CFNoInterpolSrc; settype = CFInter|CFNoInterpolSrc;
rhomass = 1.5; rhomass = 1.5;
if(!fillCells) rhomass = 0.; if(!fillCells) rhomass = 0.;
Context not available.
// only compute mass transfer when init is done // only compute mass transfer when init is done
if(mInitDone) { if(mInitDone) {
errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<< errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
" fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass ); " fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass );
mObjectMassMovnd[OId] += massCheck; mObjectMassMovnd[OId] += massCheck;
} }
} // bnd, active } // bnd, active
Context not available.
if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) { if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
//changeFlag(level, i,j,k, workSet, setflflag); //changeFlag(level, i,j,k, workSet, setflflag);
initVelocityCell(level, i,j,k, setflflag, 1., 1., mObjectSpeeds[OId] ); initVelocityCell(level, i,j,k, setflflag, 1., 1., mObjectSpeeds[OId] );
} }
//else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag); } //else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { changeFlag(level, i,j,k, workSet, RFLAG(level, i,j,k, workSet)|set2Flag); }
} }
} // second static inflow pass } // second static inflow pass
Context not available.
else if(ntype&CFMbndInflow){ else if(ntype&CFMbndInflow){
// inflow pass - add new fluid cells // inflow pass - add new fluid cells
// this is slightly different for standing inflows, // this is slightly different for standing inflows,
// as the fluid is forced to the desired velocity inside as // as the fluid is forced to the desired velocity inside as
// well... // well...
const LbmFloat iniRho = 1.0; const LbmFloat iniRho = 1.0;
const LbmVec vel(mObjectSpeeds[OId]); const LbmVec vel(mObjectSpeeds[OId]);
const LbmFloat usqr = (vel[0]*vel[0]+vel[1]*vel[1]+vel[2]*vel[2])*1.5; const LbmFloat usqr = (vel[0]*vel[0]+vel[1]*vel[1]+vel[2]*vel[2])*1.5;
USQRMAXCHECK(usqr,vel[0],vel[1],vel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); USQRMAXCHECK(usqr,vel[0],vel[1],vel[2], mMaxVlen, mMxvx,mMxvy,mMxvz);
//errMsg("LbmFsgrSolver::initMovingObstacles","id"<<OId<<" "<<obj->getName()<<" inflow "<<staticInit<<" "<<mMOIVertices.size() ); //errMsg("LbmFsgrSolver::initMovingObstacles","id"<<OId<<" "<<obj->getName()<<" inflow "<<staticInit<<" "<<mMOIVertices.size() );
for(size_t n=0; n<mMOIVertices.size(); n++) { for(size_t n=0; n<mMOIVertices.size(); n++) {
POS2GRID_CHECK(mMOIVertices,n); POS2GRID_CHECK(mMOIVertices,n);
// TODO - also reinit interface cells !? // TODO - also reinit interface cells !?
if(RFLAG(level, i,j,k, workSet)!=CFEmpty) { if(RFLAG(level, i,j,k, workSet)!=CFEmpty) {
// test prevent particle gen for inflow inter cells // test prevent particle gen for inflow inter cells
if(RFLAG(level, i,j,k, workSet)&(CFInter)) { RFLAG(level, i,j,k, workSet) &= (~CFNoBndFluid); } // remove CFNoBndFluid flag if(RFLAG(level, i,j,k, workSet)&(CFInter)) { RFLAG(level, i,j,k, workSet) &= (~CFNoBndFluid); } // remove CFNoBndFluid flag
continue; } continue; }
Context not available.
if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) { if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
forceChangeFlag(level, i,j,k, workSet, set2Flag); forceChangeFlag(level, i,j,k, workSet, set2Flag);
} else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
forceChangeFlag(level, i,j,k, workSet, forceChangeFlag(level, i,j,k, workSet,
(RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag); (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
} }
} }
Context not available.
if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) { if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
forceChangeFlag(level, i,j,k, workSet, set2Flag); forceChangeFlag(level, i,j,k, workSet, set2Flag);
} else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) { } else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
forceChangeFlag(level, i,j,k, workSet, forceChangeFlag(level, i,j,k, workSet,
(RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag); (RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
} }
} }
Context not available.
} }
} }
} // DEBUG */ } // DEBUG */
#undef POS2GRID_CHECK #undef POS2GRID_CHECK
myTime_t monend = getTime(); myTime_t monend = getTime();
if(monend-monstart>0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7); if(monend-monstart>0) debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"Total: "<<monTotal<<" Points :"<<monPoints<<" ObstInits:"<<monObsts<<" FlInits:"<<monFluids<<" Trafos:"<<monTotal<<", took "<<getTimeString(monend-monstart), 7);
Context not available.
extern int globGeoInitDebug; //solver_interface extern int globGeoInitDebug; //solver_interface
bool LbmFsgrSolver::initGeometryFlags() { bool LbmFsgrSolver::initGeometryFlags() {
int level = mMaxRefine; int level = mMaxRefine;
myTime_t geotimestart = getTime(); myTime_t geotimestart = getTime();
ntlGeometryObject *pObj; ntlGeometryObject *pObj;
ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0; ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3); debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3);
Context not available.
/* init object velocities, this has always to be called for init */ /* init object velocities, this has always to be called for init */
initGeoTree(); initGeoTree();
if(mAllfluid) { if(mAllfluid) {
freeGeoTree(); freeGeoTree();
return true; } return true; }
Context not available.
if( if(
((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) || ((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) ||
(obj->getVolumeInit()&VOLUMEINIT_SHELL) ) { (obj->getVolumeInit()&VOLUMEINIT_SHELL) ) {
if(!obj->getMeshAnimated()) { if(!obj->getMeshAnimated()) {
debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9); debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
obj->initMovingPoints(mSimulationTime, mLevel[mMaxRefine].nodeSize); obj->initMovingPoints(mSimulationTime, mLevel[mMaxRefine].nodeSize);
} }
Context not available.
// 2d display as rectangles // 2d display as rectangles
if(LBMDIM==2) { if(LBMDIM==2) {
dvec[2] = 0.0; dvec[2] = 0.0;
iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5); iniPos =(mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))+(dvec*0.5);
//if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } } //if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
} else { } else {
Context not available.
for(int j=1;j<mLevel[level].lSizey-1;j++) { for(int j=1;j<mLevel[level].lSizey-1;j++) {
for(int i=1;i<mLevel[level].lSizex-1;i++) { for(int i=1;i<mLevel[level].lSizex-1;i++) {
ntype = CFInvalid; ntype = CFInvalid;
GETPOS(i,j,k); GETPOS(i,j,k);
const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance); const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
if(inside) { if(inside) {
pObj = (*mpGiObjects)[OId]; pObj = (*mpGiObjects)[OId];
switch( pObj->getGeoInitType() ){ switch( pObj->getGeoInitType() ){
case FGI_MBNDINFLOW: case FGI_MBNDINFLOW:
if(! pObj->getIsAnimated() ) { if(! pObj->getIsAnimated() ) {
rhomass = 1.0; rhomass = 1.0;
ntype = CFFluid | CFMbndInflow; ntype = CFFluid | CFMbndInflow;
Context not available.
ntype = CFInvalid; ntype = CFInvalid;
} }
break; break;
case FGI_MBNDOUTFLOW: case FGI_MBNDOUTFLOW:
if(! pObj->getIsAnimated() ) { if(! pObj->getIsAnimated() ) {
rhomass = 0.0; rhomass = 0.0;
ntype = CFEmpty|CFMbndOutflow; ntype = CFEmpty|CFMbndOutflow;
} else { } else {
ntype = CFInvalid; ntype = CFInvalid;
} }
break; break;
case FGI_BNDNO: case FGI_BNDNO:
rhomass = BND_FILL; rhomass = BND_FILL;
ntype = CFBnd|CFBndNoslip; ntype = CFBnd|CFBndNoslip;
break; break;
case FGI_BNDPART: case FGI_BNDPART:
rhomass = BND_FILL; rhomass = BND_FILL;
ntype = CFBnd|CFBndPartslip; break; ntype = CFBnd|CFBndPartslip; break;
case FGI_BNDFREE: case FGI_BNDFREE:
rhomass = BND_FILL; rhomass = BND_FILL;
ntype = CFBnd|CFBndFreeslip; break; ntype = CFBnd|CFBndFreeslip; break;
default: // warn here? default: // warn here?
Context not available.
initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] ); initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
} }
} }
} }
// */ // */
} }
} }
} // zmax } // zmax
// */ // */
Context not available.
} }
} }
} // distance>0 } // distance>0
} }
} }
} // zmax } // zmax
// reset invalid to empty again // reset invalid to empty again
Context not available.
for(int j=1;j<mLevel[level].lSizey-1;j++) { for(int j=1;j<mLevel[level].lSizey-1;j++) {
for(int i=1;i<mLevel[level].lSizex-1;i++) { for(int i=1;i<mLevel[level].lSizex-1;i++) {
if(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFInvalid) { if(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFInvalid) {
RFLAG(level, i,j,k, mLevel[level].setOther) = RFLAG(level, i,j,k, mLevel[level].setOther) =
RFLAG(level, i,j,k, mLevel[level].setCurr) = CFEmpty; RFLAG(level, i,j,k, mLevel[level].setCurr) = CFEmpty;
} }
} }
Context not available.
} }
freeGeoTree(); freeGeoTree();
myTime_t geotimeend = getTime(); myTime_t geotimeend = getTime();
debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<<savedNodes<<") " , 10 ); debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Geometry init done ("<< getTimeString(geotimeend-geotimestart)<<","<<savedNodes<<") " , 10 );
//errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez)); //errMsg(" SAVED "," "<<savedNodes<<" of "<<(mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez));
return true; return true;
} }
Context not available.
double interfaceFill = 0.45; // filling level of interface cells double interfaceFill = 0.45; // filling level of interface cells
//interfaceFill = 1.0; // DEUG!! GEOMTEST!!!!!!!!!!!! //interfaceFill = 1.0; // DEUG!! GEOMTEST!!!!!!!!!!!!
// set interface cells // set interface cells
FSGR_FORIJK1(mMaxRefine) { FSGR_FORIJK1(mMaxRefine) {
if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) { if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) {
FORDF1 { FORDF1 {
Context not available.
} }
} }
// remove invalid interface cells // remove invalid interface cells
FSGR_FORIJK1(mMaxRefine) { FSGR_FORIJK1(mMaxRefine) {
// remove invalid interface cells // remove invalid interface cells
if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) { if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) {
int delit = 0; int delit = 0;
int NBs = 0; // neighbor flags or'ed int NBs = 0; // neighbor flags or'ed
int noEmptyNB = 1; int noEmptyNB = 1;
FORDF1 { FORDF1 {
Context not available.
// remove cells with no empty neighbors // remove cells with no empty neighbors
if(noEmptyNB) { delit = 2; } if(noEmptyNB) { delit = 2; }
// now we can remove the cell // now we can remove the cell
if(delit==1) { if(delit==1) {
initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0); initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
} }
Context not available.
interpolateCellValues(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz); interpolateCellValues(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
initVelocityCell(mMaxRefine, i,j,k, CFFluid, arho,1., LbmVec(aux,auy,auz) ); initVelocityCell(mMaxRefine, i,j,k, CFFluid, arho,1., LbmVec(aux,auy,auz) );
} }
} // interface } // interface
} // */ } // */
// another brute force init, make sure the fill values are right... // another brute force init, make sure the fill values are right...
// and make sure both sets are equal // and make sure both sets are equal
for(int lev=0; lev<=mMaxRefine; lev++) { for(int lev=0; lev<=mMaxRefine; lev++) {
FSGR_FORIJK_BOUNDS(lev) { FSGR_FORIJK_BOUNDS(lev) {
if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) { if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) {
QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL; QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL;
continue; continue;
} }
if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) { if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) {
QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0; QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0;
continue; continue;
} }
Context not available.
LbmFloat mass = 0.0; LbmFloat mass = 0.0;
//LbmFloat nbdiv; //LbmFloat nbdiv;
//FORDF0 { //FORDF0 {
for(int l=0;(l<cDirNum); l++) { for(int l=0;(l<cDirNum); l++) {
int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l]; int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){ if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){
mass += 1.0; mass += 1.0;
Context not available.
if((iindex)>1) { haveStandingFluid=(iindex); } \ if((iindex)>1) { haveStandingFluid=(iindex); } \
j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=getForZMaxBnd(); \ j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=getForZMaxBnd(); \
continue; \ continue; \
} }
int gravIndex[3] = {0,0,0}; int gravIndex[3] = {0,0,0};
int gravDir[3] = {1,1,1}; int gravDir[3] = {1,1,1};
int maxGravComp = 1; // by default y int maxGravComp = 1; // by default y
Context not available.
#define PRINTGDIRS \ #define PRINTGDIRS \
errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \ errMsg("Standing fp","X start="<<gravIMin[0]<<" end="<<gravIMax[0]<<" dir="<<gravDir[0] ); \
errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \ errMsg("Standing fp","Y start="<<gravIMin[1]<<" end="<<gravIMax[1]<<" dir="<<gravDir[1] ); \
errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] ); errMsg("Standing fp","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] );
// _PRINTGDIRS; // _PRINTGDIRS;
bool gravAbort = false; bool gravAbort = false;
Context not available.
gravAbort=false; \ gravAbort=false; \
for(gravIndex[2]= gravIMin[2]; (gravIndex[2]!=gravIMax[2])&&(!gravAbort); gravIndex[2] += gravDir[2]) \ for(gravIndex[2]= gravIMin[2]; (gravIndex[2]!=gravIMax[2])&&(!gravAbort); gravIndex[2] += gravDir[2]) \
for(gravIndex[1]= gravIMin[1]; (gravIndex[1]!=gravIMax[1])&&(!gravAbort); gravIndex[1] += gravDir[1]) \ for(gravIndex[1]= gravIMin[1]; (gravIndex[1]!=gravIMax[1])&&(!gravAbort); gravIndex[1] += gravDir[1]) \
for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort); gravIndex[0] += gravDir[0]) for(gravIndex[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort); gravIndex[0] += gravDir[0])
GRAVLOOP { GRAVLOOP {
int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2]; int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) ||
( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFBndMoving)) ) || ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFBndMoving)) ) ||
( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ) { ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ) {
int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp])); int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] ); if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
if(fluidHeight>1) { if(fluidHeight>1) {
haveStandingFluid = fluidHeight; //gravIndex[maxGravComp]; haveStandingFluid = fluidHeight; //gravIndex[maxGravComp];
gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp]; gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
} }
gravAbort = true; continue; gravAbort = true; continue;
} }
} // GRAVLOOP } // GRAVLOOP
// _PRINTGDIRS; // _PRINTGDIRS;
Context not available.
#endif // ELBEEM_PLUGIN!=1 #endif // ELBEEM_PLUGIN!=1
if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "fheight="<<fluidHeight<<" min="<<PRINT_VEC(gravIMin[0],gravIMin[1], gravIMin[2])<<" max="<<PRINT_VEC(gravIMax[0], gravIMax[1],gravIMax[2])<< if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "fheight="<<fluidHeight<<" min="<<PRINT_VEC(gravIMin[0],gravIMin[1], gravIMin[2])<<" max="<<PRINT_VEC(gravIMax[0], gravIMax[1],gravIMax[2])<<
" mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10); " mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
if(mDisableStandingFluidInit) { if(mDisableStandingFluidInit) {
debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2); debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
haveStandingFluid=0; haveStandingFluid=0;
Context not available.
// copy flags and init , as no flags will be changed during grav init // copy flags and init , as no flags will be changed during grav init
// also important for corasening later on // also important for corasening later on
const int lev = mMaxRefine; const int lev = mMaxRefine;
CellFlagType nbflag[LBM_DFNUM], nbored; CellFlagType nbflag[LBM_DFNUM], nbored;
for(int k=getForZMinBnd();k<getForZMaxBnd(mMaxRefine);++k) { for(int k=getForZMinBnd();k<getForZMaxBnd(mMaxRefine);++k) {
for(int j=0;j<mLevel[lev].lSizey-0;++j) { for(int j=0;j<mLevel[lev].lSizey-0;++j) {
for(int i=0;i<mLevel[lev].lSizex-0;++i) { for(int i=0;i<mLevel[lev].lSizex-0;++i) {
Context not available.
FORDF1 { FORDF1 {
nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l); nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
nbored |= nbflag[l]; nbored |= nbflag[l];
} }
if(nbored&CFBnd) { if(nbored&CFBnd) {
RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid); RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
} else { } else {
Context not available.
int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2]; int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
//debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 ); //debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) || if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){ ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){
//gravAbort = true; //gravAbort = true;
continue; continue;
} }
LbmFloat rho = 1.0; LbmFloat rho = 1.0;
// 1/6 velocity from denisty gradient, 1/2 for delta of two cells // 1/6 velocity from denisty gradient, 1/2 for delta of two cells
rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) * rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) *
(mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega); (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega);
if(debugStandingPreinit) if(debugStandingPreinit)
if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) { if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) {
errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK); errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK);
} }
if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) || if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
Context not available.
debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited (max-rho:"<< debMsgStd("Standing fluid preinit", DM_MSG, "Density gradient inited (max-rho:"<<
(1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega)) <<", h:"<< fluidHeight<<") ", 8); (1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega)) <<", h:"<< fluidHeight<<") ", 8);
int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) ); int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
preinitSteps = (haveStandingFluid>>2); // not much use...? preinitSteps = (haveStandingFluid>>2); // not much use...?
//preinitSteps = 0; //preinitSteps = 0;
Context not available.
for(int s=0; s<2; s++) { for(int s=0; s<2; s++) {
FSGR_FORIJK1(lev) { FSGR_FORIJK1(lev) {
if(i<(mLevel[lev].lSizex/2)) { if(i<(mLevel[lev].lSizex/2)) {
int inb = (mLevel[lev].lSizey-1-i); int inb = (mLevel[lev].lSizey-1-i);
if(lev==mMaxRefine) inb -= 1; // FSGR_SYMM_T if(lev==mMaxRefine) inb -= 1; // FSGR_SYMM_T
Context not available.
LbmFloat maxdiv =0.0; LbmFloat maxdiv =0.0;
if(erro) { if(erro) {
errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv ); errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
//if(LBMDIM==2) mPanic = true; //if(LBMDIM==2) mPanic = true;
//return false; //return false;
} else { } else {
errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv ); errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
Context not available.
}// */ }// */
void void
LbmFsgrSolver::interpolateCellValues( LbmFsgrSolver::interpolateCellValues(
int level,int ei,int ej,int ek,int workSet, int level,int ei,int ej,int ek,int workSet,
LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz) LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz)
{ {
LbmFloat avgrho = 0.0; LbmFloat avgrho = 0.0;
LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0; LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
Context not available.
if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue; if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue;
if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & (CFFluid|CFInter)) ){ if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & (CFFluid|CFInter)) ){
//((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) && //((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
//(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) ) { //(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) ) {
cellcnt += 1.0; cellcnt += 1.0;
for(int rl=0; rl< cDfNum ; ++rl) { for(int rl=0; rl< cDfNum ; ++rl) {
LbmFloat nbdf = QCELL_NB(level,ei,ej,ek, workSet,nbl, rl); LbmFloat nbdf = QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
avgux += (dfDvecX[rl]*nbdf); avgux += (dfDvecX[rl]*nbdf);
avguy += (dfDvecY[rl]*nbdf); avguy += (dfDvecY[rl]*nbdf);
avguz += (dfDvecZ[rl]*nbdf); avguz += (dfDvecZ[rl]*nbdf);
avgrho += nbdf; avgrho += nbdf;
} }
} }
Context not available.
avgux = avguy = avguz = 0.0; avgux = avguy = avguz = 0.0;
//TTT mNumProblems++; //TTT mNumProblems++;
#if ELBEEM_PLUGIN!=1 #if ELBEEM_PLUGIN!=1
//mPanic=1; //mPanic=1;
// this can happen for worst case moving obj scenarios... // this can happen for worst case moving obj scenarios...
errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek)); errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek));
#endif // ELBEEM_PLUGIN #endif // ELBEEM_PLUGIN
Context not available.
//! lbm factory functions //! lbm factory functions
LbmSolverInterface* createSolver() { return new LbmFsgrSolver(); } LbmSolverInterface* createSolver() { return new LbmFsgrSolver(); }
Context not available.