5745f99dee
a) Enable the possibility to remove the "air bubble" around submerged collision object. This feature is enabled as standard for new files. The code was found in elbeem by nudelZ, coded and provided by Nils Thürey (thanks!) b) Old baked files gets deleted if a new bake gets started (were overwritten before and resulted in weird old bake + new bake mixture) (idea by nudelZ)
2383 lines
84 KiB
C++
2383 lines
84 KiB
C++
/** \file elbeem/intern/solver_init.cpp
|
|
* \ingroup elbeem
|
|
*/
|
|
/******************************************************************************
|
|
*
|
|
* El'Beem - Free Surface Fluid Simulation with the Lattice Boltzmann Method
|
|
* Copyright 2003-2006 Nils Thuerey
|
|
*
|
|
* Standard LBM Factory implementation
|
|
*
|
|
*****************************************************************************/
|
|
|
|
|
|
#include "solver_class.h"
|
|
#include "solver_relax.h"
|
|
// for geo init FGI_ defines
|
|
#include "elbeem.h"
|
|
|
|
// helper for 2d init
|
|
#define SWAPYZ(vec) { \
|
|
const LbmFloat tmp = (vec)[2]; \
|
|
(vec)[2] = (vec)[1]; (vec)[1] = tmp; }
|
|
|
|
|
|
/*****************************************************************************/
|
|
//! common variables
|
|
|
|
/*****************************************************************************/
|
|
/*! 3D implementation D3Q19 */
|
|
#if LBMDIM==3
|
|
|
|
//! how many dimensions?
|
|
const int LbmFsgrSolver::cDimension = 3;
|
|
|
|
// Wi factors for collide step
|
|
const LbmFloat LbmFsgrSolver::cCollenZero = (1.0/3.0);
|
|
const LbmFloat LbmFsgrSolver::cCollenOne = (1.0/18.0);
|
|
const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
|
|
|
|
//! threshold value for filled/emptied cells
|
|
const LbmFloat LbmFsgrSolver::cMagicNr2 = 1.0005;
|
|
const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
|
|
const LbmFloat LbmFsgrSolver::cMagicNr = 1.010001;
|
|
const LbmFloat LbmFsgrSolver::cMagicNrNeg = -0.010001;
|
|
|
|
//! size of a single set of distribution functions
|
|
const int LbmFsgrSolver::cDfNum = 19;
|
|
//! direction vector contain vecs for all spatial dirs, even if not used for LBM model
|
|
const int LbmFsgrSolver::cDirNum = 27;
|
|
|
|
//const string LbmFsgrSolver::dfString[ cDfNum ] = {
|
|
const char* LbmFsgrSolver::dfString[ cDfNum ] = {
|
|
" C", " N"," S"," E"," W"," T"," B",
|
|
"NE","NW","SE","SW",
|
|
"NT","NB","ST","SB",
|
|
"ET","EB","WT","WB"
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfNorm[ cDfNum ] = {
|
|
cDirC, cDirN, cDirS, cDirE, cDirW, cDirT, cDirB,
|
|
cDirNE, cDirNW, cDirSE, cDirSW,
|
|
cDirNT, cDirNB, cDirST, cDirSB,
|
|
cDirET, cDirEB, cDirWT, cDirWB
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfInv[ cDfNum ] = {
|
|
cDirC, cDirS, cDirN, cDirW, cDirE, cDirB, cDirT,
|
|
cDirSW, cDirSE, cDirNW, cDirNE,
|
|
cDirSB, cDirST, cDirNB, cDirNT,
|
|
cDirWB, cDirWT, cDirEB, cDirET
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfRefX[ cDfNum ] = {
|
|
0, 0, 0, 0, 0, 0, 0,
|
|
cDirSE, cDirSW, cDirNE, cDirNW,
|
|
0, 0, 0, 0,
|
|
cDirEB, cDirET, cDirWB, cDirWT
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfRefY[ cDfNum ] = {
|
|
0, 0, 0, 0, 0, 0, 0,
|
|
cDirNW, cDirNE, cDirSW, cDirSE,
|
|
cDirNB, cDirNT, cDirSB, cDirST,
|
|
0, 0, 0, 0
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfRefZ[ cDfNum ] = {
|
|
0, 0, 0, 0, 0, 0, 0,
|
|
0, 0, 0, 0,
|
|
cDirST, cDirSB, cDirNT, cDirNB,
|
|
cDirWT, cDirWB, cDirET, cDirEB
|
|
};
|
|
|
|
// Vector Order 3D:
|
|
// 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
|
|
// 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
|
|
// 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
|
|
|
|
const int LbmFsgrSolver::dfVecX[ cDirNum ] = {
|
|
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,
|
|
};
|
|
const int LbmFsgrSolver::dfVecY[ cDirNum ] = {
|
|
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
|
|
};
|
|
const int LbmFsgrSolver::dfVecZ[ cDirNum ] = {
|
|
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 LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
|
|
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
|
|
};
|
|
const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
|
|
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
|
|
};
|
|
const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
|
|
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
|
|
};
|
|
|
|
/* principal directions */
|
|
const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = {
|
|
1,-1, 0,0, 0,0
|
|
};
|
|
const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = {
|
|
0,0, 1,-1, 0,0
|
|
};
|
|
const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = {
|
|
0,0, 0,0, 1,-1
|
|
};
|
|
|
|
/*! arrays for les model coefficients, inited in lbmsolver constructor */
|
|
LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
|
|
LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
|
|
|
|
|
|
const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= {
|
|
cCollenZero,
|
|
cCollenOne, cCollenOne, cCollenOne,
|
|
cCollenOne, cCollenOne, cCollenOne,
|
|
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo,
|
|
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo,
|
|
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
|
|
};
|
|
|
|
/* precalculated equilibrium dfs, inited in lbmsolver constructor */
|
|
LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
|
|
|
|
#else // end LBMDIM==3 , LBMDIM==2
|
|
|
|
/*****************************************************************************/
|
|
/*! 2D implementation D2Q9 */
|
|
|
|
//! how many dimensions?
|
|
const int LbmFsgrSolver::cDimension = 2;
|
|
|
|
//! Wi factors for collide step
|
|
const LbmFloat LbmFsgrSolver::cCollenZero = (4.0/9.0);
|
|
const LbmFloat LbmFsgrSolver::cCollenOne = (1.0/9.0);
|
|
const LbmFloat LbmFsgrSolver::cCollenSqrtTwo = (1.0/36.0);
|
|
|
|
//! threshold value for filled/emptied cells
|
|
const LbmFloat LbmFsgrSolver::cMagicNr2 = 1.0005;
|
|
const LbmFloat LbmFsgrSolver::cMagicNr2Neg = -0.0005;
|
|
const LbmFloat LbmFsgrSolver::cMagicNr = 1.010001;
|
|
const LbmFloat LbmFsgrSolver::cMagicNrNeg = -0.010001;
|
|
|
|
//! size of a single set of distribution functions
|
|
const int LbmFsgrSolver::cDfNum = 9;
|
|
const int LbmFsgrSolver::cDirNum = 9;
|
|
|
|
//const string LbmFsgrSolver::dfString[ cDfNum ] = {
|
|
const char* LbmFsgrSolver::dfString[ cDfNum ] = {
|
|
" C",
|
|
" N", " S", " E", " W",
|
|
"NE", "NW", "SE","SW"
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfNorm[ cDfNum ] = {
|
|
cDirC,
|
|
cDirN, cDirS, cDirE, cDirW,
|
|
cDirNE, cDirNW, cDirSE, cDirSW
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfInv[ cDfNum ] = {
|
|
cDirC,
|
|
cDirS, cDirN, cDirW, cDirE,
|
|
cDirSW, cDirSE, cDirNW, cDirNE
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfRefX[ cDfNum ] = {
|
|
0,
|
|
0, 0, 0, 0,
|
|
cDirSE, cDirSW, cDirNE, cDirNW
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfRefY[ cDfNum ] = {
|
|
0,
|
|
0, 0, 0, 0,
|
|
cDirNW, cDirNE, cDirSW, cDirSE
|
|
};
|
|
|
|
const int LbmFsgrSolver::dfRefZ[ cDfNum ] = {
|
|
0, 0, 0, 0, 0,
|
|
0, 0, 0, 0
|
|
};
|
|
|
|
// Vector Order 2D:
|
|
// 0 1 2 3 4 5 6 7 8
|
|
// 0, 0,0, 1,-1, 1,-1,1,-1
|
|
// 0, 1,-1, 0,0, 1,1,-1,-1
|
|
|
|
const int LbmFsgrSolver::dfVecX[ cDirNum ] = {
|
|
0,
|
|
0,0, 1,-1,
|
|
1,-1,1,-1
|
|
};
|
|
const int LbmFsgrSolver::dfVecY[ cDirNum ] = {
|
|
0,
|
|
1,-1, 0,0,
|
|
1,1,-1,-1
|
|
};
|
|
const int LbmFsgrSolver::dfVecZ[ cDirNum ] = {
|
|
0, 0,0,0,0, 0,0,0,0
|
|
};
|
|
|
|
const LbmFloat LbmFsgrSolver::dfDvecX[ cDirNum ] = {
|
|
0,
|
|
0,0, 1,-1,
|
|
1,-1,1,-1
|
|
};
|
|
const LbmFloat LbmFsgrSolver::dfDvecY[ cDirNum ] = {
|
|
0,
|
|
1,-1, 0,0,
|
|
1,1,-1,-1
|
|
};
|
|
const LbmFloat LbmFsgrSolver::dfDvecZ[ cDirNum ] = {
|
|
0, 0,0,0,0, 0,0,0,0
|
|
};
|
|
|
|
const int LbmFsgrSolver::princDirX[ 2*LbmFsgrSolver::cDimension ] = {
|
|
1,-1, 0,0
|
|
};
|
|
const int LbmFsgrSolver::princDirY[ 2*LbmFsgrSolver::cDimension ] = {
|
|
0,0, 1,-1
|
|
};
|
|
const int LbmFsgrSolver::princDirZ[ 2*LbmFsgrSolver::cDimension ] = {
|
|
0,0, 0,0
|
|
};
|
|
|
|
|
|
/*! arrays for les model coefficients, inited in lbmsolver constructor */
|
|
LbmFloat LbmFsgrSolver::lesCoeffDiag[ (cDimension-1)*(cDimension-1) ][ cDirNum ];
|
|
LbmFloat LbmFsgrSolver::lesCoeffOffdiag[ cDimension ][ cDirNum ];
|
|
|
|
|
|
const LbmFloat LbmFsgrSolver::dfLength[ cDfNum ]= {
|
|
cCollenZero,
|
|
cCollenOne, cCollenOne, cCollenOne, cCollenOne,
|
|
cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo, cCollenSqrtTwo
|
|
};
|
|
|
|
/* precalculated equilibrium dfs, inited in lbmsolver constructor */
|
|
LbmFloat LbmFsgrSolver::dfEquil[ dTotalNum ];
|
|
|
|
// D2Q9 end
|
|
#endif // LBMDIM==2
|
|
|
|
|
|
// required globals
|
|
extern bool glob_mpactive;
|
|
extern int glob_mpnum, glob_mpindex;
|
|
|
|
|
|
/******************************************************************************
|
|
* Lbm Constructor
|
|
*****************************************************************************/
|
|
LbmFsgrSolver::LbmFsgrSolver() :
|
|
//D(),
|
|
mCurrentMass(0.0), mCurrentVolume(0.0),
|
|
mNumProblems(0),
|
|
mAvgMLSUPS(0.0), mAvgMLSUPSCnt(0.0),
|
|
mpPreviewSurface(NULL),
|
|
mTimeAdap(true), mForceTimeStepReduce(false),
|
|
mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
|
|
mInitSurfaceSmoothing(0), mFsSurfGenSetting(0),
|
|
mTimestepReduceLock(0),
|
|
mTimeSwitchCounts(0), mTimeMaxvelStepCnt(0),
|
|
mSimulationTime(0.0), mLastSimTime(0.0),
|
|
mMinTimestep(0.0), mMaxTimestep(0.0),
|
|
mMaxNoCells(0), mMinNoCells(0), mAvgNumUsedCells(0),
|
|
mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
|
|
mMOIVertices(), mMOIVerticesOld(), mMOINormals(),
|
|
mIsoWeightMethod(1),
|
|
mMaxRefine(1),
|
|
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
|
|
mDebugOmegaRet(0.0),
|
|
mLastOmega(1e10), mLastGravity(1e10),
|
|
mNumInvIfTotal(0), mNumFsgrChanges(0),
|
|
mDisableStandingFluidInit(0),
|
|
mInit2dYZ(false),
|
|
mForceTadapRefine(-1), mCutoff(-1)
|
|
{
|
|
mpControl = new LbmControlData();
|
|
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
mpTest = new LbmTestdata();
|
|
mMpNum = mMpIndex = 0;
|
|
mOrgSizeX = 0;
|
|
mOrgStartX = 0.;
|
|
mOrgEndX = 0.;
|
|
#endif // LBM_INCLUDE_TESTSOLVERS!=1
|
|
mpIso = new IsoSurface( mIsoValue );
|
|
|
|
// init equilibrium dist. func
|
|
LbmFloat rho=1.0;
|
|
FORDF0 {
|
|
dfEquil[l] = this->getCollideEq( l,rho, 0.0, 0.0, 0.0);
|
|
}
|
|
dfEquil[dMass] = 1.;
|
|
dfEquil[dFfrac] = 1.;
|
|
dfEquil[dFlux] = FLUX_INIT;
|
|
|
|
// init LES
|
|
int odm = 0;
|
|
for(int m=0; m<LBMDIM; m++) {
|
|
for(int l=0; l<cDfNum; l++) {
|
|
this->lesCoeffDiag[m][l] =
|
|
this->lesCoeffOffdiag[m][l] = 0.0;
|
|
}
|
|
}
|
|
for(int m=0; m<LBMDIM; m++) {
|
|
for(int n=0; n<LBMDIM; n++) {
|
|
for(int l=1; l<cDfNum; l++) {
|
|
LbmFloat em;
|
|
switch(m) {
|
|
case 0: em = dfDvecX[l]; break;
|
|
case 1: em = dfDvecY[l]; break;
|
|
case 2: em = dfDvecZ[l]; break;
|
|
default: em = -1.0; errFatal("SMAGO1","err m="<<m, SIMWORLD_GENERICERROR);
|
|
}
|
|
LbmFloat en;
|
|
switch(n) {
|
|
case 0: en = dfDvecX[l]; break;
|
|
case 1: en = dfDvecY[l]; break;
|
|
case 2: en = dfDvecZ[l]; break;
|
|
default: en = -1.0; errFatal("SMAGO2","err n="<<n, SIMWORLD_GENERICERROR);
|
|
}
|
|
const LbmFloat coeff = em*en;
|
|
if(m==n) {
|
|
this->lesCoeffDiag[m][l] = coeff;
|
|
} else {
|
|
if(m>n) {
|
|
this->lesCoeffOffdiag[odm][l] = coeff;
|
|
}
|
|
}
|
|
}
|
|
|
|
if(m==n) {
|
|
} else {
|
|
if(m>n) odm++;
|
|
}
|
|
}
|
|
}
|
|
|
|
mDvecNrm[0] = LbmVec(0.0);
|
|
FORDF1 {
|
|
mDvecNrm[l] = getNormalized(
|
|
LbmVec(dfDvecX[dfInv[l]], dfDvecY[dfInv[l]], dfDvecZ[dfInv[l]] )
|
|
) * -1.0;
|
|
}
|
|
|
|
// calculate gauss weights for restriction
|
|
//LbmFloat mGaussw[27];
|
|
LbmFloat totGaussw = 0.0;
|
|
const LbmFloat alpha = 1.0;
|
|
const LbmFloat gw = sqrt(2.0*LBMDIM);
|
|
#if ELBEEM_PLUGIN!=1
|
|
errMsg("coarseRestrictFromFine", "TCRFF_DFDEBUG2 test df/dir num!");
|
|
#endif
|
|
for(int n=0;(n<cDirNum); n++) { mGaussw[n] = 0.0; }
|
|
//for(int n=0;(n<cDirNum); n++) {
|
|
for(int n=0;(n<cDfNum); n++) {
|
|
const LbmFloat d = norm(LbmVec(dfVecX[n], dfVecY[n], dfVecZ[n]));
|
|
LbmFloat w = expf( -alpha*d*d ) - expf( -alpha*gw*gw );
|
|
mGaussw[n] = w;
|
|
totGaussw += w;
|
|
}
|
|
for(int n=0;(n<cDirNum); n++) {
|
|
mGaussw[n] = mGaussw[n]/totGaussw;
|
|
}
|
|
|
|
}
|
|
|
|
/*****************************************************************************/
|
|
/* Destructor */
|
|
/*****************************************************************************/
|
|
LbmFsgrSolver::~LbmFsgrSolver()
|
|
{
|
|
if(!mInitDone){ debMsgStd("LbmFsgrSolver::LbmFsgrSolver",DM_MSG,"not inited...",0); return; }
|
|
#if COMPRESSGRIDS==1
|
|
delete [] mLevel[mMaxRefine].mprsCells[1];
|
|
mLevel[mMaxRefine].mprsCells[0] = mLevel[mMaxRefine].mprsCells[1] = NULL;
|
|
#endif // COMPRESSGRIDS==1
|
|
|
|
for(int i=0; i<=mMaxRefine; i++) {
|
|
for(int s=0; s<2; s++) {
|
|
if(mLevel[i].mprsCells[s]) delete [] mLevel[i].mprsCells[s];
|
|
if(mLevel[i].mprsFlags[s]) delete [] mLevel[i].mprsFlags[s];
|
|
}
|
|
}
|
|
delete mpIso;
|
|
if(mpPreviewSurface) delete mpPreviewSurface;
|
|
// cleanup done during scene deletion...
|
|
|
|
if(mpControl) delete mpControl;
|
|
|
|
// always output performance estimate
|
|
debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG," Avg. MLSUPS:"<<(mAvgMLSUPS/mAvgMLSUPSCnt), 5);
|
|
if(!mSilent) debMsgStd("LbmFsgrSolver::~LbmFsgrSolver",DM_MSG,"Deleted...",10);
|
|
}
|
|
|
|
|
|
|
|
|
|
/******************************************************************************
|
|
* initilize variables fom attribute list
|
|
*****************************************************************************/
|
|
void LbmFsgrSolver::parseAttrList()
|
|
{
|
|
LbmSolverInterface::parseStdAttrList();
|
|
|
|
string matIso("default");
|
|
matIso = mpSifAttrs->readString("material_surf", matIso, "SimulationLbm","mpIso->material", false );
|
|
mpIso->setMaterialName( matIso );
|
|
mOutputSurfacePreview = mpSifAttrs->readInt("surfacepreview", mOutputSurfacePreview, "SimulationLbm","mOutputSurfacePreview", false );
|
|
mTimeAdap = mpSifAttrs->readBool("timeadap", mTimeAdap, "SimulationLbm","mTimeAdap", false );
|
|
mDomainBound = mpSifAttrs->readString("domainbound", mDomainBound, "SimulationLbm","mDomainBound", false );
|
|
mDomainPartSlipValue = mpSifAttrs->readFloat("domainpartslip", mDomainPartSlipValue, "SimulationLbm","mDomainPartSlipValue", false );
|
|
|
|
mIsoWeightMethod= mpSifAttrs->readInt("isoweightmethod", mIsoWeightMethod, "SimulationLbm","mIsoWeightMethod", false );
|
|
mInitSurfaceSmoothing = mpSifAttrs->readInt("initsurfsmooth", mInitSurfaceSmoothing, "SimulationLbm","mInitSurfaceSmoothing", false );
|
|
mSmoothSurface = mpSifAttrs->readFloat("smoothsurface", mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
|
|
mSmoothNormals = mpSifAttrs->readFloat("smoothnormals", mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
|
|
mFsSurfGenSetting = mpSifAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );
|
|
|
|
// refinement
|
|
mMaxRefine = mRefinementDesired;
|
|
mMaxRefine = mpSifAttrs->readInt("maxrefine", mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
|
|
if(mMaxRefine<0) mMaxRefine=0;
|
|
if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1;
|
|
mDisableStandingFluidInit = mpSifAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
|
|
mInit2dYZ = mpSifAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
|
|
mForceTadapRefine = mpSifAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
|
|
|
|
// demo mode settings
|
|
mFVHeight = mpSifAttrs->readFloat("fvolheight", mFVHeight, "LbmFsgrSolver","mFVHeight", false );
|
|
// FIXME check needed?
|
|
mFVArea = mpSifAttrs->readFloat("fvolarea", mFVArea, "LbmFsgrSolver","mFArea", false );
|
|
|
|
// debugging - skip some time...
|
|
double starttimeskip = 0.;
|
|
starttimeskip = mpSifAttrs->readFloat("forcestarttimeskip", starttimeskip, "LbmFsgrSolver","starttimeskip", false );
|
|
mSimulationTime += starttimeskip;
|
|
if(starttimeskip>0.) debMsgStd("LbmFsgrSolver::parseStdAttrList",DM_NOTIFY,"Used starttimeskip="<<starttimeskip<<", t="<<mSimulationTime, 1);
|
|
|
|
mpControl->parseControldataAttrList(mpSifAttrs);
|
|
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
mUseTestdata = 0;
|
|
mUseTestdata = mpSifAttrs->readBool("use_testdata", mUseTestdata,"LbmFsgrSolver", "mUseTestdata", false);
|
|
mpTest->parseTestdataAttrList(mpSifAttrs);
|
|
#ifdef ELBEEM_PLUGIN
|
|
mUseTestdata=1; // DEBUG
|
|
#endif // ELBEEM_PLUGIN
|
|
|
|
mMpNum = mpSifAttrs->readInt("mpnum", mMpNum ,"LbmFsgrSolver", "mMpNum", false);
|
|
mMpIndex = mpSifAttrs->readInt("mpindex", mMpIndex ,"LbmFsgrSolver", "mMpIndex", false);
|
|
if(glob_mpactive) {
|
|
// used instead...
|
|
mMpNum = glob_mpnum;
|
|
mMpIndex = glob_mpindex;
|
|
} else {
|
|
glob_mpnum = mMpNum;
|
|
glob_mpindex = 0;
|
|
}
|
|
errMsg("LbmFsgrSolver::parseAttrList"," mpactive:"<<glob_mpactive<<", "<<glob_mpindex<<"/"<<glob_mpnum);
|
|
if(mMpNum>0) {
|
|
mUseTestdata=1; // needed in this case...
|
|
}
|
|
|
|
errMsg("LbmFsgrSolver::LBM_INCLUDE_TESTSOLVERS","Active, mUseTestdata:"<<mUseTestdata<<" ");
|
|
#else // LBM_INCLUDE_TESTSOLVERS!=1
|
|
// not testsolvers, off by default
|
|
mUseTestdata = 0;
|
|
if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
|
|
#endif // LBM_INCLUDE_TESTSOLVERS!=1
|
|
|
|
mInitialCsmago = mpSifAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
|
|
// deprecated!
|
|
float mInitialCsmagoCoarse = 0.0;
|
|
mInitialCsmagoCoarse = mpSifAttrs->readFloat("csmago_coarse", mInitialCsmagoCoarse, "SimulationLbm","mInitialCsmagoCoarse", false );
|
|
#if USE_LES==1
|
|
#else // USE_LES==1
|
|
debMsgStd("LbmFsgrSolver", DM_WARNING, "LES model switched off!",2);
|
|
mInitialCsmago = 0.0;
|
|
#endif // USE_LES==1
|
|
}
|
|
|
|
|
|
/******************************************************************************
|
|
* (part of enabling chapter 6 of "Free Surface Flows with Moving and Deforming Objects for LBM")
|
|
*****************************************************************************/
|
|
void LbmFsgrSolver::setSurfGenSettings(short value)
|
|
{
|
|
mFsSurfGenSetting = value;
|
|
}
|
|
|
|
|
|
/******************************************************************************
|
|
* Initialize omegas and forces on all levels (for init/timestep change)
|
|
*****************************************************************************/
|
|
void LbmFsgrSolver::initLevelOmegas()
|
|
{
|
|
// no explicit settings
|
|
mOmega = mpParam->calculateOmega(mSimulationTime);
|
|
mGravity = vec2L( mpParam->calculateGravity(mSimulationTime) );
|
|
mSurfaceTension = 0.; //mpParam->calculateSurfaceTension(); // unused
|
|
if(mInit2dYZ) { SWAPYZ(mGravity); }
|
|
|
|
// check if last init was ok
|
|
LbmFloat gravDelta = norm(mGravity-mLastGravity);
|
|
//errMsg("ChannelAnimDebug","t:"<<mSimulationTime<<" om:"<<mOmega<<" - lom:"<<mLastOmega<<" gv:"<<mGravity<<" - "<<mLastGravity<<" , "<<gravDelta );
|
|
if((mOmega == mLastOmega) && (gravDelta<=0.0)) return;
|
|
|
|
if(mInitialCsmago<=0.0) {
|
|
if(OPT3D==1) {
|
|
errFatal("LbmFsgrSolver::initLevelOmegas","Csmago-LES = 0 not supported for optimized 3D version...",SIMWORLD_INITERROR);
|
|
return;
|
|
}
|
|
}
|
|
|
|
LbmFloat fineCsmago = mInitialCsmago;
|
|
LbmFloat coarseCsmago = mInitialCsmago;
|
|
LbmFloat maxFineCsmago1 = 0.026;
|
|
LbmFloat maxCoarseCsmago1 = 0.029; // try stabilizing
|
|
LbmFloat maxFineCsmago2 = 0.028;
|
|
LbmFloat maxCoarseCsmago2 = 0.032; // try stabilizing some more
|
|
if((mMaxRefine==1)&&(mInitialCsmago<maxFineCsmago1)) {
|
|
fineCsmago = maxFineCsmago1;
|
|
coarseCsmago = maxCoarseCsmago1;
|
|
}
|
|
if((mMaxRefine>1)&&(mInitialCsmago<maxFineCsmago2)) {
|
|
fineCsmago = maxFineCsmago2;
|
|
coarseCsmago = maxCoarseCsmago2;
|
|
}
|
|
|
|
|
|
// use Tau instead of Omega for calculations
|
|
{ // init base level
|
|
int i = mMaxRefine;
|
|
mLevel[i].omega = mOmega;
|
|
mLevel[i].timestep = mpParam->getTimestep();
|
|
mLevel[i].lcsmago = fineCsmago; //CSMAGO_INITIAL;
|
|
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);
|
|
}
|
|
|
|
// init all sub levels
|
|
for(int i=mMaxRefine-1; i>=0; i--) {
|
|
//mLevel[i].omega = 2.0 * (mLevel[i+1].omega-0.5) + 0.5;
|
|
double nomega = 0.5 * ( (1.0/(double)mLevel[i+1].omega) -0.5) + 0.5;
|
|
nomega = 1.0/nomega;
|
|
mLevel[i].omega = (LbmFloat)nomega;
|
|
mLevel[i].timestep = 2.0 * mLevel[i+1].timestep;
|
|
mLevel[i].lcsmago = coarseCsmago;
|
|
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);
|
|
}
|
|
|
|
// for lbgk
|
|
mLevel[ mMaxRefine ].gravity = mGravity / mLevel[ mMaxRefine ].omega;
|
|
for(int i=mMaxRefine-1; i>=0; i--) {
|
|
// should be the same on all levels...
|
|
// for lbgk
|
|
mLevel[i].gravity = (mLevel[i+1].gravity * mLevel[i+1].omega) * 2.0 / mLevel[i].omega;
|
|
}
|
|
|
|
mLastOmega = mOmega;
|
|
mLastGravity = mGravity;
|
|
// debug? invalidate old values...
|
|
mGravity = -100.0;
|
|
mOmega = -100.0;
|
|
|
|
for(int i=0; i<=mMaxRefine; i++) {
|
|
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
|
|
<<" omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity<< ", "
|
|
<<" cmsagp:"<<mLevel[i].lcsmago<<", "
|
|
<< " ss"<<mLevel[i].timestep<<" ns"<<mLevel[i].nodeSize<<" cs"<<mLevel[i].simCellSize );
|
|
} else {
|
|
if(!mInitDone) {
|
|
debMsgStd("LbmFsgrSolver", DM_MSG, "Level init "<<i<<" - sizes:"<<mLevel[i].lSizex<<","<<mLevel[i].lSizey<<","<<mLevel[i].lSizez<<" "
|
|
<<"omega:"<<mLevel[i].omega<<" grav:"<<mLevel[i].gravity , 5);
|
|
}
|
|
}
|
|
}
|
|
if(mMaxRefine>0) {
|
|
mDfScaleUp = (mLevel[0 ].timestep/mLevel[0+1].timestep)* (1.0/mLevel[0 ].omega-1.0)/ (1.0/mLevel[0+1].omega-1.0); // yu
|
|
mDfScaleDown = (mLevel[0+1].timestep/mLevel[0 ].timestep)* (1.0/mLevel[0+1].omega-1.0)/ (1.0/mLevel[0 ].omega-1.0); // yu
|
|
}
|
|
}
|
|
|
|
|
|
/******************************************************************************
|
|
* Init Solver (values should be read from config file)
|
|
*****************************************************************************/
|
|
|
|
/*! finish the init with config file values (allocate arrays...) */
|
|
bool LbmFsgrSolver::initializeSolverMemory()
|
|
{
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init start... "<<mInitDone<<" "<<(void*)this,1);
|
|
|
|
// init cppf stage
|
|
if(mCppfStage>0) {
|
|
mSizex *= mCppfStage;
|
|
mSizey *= mCppfStage;
|
|
mSizez *= mCppfStage;
|
|
}
|
|
if(mFsSurfGenSetting==-1) {
|
|
// all on
|
|
mFsSurfGenSetting =
|
|
fssgNormal | fssgNoNorth | fssgNoSouth | fssgNoEast |
|
|
fssgNoWest | fssgNoTop | fssgNoBottom | fssgNoObs ;
|
|
}
|
|
|
|
// size inits to force cubic cells and mult4 level dimensions
|
|
// and make sure we dont allocate too much...
|
|
bool memOk = false;
|
|
int orgSx = mSizex;
|
|
int orgSy = mSizey;
|
|
int orgSz = mSizez;
|
|
double sizeReduction = 1.0;
|
|
double memEstFromFunc = -1.0;
|
|
double memEstFine = -1.0;
|
|
string memreqStr("");
|
|
bool firstMInit = true;
|
|
int minitTries=0;
|
|
while(!memOk) {
|
|
minitTries++;
|
|
initGridSizes( mSizex, mSizey, mSizez,
|
|
mvGeoStart, mvGeoEnd, mMaxRefine, PARALLEL);
|
|
|
|
// MPT
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
if(firstMInit) {
|
|
mrSetup();
|
|
}
|
|
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
|
firstMInit=false;
|
|
|
|
calculateMemreqEstimate( mSizex, mSizey, mSizez,
|
|
mMaxRefine, mFarFieldSize, &memEstFromFunc, &memEstFine, &memreqStr );
|
|
|
|
double memLimit;
|
|
string memLimStr("-");
|
|
if(sizeof(void*)==4) {
|
|
// 32bit system, limit to 2GB
|
|
memLimit = 2.0* 1024.0*1024.0*1024.0;
|
|
memLimStr = string("2GB");
|
|
} else {
|
|
// 64bit, just take 16GB as limit for now...
|
|
memLimit = 16.0* 1024.0*1024.0*1024.0;
|
|
memLimStr = string("16GB");
|
|
}
|
|
|
|
// restrict max. chunk of 1 mem block to 1GB for windos
|
|
bool memBlockAllocProblem = false;
|
|
double maxDefaultMemChunk = 2.*1024.*1024.*1024.;
|
|
//std::cerr<<" memEstFine "<< memEstFine <<" maxWin:" <<maxWinMemChunk <<" maxMac:" <<maxMacMemChunk ; // DEBUG
|
|
#ifdef WIN32
|
|
double maxWinMemChunk = 1100.*1024.*1024.;
|
|
if(sizeof(void *)==4 && memEstFine>maxWinMemChunk) {
|
|
memBlockAllocProblem = true;
|
|
}
|
|
#endif // WIN32
|
|
#ifdef __APPLE__
|
|
double maxMacMemChunk = 1200.*1024.*1024.;
|
|
if(memEstFine> maxMacMemChunk) {
|
|
memBlockAllocProblem = true;
|
|
}
|
|
#endif // Mac
|
|
if(sizeof(void*)==4 && memEstFine>maxDefaultMemChunk) {
|
|
// max memory chunk for 32bit systems 2gig
|
|
memBlockAllocProblem = true;
|
|
}
|
|
|
|
if(memEstFromFunc>memLimit || memBlockAllocProblem) {
|
|
sizeReduction *= 0.9;
|
|
mSizex = (int)(orgSx * sizeReduction);
|
|
mSizey = (int)(orgSy * sizeReduction);
|
|
mSizez = (int)(orgSz * sizeReduction);
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_WARNING,"initGridSizes: memory limit exceeded "<<
|
|
//memEstFromFunc<<"/"<<memLimit<<", "<<
|
|
//memEstFine<<"/"<<maxWinMemChunk<<", "<<
|
|
memreqStr<<"/"<<memLimStr<<", "<<
|
|
"retrying: "<<PRINT_VEC(mSizex,mSizey,mSizez)<<" org:"<<PRINT_VEC(orgSx,orgSy,orgSz)
|
|
, 3 );
|
|
} else {
|
|
memOk = true;
|
|
}
|
|
}
|
|
|
|
mPreviewFactor = (LbmFloat)mOutputSurfacePreview / (LbmFloat)mSizex;
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"initGridSizes: Final domain size X:"<<mSizex<<" Y:"<<mSizey<<" Z:"<<mSizez<<
|
|
", Domain: "<<mvGeoStart<<":"<<mvGeoEnd<<", "<<(mvGeoEnd-mvGeoStart)<<
|
|
", PointerSize: "<< sizeof(void*) << ", IntSize: "<< sizeof(int) <<
|
|
", est. Mem.Req.: "<<memreqStr ,2);
|
|
mpParam->setSize(mSizex, mSizey, mSizez);
|
|
if((minitTries>1)&&(glob_mpnum)) { errMsg("LbmFsgrSolver::initialize","Warning!!!!!!!!!!!!!!! Original gridsize changed........."); }
|
|
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Definitions: "
|
|
<<"LBM_EPSILON="<<LBM_EPSILON <<" "
|
|
<<"FSGR_STRICT_DEBUG="<<FSGR_STRICT_DEBUG <<" "
|
|
<<"OPT3D="<<OPT3D <<" "
|
|
<<"COMPRESSGRIDS="<<COMPRESSGRIDS<<" "
|
|
<<"MASS_INVALID="<<MASS_INVALID <<" "
|
|
<<"FSGR_LISTTRICK="<<FSGR_LISTTRICK <<" "
|
|
<<"FSGR_LISTTTHRESHEMPTY="<<FSGR_LISTTTHRESHEMPTY <<" "
|
|
<<"FSGR_LISTTTHRESHFULL="<<FSGR_LISTTTHRESHFULL <<" "
|
|
<<"FSGR_MAGICNR="<<FSGR_MAGICNR <<" "
|
|
<<"USE_LES="<<USE_LES <<" "
|
|
,10);
|
|
|
|
// perform 2D corrections...
|
|
if(LBMDIM == 2) mSizez = 1;
|
|
|
|
mpParam->setSimulationMaxSpeed(0.0);
|
|
if(mFVHeight>0.0) mpParam->setFluidVolumeHeight(mFVHeight);
|
|
mpParam->setTadapLevels( mMaxRefine+1 );
|
|
|
|
if(mForceTadapRefine>mMaxRefine) {
|
|
mpParam->setTadapLevels( mForceTadapRefine+1 );
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Forcing a t-adap refine level of "<<mForceTadapRefine, 6);
|
|
}
|
|
|
|
if(!mpParam->calculateAllMissingValues(mSimulationTime, false)) {
|
|
errFatal("LbmFsgrSolver::initialize","Fatal: failed to init parameters! Aborting...",SIMWORLD_INITERROR);
|
|
return false;
|
|
}
|
|
|
|
|
|
// init vectors
|
|
for(int i=0; i<=mMaxRefine; i++) {
|
|
mLevel[i].id = i;
|
|
mLevel[i].nodeSize = 0.0;
|
|
mLevel[i].simCellSize = 0.0;
|
|
mLevel[i].omega = 0.0;
|
|
mLevel[i].time = 0.0;
|
|
mLevel[i].timestep = 1.0;
|
|
mLevel[i].gravity = LbmVec(0.0);
|
|
mLevel[i].mprsCells[0] = NULL;
|
|
mLevel[i].mprsCells[1] = NULL;
|
|
mLevel[i].mprsFlags[0] = NULL;
|
|
mLevel[i].mprsFlags[1] = NULL;
|
|
|
|
mLevel[i].avgOmega = 0.0;
|
|
mLevel[i].avgOmegaCnt = 0.0;
|
|
}
|
|
|
|
// init sizes
|
|
mLevel[mMaxRefine].lSizex = mSizex;
|
|
mLevel[mMaxRefine].lSizey = mSizey;
|
|
mLevel[mMaxRefine].lSizez = mSizez;
|
|
for(int i=mMaxRefine-1; i>=0; i--) {
|
|
mLevel[i].lSizex = mLevel[i+1].lSizex/2;
|
|
mLevel[i].lSizey = mLevel[i+1].lSizey/2;
|
|
mLevel[i].lSizez = mLevel[i+1].lSizez/2;
|
|
}
|
|
|
|
// safety check
|
|
if(sizeof(CellFlagType) != CellFlagTypeSize) {
|
|
errFatal("LbmFsgrSolver::initialize","Fatal Error: CellFlagType has wrong size! Is:"<<sizeof(CellFlagType)<<", should be:"<<CellFlagTypeSize, SIMWORLD_GENERICERROR);
|
|
return false;
|
|
}
|
|
|
|
double ownMemCheck = 0.0;
|
|
mLevel[ mMaxRefine ].nodeSize = ((mvGeoEnd[0]-mvGeoStart[0]) / (LbmFloat)(mSizex));
|
|
mLevel[ mMaxRefine ].simCellSize = mpParam->getCellSize();
|
|
mLevel[ mMaxRefine ].lcellfactor = 1.0;
|
|
LONGINT rcellSize = ((mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*mLevel[mMaxRefine].lSizez) *dTotalNum);
|
|
|
|
#if COMPRESSGRIDS==0
|
|
mLevel[ mMaxRefine ].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
|
|
mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
|
|
ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
|
|
#else // COMPRESSGRIDS==0
|
|
LONGINT compressOffset = (mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum*2);
|
|
// D int tmp = ( (rcellSize +compressOffset +4)/(1024*1024) )*4;
|
|
// D printf("Debug MEMMMM excee: %d\n", tmp);
|
|
mLevel[ mMaxRefine ].mprsCells[1] = new LbmFloat[ rcellSize +compressOffset +4 ];
|
|
mLevel[ mMaxRefine ].mprsCells[0] = mLevel[ mMaxRefine ].mprsCells[1]+compressOffset;
|
|
ownMemCheck += sizeof(LbmFloat) * (rcellSize +compressOffset +4);
|
|
#endif // COMPRESSGRIDS==0
|
|
|
|
if(!mLevel[ mMaxRefine ].mprsCells[1] || !mLevel[ mMaxRefine ].mprsCells[0]) {
|
|
errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (1)! Aborting...",SIMWORLD_INITERROR);
|
|
return false;
|
|
}
|
|
|
|
// +4 for safety ?
|
|
mLevel[ mMaxRefine ].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
|
|
mLevel[ mMaxRefine ].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
|
|
ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
|
|
if(!mLevel[ mMaxRefine ].mprsFlags[1] || !mLevel[ mMaxRefine ].mprsFlags[0]) {
|
|
errFatal("LbmFsgrSolver::initialize","Fatal: Couldnt allocate memory (2)! Aborting...",SIMWORLD_INITERROR);
|
|
|
|
#if COMPRESSGRIDS==0
|
|
delete[] mLevel[ mMaxRefine ].mprsCells[0];
|
|
delete[] mLevel[ mMaxRefine ].mprsCells[1];
|
|
#else // COMPRESSGRIDS==0
|
|
delete[] mLevel[ mMaxRefine ].mprsCells[1];
|
|
#endif // COMPRESSGRIDS==0
|
|
return false;
|
|
}
|
|
|
|
LbmFloat lcfdimFac = 8.0;
|
|
if(LBMDIM==2) lcfdimFac = 4.0;
|
|
for(int i=mMaxRefine-1; i>=0; i--) {
|
|
mLevel[i].nodeSize = 2.0 * mLevel[i+1].nodeSize;
|
|
mLevel[i].simCellSize = 2.0 * mLevel[i+1].simCellSize;
|
|
mLevel[i].lcellfactor = mLevel[i+1].lcellfactor * lcfdimFac;
|
|
|
|
if(LBMDIM==2){ mLevel[i].lSizez = 1; } // 2D
|
|
rcellSize = ((mLevel[i].lSizex*mLevel[i].lSizey*mLevel[i].lSizez) *dTotalNum);
|
|
mLevel[i].mprsFlags[0] = new CellFlagType[ rcellSize/dTotalNum +4 ];
|
|
mLevel[i].mprsFlags[1] = new CellFlagType[ rcellSize/dTotalNum +4 ];
|
|
ownMemCheck += 2 * sizeof(CellFlagType) * (rcellSize/dTotalNum +4);
|
|
mLevel[i].mprsCells[0] = new LbmFloat[ rcellSize +4 ];
|
|
mLevel[i].mprsCells[1] = new LbmFloat[ rcellSize +4 ];
|
|
ownMemCheck += 2 * sizeof(LbmFloat) * (rcellSize+4);
|
|
}
|
|
|
|
// isosurface memory, use orig res values
|
|
if(mFarFieldSize>0.) {
|
|
ownMemCheck += (double)( (3*sizeof(int)+sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
|
|
} else {
|
|
// ignore 3 int slices...
|
|
ownMemCheck += (double)( ( sizeof(float)) * ((mSizex+2)*(mSizey+2)*(mSizez+2)) );
|
|
}
|
|
|
|
// sanity check
|
|
#if ELBEEM_PLUGIN!=1
|
|
if(ABS(1.0-ownMemCheck/memEstFromFunc)>0.01) {
|
|
errMsg("LbmFsgrSolver::initialize","Sanity Error - memory estimate is off! real:"<<ownMemCheck<<" vs. estimate:"<<memEstFromFunc );
|
|
}
|
|
#endif // ELBEEM_PLUGIN!=1
|
|
|
|
// init sizes for _all_ levels
|
|
for(int i=mMaxRefine; i>=0; i--) {
|
|
mLevel[i].lOffsx = mLevel[i].lSizex;
|
|
mLevel[i].lOffsy = mLevel[i].lOffsx*mLevel[i].lSizey;
|
|
mLevel[i].lOffsz = mLevel[i].lOffsy*mLevel[i].lSizez;
|
|
mLevel[i].setCurr = 0;
|
|
mLevel[i].setOther = 1;
|
|
mLevel[i].lsteps = 0;
|
|
mLevel[i].lmass = 0.0;
|
|
mLevel[i].lvolume = 0.0;
|
|
}
|
|
|
|
// calc omega, force for all levels
|
|
initLevelOmegas();
|
|
mMinTimestep = mpParam->getTimestep();
|
|
mMaxTimestep = mpParam->getTimestep();
|
|
|
|
// init isosurf
|
|
mpIso->setIsolevel( mIsoValue );
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
if(mUseTestdata) {
|
|
mpTest->setMaterialName( mpIso->getMaterialName() );
|
|
delete mpIso;
|
|
mpIso = mpTest;
|
|
if(mpTest->mFarfMode>0) { // 3d off
|
|
mpTest->setIsolevel(-100.0);
|
|
} else {
|
|
mpTest->setIsolevel( mIsoValue );
|
|
}
|
|
}
|
|
#endif // LBM_INCLUDE_TESTSOLVERS!=1
|
|
// approximate feature size with mesh resolution
|
|
float featureSize = mLevel[ mMaxRefine ].nodeSize*0.5;
|
|
// smooth vars defined in solver_interface, set by simulation object
|
|
// reset for invalid values...
|
|
if((mSmoothSurface<0.)||(mSmoothSurface>50.)) mSmoothSurface = 1.;
|
|
if((mSmoothNormals<0.)||(mSmoothNormals>50.)) mSmoothNormals = 1.;
|
|
mpIso->setSmoothSurface( mSmoothSurface * featureSize );
|
|
mpIso->setSmoothNormals( mSmoothNormals * featureSize );
|
|
|
|
// init iso weight values mIsoWeightMethod
|
|
int wcnt = 0;
|
|
float totw = 0.0;
|
|
for(int ak=-1;ak<=1;ak++)
|
|
for(int aj=-1;aj<=1;aj++)
|
|
for(int ai=-1;ai<=1;ai++) {
|
|
switch(mIsoWeightMethod) {
|
|
case 1: // light smoothing
|
|
mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
|
|
break;
|
|
case 2: // very light smoothing
|
|
mIsoWeight[wcnt] = sqrt(3.0) - sqrt( (LbmFloat)(ak*ak + aj*aj + ai*ai) );
|
|
mIsoWeight[wcnt] *= mIsoWeight[wcnt];
|
|
break;
|
|
case 3: // no smoothing
|
|
if(ai==0 && aj==0 && ak==0) mIsoWeight[wcnt] = 1.0;
|
|
else mIsoWeight[wcnt] = 0.0;
|
|
break;
|
|
default: // strong smoothing (=0)
|
|
mIsoWeight[wcnt] = 1.0;
|
|
break;
|
|
}
|
|
totw += mIsoWeight[wcnt];
|
|
wcnt++;
|
|
}
|
|
for(int i=0; i<27; i++) mIsoWeight[i] /= totw;
|
|
|
|
LbmVec isostart = vec2L(mvGeoStart);
|
|
LbmVec isoend = vec2L(mvGeoEnd);
|
|
int twodOff = 0; // 2d slices
|
|
if(LBMDIM==2) {
|
|
LbmFloat sn,se;
|
|
sn = isostart[2]+(isoend[2]-isostart[2])*0.5 - ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
|
|
se = isostart[2]+(isoend[2]-isostart[2])*0.5 + ((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5;
|
|
isostart[2] = sn;
|
|
isoend[2] = se;
|
|
twodOff = 2;
|
|
}
|
|
int isosx = mSizex+2;
|
|
int isosy = mSizey+2;
|
|
int isosz = mSizez+2+twodOff;
|
|
|
|
// MPT
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
//if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
|
|
if( (mMpNum>0) && (mMpIndex==0) ) {
|
|
//? mpindex==0
|
|
// restore original value for node0
|
|
isosx = mOrgSizeX + 2;
|
|
isostart[0] = mOrgStartX;
|
|
isoend[0] = mOrgEndX;
|
|
}
|
|
errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMin.mSrcx,mpTest->mGCMin.mSrcy,mpTest->mGCMin.mSrcz)<<" dst"
|
|
<< PRINT_VEC(mpTest->mGCMin.mDstx,mpTest->mGCMin.mDsty,mpTest->mGCMin.mDstz)<<" consize"
|
|
<< PRINT_VEC(mpTest->mGCMin.mConSizex,mpTest->mGCMin.mConSizey,mpTest->mGCMin.mConSizez)<<" ");
|
|
errMsg("LbmFsgrSolver::initialize", "MPT: gcon "<<mMpNum<<","<<mMpIndex<<" src"<< PRINT_VEC(mpTest->mGCMax.mSrcx,mpTest->mGCMax.mSrcy,mpTest->mGCMax.mSrcz)<<" dst"
|
|
<< PRINT_VEC(mpTest->mGCMax.mDstx,mpTest->mGCMax.mDsty,mpTest->mGCMax.mDstz)<<" consize"
|
|
<< PRINT_VEC(mpTest->mGCMax.mConSizex,mpTest->mGCMax.mConSizey,mpTest->mGCMax.mConSizez)<<" ");
|
|
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
|
|
|
errMsg(" SETISO ", "iso "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(mSizex+1.0))*0.5)<<" "<<(LbmFloat)(mSizex+1.0) );
|
|
errMsg("LbmFsgrSolver::initialize", "MPT: geo "<< mvGeoStart<<","<<mvGeoEnd<<
|
|
" grid:"<<PRINT_VEC(mSizex,mSizey,mSizez)<<",iso:"<< PRINT_VEC(isosx,isosy,isosz) );
|
|
mpIso->setStart( vec2G(isostart) );
|
|
mpIso->setEnd( vec2G(isoend) );
|
|
LbmVec isodist = isoend-isostart;
|
|
|
|
int isosubs = mIsoSubdivs;
|
|
if(mFarFieldSize>1.) {
|
|
errMsg("LbmFsgrSolver::initialize","Warning - resetting isosubdivs, using fulledge!");
|
|
isosubs = 1;
|
|
mpIso->setUseFulledgeArrays(true);
|
|
}
|
|
mpIso->setSubdivs(isosubs);
|
|
|
|
mpIso->initializeIsosurface( isosx,isosy,isosz, vec2G(isodist) );
|
|
|
|
// reset iso field
|
|
for(int ak=0;ak<isosz;ak++)
|
|
for(int aj=0;aj<isosy;aj++)
|
|
for(int ai=0;ai<isosx;ai++) { *mpIso->getData(ai,aj,ak) = 0.0; }
|
|
|
|
|
|
/* init array (set all invalid first) */
|
|
preinitGrids();
|
|
for(int lev=0; lev<=mMaxRefine; lev++) {
|
|
FSGR_FORIJK_BOUNDS(lev) {
|
|
RFLAG(lev,i,j,k,0) = 0, RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
|
|
if(!mAllfluid) {
|
|
initEmptyCell(lev, i,j,k, CFEmpty, -1.0, -1.0);
|
|
} else {
|
|
initEmptyCell(lev, i,j,k, CFFluid, 1.0, 1.0);
|
|
}
|
|
}
|
|
}
|
|
|
|
|
|
if(LBMDIM==2) {
|
|
if(mOutputSurfacePreview) {
|
|
errMsg("LbmFsgrSolver::init","No preview in 2D allowed!");
|
|
mOutputSurfacePreview = 0; }
|
|
}
|
|
if((glob_mpactive) && (glob_mpindex>0)) {
|
|
mOutputSurfacePreview = 0;
|
|
}
|
|
|
|
#if LBM_USE_GUI==1
|
|
if(mOutputSurfacePreview) {
|
|
errMsg("LbmFsgrSolver::init","No preview in GUI mode... mOutputSurfacePreview=0");
|
|
mOutputSurfacePreview = 0; }
|
|
#endif // LBM_USE_GUI==1
|
|
if(mOutputSurfacePreview) {
|
|
// same as normal one, but use reduced size
|
|
mpPreviewSurface = new IsoSurface( mIsoValue );
|
|
mpPreviewSurface->setMaterialName( mpPreviewSurface->getMaterialName() );
|
|
mpPreviewSurface->setIsolevel( mIsoValue );
|
|
// usually dont display for rendering
|
|
mpPreviewSurface->setVisible( false );
|
|
|
|
mpPreviewSurface->setStart( vec2G(isostart) );
|
|
mpPreviewSurface->setEnd( vec2G(isoend) );
|
|
LbmVec pisodist = isoend-isostart;
|
|
LbmFloat pfac = mPreviewFactor;
|
|
mpPreviewSurface->initializeIsosurface( (int)(pfac*mSizex)+2, (int)(pfac*mSizey)+2, (int)(pfac*mSizez)+2, vec2G(pisodist) );
|
|
//mpPreviewSurface->setName( getName() + "preview" );
|
|
mpPreviewSurface->setName( "preview" );
|
|
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Preview with sizes "<<(pfac*mSizex)<<","<<(pfac*mSizey)<<","<<(pfac*mSizez)<<" enabled",10);
|
|
}
|
|
|
|
// init defaults
|
|
mAvgNumUsedCells = 0;
|
|
mFixMass= 0.0;
|
|
return true;
|
|
}
|
|
|
|
/*! init solver arrays */
|
|
bool LbmFsgrSolver::initializeSolverGrids() {
|
|
/* init boundaries */
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Boundary init...",10);
|
|
// init obstacles, and reinit time step size
|
|
initGeometryFlags();
|
|
mLastSimTime = -1.0;
|
|
// TODO check for invalid cells? nitGenericTestCases();
|
|
|
|
// new - init noslip 1 everywhere...
|
|
// half fill boundary cells?
|
|
|
|
CellFlagType domainBoundType = CFInvalid;
|
|
// TODO use normal object types instad...
|
|
if(mDomainBound.find(string("free")) != string::npos) {
|
|
domainBoundType = CFBnd | CFBndFreeslip;
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: FreeSlip, value:"<<mDomainBound,10);
|
|
} else if(mDomainBound.find(string("part")) != string::npos) {
|
|
domainBoundType = CFBnd | CFBndPartslip; // part slip type
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: PartSlip ("<<mDomainPartSlipValue<<"), value:"<<mDomainBound,10);
|
|
} else {
|
|
domainBoundType = CFBnd | CFBndNoslip;
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Domain Boundary Type: NoSlip, value:"<<mDomainBound,10);
|
|
}
|
|
|
|
// use ar[numobjs] as entry for domain (used e.g. for mDomainPartSlipValue in mObjectPartslips)
|
|
int domainobj = (int)(mpGiObjects->size());
|
|
domainBoundType |= (domainobj<<24);
|
|
//for(int i=0; i<(int)(domainobj+0); i++) {
|
|
//errMsg("GEOIN","i"<<i<<" "<<(*mpGiObjects)[i]->getName());
|
|
//if((*mpGiObjects)[i] == mpIso) { //check...
|
|
//}
|
|
//}
|
|
//errMsg("GEOIN"," dm "<<(domainBoundType>>24));
|
|
|
|
for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
|
|
for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
|
|
initEmptyCell(mMaxRefine, i,0,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 j=0;j<mLevel[mMaxRefine].lSizey;j++) {
|
|
initEmptyCell(mMaxRefine, 0,j,k, domainBoundType, 0.0, BND_FILL);
|
|
initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-1,j,k, domainBoundType, 0.0, BND_FILL);
|
|
// DEBUG BORDER!
|
|
//initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL);
|
|
}
|
|
|
|
if(LBMDIM == 3) {
|
|
// only for 3D
|
|
for(int j=0;j<mLevel[mMaxRefine].lSizey;j++)
|
|
for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
|
|
initEmptyCell(mMaxRefine, i,j,0, domainBoundType, 0.0, BND_FILL);
|
|
initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-1, domainBoundType, 0.0, BND_FILL);
|
|
}
|
|
}
|
|
|
|
// TEST!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!11
|
|
/*for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
|
|
for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
|
|
initEmptyCell(mMaxRefine, mLevel[mMaxRefine].lSizex-2,j,k, domainBoundType, 0.0, BND_FILL);
|
|
}
|
|
for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
|
|
for(int i=0;i<mLevel[mMaxRefine].lSizex;i++) {
|
|
initEmptyCell(mMaxRefine, i,1,k, domainBoundType, 0.0, BND_FILL);
|
|
}
|
|
// */
|
|
|
|
/*for(int ii=0; ii<(int)po w_change?(2.0,mMaxRefine)-1; ii++) {
|
|
errMsg("BNDTESTSYMM","set "<<mLevel[mMaxRefine].lSizex-2-ii );
|
|
for(int k=0;k<mLevel[mMaxRefine].lSizez;k++)
|
|
for(int j=0;j<mLevel[mMaxRefine].lSizey;j++) {
|
|
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 i=0;i<mLevel[mMaxRefine].lSizex;i++) {
|
|
initEmptyCell(mMaxRefine, i,j,mLevel[mMaxRefine].lSizez-2-ii, domainBoundType, 0.0, BND_FILL); // SYMM!? 3D?
|
|
}
|
|
}
|
|
// Symmetry tests */
|
|
// vortt
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
if(( strstr( this->getName().c_str(), "vorttfluid" ) != NULL) && (LBMDIM==2)) {
|
|
errMsg("VORTT","init");
|
|
int level=mMaxRefine;
|
|
int cx = mLevel[level].lSizex/2;
|
|
int cyo = mLevel[level].lSizey/2;
|
|
int sx = mLevel[level].lSizex/8;
|
|
int sy = mLevel[level].lSizey/8;
|
|
LbmFloat rho = 1.;
|
|
LbmFloat rhomass = 1.;
|
|
LbmFloat uFactor = 0.15;
|
|
LbmFloat vdist = 1.0;
|
|
|
|
int cy1=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=1;j<mLevel[level].lSizey-1;j++)
|
|
for(int i=1;i<mLevel[level].lSizex-1;i++) {
|
|
LbmFloat d1 = norm(LbmVec(cx,cy1,0.)-LbmVec(i,j,0));
|
|
LbmFloat d2 = norm(LbmVec(cx,cy2,0.)-LbmVec(i,j,0));
|
|
bool in1 = (d1<=(LbmFloat)(sx));
|
|
bool in2 = (d2<=(LbmFloat)(sx));
|
|
LbmVec uvec(0.);
|
|
LbmVec v1 = getNormalized( cross( LbmVec(cx,cy1,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )* uFactor;
|
|
LbmVec v2 = getNormalized( cross( LbmVec(cx,cy2,0.)-LbmVec(i,j,0), LbmVec(0.,0.,1.)) )* uFactor;
|
|
LbmFloat w1=1., w2=1.;
|
|
if(!in1) w1=(LbmFloat)(sx)/(1.5*d1);
|
|
if(!in2) w2=(LbmFloat)(sx)/(1.5*d2);
|
|
if(!in1) w1=0.; if(!in2) w2=0.; // sharp falloff
|
|
uvec += v1*w1;
|
|
uvec += v2*w2;
|
|
initVelocityCell(level, i,j,0, CFFluid, rho, rhomass, uvec );
|
|
//errMsg("VORTT","init uvec"<<uvec);
|
|
}
|
|
|
|
}
|
|
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
|
|
|
//if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal1, causing panic, aborting..."); return false; }
|
|
|
|
// prepare interface cells
|
|
initFreeSurfaces();
|
|
initStandingFluidGradient();
|
|
|
|
// perform first step to init initial mass
|
|
mInitialMass = 0.0;
|
|
int inmCellCnt = 0;
|
|
FSGR_FORIJK1(mMaxRefine) {
|
|
if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid) {
|
|
LbmFloat fluidRho = QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, 0);
|
|
FORDF1 { fluidRho += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, l); }
|
|
mInitialMass += fluidRho;
|
|
inmCellCnt ++;
|
|
} else if( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) {
|
|
mInitialMass += QCELL(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, dMass);
|
|
inmCellCnt ++;
|
|
}
|
|
}
|
|
mCurrentVolume = mCurrentMass = mInitialMass;
|
|
|
|
ParamVec cspv = mpParam->calculateCellSize();
|
|
if(LBMDIM==2) cspv[2] = 1.0;
|
|
inmCellCnt = 1;
|
|
double nrmMass = (double)mInitialMass / (double)(inmCellCnt) *cspv[0]*cspv[1]*cspv[2] * 1000.0;
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Initial Mass:"<<mInitialMass<<" normalized:"<<nrmMass, 3);
|
|
mInitialMass = 0.0; // reset, and use actual value after first step
|
|
|
|
//mStartSymm = false;
|
|
#if ELBEEM_PLUGIN!=1
|
|
if((LBMDIM==2)&&(mSizex<200)) {
|
|
if(!checkSymmetry("init")) {
|
|
errMsg("LbmFsgrSolver::initialize","Unsymmetric init...");
|
|
} else {
|
|
errMsg("LbmFsgrSolver::initialize","Symmetric init!");
|
|
}
|
|
}
|
|
#endif // ELBEEM_PLUGIN!=1
|
|
return true;
|
|
}
|
|
|
|
|
|
/*! prepare actual simulation start, setup viz etc */
|
|
bool LbmFsgrSolver::initializeSolverPostinit() {
|
|
// coarsen region
|
|
myTime_t fsgrtstart = getTime();
|
|
for(int lev=mMaxRefine-1; lev>=0; lev--) {
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Coarsening level "<<lev<<".",8);
|
|
adaptGrid(lev);
|
|
coarseRestrictFromFine(lev);
|
|
adaptGrid(lev);
|
|
coarseRestrictFromFine(lev);
|
|
}
|
|
markedClearList();
|
|
myTime_t fsgrtend = getTime();
|
|
if(!mSilent){ debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"FSGR init done ("<< getTimeString(fsgrtend-fsgrtstart)<<"), changes:"<<mNumFsgrChanges , 10 ); }
|
|
mNumFsgrChanges = 0;
|
|
|
|
for(int l=0; l<cDirNum; l++) {
|
|
LbmFloat area = 0.5 * 0.5 *0.5;
|
|
if(LBMDIM==2) area = 0.5 * 0.5;
|
|
|
|
if(dfVecX[l]!=0) area *= 0.5;
|
|
if(dfVecY[l]!=0) area *= 0.5;
|
|
if(dfVecZ[l]!=0) area *= 0.5;
|
|
mFsgrCellArea[l] = area;
|
|
} // l
|
|
|
|
// make sure both sets are ok
|
|
// copy from other to curr
|
|
for(int lev=0; lev<=mMaxRefine; lev++) {
|
|
FSGR_FORIJK_BOUNDS(lev) {
|
|
RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
|
|
} } // first copy flags */
|
|
|
|
|
|
// old mpPreviewSurface init
|
|
//if(getGlobalBakeState()<0) { CAUSE_PANIC; errMsg("LbmFsgrSolver::initialize","Got abort signal2, causing panic, aborting..."); return false; }
|
|
// make sure fill fracs are right for first surface generation
|
|
stepMain();
|
|
|
|
// prepare once...
|
|
mpIso->setParticles(mpParticles, mPartDropMassSub);
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Iso Settings, subdivs="<<mpIso->getSubdivs()<<", partsize="<<mPartDropMassSub, 9);
|
|
prepareVisualization();
|
|
// copy again for stats counting
|
|
for(int lev=0; lev<=mMaxRefine; lev++) {
|
|
FSGR_FORIJK_BOUNDS(lev) {
|
|
RFLAG(lev, i,j,k,mLevel[lev].setOther) = RFLAG(lev, i,j,k,mLevel[lev].setCurr);
|
|
} } // first copy flags */
|
|
|
|
|
|
// now really done...
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"SurfaceGen: SmsOrg("<<mSmoothSurface<<","<<mSmoothNormals<< /*","<<featureSize<<*/ "), Iso("<<mpIso->getSmoothSurface()<<","<<mpIso->getSmoothNormals()<<") ",10);
|
|
debMsgStd("LbmFsgrSolver::initialize",DM_MSG,"Init done ... ",10);
|
|
mInitDone = 1;
|
|
|
|
// init fluid control
|
|
initCpdata();
|
|
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
initTestdata();
|
|
#endif // ELBEEM_PLUGIN!=1
|
|
// not inited? dont use...
|
|
if(mCutoff<0) mCutoff=0;
|
|
|
|
initParticles();
|
|
return true;
|
|
}
|
|
|
|
|
|
|
|
// macros for mov obj init
|
|
#if LBMDIM==2
|
|
|
|
#define POS2GRID_CHECK(vec,n) \
|
|
monTotal++;\
|
|
int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
|
|
if(k!=0) continue; \
|
|
const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
|
|
if(i<=0) continue; \
|
|
if(i>=mLevel[level].lSizex-1) continue; \
|
|
const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
|
|
if(j<=0) continue; \
|
|
if(j>=mLevel[level].lSizey-1) continue; \
|
|
|
|
#else // LBMDIM -> 3
|
|
#define POS2GRID_CHECK(vec,n) \
|
|
monTotal++;\
|
|
const int i=(int)( ((vec)[n][0]-iniPos[0])/dvec[0] +0.0); \
|
|
if(i<=0) continue; \
|
|
if(i>=mLevel[level].lSizex-1) continue; \
|
|
const int j=(int)( ((vec)[n][1]-iniPos[1])/dvec[1] +0.0); \
|
|
if(j<=0) continue; \
|
|
if(j>=mLevel[level].lSizey-1) continue; \
|
|
const int k=(int)( ((vec)[n][2]-iniPos[2])/dvec[2] +0.0); \
|
|
if(k<=0) continue; \
|
|
if(k>=mLevel[level].lSizez-1) continue; \
|
|
|
|
#endif // LBMDIM
|
|
|
|
// calculate object velocity from vert arrays in objvel vec
|
|
#define OBJVEL_CALC \
|
|
LbmVec objvel = vec2L((mMOIVertices[n]-mMOIVerticesOld[n]) /dvec); { \
|
|
const LbmFloat usqr = (objvel[0]*objvel[0]+objvel[1]*objvel[1]+objvel[2]*objvel[2])*1.5; \
|
|
USQRMAXCHECK(usqr, objvel[0],objvel[1],objvel[2], mMaxVlen, mMxvx,mMxvy,mMxvz); \
|
|
if(usqr>maxusqr) { \
|
|
/* cutoff at maxVelVal */ \
|
|
for(int jj=0; jj<3; jj++) { \
|
|
if(objvel[jj]>0.) objvel[jj] = maxVelVal; \
|
|
if(objvel[jj]<0.) objvel[jj] = -maxVelVal; \
|
|
} \
|
|
} } \
|
|
if(ntype&(CFBndFreeslip)) { \
|
|
const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
|
|
const LbmVec oldov=objvel; /*DEBUG*/ \
|
|
objvel = vec2L((*pNormals)[n]) *dp; \
|
|
/* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
|
|
} \
|
|
else if(ntype&(CFBndPartslip)) { \
|
|
const LbmFloat dp=dot(objvel, vec2L((*pNormals)[n]) ); \
|
|
const LbmVec oldov=objvel; /*DEBUG*/ \
|
|
/* if((j==24)&&(n%5==2)) errMsg("FSBT","n"<<n<<" v"<<objvel<<" nn"<<(*pNormals)[n]<<" dp"<<dp<<" oldov"<<oldov ); */ \
|
|
const LbmFloat partv = mObjectPartslips[OId]; \
|
|
/*errMsg("PARTSLIP_DEBUG","l="<<l<<" ccel="<<RAC(ccel, dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
|
|
/* m[l] = (RAC(ccel, dfInv[l] ) ) * partv + newval * (1.-partv); part slip */ \
|
|
objvel = objvel*partv + vec2L((*pNormals)[n]) *dp*(1.-partv); \
|
|
}
|
|
|
|
#define TTT \
|
|
|
|
|
|
/*****************************************************************************/
|
|
//! init moving obstacles for next sim step sim
|
|
/*****************************************************************************/
|
|
void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
|
myTime_t monstart = getTime();
|
|
|
|
// movobj init
|
|
const int level = mMaxRefine;
|
|
const int workSet = mLevel[level].setCurr;
|
|
const int otherSet = mLevel[level].setOther;
|
|
LbmFloat sourceTime = mSimulationTime; // should be equal to mLastSimTime!
|
|
// for debugging - check targetTime check during DEFAULT STREAM
|
|
LbmFloat targetTime = mSimulationTime + mpParam->getTimestep();
|
|
if(mLastSimTime == targetTime) {
|
|
debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"Called for same time! (t="<<mSimulationTime<<" , targett="<<targetTime<<")", 1);
|
|
return;
|
|
}
|
|
//debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_WARNING,"time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime,10);
|
|
//if(mSimulationTime!=mLastSimTime) errMsg("LbmFsgrSolver::initMovingObstacles","time: "<<mSimulationTime<<" lasttt:"<<mLastSimTime);
|
|
|
|
const LbmFloat maxVelVal = 0.1666;
|
|
const LbmFloat maxusqr = maxVelVal*maxVelVal*3. *1.5;
|
|
|
|
LbmFloat rhomass = 0.0;
|
|
CellFlagType otype = CFInvalid; // verify type of last step, might be non moving obs!
|
|
CellFlagType ntype = CFInvalid;
|
|
// WARNING - copied from geo init!
|
|
int numobjs = (int)(mpGiObjects->size());
|
|
ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
|
|
// 2d display as rectangles
|
|
ntlVec3Gfx iniPos(0.0);
|
|
if(LBMDIM==2) {
|
|
dvec[2] = 1.0;
|
|
iniPos = (mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (mvGeoEnd[2]-mvGeoStart[2])*0.5 ))-(dvec*0.0);
|
|
} else {
|
|
iniPos = (mvGeoStart + ntlVec3Gfx( 0.0 ))-(dvec*0.0);
|
|
}
|
|
|
|
if( (int)mObjectMassMovnd.size() < numobjs) {
|
|
for(int i=mObjectMassMovnd.size(); i<numobjs; i++) {
|
|
mObjectMassMovnd.push_back(0.);
|
|
}
|
|
}
|
|
|
|
// stats
|
|
int monPoints=0, monObsts=0, monFluids=0, monTotal=0, monTrafo=0;
|
|
int nbored;
|
|
for(int OId=0; OId<numobjs; OId++) {
|
|
ntlGeometryObject *obj = (*mpGiObjects)[OId];
|
|
bool skip = false;
|
|
if(obj->getGeoInitId() != mLbmInitId) skip=true;
|
|
if( (!staticInit) && (!obj->getIsAnimated()) ) skip=true;
|
|
if( ( staticInit) && ( obj->getIsAnimated()) ) skip=true;
|
|
if(skip) continue;
|
|
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) ||
|
|
((obj->getGeoInitType()&FGI_FLUID) && staticInit) ) {
|
|
|
|
otype = ntype = CFInvalid;
|
|
switch(obj->getGeoInitType()) {
|
|
/* case FGI_BNDPART: // old, use noslip for moving part/free objs
|
|
case FGI_BNDFREE:
|
|
if(!staticInit) {
|
|
errMsg("LbmFsgrSolver::initMovingObstacles","Warning - moving free/part slip objects NYI "<<obj->getName() );
|
|
otype = ntype = CFBnd|CFBndNoslip;
|
|
} else {
|
|
if(obj->getGeoInitType()==FGI_BNDPART) otype = ntype = CFBnd|CFBndPartslip;
|
|
if(obj->getGeoInitType()==FGI_BNDFREE) otype = ntype = CFBnd|CFBndFreeslip;
|
|
}
|
|
break;
|
|
// off */
|
|
case FGI_BNDPART: rhomass = BND_FILL;
|
|
otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
|
|
break;
|
|
case FGI_BNDFREE: rhomass = BND_FILL;
|
|
otype = ntype = CFBnd|CFBndFreeslip|(OId<<24);
|
|
break;
|
|
// off */
|
|
case FGI_BNDNO: rhomass = BND_FILL;
|
|
otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
|
|
break;
|
|
case FGI_FLUID:
|
|
otype = ntype = CFFluid;
|
|
break;
|
|
case FGI_MBNDINFLOW:
|
|
otype = ntype = CFMbndInflow;
|
|
break;
|
|
case FGI_MBNDOUTFLOW:
|
|
otype = ntype = CFMbndOutflow;
|
|
break;
|
|
}
|
|
int wasActive = ((obj->getGeoActive(sourceTime)>0.)? 1:0);
|
|
int active = ((obj->getGeoActive(targetTime)>0.)? 1:0);
|
|
//errMsg("GEOACTT"," obj "<<obj->getName()<<" a:"<<active<<","<<wasActive<<" s"<<sourceTime<<" t"<<targetTime <<" v"<<mObjectSpeeds[OId] );
|
|
// skip inactive in/out flows
|
|
if(ntype==CFInvalid){ errMsg("LbmFsgrSolver::initMovingObstacles","Invalid obj type "<<obj->getGeoInitType()); continue; }
|
|
if((!active) && (otype&(CFMbndOutflow|CFMbndInflow)) ) continue;
|
|
|
|
// copied from recalculateObjectSpeeds
|
|
mObjectSpeeds[OId] = vec2L(mpParam->calculateLattVelocityFromRw( vec2P( (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
|
|
debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) ,10 );
|
|
|
|
//vector<ntlVec3Gfx> tNormals;
|
|
vector<ntlVec3Gfx> *pNormals = NULL;
|
|
mMOINormals.clear();
|
|
if(ntype&(CFBndFreeslip|CFBndPartslip)) { pNormals = &mMOINormals; }
|
|
|
|
mMOIVertices.clear();
|
|
if(obj->getMeshAnimated()) {
|
|
// do two full update
|
|
// TODO tNormals handling!?
|
|
mMOIVerticesOld.clear();
|
|
obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals, mLevel[mMaxRefine].nodeSize, mvGeoStart, mvGeoEnd);
|
|
monTrafo += mMOIVerticesOld.size();
|
|
obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
|
|
monTrafo += mMOIVertices.size();
|
|
obj->applyTransformation(targetTime, &mMOIVertices,NULL /* no old normals needed */, 0, mMOIVertices.size(), false );
|
|
} else {
|
|
// only do transform update
|
|
obj->getMovingPoints(mMOIVertices,pNormals); // mMOIVertices = mCachedMovPoints
|
|
mMOIVerticesOld = mMOIVertices;
|
|
// WARNING - assumes mSimulationTime is global!?
|
|
obj->applyTransformation(targetTime, &mMOIVertices,pNormals, 0, mMOIVertices.size(), false );
|
|
monTrafo += mMOIVertices.size();
|
|
|
|
// correct flags from last position, but extrapolate
|
|
// velocity to next timestep
|
|
obj->applyTransformation(sourceTime, &mMOIVerticesOld, NULL /* no old normals needed */, 0, mMOIVerticesOld.size(), false );
|
|
monTrafo += mMOIVerticesOld.size();
|
|
}
|
|
|
|
// object types
|
|
if(ntype&CFBnd){
|
|
|
|
// check if object is moving at all
|
|
if(obj->getIsAnimated()) {
|
|
ntlVec3Gfx objMaxVel = obj->calculateMaxVel(sourceTime,targetTime);
|
|
// FIXME?
|
|
if(normNoSqrt(objMaxVel)>0.0) { ntype |= CFBndMoving; }
|
|
// get old type - CHECK FIXME , timestep could have changed - cause trouble?
|
|
ntlVec3Gfx oldobjMaxVel = obj->calculateMaxVel(sourceTime - mpParam->getTimestep(),sourceTime);
|
|
if(normNoSqrt(oldobjMaxVel)>0.0) { otype |= CFBndMoving; }
|
|
}
|
|
if(obj->getMeshAnimated()) { ntype |= CFBndMoving; otype |= CFBndMoving; }
|
|
CellFlagType rflagnb[27];
|
|
LbmFloat massCheck = 0.;
|
|
int massReinits=0;
|
|
bool fillCells = (mObjectMassMovnd[OId]<=-1.);
|
|
LbmFloat impactCorrFactor = obj->getGeoImpactFactor(targetTime);
|
|
|
|
|
|
// first pass - set new obs. cells
|
|
if(active) {
|
|
for(size_t n=0; n<mMOIVertices.size(); n++) {
|
|
//errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
|
|
POS2GRID_CHECK(mMOIVertices,n);
|
|
//{ errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<", t="<<targetTime); }
|
|
if(QCELL(level, i,j,k, workSet, dFlux)==targetTime) continue;
|
|
monPoints++;
|
|
|
|
// check mass
|
|
if(RFLAG(level, i,j,k, workSet)&(CFFluid)) {
|
|
FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); }
|
|
massReinits++;
|
|
}
|
|
else if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
|
|
massCheck -= QCELL(level, i,j,k, workSet, dMass);
|
|
massReinits++;
|
|
}
|
|
|
|
RFLAG(level, i,j,k, workSet) = ntype;
|
|
FORDF1 {
|
|
//CellFlagType flag = RFLAG_NB(level, i,j,k,workSet,l);
|
|
rflagnb[l] = RFLAG_NB(level, i,j,k,workSet,l);
|
|
if(rflagnb[l]&(CFFluid|CFInter)) {
|
|
rflagnb[l] &= (~CFNoBndFluid); // remove CFNoBndFluid flag
|
|
RFLAG_NB(level, i,j,k,workSet,l) &= rflagnb[l];
|
|
}
|
|
}
|
|
LbmFloat *dstCell = RACPNT(level, i,j,k,workSet);
|
|
RAC(dstCell,0) = 0.0;
|
|
if(ntype&CFBndMoving) {
|
|
OBJVEL_CALC;
|
|
|
|
// compute fluid acceleration
|
|
FORDF1 {
|
|
if(rflagnb[l]&(CFFluid|CFInter)) {
|
|
const LbmFloat ux = dfDvecX[l]*objvel[0];
|
|
const LbmFloat uy = dfDvecY[l]*objvel[1];
|
|
const LbmFloat uz = dfDvecZ[l]*objvel[2];
|
|
|
|
LbmFloat factor = 2. * dfLength[l] * 3.0 * (ux+uy+uz); //
|
|
if(ntype&(CFBndFreeslip|CFBndPartslip)) {
|
|
// missing, diag mass checks...
|
|
//if(l<=LBMDIM*2) factor *= 1.4142;
|
|
factor *= 2.0; // TODO, optimize
|
|
} else {
|
|
factor *= 1.2; // TODO, optimize
|
|
}
|
|
factor *= impactCorrFactor; // use manual tweaking channel
|
|
RAC(dstCell,l) = factor;
|
|
massCheck += factor;
|
|
} else {
|
|
RAC(dstCell,l) = 0.;
|
|
}
|
|
}
|
|
|
|
#if NEWDIRVELMOTEST==1
|
|
FORDF1 { RAC(dstCell,l) = 0.; }
|
|
RAC(dstCell,dMass) = objvel[0];
|
|
RAC(dstCell,dFfrac) = objvel[1];
|
|
RAC(dstCell,dC) = objvel[2];
|
|
#endif // NEWDIRVELMOTEST==1
|
|
} else {
|
|
FORDF1 { RAC(dstCell,l) = 0.0; }
|
|
}
|
|
RAC(dstCell, dFlux) = targetTime;
|
|
//errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK" dflt="<<RAC(dstCell, dFlux) );
|
|
monObsts++;
|
|
} // points
|
|
} // bnd, is active?
|
|
|
|
// second pass, remove old ones
|
|
// warning - initEmptyCell et. al dont overwrite OId or persist flags...
|
|
if(wasActive) {
|
|
for(size_t n=0; n<mMOIVerticesOld.size(); n++) {
|
|
POS2GRID_CHECK(mMOIVerticesOld,n);
|
|
monPoints++;
|
|
if((RFLAG(level, i,j,k, workSet) == otype) &&
|
|
(QCELL(level, i,j,k, workSet, dFlux) != targetTime)) {
|
|
// from mainloop
|
|
nbored = 0;
|
|
// TODO: optimize for OPT3D==0
|
|
FORDF1 {
|
|
//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
|
|
nbored |= rflagnb[l];
|
|
}
|
|
CellFlagType settype = CFInvalid;
|
|
if(nbored&CFFluid) {
|
|
settype = CFInter|CFNoInterpolSrc;
|
|
rhomass = 1.5;
|
|
if(!fillCells) rhomass = 0.;
|
|
|
|
OBJVEL_CALC;
|
|
if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }
|
|
|
|
// new interpolate values
|
|
LbmFloat avgrho = 0.0;
|
|
LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
|
|
interpolateCellValues(level,i,j,k,workSet, avgrho,avgux,avguy,avguz);
|
|
initVelocityCell(level, i,j,k, settype, avgrho, rhomass, LbmVec(avgux,avguy,avguz) );
|
|
//errMsg("NMOCIT"," at "<<PRINT_IJK<<" "<<avgrho<<" "<<norm(LbmVec(avgux,avguy,avguz))<<" "<<LbmVec(avgux,avguy,avguz) );
|
|
massCheck += rhomass;
|
|
} else {
|
|
settype = CFEmpty; rhomass = 0.0;
|
|
initEmptyCell(level, i,j,k, settype, 1., rhomass );
|
|
}
|
|
monFluids++;
|
|
massReinits++;
|
|
} // flag & simtime
|
|
}
|
|
} // wasactive
|
|
|
|
// only compute mass transfer when init is done
|
|
if(mInitDone) {
|
|
errMsg("initMov","dccd\n\nMassd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
|
|
" fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" inim:"<<mInitialMass );
|
|
mObjectMassMovnd[OId] += massCheck;
|
|
}
|
|
} // bnd, active
|
|
|
|
else if(ntype&CFFluid){
|
|
// second static init pass
|
|
if(staticInit) {
|
|
//debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" verts"<<mMOIVertices.size() ,9);
|
|
CellFlagType setflflag = CFFluid; //|(OId<<24);
|
|
for(size_t n=0; n<mMOIVertices.size(); n++) {
|
|
POS2GRID_CHECK(mMOIVertices,n);
|
|
if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
|
|
if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
|
|
//changeFlag(level, i,j,k, workSet, setflflag);
|
|
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); }
|
|
}
|
|
} // second static inflow pass
|
|
} // fluid
|
|
|
|
else if(ntype&CFMbndInflow){
|
|
// inflow pass - add new fluid cells
|
|
// this is slightly different for standing inflows,
|
|
// as the fluid is forced to the desired velocity inside as
|
|
// well...
|
|
const LbmFloat iniRho = 1.0;
|
|
const LbmVec vel(mObjectSpeeds[OId]);
|
|
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);
|
|
//errMsg("LbmFsgrSolver::initMovingObstacles","id"<<OId<<" "<<obj->getName()<<" inflow "<<staticInit<<" "<<mMOIVertices.size() );
|
|
|
|
for(size_t n=0; n<mMOIVertices.size(); n++) {
|
|
POS2GRID_CHECK(mMOIVertices,n);
|
|
// TODO - also reinit interface cells !?
|
|
if(RFLAG(level, i,j,k, workSet)!=CFEmpty) {
|
|
// 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
|
|
continue; }
|
|
monFluids++;
|
|
|
|
// TODO add OPT3D treatment
|
|
LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
|
|
FORDF0 { RAC(tcel, l) = this->getCollideEq(l, iniRho,vel[0],vel[1],vel[2]); }
|
|
RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
|
|
RAC(tcel, dFlux) = FLUX_INIT;
|
|
CellFlagType setFlag = CFInter;
|
|
changeFlag(level, i,j,k, workSet, setFlag);
|
|
mInitialMass += iniRho;
|
|
}
|
|
// second static init pass
|
|
if(staticInit) {
|
|
CellFlagType set2Flag = CFMbndInflow|(OId<<24);
|
|
for(size_t n=0; n<mMOIVertices.size(); n++) {
|
|
POS2GRID_CHECK(mMOIVertices,n);
|
|
if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
|
|
if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
|
|
forceChangeFlag(level, i,j,k, workSet, set2Flag);
|
|
} else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
|
|
forceChangeFlag(level, i,j,k, workSet,
|
|
(RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
|
|
}
|
|
}
|
|
} // second static inflow pass
|
|
|
|
} // inflow
|
|
|
|
else if(ntype&CFMbndOutflow){
|
|
const LbmFloat iniRho = 0.0;
|
|
for(size_t n=0; n<mMOIVertices.size(); n++) {
|
|
POS2GRID_CHECK(mMOIVertices,n);
|
|
// FIXME check fluid/inter cells for non-static!?
|
|
if(!(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter))) {
|
|
if((staticInit)&&(RFLAG(level, i,j,k, workSet)==CFEmpty)) {
|
|
forceChangeFlag(level, i,j,k, workSet, CFMbndOutflow); }
|
|
continue;
|
|
}
|
|
monFluids++;
|
|
// interface cell - might be double empty? FIXME check?
|
|
|
|
// remove fluid cells, shouldnt be here anyway
|
|
LbmFloat *tcel = RACPNT(level, i,j,k,workSet);
|
|
LbmFloat fluidRho = RAC(tcel,0); FORDF1 { fluidRho += RAC(tcel,l); }
|
|
mInitialMass -= fluidRho;
|
|
RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
|
|
RAC(tcel, dFlux) = FLUX_INIT;
|
|
CellFlagType setFlag = CFInter;
|
|
changeFlag(level, i,j,k, workSet, setFlag);
|
|
|
|
// same as ifemptied for if below
|
|
LbmPoint emptyp;
|
|
emptyp.x = i; emptyp.y = j; emptyp.z = k;
|
|
mListEmpty.push_back( emptyp );
|
|
//calcCellsEmptied++;
|
|
} // points
|
|
// second static init pass
|
|
if(staticInit) {
|
|
CellFlagType set2Flag = CFMbndOutflow|(OId<<24);
|
|
for(size_t n=0; n<mMOIVertices.size(); n++) {
|
|
POS2GRID_CHECK(mMOIVertices,n);
|
|
if(RFLAG(level, i,j,k, workSet)&(CFMbndInflow|CFMbndOutflow)){ continue; }
|
|
if(RFLAG(level, i,j,k, workSet)&(CFEmpty)) {
|
|
forceChangeFlag(level, i,j,k, workSet, set2Flag);
|
|
} else if(RFLAG(level, i,j,k, workSet)&(CFFluid|CFInter)) {
|
|
forceChangeFlag(level, i,j,k, workSet,
|
|
(RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
|
|
}
|
|
}
|
|
} // second static outflow pass
|
|
} // outflow
|
|
|
|
} // allbound check
|
|
} // OId
|
|
|
|
|
|
/* { // DEBUG check
|
|
int workSet = mLevel[level].setCurr;
|
|
FSGR_FORIJK1(level) {
|
|
if( (RFLAG(level,i,j,k, workSet) & CFBndMoving) ) {
|
|
if(QCELL(level, i,j,k, workSet, dFlux)!=targetTime) {
|
|
errMsg("lastt"," old val!? at "<<PRINT_IJK<<" t="<<QCELL(level, i,j,k, workSet, dFlux)<<" target="<<targetTime);
|
|
}
|
|
}
|
|
}
|
|
} // DEBUG */
|
|
|
|
#undef POS2GRID_CHECK
|
|
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);
|
|
mLastSimTime = targetTime;
|
|
}
|
|
|
|
|
|
// geoinit position
|
|
|
|
#define GETPOS(i,j,k) \
|
|
ntlVec3Gfx ggpos = \
|
|
ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
|
|
iniPos[1]+ dvec[1]*(gfxReal)(j), \
|
|
iniPos[2]+ dvec[2]*(gfxReal)(k) ); \
|
|
if((LBMDIM==2)&&(mInit2dYZ)) { SWAPYZ(ggpos); }
|
|
|
|
/*****************************************************************************/
|
|
/*! perform geometry init (if switched on) */
|
|
/*****************************************************************************/
|
|
extern int globGeoInitDebug; //solver_interface
|
|
bool LbmFsgrSolver::initGeometryFlags() {
|
|
int level = mMaxRefine;
|
|
myTime_t geotimestart = getTime();
|
|
ntlGeometryObject *pObj;
|
|
ntlVec3Gfx dvec = ntlVec3Gfx(mLevel[level].nodeSize); //dvec*1.0;
|
|
debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing geometry init ("<< mLbmInitId <<") v"<<dvec,3);
|
|
// WARNING - copied to movobj init!
|
|
|
|
/* init object velocities, this has always to be called for init */
|
|
initGeoTree();
|
|
if(mAllfluid) {
|
|
freeGeoTree();
|
|
return true; }
|
|
|
|
// make sure moving obstacles are inited correctly
|
|
// preinit moving obj points
|
|
int numobjs = (int)(mpGiObjects->size());
|
|
for(int o=0; o<numobjs; o++) {
|
|
ntlGeometryObject *obj = (*mpGiObjects)[o];
|
|
//debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
|
|
if(
|
|
((obj->getGeoInitType()&FGI_ALLBOUNDS) && (obj->getIsAnimated())) ||
|
|
(obj->getVolumeInit()&VOLUMEINIT_SHELL) ) {
|
|
if(!obj->getMeshAnimated()) {
|
|
debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG," obj "<<obj->getName()<<" type "<<obj->getGeoInitType()<<" anim"<<obj->getIsAnimated()<<" "<<obj->getVolumeInit() ,9);
|
|
obj->initMovingPoints(mSimulationTime, mLevel[mMaxRefine].nodeSize);
|
|
}
|
|
}
|
|
}
|
|
|
|
// max speed init
|
|
ntlVec3Gfx maxMovobjVelRw = getGeoMaxMovementVelocity( mSimulationTime, mpParam->getTimestep() );
|
|
ntlVec3Gfx maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P( maxMovobjVelRw )) );
|
|
mpParam->setSimulationMaxSpeed( norm(maxMovobjVel) + norm(mLevel[level].gravity) );
|
|
LbmFloat allowMax = mpParam->getTadapMaxSpeed(); // maximum allowed velocity
|
|
debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Maximum Velocity from geo init="<< maxMovobjVel <<" from mov. obj.="<<maxMovobjVelRw<<" , allowed Max="<<allowMax ,5);
|
|
if(mpParam->getSimulationMaxSpeed() > allowMax) {
|
|
// similar to adaptTimestep();
|
|
LbmFloat nextmax = mpParam->getSimulationMaxSpeed();
|
|
LbmFloat newdt = mpParam->getTimestep() * (allowMax / nextmax); // newtr
|
|
debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"Performing reparametrization, newdt="<< newdt<<" prevdt="<< mpParam->getTimestep() <<" ",5);
|
|
mpParam->setDesiredTimestep( newdt );
|
|
mpParam->calculateAllMissingValues( mSimulationTime, mSilent );
|
|
maxMovobjVel = vec2G( mpParam->calculateLattVelocityFromRw( vec2P(getGeoMaxMovementVelocity(
|
|
mSimulationTime, mpParam->getTimestep() )) ));
|
|
debMsgStd("LbmFsgrSolver::initGeometryFlags",DM_MSG,"New maximum Velocity from geo init="<< maxMovobjVel,5);
|
|
}
|
|
recalculateObjectSpeeds();
|
|
// */
|
|
|
|
// init obstacles for first time step (requires obj speeds)
|
|
initMovingObstacles(true);
|
|
|
|
/* set interface cells */
|
|
ntlVec3Gfx pos,iniPos; // position of current cell
|
|
LbmFloat rhomass = 0.0;
|
|
CellFlagType ntype = CFInvalid;
|
|
int savedNodes = 0;
|
|
int OId = -1;
|
|
gfxReal distance;
|
|
|
|
// 2d display as rectangles
|
|
if(LBMDIM==2) {
|
|
dvec[2] = 0.0;
|
|
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 ); } }
|
|
} else {
|
|
iniPos =(mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
|
|
}
|
|
|
|
|
|
// first init boundary conditions
|
|
// invalid cells are set to empty afterwards
|
|
// TODO use floop macros!?
|
|
for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
|
|
for(int j=1;j<mLevel[level].lSizey-1;j++) {
|
|
for(int i=1;i<mLevel[level].lSizex-1;i++) {
|
|
ntype = CFInvalid;
|
|
|
|
GETPOS(i,j,k);
|
|
const bool inside = geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
|
|
if(inside) {
|
|
pObj = (*mpGiObjects)[OId];
|
|
switch( pObj->getGeoInitType() ){
|
|
case FGI_MBNDINFLOW:
|
|
if(! pObj->getIsAnimated() ) {
|
|
rhomass = 1.0;
|
|
ntype = CFFluid | CFMbndInflow;
|
|
} else {
|
|
ntype = CFInvalid;
|
|
}
|
|
break;
|
|
case FGI_MBNDOUTFLOW:
|
|
if(! pObj->getIsAnimated() ) {
|
|
rhomass = 0.0;
|
|
ntype = CFEmpty|CFMbndOutflow;
|
|
} else {
|
|
ntype = CFInvalid;
|
|
}
|
|
break;
|
|
case FGI_BNDNO:
|
|
rhomass = BND_FILL;
|
|
ntype = CFBnd|CFBndNoslip;
|
|
break;
|
|
case FGI_BNDPART:
|
|
rhomass = BND_FILL;
|
|
ntype = CFBnd|CFBndPartslip; break;
|
|
case FGI_BNDFREE:
|
|
rhomass = BND_FILL;
|
|
ntype = CFBnd|CFBndFreeslip; break;
|
|
default: // warn here?
|
|
rhomass = BND_FILL;
|
|
ntype = CFBnd|CFBndNoslip; break;
|
|
}
|
|
}
|
|
if(ntype != CFInvalid) {
|
|
// initDefaultCell
|
|
if((ntype & CFMbndInflow) || (ntype & CFMbndOutflow) ) { }
|
|
ntype |= (OId<<24); // NEWTEST2
|
|
initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
|
|
}
|
|
|
|
// walk along x until hit for following inits
|
|
if(distance<=-1.0) { distance = 100.0; } // FIXME dangerous
|
|
if(distance>=0.0) {
|
|
gfxReal dcnt=dvec[0];
|
|
while(( dcnt< distance-dvec[0] )&&(i+1<mLevel[level].lSizex-1)) {
|
|
dcnt += dvec[0]; i++;
|
|
savedNodes++;
|
|
if(ntype != CFInvalid) {
|
|
// rho,mass,OId are still inited from above
|
|
initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
|
|
}
|
|
}
|
|
}
|
|
// */
|
|
|
|
}
|
|
}
|
|
} // zmax
|
|
// */
|
|
|
|
// now init fluid layer
|
|
for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
|
|
for(int j=1;j<mLevel[level].lSizey-1;j++) {
|
|
for(int i=1;i<mLevel[level].lSizex-1;i++) {
|
|
if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
|
|
ntype = CFInvalid;
|
|
int inits = 0;
|
|
GETPOS(i,j,k);
|
|
const bool inside = geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
|
|
if(inside) {
|
|
ntype = CFFluid;
|
|
}
|
|
if(ntype != CFInvalid) {
|
|
// initDefaultCell
|
|
rhomass = 1.0;
|
|
initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
|
|
inits++;
|
|
}
|
|
|
|
// walk along x until hit for following inits
|
|
if(distance<=-1.0) { distance = 100.0; }
|
|
if(distance>=0.0) {
|
|
gfxReal dcnt=dvec[0];
|
|
while((dcnt< distance )&&(i+1<mLevel[level].lSizex-1)) {
|
|
dcnt += dvec[0]; i++;
|
|
savedNodes++;
|
|
if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
|
|
if(ntype != CFInvalid) {
|
|
// rhomass are still inited from above
|
|
initVelocityCell(level, i,j,k, ntype, rhomass, rhomass, mObjectSpeeds[OId] );
|
|
inits++;
|
|
}
|
|
}
|
|
} // distance>0
|
|
|
|
}
|
|
}
|
|
} // zmax
|
|
|
|
// reset invalid to empty again
|
|
for(int k= getForZMin1(); k< getForZMax1(level); ++k) {
|
|
for(int j=1;j<mLevel[level].lSizey-1;j++) {
|
|
for(int i=1;i<mLevel[level].lSizex-1;i++) {
|
|
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].setCurr) = CFEmpty;
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
freeGeoTree();
|
|
myTime_t geotimeend = getTime();
|
|
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));
|
|
return true;
|
|
}
|
|
|
|
#undef GETPOS
|
|
|
|
/*****************************************************************************/
|
|
/* init part for all freesurface testcases */
|
|
void LbmFsgrSolver::initFreeSurfaces() {
|
|
double interfaceFill = 0.45; // filling level of interface cells
|
|
//interfaceFill = 1.0; // DEUG!! GEOMTEST!!!!!!!!!!!!
|
|
|
|
// set interface cells
|
|
FSGR_FORIJK1(mMaxRefine) {
|
|
if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFFluid )) {
|
|
FORDF1 {
|
|
int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
|
|
if( ( RFLAG(mMaxRefine, ni, nj, nk, mLevel[mMaxRefine].setCurr)& CFEmpty ) ) {
|
|
LbmFloat arho=0., aux=0., auy=0., auz=0.;
|
|
interpolateCellValues(mMaxRefine, ni,nj,nk, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
|
|
//errMsg("TINI"," "<<PRINT_VEC(ni,nj,nk)<<" v"<<LbmVec(aux,auy,auz) );
|
|
// unnecessary? initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
|
|
initVelocityCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill, LbmVec(aux,auy,auz) );
|
|
}
|
|
}
|
|
}
|
|
}
|
|
|
|
// remove invalid interface cells
|
|
FSGR_FORIJK1(mMaxRefine) {
|
|
// remove invalid interface cells
|
|
if( ( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr)& CFInter) ) {
|
|
int delit = 0;
|
|
int NBs = 0; // neighbor flags or'ed
|
|
int noEmptyNB = 1;
|
|
|
|
FORDF1 {
|
|
if( ( RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr,l )& CFEmpty ) ) {
|
|
noEmptyNB = 0;
|
|
}
|
|
NBs |= RFLAG_NBINV(mMaxRefine, i, j, k, mLevel[mMaxRefine].setCurr, l);
|
|
}
|
|
// remove cells with no fluid or interface neighbors
|
|
if((NBs & CFFluid)==0) { delit = 1; }
|
|
if((NBs & CFInter)==0) { delit = 1; }
|
|
// remove cells with no empty neighbors
|
|
if(noEmptyNB) { delit = 2; }
|
|
|
|
// now we can remove the cell
|
|
if(delit==1) {
|
|
initEmptyCell(mMaxRefine, i,j,k, CFEmpty, 1.0, 0.0);
|
|
}
|
|
if(delit==2) {
|
|
//initEmptyCell(mMaxRefine, i,j,k, CFFluid, 1.0, 1.0);
|
|
LbmFloat arho=0., aux=0., auy=0., auz=0.;
|
|
interpolateCellValues(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr, arho,aux,auy,auz);
|
|
initVelocityCell(mMaxRefine, i,j,k, CFFluid, arho,1., LbmVec(aux,auy,auz) );
|
|
}
|
|
} // interface
|
|
} // */
|
|
|
|
// another brute force init, make sure the fill values are right...
|
|
// and make sure both sets are equal
|
|
for(int lev=0; lev<=mMaxRefine; lev++) {
|
|
FSGR_FORIJK_BOUNDS(lev) {
|
|
if( (RFLAG(lev, i,j,k,0) & (CFBnd)) ) {
|
|
QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = BND_FILL;
|
|
continue;
|
|
}
|
|
if( (RFLAG(lev, i,j,k,0) & (CFEmpty)) ) {
|
|
QCELL(lev, i,j,k,mLevel[mMaxRefine].setCurr, dFfrac) = 0.0;
|
|
continue;
|
|
}
|
|
} }
|
|
|
|
// ----------------------------------------------------------------------
|
|
// smoother surface...
|
|
if(mInitSurfaceSmoothing>0) {
|
|
debMsgStd("Surface Smoothing init", DM_MSG, "Performing "<<(mInitSurfaceSmoothing)<<" smoothing timestep ",10);
|
|
#if COMPRESSGRIDS==1
|
|
//errFatal("NYI","COMPRESSGRIDS mInitSurfaceSmoothing",SIMWORLD_INITERROR); return;
|
|
#endif // COMPRESSGRIDS==0
|
|
}
|
|
for(int s=0; s<mInitSurfaceSmoothing; s++) {
|
|
//SGR_FORIJK1(mMaxRefine) {
|
|
|
|
int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
|
|
int lev = mMaxRefine;
|
|
#if COMPRESSGRIDS==0
|
|
for(int k=kstart;k<kend;++k) {
|
|
#else // COMPRESSGRIDS==0
|
|
int kdir = 1; // COMPRT ON
|
|
if(mLevel[lev].setCurr==1) {
|
|
kdir = -1;
|
|
int temp = kend;
|
|
kend = kstart-1;
|
|
kstart = temp-1;
|
|
} // COMPRT
|
|
for(int k=kstart;k!=kend;k+=kdir) {
|
|
#endif // COMPRESSGRIDS==0
|
|
for(int j=1;j<mLevel[lev].lSizey-1;++j) {
|
|
for(int i=1;i<mLevel[lev].lSizex-1;++i) {
|
|
if( ( RFLAG(lev, i,j,k, mLevel[lev].setCurr)& CFInter) ) {
|
|
LbmFloat mass = 0.0;
|
|
//LbmFloat nbdiv;
|
|
//FORDF0 {
|
|
for(int l=0;(l<cDirNum); l++) {
|
|
int ni=i+dfVecX[l], nj=j+dfVecY[l], nk=k+dfVecZ[l];
|
|
if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFFluid ){
|
|
mass += 1.0;
|
|
}
|
|
if( RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr) & CFInter ){
|
|
mass += QCELL(lev, ni,nj,nk, mLevel[lev].setCurr, dMass);
|
|
}
|
|
//nbdiv+=1.0;
|
|
}
|
|
|
|
//errMsg(" I ", PRINT_IJK<<" m"<<mass );
|
|
QCELL(lev, i,j,k, mLevel[lev].setOther, dMass) = (mass/ ((LbmFloat)cDirNum) );
|
|
QCELL(lev, i,j,k, mLevel[lev].setOther, dFfrac) = QCELL(lev, i,j,k, mLevel[lev].setOther, dMass);
|
|
}
|
|
}}}
|
|
|
|
mLevel[lev].setOther = mLevel[lev].setCurr;
|
|
mLevel[lev].setCurr ^= 1;
|
|
}
|
|
// copy back...?
|
|
}
|
|
|
|
/*****************************************************************************/
|
|
/* init part for all freesurface testcases */
|
|
void LbmFsgrSolver::initStandingFluidGradient() {
|
|
// ----------------------------------------------------------------------
|
|
// standing fluid preinit
|
|
const int debugStandingPreinit = 0;
|
|
int haveStandingFluid = 0;
|
|
|
|
#define STANDFLAGCHECK(iindex) \
|
|
if( ( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFInter)) ) || \
|
|
( (RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) & (CFEmpty)) ) ){ \
|
|
if((iindex)>1) { haveStandingFluid=(iindex); } \
|
|
j = mLevel[mMaxRefine].lSizey; i=mLevel[mMaxRefine].lSizex; k=getForZMaxBnd(); \
|
|
continue; \
|
|
}
|
|
int gravIndex[3] = {0,0,0};
|
|
int gravDir[3] = {1,1,1};
|
|
int maxGravComp = 1; // by default y
|
|
int gravComp1 = 0; // by default y
|
|
int gravComp2 = 2; // by default y
|
|
if( ABS(mLevel[mMaxRefine].gravity[0]) > ABS(mLevel[mMaxRefine].gravity[1]) ){ maxGravComp = 0; gravComp1=1; gravComp2=2; }
|
|
if( ABS(mLevel[mMaxRefine].gravity[2]) > ABS(mLevel[mMaxRefine].gravity[0]) ){ maxGravComp = 2; gravComp1=0; gravComp2=1; }
|
|
|
|
int gravIMin[3] = { 0 , 0 , 0 };
|
|
int gravIMax[3] = {
|
|
mLevel[mMaxRefine].lSizex + 0,
|
|
mLevel[mMaxRefine].lSizey + 0,
|
|
mLevel[mMaxRefine].lSizez + 0 };
|
|
if(LBMDIM==2) gravIMax[2] = 1;
|
|
|
|
//int gravDir = 1;
|
|
if( mLevel[mMaxRefine].gravity[maxGravComp] > 0.0 ) {
|
|
// swap directions
|
|
int i=maxGravComp;
|
|
int tmp = gravIMin[i];
|
|
gravIMin[i] = gravIMax[i] - 1;
|
|
gravIMax[i] = tmp - 1;
|
|
gravDir[i] = -1;
|
|
}
|
|
#define PRINTGDIRS \
|
|
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","Z start="<<gravIMin[2]<<" end="<<gravIMax[2]<<" dir="<<gravDir[2] );
|
|
// _PRINTGDIRS;
|
|
|
|
bool gravAbort = false;
|
|
#define GRAVLOOP \
|
|
gravAbort=false; \
|
|
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[0]= gravIMin[0]; (gravIndex[0]!=gravIMax[0])&&(!gravAbort); gravIndex[0] += gravDir[0])
|
|
|
|
GRAVLOOP {
|
|
int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
|
|
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) & (CFEmpty)) ) ) {
|
|
int fluidHeight = (ABS(gravIndex[maxGravComp] - gravIMin[maxGravComp]));
|
|
if(debugStandingPreinit) errMsg("Standing fp","fh="<<fluidHeight<<" gmax="<<gravIMax[maxGravComp]<<" gi="<<gravIndex[maxGravComp] );
|
|
if(fluidHeight>1) {
|
|
haveStandingFluid = fluidHeight; //gravIndex[maxGravComp];
|
|
gravIMax[maxGravComp] = gravIndex[maxGravComp] + gravDir[maxGravComp];
|
|
}
|
|
gravAbort = true; continue;
|
|
}
|
|
} // GRAVLOOP
|
|
// _PRINTGDIRS;
|
|
|
|
LbmFloat fluidHeight;
|
|
fluidHeight = (LbmFloat)(ABS(gravIMax[maxGravComp]-gravIMin[maxGravComp]));
|
|
#if LBM_INCLUDE_TESTSOLVERS==1
|
|
if(mUseTestdata) mpTest->mFluidHeight = (int)fluidHeight;
|
|
#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])<<
|
|
" mgc="<<maxGravComp<<" mc1="<<gravComp1<<" mc2="<<gravComp2<<" dir="<<gravDir[maxGravComp]<<" have="<<haveStandingFluid ,10);
|
|
|
|
if(mDisableStandingFluidInit) {
|
|
debMsgStd("Standing fluid preinit", DM_MSG, "Should be performed - but skipped due to mDisableStandingFluidInit flag set!", 2);
|
|
haveStandingFluid=0;
|
|
}
|
|
|
|
// copy flags and init , as no flags will be changed during grav init
|
|
// also important for corasening later on
|
|
const int lev = mMaxRefine;
|
|
CellFlagType nbflag[LBM_DFNUM], nbored;
|
|
for(int k=getForZMinBnd();k<getForZMaxBnd(mMaxRefine);++k) {
|
|
for(int j=0;j<mLevel[lev].lSizey-0;++j) {
|
|
for(int i=0;i<mLevel[lev].lSizex-0;++i) {
|
|
if( (RFLAG(lev, i,j,k,SRCS(lev)) & (CFFluid)) ) {
|
|
nbored = 0;
|
|
FORDF1 {
|
|
nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
|
|
nbored |= nbflag[l];
|
|
}
|
|
if(nbored&CFBnd) {
|
|
RFLAG(lev, i,j,k,SRCS(lev)) &= (~CFNoBndFluid);
|
|
} else {
|
|
RFLAG(lev, i,j,k,SRCS(lev)) |= CFNoBndFluid;
|
|
}
|
|
}
|
|
RFLAG(lev, i,j,k,TSET(lev)) = RFLAG(lev, i,j,k,SRCS(lev));
|
|
} } }
|
|
|
|
if(haveStandingFluid) {
|
|
int rhoworkSet = mLevel[lev].setCurr;
|
|
myTime_t timestart = getTime(); // FIXME use user time here?
|
|
|
|
GRAVLOOP {
|
|
int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
|
|
//debMsgStd("Standing fluid preinit", DM_MSG, " init check "<<PRINT_IJK<<" "<< haveStandingFluid, 1 );
|
|
if( ( (RFLAG(lev, i,j,k,rhoworkSet) & (CFInter)) ) ||
|
|
( (RFLAG(lev, i,j,k,rhoworkSet) & (CFEmpty)) ) ){
|
|
//gravAbort = true;
|
|
continue;
|
|
}
|
|
|
|
LbmFloat rho = 1.0;
|
|
// 1/6 velocity from denisty gradient, 1/2 for delta of two cells
|
|
rho += 1.0* (fluidHeight-gravIndex[maxGravComp]) *
|
|
(mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega);
|
|
if(debugStandingPreinit)
|
|
if((gravIndex[gravComp1]==gravIMin[gravComp1]) && (gravIndex[gravComp2]==gravIMin[gravComp2])) {
|
|
errMsg("Standing fp","gi="<<gravIndex[maxGravComp]<<" rho="<<rho<<" at "<<PRINT_IJK);
|
|
}
|
|
|
|
if( (RFLAG(lev, i,j,k, rhoworkSet) & CFFluid) ||
|
|
(RFLAG(lev, i,j,k, rhoworkSet) & CFInter) ) {
|
|
FORDF0 { QCELL(lev, i,j,k, rhoworkSet, l) *= rho; }
|
|
QCELL(lev, i,j,k, rhoworkSet, dMass) *= rho;
|
|
}
|
|
|
|
} // GRAVLOOP
|
|
|
|
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);
|
|
|
|
int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
|
|
preinitSteps = (haveStandingFluid>>2); // not much use...?
|
|
//preinitSteps = 0;
|
|
debMsgStd("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
|
|
for(int s=0; s<preinitSteps; s++) {
|
|
// in solver main cpp
|
|
standingFluidPreinit();
|
|
}
|
|
|
|
myTime_t timeend = getTime();
|
|
debMsgStd("Standing fluid preinit", DM_MSG, " done, "<<getTimeString(timeend-timestart), 9);
|
|
#undef NBFLAG
|
|
}
|
|
}
|
|
|
|
|
|
|
|
bool LbmFsgrSolver::checkSymmetry(string idstring)
|
|
{
|
|
bool erro = false;
|
|
bool symm = true;
|
|
int msgs = 0;
|
|
const int maxMsgs = 10;
|
|
const bool markCells = false;
|
|
|
|
//for(int lev=0; lev<=mMaxRefine; lev++) {
|
|
{ int lev = mMaxRefine;
|
|
|
|
// no point if not symm.
|
|
if( (mLevel[lev].lSizex==mLevel[lev].lSizey) && (mLevel[lev].lSizex==mLevel[lev].lSizez)) {
|
|
// ok
|
|
} else {
|
|
return false;
|
|
}
|
|
|
|
for(int s=0; s<2; s++) {
|
|
FSGR_FORIJK1(lev) {
|
|
if(i<(mLevel[lev].lSizex/2)) {
|
|
int inb = (mLevel[lev].lSizey-1-i);
|
|
|
|
if(lev==mMaxRefine) inb -= 1; // FSGR_SYMM_T
|
|
|
|
if( RFLAG(lev, i,j,k,s) != RFLAG(lev, inb,j,k,s) ) { erro = true;
|
|
if(LBMDIM==2) {
|
|
if(msgs<maxMsgs) { msgs++;
|
|
errMsg("EFLAG", PRINT_IJK<<"s"<<s<<" flag "<<RFLAG(lev, i,j,k,s)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" flag "<<RFLAG(lev, inb,j,k,s) );
|
|
}
|
|
}
|
|
if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
|
|
symm = false;
|
|
}
|
|
if( LBM_FLOATNEQ(QCELL(lev, i,j,k,s, dMass), QCELL(lev, inb,j,k,s, dMass)) ) { erro = true;
|
|
if(LBMDIM==2) {
|
|
if(msgs<maxMsgs) { msgs++;
|
|
//debMsgDirect(" mass1 "<<QCELL(lev, i,j,k,s, dMass)<<" mass2 "<<QCELL(lev, inb,j,k,s, dMass) <<std::endl);
|
|
errMsg("EMASS", PRINT_IJK<<"s"<<s<<" mass "<<QCELL(lev, i,j,k,s, dMass)<<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" mass "<<QCELL(lev, inb,j,k,s, dMass) );
|
|
}
|
|
}
|
|
if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
|
|
symm = false;
|
|
}
|
|
|
|
LbmFloat nbrho = QCELL(lev, i,j,k, s, dC);
|
|
FORDF1 { nbrho += QCELL(lev, i,j,k, s, l); }
|
|
LbmFloat otrho = QCELL(lev, inb,j,k, s, dC);
|
|
FORDF1 { otrho += QCELL(lev, inb,j,k, s, l); }
|
|
if( LBM_FLOATNEQ(nbrho, otrho) ) { erro = true;
|
|
if(LBMDIM==2) {
|
|
if(msgs<maxMsgs) { msgs++;
|
|
//debMsgDirect(" rho 1 "<<nbrho <<" rho 2 "<<otrho <<std::endl);
|
|
errMsg("ERHO ", PRINT_IJK<<"s"<<s<<" rho "<<nbrho <<" , at "<<PRINT_VEC(inb,j,k)<<"s"<<s<<" rho "<<otrho );
|
|
}
|
|
}
|
|
if(markCells){ debugMarkCell(lev, i,j,k); debugMarkCell(lev, inb,j,k); }
|
|
symm = false;
|
|
}
|
|
}
|
|
} }
|
|
} // lev
|
|
LbmFloat maxdiv =0.0;
|
|
if(erro) {
|
|
errMsg("SymCheck Failed!", idstring<<" rho maxdiv:"<< maxdiv );
|
|
//if(LBMDIM==2) mPanic = true;
|
|
//return false;
|
|
} else {
|
|
errMsg("SymCheck OK!", idstring<<" rho maxdiv:"<< maxdiv );
|
|
}
|
|
// all ok...
|
|
return symm;
|
|
}// */
|
|
|
|
|
|
void
|
|
LbmFsgrSolver::interpolateCellValues(
|
|
int level,int ei,int ej,int ek,int workSet,
|
|
LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz)
|
|
{
|
|
LbmFloat avgrho = 0.0;
|
|
LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
|
|
LbmFloat cellcnt = 0.0;
|
|
|
|
for(int nbl=1; nbl< cDfNum ; ++nbl) {
|
|
if(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) continue;
|
|
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) & CFInter) ) {
|
|
cellcnt += 1.0;
|
|
for(int rl=0; rl< cDfNum ; ++rl) {
|
|
LbmFloat nbdf = QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
|
|
avgux += (dfDvecX[rl]*nbdf);
|
|
avguy += (dfDvecY[rl]*nbdf);
|
|
avguz += (dfDvecZ[rl]*nbdf);
|
|
avgrho += nbdf;
|
|
}
|
|
}
|
|
}
|
|
|
|
if(cellcnt<=0.0) {
|
|
// no nbs? just use eq.
|
|
avgrho = 1.0;
|
|
avgux = avguy = avguz = 0.0;
|
|
//TTT mNumProblems++;
|
|
#if ELBEEM_PLUGIN!=1
|
|
//mPanic=1;
|
|
// this can happen for worst case moving obj scenarios...
|
|
errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek));
|
|
#endif // ELBEEM_PLUGIN
|
|
} else {
|
|
// init speed
|
|
avgux /= cellcnt; avguy /= cellcnt; avguz /= cellcnt;
|
|
avgrho /= cellcnt;
|
|
}
|
|
|
|
retrho = avgrho;
|
|
retux = avgux;
|
|
retuy = avguy;
|
|
retuz = avguz;
|
|
} // interpolateCellValues
|
|
|
|
|
|
/******************************************************************************
|
|
* instantiation
|
|
*****************************************************************************/
|
|
|
|
//! lbm factory functions
|
|
LbmSolverInterface* createSolver() { return new LbmFsgrSolver(); }
|
|
|
|
|