elbeem update:
- slightly improved and optimized handling of moving obstacles - cleanup of debug messages and minor fixes
This commit is contained in:
Nils Thuerey
parent
d6d2d6f063
commit
b6257305c9
@ -78,6 +78,7 @@ typedef struct elbeemSimulationSettings {
|
||||
short generateVertexVectors;
|
||||
/* strength of surface smoothing */
|
||||
float surfaceSmoothing;
|
||||
// TODO add surf gen flags
|
||||
|
||||
/* global transformation to apply to fluidsim mesh */
|
||||
float surfaceTrafo[4*4];
|
||||
|
||||
@ -86,6 +86,12 @@ void IsoSurface::initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx e
|
||||
|
||||
|
||||
|
||||
/*! Reset all values */
|
||||
void IsoSurface::resetAll(gfxReal val) {
|
||||
int nodes = mSizez*mSizey*mSizex;
|
||||
for(int i=0;i<nodes;i++) { mpData[i] = val; }
|
||||
}
|
||||
|
||||
|
||||
/******************************************************************************
|
||||
* Destructor
|
||||
@ -174,6 +180,10 @@ void IsoSurface::triangulate( void )
|
||||
value[6] = *getData(i+1,j+1,k+1);
|
||||
value[7] = *getData(i ,j+1,k+1);
|
||||
|
||||
/*int bndskip = 0; // BNDOFFT
|
||||
for(int s=0; s<8; s++) if(value[s]==-76.) bndskip++;
|
||||
if(bndskip>0) continue; // */
|
||||
|
||||
// check intersections of isosurface with edges, and calculate cubie index
|
||||
cubeIndex = 0;
|
||||
if (value[0] < mIsoValue) cubeIndex |= 1;
|
||||
|
||||
@ -47,6 +47,9 @@ class IsoSurface :
|
||||
/*! Init ararys etc. */
|
||||
virtual void initializeIsosurface(int setx, int sety, int setz, ntlVec3Gfx extent);
|
||||
|
||||
/*! Reset all values */
|
||||
void resetAll(gfxReal val);
|
||||
|
||||
/*! triangulate the scalar field given by pointer*/
|
||||
void triangulate( void );
|
||||
|
||||
|
||||
105
intern/elbeem/intern/loop_tools.h
Normal file
105
intern/elbeem/intern/loop_tools.h
Normal file
@ -0,0 +1,105 @@
|
||||
|
||||
// advance pointer in main loop
|
||||
#define ADVANCE_POINTERS(p) \
|
||||
ccel += (QCELLSTEP*(p)); \
|
||||
tcel += (QCELLSTEP*(p)); \
|
||||
pFlagSrc+= (p); \
|
||||
pFlagDst+= (p); \
|
||||
i+= (p);
|
||||
|
||||
#define MAX_CALC_ARR 4
|
||||
|
||||
// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
||||
// init region vars
|
||||
#define GRID_REGION_INIT() \
|
||||
const int istart = -1+gridLoopBound; \
|
||||
const int iend = mLevel[mMaxRefine].lSizex-1-gridLoopBound; \
|
||||
LbmFloat calcCurrentMass=0; \
|
||||
LbmFloat calcCurrentVolume=0; \
|
||||
int calcCellsFilled=0; \
|
||||
int calcCellsEmptied=0; \
|
||||
int calcNumUsedCells=0; \
|
||||
|
||||
|
||||
|
||||
|
||||
// -----------------------------------------------------------------------------------
|
||||
// serial stuff
|
||||
#if PARALLEL!=1
|
||||
|
||||
#define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
|
||||
#define LIST_EMPTY(x) mListEmpty.push_back( x );
|
||||
#define LIST_FULL(x) mListFull.push_back( x );
|
||||
|
||||
// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
||||
#define GRID_REGION_START() \
|
||||
{ /* main_region */ \
|
||||
int kstart=getForZMinBnd(), kend=getForZMaxBnd(mMaxRefine); \
|
||||
if(gridLoopBound>0){ kstart=getForZMin1(), kend=getForZMax1(mMaxRefine); } \
|
||||
int kdir = 1; \
|
||||
int jstart = gridLoopBound; \
|
||||
int jend = mLevel[mMaxRefine].lSizey-gridLoopBound; \
|
||||
const int id=0; \
|
||||
LbmFloat *ccel = NULL, *tcel = NULL; \
|
||||
CellFlagType *pFlagSrc=NULL, *pFlagDst=NULL; \
|
||||
if(mLevel[mMaxRefine].setCurr==1) { \
|
||||
kdir = -1; \
|
||||
int temp = kend; \
|
||||
kend = kstart-1; \
|
||||
kstart = temp-1; \
|
||||
temp = id; /* dummy remove warning */ \
|
||||
} \
|
||||
|
||||
|
||||
|
||||
|
||||
#define unused_GRID_REGION_END() \
|
||||
} /* main_region */ \
|
||||
// end unusedGRID_REGION_END
|
||||
|
||||
|
||||
// -----------------------------------------------------------------------------------
|
||||
#else // PARALLEL==1
|
||||
|
||||
#include "paraloop.h"
|
||||
|
||||
#endif // PARALLEL==1
|
||||
|
||||
|
||||
// -----------------------------------------------------------------------------------
|
||||
|
||||
// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
||||
#define GRID_LOOP_START() \
|
||||
for(int k=kstart;k!=kend;k+=kdir) { \
|
||||
pFlagSrc = &RFLAG(lev, istart, jstart, k, SRCS(lev)); \
|
||||
pFlagDst = &RFLAG(lev, istart, jstart, k, TSET(lev)); \
|
||||
ccel = RACPNT(lev, istart, jstart, k, SRCS(lev)); \
|
||||
tcel = RACPNT(lev, istart, jstart, k, TSET(lev)); \
|
||||
for(int j=jstart;j!=jend;++j) { \
|
||||
/* for(int i=0;i<mLevel[lev].lSizex-2; ) { */ \
|
||||
for(int i=istart;i!=iend; ) { \
|
||||
ADVANCE_POINTERS(1); \
|
||||
|
||||
|
||||
|
||||
|
||||
// >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
|
||||
#define GRID_LOOPREG_END() \
|
||||
\
|
||||
} /* i */ \
|
||||
int i=0; \
|
||||
ADVANCE_POINTERS(2*gridLoopBound); \
|
||||
} /* j */ \
|
||||
} /* all cell loop k,j,i */ \
|
||||
if(doReduce) { } /* dummy remove warning */ \
|
||||
} /* main_region */ \
|
||||
\
|
||||
|
||||
|
||||
|
||||
// old loop for COMPRESSGRIDS==0
|
||||
#define old__GRID_LOOP_START() \
|
||||
for(int k=kstart;k<kend;++k) { \
|
||||
for(int j=1;j<mLevel[lev].lSizey-1;++j) { \
|
||||
for(int i=0;i<mLevel[lev].lSizex-2; ) {
|
||||
|
||||
@ -33,9 +33,9 @@ ntlGeometryObject::ntlGeometryObject() :
|
||||
mInitialPos(0.),
|
||||
mcTrans(0.), mcRot(0.), mcScale(1.),
|
||||
mIsAnimated(false),
|
||||
mMovPoints(), //mMovNormals(),
|
||||
mMovPoints(), mMovNormals(),
|
||||
mHaveCachedMov(false),
|
||||
mCachedMovPoints(), //mCachedMovNormals(),
|
||||
mCachedMovPoints(), mCachedMovNormals(),
|
||||
mTriangleDivs1(), mTriangleDivs2(), mTriangleDivs3(),
|
||||
mMovPntsInited(-100.0), mMaxMovPnt(-1),
|
||||
mcGeoActive(1.)
|
||||
@ -169,14 +169,6 @@ void ntlGeometryObject::initialize(ntlRenderGlobals *glob)
|
||||
// always use channel
|
||||
if(!mcGeoActive.isInited()) { mcGeoActive = AnimChannel<double>(geoactive); }
|
||||
|
||||
//if( (mcTrans.accessValues().size()>1) // VALIDATE
|
||||
//|| (mcRot.accessValues().size()>1)
|
||||
//|| (mcScale.accessValues().size()>1)
|
||||
//|| (mcGeoActive.accessValues().size()>1)
|
||||
//|| (mcInitialVelocity.accessValues().size()>1)
|
||||
//) {
|
||||
//mIsAnimated = true;
|
||||
//}
|
||||
checkIsAnimated();
|
||||
|
||||
mIsInitialized = true;
|
||||
@ -340,14 +332,6 @@ void ntlGeometryObject::initChannels(
|
||||
if((act)&&(nAct>0)) { ADD_CHANNEL_FLOAT(mcGeoActive, nAct, act); }
|
||||
if((ivel)&&(nIvel>0)) { ADD_CHANNEL_VEC(mcInitialVelocity, nIvel, ivel); }
|
||||
|
||||
//if( (mcTrans.accessValues().size()>1) // VALIDATE
|
||||
//|| (mcRot.accessValues().size()>1)
|
||||
//|| (mcScale.accessValues().size()>1)
|
||||
//|| (mcGeoActive.accessValues().size()>1)
|
||||
//|| (mcInitialVelocity.accessValues().size()>1)
|
||||
//) {
|
||||
//mIsAnimated = true;
|
||||
//}
|
||||
checkIsAnimated();
|
||||
if(debugInitc) {
|
||||
debMsgStd("ntlGeometryObject::initChannels",DM_MSG,getName()<<
|
||||
@ -422,6 +406,9 @@ void ntlGeometryObject::calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTr
|
||||
mTriangleDivs1.resize( tris.size() );
|
||||
mTriangleDivs2.resize( tris.size() );
|
||||
mTriangleDivs3.resize( tris.size() );
|
||||
|
||||
//fsTri *= 2.; // DEBUG! , wrong init!
|
||||
|
||||
for(size_t i=0; i<tris.size(); i++) {
|
||||
const ntlVec3Gfx p0 = verts[ tris[i].getPoints()[0] ];
|
||||
const ntlVec3Gfx p1 = verts[ tris[i].getPoints()[1] ];
|
||||
@ -432,17 +419,7 @@ void ntlGeometryObject::calcTriangleDivs(vector<ntlVec3Gfx> &verts, vector<ntlTr
|
||||
int divs1=0, divs2=0, divs3=0;
|
||||
if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
|
||||
if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
|
||||
//if(normNoSqrt(side3) > fsTri*fsTri) { divs3 = (int)(norm(side3)/fsTri); }
|
||||
/*if(getMeshAnimated()) {
|
||||
vector<ntlSetVec3f> *verts =mcAniVerts.accessValues();
|
||||
for(int s=0; s<verts->size(); s++) {
|
||||
int tdivs1=0, tdivs2=0, tdivs3=0;
|
||||
if(normNoSqrt(side1) > fsTri*fsTri) { tdivs1 = (int)(norm(side1)/fsTri); }
|
||||
if(normNoSqrt(side2) > fsTri*fsTri) { tdivs2 = (int)(norm(side2)/fsTri); }
|
||||
if(tdivs1>divs1) divs1=tdivs1;
|
||||
if(tdivs2>divs2) divs2=tdivs2;
|
||||
}
|
||||
}*/
|
||||
|
||||
mTriangleDivs1[i] = divs1;
|
||||
mTriangleDivs2[i] = divs2;
|
||||
mTriangleDivs3[i] = divs3;
|
||||
@ -454,15 +431,14 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
|
||||
if((mMovPntsInited==featureSize)&&(!getMeshAnimated())) return;
|
||||
const bool debugMoinit=false;
|
||||
|
||||
//vector<ntlVec3Gfx> movNormals;
|
||||
vector<ntlTriangle> triangles;
|
||||
vector<ntlVec3Gfx> vertices;
|
||||
vector<ntlVec3Gfx> vnormals;
|
||||
int objectId = 1;
|
||||
this->getTriangles(time, &triangles,&vertices,&vnormals,objectId);
|
||||
|
||||
mMovPoints.clear(); //= vertices;
|
||||
//movNormals.clear(); //= vnormals;
|
||||
mMovPoints.clear();
|
||||
mMovNormals.clear();
|
||||
if(debugMoinit) errMsg("ntlGeometryObject::initMovingPoints","Object "<<getName()<<" has v:"<<vertices.size()<<" t:"<<triangles.size() );
|
||||
// no points?
|
||||
if(vertices.size()<1) {
|
||||
@ -523,7 +499,7 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
|
||||
if(dot(mInitialVelocity,n)<0.0) continue;
|
||||
}
|
||||
mMovPoints.push_back(p);
|
||||
//movNormals.push_back(n);
|
||||
mMovNormals.push_back(n);
|
||||
//errMsg("ntlGeometryObject::initMovingPoints","std"<<i<<" p"<<p<<" n"<<n<<" ");
|
||||
}
|
||||
// init points & refine...
|
||||
@ -533,9 +509,9 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
|
||||
const ntlVec3Gfx side1 = vertices[ trips[1] ] - p0;
|
||||
const ntlVec3Gfx side2 = vertices[ trips[2] ] - p0;
|
||||
int divs1=mTriangleDivs1[i], divs2=mTriangleDivs2[i];
|
||||
//if(normNoSqrt(side1) > fsTri*fsTri) { divs1 = (int)(norm(side1)/fsTri); }
|
||||
//if(normNoSqrt(side2) > fsTri*fsTri) { divs2 = (int)(norm(side2)/fsTri); }
|
||||
const ntlVec3Gfx trinorm = getNormalized(cross(side1,side2))*0.5*featureSize;
|
||||
|
||||
const ntlVec3Gfx trinormOrg = getNormalized(cross(side1,side2));
|
||||
const ntlVec3Gfx trinorm = trinormOrg*0.25*featureSize;
|
||||
if(discardInflowBack) {
|
||||
if(dot(mInitialVelocity,trinorm)<0.0) continue;
|
||||
}
|
||||
@ -557,25 +533,16 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
|
||||
//normalize(n);
|
||||
// discard inflow backsides
|
||||
|
||||
mMovPoints.push_back(p);
|
||||
mMovPoints.push_back(p + trinorm); // NEW!?
|
||||
mMovPoints.push_back(p - trinorm);
|
||||
//movNormals.push_back(n);
|
||||
//movNormals.push_back(trinorm);
|
||||
mMovNormals.push_back(trinormOrg);
|
||||
mMovNormals.push_back(trinormOrg);
|
||||
//errMsg("TRINORM","p"<<p<<" n"<<n<<" trin"<<trinorm);
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// duplicate insides
|
||||
//size_t mpsize = mMovPoints.size();
|
||||
//for(size_t i=0; i<mpsize; i++) {
|
||||
//normalize(vnormals[i]);
|
||||
//errMsg("TTAT"," moved:"<<(mMovPoints[i] - mMovPoints[i]*featureSize)<<" org"<<mMovPoints[i]<<" norm"<<mMovPoints[i]<<" fs"<<featureSize);
|
||||
//mMovPoints.push_back(mMovPoints[i] - movNormals[i]*0.5*featureSize);
|
||||
//movNormals.push_back(movNormals[i]);
|
||||
//}
|
||||
|
||||
// find max point
|
||||
mMaxMovPnt = 0;
|
||||
gfxReal dist = normNoSqrt(mMovPoints[0]);
|
||||
@ -594,7 +561,7 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
|
||||
// also do trafo...
|
||||
} else {
|
||||
mCachedMovPoints = mMovPoints;
|
||||
//mCachedMovNormals = movNormals;
|
||||
mCachedMovNormals = mMovNormals;
|
||||
//applyTransformation(time, &mCachedMovPoints, &mCachedMovNormals, 0, mCachedMovPoints.size(), true);
|
||||
applyTransformation(time, &mCachedMovPoints, NULL, 0, mCachedMovPoints.size(), true);
|
||||
mHaveCachedMov = true;
|
||||
@ -610,31 +577,27 @@ void ntlGeometryObject::initMovingPoints(double time, gfxReal featureSize) {
|
||||
void ntlGeometryObject::initMovingPointsAnim(
|
||||
double srctime, vector<ntlVec3Gfx> &srcmovPoints,
|
||||
double dsttime, vector<ntlVec3Gfx> &dstmovPoints,
|
||||
vector<ntlVec3Gfx> *dstmovNormals,
|
||||
gfxReal featureSize,
|
||||
ntlVec3Gfx geostart, ntlVec3Gfx geoend
|
||||
) {
|
||||
const bool debugMoinit=false;
|
||||
|
||||
//vector<ntlVec3Gfx> srcmovNormals;
|
||||
//vector<ntlVec3Gfx> dstmovNormals;
|
||||
|
||||
vector<ntlTriangle> srctriangles;
|
||||
vector<ntlVec3Gfx> srcvertices;
|
||||
vector<ntlVec3Gfx> unused_normals;
|
||||
vector<ntlTriangle> dsttriangles;
|
||||
vector<ntlVec3Gfx> dstvertices;
|
||||
//vector<ntlVec3Gfx> dstnormals;
|
||||
vector<ntlVec3Gfx> dstnormals;
|
||||
int objectId = 1;
|
||||
// TODO optimize? , get rid of normals?
|
||||
unused_normals.clear();
|
||||
this->getTriangles(srctime, &srctriangles,&srcvertices,&unused_normals,objectId);
|
||||
unused_normals.clear();
|
||||
this->getTriangles(dsttime, &dsttriangles,&dstvertices,&unused_normals,objectId);
|
||||
this->getTriangles(dsttime, &dsttriangles,&dstvertices,&dstnormals,objectId);
|
||||
|
||||
srcmovPoints.clear();
|
||||
dstmovPoints.clear();
|
||||
//srcmovNormals.clear();
|
||||
//dstmovNormals.clear();
|
||||
if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has srcv:"<<srcvertices.size()<<" srct:"<<srctriangles.size() );
|
||||
if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","Object "<<getName()<<" has dstv:"<<dstvertices.size()<<" dstt:"<<dsttriangles.size() );
|
||||
// no points?
|
||||
@ -668,7 +631,7 @@ void ntlGeometryObject::initMovingPointsAnim(
|
||||
}
|
||||
for(size_t i=0; i<dstvertices.size(); i++) {
|
||||
dstmovPoints.push_back(dstvertices[i]);
|
||||
//dstmovNormals.push_back(dstnormals[i]);
|
||||
if(dstmovNormals) (*dstmovNormals).push_back(dstnormals[i]);
|
||||
}
|
||||
if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<" dst:"<<dstmovPoints.size()<<" " );
|
||||
// init points & refine...
|
||||
@ -684,8 +647,9 @@ void ntlGeometryObject::initMovingPointsAnim(
|
||||
const ntlVec3Gfx dstp0 = dstvertices[ dsttrips[0] ];
|
||||
const ntlVec3Gfx dstside1 = dstvertices[ dsttrips[1] ] - dstp0;
|
||||
const ntlVec3Gfx dstside2 = dstvertices[ dsttrips[2] ] - dstp0;
|
||||
const ntlVec3Gfx src_trinorm = getNormalized(cross(srcside1,srcside2))*0.5*featureSize;
|
||||
const ntlVec3Gfx dst_trinorm = getNormalized(cross(dstside1,dstside2))*0.5*featureSize;
|
||||
const ntlVec3Gfx src_trinorm = getNormalized(cross(srcside1,srcside2))*0.25*featureSize;
|
||||
const ntlVec3Gfx dst_trinormOrg = getNormalized(cross(dstside1,dstside2));
|
||||
const ntlVec3Gfx dst_trinorm = dst_trinormOrg *0.25*featureSize;
|
||||
//errMsg("ntlGeometryObject::initMovingPointsAnim","Tri1 "<<srcvertices[srctrips[0]]<<","<<srcvertices[srctrips[1]]<<","<<srcvertices[srctrips[2]]<<" "<<divs1<<","<<divs2 );
|
||||
for(int u=0; u<=divs1; u++) {
|
||||
for(int v=0; v<=divs2; v++) {
|
||||
@ -709,60 +673,19 @@ void ntlGeometryObject::initMovingPointsAnim(
|
||||
if((srcp[1]>geoend[1] ) && (dstp[1]>geoend[1] )) continue;
|
||||
if((srcp[2]>geoend[2] ) && (dstp[2]>geoend[2] )) continue;
|
||||
|
||||
//ntlVec3Gfx srcn =
|
||||
//srcnormals[ srctriangles[i].getPoints()[0] ] * (1.0-uf-vf)+
|
||||
//srcnormals[ srctriangles[i].getPoints()[1] ] * uf +
|
||||
//srcnormals[ srctriangles[i].getPoints()[2] ] * vf;
|
||||
//normalize(srcn);
|
||||
srcmovPoints.push_back(srcp);
|
||||
srcmovPoints.push_back(srcp+src_trinorm); // SURFENHTEST
|
||||
srcmovPoints.push_back(srcp-src_trinorm);
|
||||
//srcmovNormals.push_back(srcn);
|
||||
|
||||
//ntlVec3Gfx dstn =
|
||||
//dstnormals[ dsttriangles[i].getPoints()[0] ] * (1.0-uf-vf)+
|
||||
//dstnormals[ dsttriangles[i].getPoints()[1] ] * uf +
|
||||
//dstnormals[ dsttriangles[i].getPoints()[2] ] * vf;
|
||||
//normalize(dstn);
|
||||
dstmovPoints.push_back(dstp);
|
||||
dstmovPoints.push_back(dstp+dst_trinorm); // SURFENHTEST
|
||||
dstmovPoints.push_back(dstp-dst_trinorm);
|
||||
//dstmovNormals.push_back(dstn);
|
||||
|
||||
/*if(debugMoinit && (i>=0)) errMsg("ntlGeometryObject::initMovingPointsAnim"," "<<
|
||||
//srcmovPoints[ srcmovPoints.size()-1 ]<<","<<
|
||||
//srcmovNormals[ srcmovNormals.size()-1 ]<<" "<<
|
||||
//dstmovPoints[ dstmovPoints.size()-1 ]<<","<<
|
||||
//dstmovNormals[ dstmovNormals.size()-1 ]<<" "
|
||||
(srcmovNormals[ srcmovPoints.size()-1 ]-
|
||||
dstmovNormals[ dstmovPoints.size()-1 ])<<" "
|
||||
);
|
||||
// */
|
||||
if(dstmovNormals) {
|
||||
(*dstmovNormals).push_back(dst_trinormOrg);
|
||||
(*dstmovNormals).push_back(dst_trinormOrg); }
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
/*if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<","<<srcmovNormals.size()<<" dst:"<<dstmovPoints.size()<<","<<dstmovNormals.size() );
|
||||
// duplicate insides
|
||||
size_t mpsize = srcmovPoints.size();
|
||||
for(size_t i=0; i<mpsize; i++) {
|
||||
srcmovPoints.push_back(srcmovPoints[i] - srcmovNormals[i]*1.0*featureSize);
|
||||
//? srcnormals.push_back(srcnormals[i]);
|
||||
}
|
||||
mpsize = dstmovPoints.size();
|
||||
for(size_t i=0; i<mpsize; i++) {
|
||||
dstmovPoints.push_back(dstmovPoints[i] - dstmovNormals[i]*1.0*featureSize);
|
||||
//? dstnormals.push_back(dstnormals[i]);
|
||||
}
|
||||
// */
|
||||
|
||||
/*if(debugMoinit) errMsg("ntlGeometryObject::initMovingPointsAnim","stats src:"<<srcmovPoints.size()<<","<<srcmovNormals.size()<<" dst:"<<dstmovPoints.size()<<","<<dstmovNormals.size() );
|
||||
for(size_t i=0; i<srcmovPoints.size(); i++) {
|
||||
ntlVec3Gfx p1 = srcmovPoints[i];
|
||||
ntlVec3Gfx p2 = dstmovPoints[i];
|
||||
gfxReal len = norm(p1-p2);
|
||||
if(len>0.01) { errMsg("ntlGeometryObject::initMovingPointsAnim"," i"<<i<<" "<< p1<<" "<<p2<<" "<<len); }
|
||||
} // */
|
||||
|
||||
// find max point not necessary
|
||||
debMsgStd("ntlGeometryObject::initMovingPointsAnim",DM_MSG,"Object "<<getName()<<" inited v:"<<srcvertices.size()<<"->"<<srcmovPoints.size()<<","<<dstmovPoints.size() , 5);
|
||||
}
|
||||
@ -772,8 +695,8 @@ void ntlGeometryObject::getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3G
|
||||
if(mHaveCachedMov) {
|
||||
ret = mCachedMovPoints;
|
||||
if(norms) {
|
||||
//*norms = mCachedMovNormals;
|
||||
errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
|
||||
*norms = mCachedMovNormals;
|
||||
//errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
|
||||
}
|
||||
//errMsg ("ntlGeometryObject::getMovingPoints","Object "<<getName()<<" used cached points "); // DEBUG
|
||||
return;
|
||||
@ -781,8 +704,8 @@ void ntlGeometryObject::getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3G
|
||||
|
||||
ret = mMovPoints;
|
||||
if(norms) {
|
||||
errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
|
||||
//*norms = mMovNormals;
|
||||
//errMsg("ntlGeometryObject","getMovingPoints - Normals currently unused!");
|
||||
*norms = mMovNormals;
|
||||
}
|
||||
}
|
||||
|
||||
|
||||
@ -118,6 +118,7 @@ class ntlGeometryObject : public ntlGeometryClass
|
||||
void initMovingPointsAnim(
|
||||
double srctime, vector<ntlVec3Gfx> &srcpoints,
|
||||
double dsttime, vector<ntlVec3Gfx> &dstpoints,
|
||||
vector<ntlVec3Gfx> *dstnormals,
|
||||
gfxReal featureSize, ntlVec3Gfx geostart, ntlVec3Gfx geoend );
|
||||
/*! Prepare points for moving objects (copy into ret) */
|
||||
void getMovingPoints(vector<ntlVec3Gfx> &ret, vector<ntlVec3Gfx> *norms = NULL);
|
||||
@ -181,11 +182,11 @@ class ntlGeometryObject : public ntlGeometryClass
|
||||
|
||||
/*! moving point/normal storage */
|
||||
vector<ntlVec3Gfx> mMovPoints;
|
||||
// unused vector<ntlVec3Gfx> mMovNormals;
|
||||
vector<ntlVec3Gfx> mMovNormals;
|
||||
/*! cached points for non moving objects/timeslots */
|
||||
bool mHaveCachedMov;
|
||||
vector<ntlVec3Gfx> mCachedMovPoints;
|
||||
// unused vector<ntlVec3Gfx> mCachedMovNormals;
|
||||
vector<ntlVec3Gfx> mCachedMovNormals;
|
||||
/*! precomputed triangle divisions */
|
||||
vector<int> mTriangleDivs1,mTriangleDivs2,mTriangleDivs3;
|
||||
/*! inited? */
|
||||
|
||||
@ -439,7 +439,7 @@ int SimulationObject::checkCallerStatus(int status, int frame) {
|
||||
}
|
||||
}
|
||||
|
||||
errMsg("SimulationObject::checkCallerStatus","s="<<status<<",f="<<frame<<" "<<this->getName()<<" ret="<<ret);
|
||||
//debMsgStd("SimulationObject::checkCallerStatus",DM_MSG, "s="<<status<<",f="<<frame<<" "<<this->getName()<<" ret="<<ret);
|
||||
if(isSimworldOk()) return 0;
|
||||
return 1;
|
||||
}
|
||||
|
||||
@ -21,10 +21,6 @@
|
||||
|
||||
void LbmFsgrSolver::coarseCalculateFluxareas(int lev)
|
||||
{
|
||||
#if LBM_NOADCOARSENING==1
|
||||
if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
|
||||
lev =0; // get rid of warnings...
|
||||
#else
|
||||
FSGR_FORIJK_BOUNDS(lev) {
|
||||
if( RFLAG(lev, i,j,k,mLevel[lev].setCurr) & CFFluid) {
|
||||
if( RFLAG(lev+1, i*2,j*2,k*2,mLevel[lev+1].setCurr) & CFGrFromCoarse) {
|
||||
@ -50,15 +46,10 @@ void LbmFsgrSolver::coarseCalculateFluxareas(int lev)
|
||||
}
|
||||
} // } TEST DEBUG
|
||||
if(!this->mSilent){ debMsgStd("coarseCalculateFluxareas",DM_MSG,"level "<<lev<<" calculated", 7); }
|
||||
#endif //! LBM_NOADCOARSENING==1
|
||||
}
|
||||
|
||||
void LbmFsgrSolver::coarseAdvance(int lev)
|
||||
{
|
||||
#if LBM_NOADCOARSENING==1
|
||||
if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
|
||||
lev =0; // get rid of warnings...
|
||||
#else
|
||||
LbmFloat calcCurrentMass = 0.0;
|
||||
LbmFloat calcCurrentVolume = 0.0;
|
||||
|
||||
@ -177,10 +168,9 @@ void LbmFsgrSolver::coarseAdvance(int lev)
|
||||
mLevel[lev].lmass = calcCurrentMass * mLevel[lev].lcellfactor;
|
||||
mLevel[lev].lvolume = calcCurrentVolume * mLevel[lev].lcellfactor;
|
||||
#if ELBEEM_PLUGIN!=1
|
||||
errMsg("DFINI", " m l"<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<" lcf="<< mLevel[lev].lcellfactor );
|
||||
errMsg("DFINI", " v l"<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<" lcf="<< mLevel[lev].lcellfactor );
|
||||
debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "mass: lev="<<lev<<" m="<<mLevel[lev].lmass<<" c="<<calcCurrentMass<<" lcf="<< mLevel[lev].lcellfactor, 8 );
|
||||
debMsgStd("LbmFsgrSolver::coarseAdvance",DM_NOTIFY, "volume: lev="<<lev<<" v="<<mLevel[lev].lvolume<<" c="<<calcCurrentVolume<<" lcf="<< mLevel[lev].lcellfactor, 8 );
|
||||
#endif // ELBEEM_PLUGIN
|
||||
#endif //! LBM_NOADCOARSENING==1
|
||||
}
|
||||
|
||||
/*****************************************************************************/
|
||||
@ -191,10 +181,6 @@ void LbmFsgrSolver::coarseAdvance(int lev)
|
||||
// get dfs from level (lev+1) to (lev) coarse border nodes
|
||||
void LbmFsgrSolver::coarseRestrictFromFine(int lev)
|
||||
{
|
||||
#if LBM_NOADCOARSENING==1
|
||||
if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
|
||||
lev =0; // get rid of warnings...
|
||||
#else
|
||||
if((lev<0) || ((lev+1)>mMaxRefine)) return;
|
||||
#if FSGR_STRICT_DEBUG==1
|
||||
// reset all unused cell values to invalid
|
||||
@ -245,15 +231,9 @@ void LbmFsgrSolver::coarseRestrictFromFine(int lev)
|
||||
} // & fluid
|
||||
}}}
|
||||
if(!this->mSilent){ errMsg("coarseRestrictFromFine"," from l"<<(lev+1)<<",s"<<mLevel[lev+1].setCurr<<" to l"<<lev<<",s"<<mLevel[lev].setCurr); }
|
||||
#endif //! LBM_NOADCOARSENING==1
|
||||
}
|
||||
|
||||
bool LbmFsgrSolver::adaptGrid(int lev) {
|
||||
#if LBM_NOADCOARSENING==1
|
||||
if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
|
||||
lev =0; // get rid of warnings...
|
||||
return false;
|
||||
#else
|
||||
if((lev<0) || ((lev+1)>mMaxRefine)) return false;
|
||||
bool change = false;
|
||||
{ // refinement, PASS 1-3
|
||||
@ -706,7 +686,6 @@ bool LbmFsgrSolver::adaptGrid(int lev) {
|
||||
|
||||
if(!this->mSilent){ errMsg("adaptGrid"," for l"<<lev<<" done " ); }
|
||||
return change;
|
||||
#endif //! LBM_NOADCOARSENING==1
|
||||
}
|
||||
|
||||
/*****************************************************************************/
|
||||
@ -715,10 +694,6 @@ bool LbmFsgrSolver::adaptGrid(int lev) {
|
||||
|
||||
void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, int dstSet)
|
||||
{
|
||||
#if LBM_NOADCOARSENING==1
|
||||
if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
|
||||
i=j=k=srcSet=dstSet=lev =0; // get rid of warnings...
|
||||
#else
|
||||
LbmFloat *ccel = RACPNT(lev+1, 2*i,2*j,2*k,srcSet);
|
||||
LbmFloat *tcel = RACPNT(lev , i,j,k ,dstSet);
|
||||
|
||||
@ -862,15 +837,9 @@ void LbmFsgrSolver::coarseRestrictCell(int lev, int i,int j,int k, int srcSet, i
|
||||
RAC(tcel, dWT) = (lcsmeq[dWT] + (MSRC_WT-lcsmeq[dWT] )*lcsmdfscale);
|
||||
RAC(tcel, dWB) = (lcsmeq[dWB] + (MSRC_WB-lcsmeq[dWB] )*lcsmdfscale);
|
||||
# endif // OPT3D==0
|
||||
#endif //! LBM_NOADCOARSENING==1
|
||||
}
|
||||
|
||||
void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int dstSet, LbmFloat t, CellFlagType flagSet, bool markNbs) {
|
||||
#if LBM_NOADCOARSENING==1
|
||||
if(mMaxRefine>0) errMsg("LbmFsgrSolver","Adaptive Coarsening not compiled, but refinement switched on ("<<mMaxRefine<<")!");
|
||||
i=j=k=dstSet=lev =0; // get rid of warnings...
|
||||
t=0.0; flagSet=0; markNbs=false;
|
||||
#else
|
||||
LbmFloat rho=0.0, ux=0.0, uy=0.0, uz=0.0;
|
||||
LbmFloat intDf[19] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
|
||||
|
||||
@ -952,7 +921,6 @@ void LbmFsgrSolver::interpolateCellFromCoarse(int lev, int i, int j,int k, int d
|
||||
|
||||
IDF_WRITEBACK;
|
||||
return;
|
||||
#endif //! LBM_NOADCOARSENING==1
|
||||
}
|
||||
|
||||
|
||||
@ -1008,8 +976,9 @@ void LbmFsgrSolver::adaptTimestep() {
|
||||
|
||||
LbmFloat dtdiff = fabs(newdt - this->mpParam->getTimestep());
|
||||
if(!this->mSilent) {
|
||||
debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt<<" max"<<this->mpParam->getMaxTimestep()<<" min"<<this->mpParam->getMinTimestep()<<" diff"<<dtdiff<<
|
||||
" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep) , 10); }
|
||||
debMsgStd("LbmFsgrSolver::TAdp",DM_MSG, "new"<<newdt
|
||||
<<" max"<<this->mpParam->getMaxTimestep()<<" min"<<this->mpParam->getMinTimestep()<<" diff"<<dtdiff
|
||||
<<" simt:"<<mSimulationTime<<" minsteps:"<<(mSimulationTime/mMaxTimestep)<<" maxsteps:"<<(mSimulationTime/mMinTimestep) , 10); }
|
||||
|
||||
// in range, and more than X% change?
|
||||
//if( newdt < this->mpParam->getTimestep() ) // DEBUG
|
||||
@ -1025,7 +994,8 @@ void LbmFsgrSolver::adaptTimestep() {
|
||||
rescale = true;
|
||||
if(!this->mSilent) {
|
||||
debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"\n\n\n\n",10);
|
||||
debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: new="<<newdt<<" old="<<this->mpParam->getTimestep()<<" maxSpeed:"<<this->mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<this->mStepCnt, 10 );
|
||||
debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: new="<<newdt<<" old="<<this->mpParam->getTimestep()
|
||||
<<" maxSpeed:"<<this->mpParam->getSimulationMaxSpeed()<<" next:"<<nextmax<<" step:"<<this->mStepCnt, 10 );
|
||||
debMsgStd("LbmFsgrSolver::TAdp",DM_NOTIFY,"Timestep change: "<<
|
||||
"rhoAvg="<<rhoAvg<<" cMass="<<mCurrentMass<<" cVol="<<mCurrentVolume,10);
|
||||
}
|
||||
|
||||
@ -106,6 +106,10 @@
|
||||
#endif
|
||||
#endif
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
#include "solver_test.h"
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
|
||||
/*****************************************************************************/
|
||||
/*! cell access classes */
|
||||
class UniformFsgrCellIdentifier :
|
||||
@ -216,10 +220,10 @@ class LbmFsgrSolver :
|
||||
//! notify object that dump is in progress (e.g. for field dump)
|
||||
virtual void notifySolverOfDump(int dumptype, int frameNr,char *frameNrStr,string outfilename);
|
||||
|
||||
#if LBM_USE_GUI==1
|
||||
# if LBM_USE_GUI==1
|
||||
//! show simulation info (implement LbmSolverInterface pure virtual func)
|
||||
virtual void debugDisplay(int set);
|
||||
#endif
|
||||
# endif
|
||||
|
||||
// implement CellIterator<UniformFsgrCellIdentifier> interface
|
||||
typedef UniformFsgrCellIdentifier stdCellId;
|
||||
@ -254,6 +258,7 @@ class LbmFsgrSolver :
|
||||
LBM_INLINED void initEmptyCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass);
|
||||
LBM_INLINED void initVelocityCell(int level, int i,int j,int k, CellFlagType flag, LbmFloat rho, LbmFloat mass, LbmVec vel);
|
||||
LBM_INLINED void changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
|
||||
LBM_INLINED void forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag);
|
||||
//! interpolate velocity and density at a given position
|
||||
void interpolateCellValues(int level,int ei,int ej,int ek,int workSet, LbmFloat &retrho, LbmFloat &retux, LbmFloat &retuy, LbmFloat &retuz);
|
||||
|
||||
@ -279,11 +284,11 @@ class LbmFsgrSolver :
|
||||
virtual vector<ntlGeometryObject*> getDebugObjects();
|
||||
|
||||
// gui/output debugging functions
|
||||
#if LBM_USE_GUI==1
|
||||
# if LBM_USE_GUI==1
|
||||
virtual void debugDisplayNode(int dispset, CellIdentifierInterface* cell );
|
||||
virtual void lbmDebugDisplay(int dispset);
|
||||
virtual void lbmMarkedCellDisplay();
|
||||
#endif // LBM_USE_GUI==1
|
||||
# endif // LBM_USE_GUI==1
|
||||
virtual void debugPrintNodeInfo(CellIdentifierInterface* cell, int forceSet=-1);
|
||||
|
||||
//! for raytracing, preprocess
|
||||
@ -299,7 +304,9 @@ class LbmFsgrSolver :
|
||||
// debugging use CellIterator interface to mark cell
|
||||
void debugMarkCellCall(int level, int vi,int vj,int vk);
|
||||
|
||||
// loop over grid, stream&collide update
|
||||
void mainLoop(int lev);
|
||||
// change time step size
|
||||
void adaptTimestep();
|
||||
//! init mObjectSpeeds for current parametrization
|
||||
void recalculateObjectSpeeds();
|
||||
@ -321,6 +328,11 @@ class LbmFsgrSolver :
|
||||
LBM_INLINED int getForZMaxBnd(int lev);
|
||||
LBM_INLINED int getForZMax1(int lev);
|
||||
|
||||
// touch grid and flags once
|
||||
void preinitGrids();
|
||||
// one relaxation step for standing fluid
|
||||
void standingFluidPreinit();
|
||||
|
||||
|
||||
// member vars
|
||||
|
||||
@ -355,6 +367,20 @@ class LbmFsgrSolver :
|
||||
LbmFloat mGfxEndTime;
|
||||
//! smoother surface initialization?
|
||||
int mInitSurfaceSmoothing;
|
||||
//! surface generation settings, default is all off (=0)
|
||||
// each flag switches side on off, fssgNoObs is for obstacle sides
|
||||
// -1 equals all on
|
||||
typedef enum {
|
||||
fssgNormal = 0,
|
||||
fssgNoNorth = 1,
|
||||
fssgNoSouth = 2,
|
||||
fssgNoEast = 4,
|
||||
fssgNoWest = 8,
|
||||
fssgNoTop = 16,
|
||||
fssgNoBottom = 32,
|
||||
fssgNoObs = 64
|
||||
} fsSurfaceGen;
|
||||
int mFsSurfGenSetting;
|
||||
|
||||
//! lock time step down switching
|
||||
int mTimestepReduceLock;
|
||||
@ -392,10 +418,6 @@ class LbmFsgrSolver :
|
||||
int mMaxNoCells, mMinNoCells;
|
||||
LONGINT mAvgNumUsedCells;
|
||||
|
||||
//! for interactive - how to drop drops?
|
||||
int mDropMode;
|
||||
LbmFloat mDropSize;
|
||||
LbmVec mDropSpeed;
|
||||
//! precalculated objects speeds for current parametrization
|
||||
vector<LbmVec> mObjectSpeeds;
|
||||
//! partslip bc. values for obstacle boundary conditions
|
||||
@ -406,6 +428,7 @@ class LbmFsgrSolver :
|
||||
//! permanent movobj vert storage
|
||||
vector<ntlVec3Gfx> mMOIVertices;
|
||||
vector<ntlVec3Gfx> mMOIVerticesOld;
|
||||
vector<ntlVec3Gfx> mMOINormals;
|
||||
|
||||
//! get isofield weights
|
||||
int mIsoWeightMethod;
|
||||
@ -443,6 +466,8 @@ class LbmFsgrSolver :
|
||||
|
||||
//! debug function to disable standing f init
|
||||
int mDisableStandingFluidInit;
|
||||
//! init 2d with skipped Y/Z coords
|
||||
bool mInit2dYZ;
|
||||
//! debug function to force tadap syncing
|
||||
int mForceTadapRefine;
|
||||
//! border cutoff value
|
||||
@ -463,14 +488,31 @@ class LbmFsgrSolver :
|
||||
# endif // FSGR_STRICT_DEBUG==1
|
||||
|
||||
bool mUseTestdata;
|
||||
# if LBM_INCLUDE_TESTSOLVERS==1
|
||||
// test functions
|
||||
LbmTestdata *mpTest;
|
||||
void initTestdata();
|
||||
void destroyTestdata();
|
||||
void handleTestdata();
|
||||
void set3dHeight(int ,int );
|
||||
|
||||
public: // former LbmModelLBGK functions
|
||||
void initCpdata();
|
||||
void handleCpdata();
|
||||
void cpDebugDisplay(int dispset);
|
||||
|
||||
void testXchng();
|
||||
public:
|
||||
// needed for testdata
|
||||
void find3dHeight(int i,int j, LbmFloat prev, LbmFloat &ret, LbmFloat *retux, LbmFloat *retuy, LbmFloat *retuz);
|
||||
# endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
|
||||
// former LbmModelLBGK functions
|
||||
// relaxation funtions - implemented together with relax macros
|
||||
static inline LbmFloat getVelVecLen(int l, LbmFloat ux,LbmFloat uy,LbmFloat uz);
|
||||
static inline LbmFloat getCollideEq(int l, LbmFloat rho, LbmFloat ux, LbmFloat uy, LbmFloat uz);
|
||||
inline LbmFloat getLesNoneqTensorCoeff( LbmFloat df[], LbmFloat feq[] );
|
||||
inline LbmFloat getLesOmega(LbmFloat omega, LbmFloat csmago, LbmFloat Qo);
|
||||
inline void collideArrays( int i, int j, int k, // position - more for debugging
|
||||
inline void collideArrays( int lev, int i, int j, int k, // position - more for debugging
|
||||
LbmFloat df[], LbmFloat &outrho, // out only!
|
||||
// velocity modifiers (returns actual velocity!)
|
||||
LbmFloat &mux, LbmFloat &muy, LbmFloat &muz,
|
||||
@ -479,12 +521,10 @@ class LbmFsgrSolver :
|
||||
|
||||
|
||||
// former LBM models
|
||||
public:
|
||||
//! shorten static const definitions
|
||||
# define STCON static const
|
||||
|
||||
//! shorten static const definitions
|
||||
#define STCON static const
|
||||
|
||||
#if LBMDIM==3
|
||||
# if LBMDIM==3
|
||||
|
||||
//! id string of solver
|
||||
virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D3Q19]"); }
|
||||
@ -581,7 +621,7 @@ class LbmFsgrSolver :
|
||||
static LbmFloat lesCoeffDiag[ (3-1)*(3-1) ][ 27 ];
|
||||
static LbmFloat lesCoeffOffdiag[ 3 ][ 27 ];
|
||||
|
||||
#else // end LBMDIM==3 , LBMDIM==2
|
||||
# else // end LBMDIM==3 , LBMDIM==2
|
||||
|
||||
//! id string of solver
|
||||
virtual string getIdString() { return string("FreeSurfaceFsgrSolver[BGK_D2Q9]"); }
|
||||
@ -667,7 +707,7 @@ class LbmFsgrSolver :
|
||||
static LbmFloat lesCoeffDiag[ (2-1)*(2-1) ][ 9 ];
|
||||
static LbmFloat lesCoeffOffdiag[ 2 ][ 9 ];
|
||||
|
||||
#endif // LBMDIM==2
|
||||
# endif // LBMDIM==2
|
||||
};
|
||||
|
||||
#undef STCON
|
||||
@ -691,7 +731,8 @@ class LbmFsgrSolver :
|
||||
|
||||
// flag array defines -----------------------------------------------------------------------------------------------
|
||||
|
||||
// lbm testsolver get index define
|
||||
// lbm testsolver get index define, note - ignores is (set) as flag
|
||||
// array is only a single entry
|
||||
#define _LBMGI(level, ii,ij,ik, is) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
|
||||
|
||||
//! flag array acces macro
|
||||
@ -701,7 +742,6 @@ class LbmFsgrSolver :
|
||||
|
||||
// array handling -----------------------------------------------------------------------------------------------
|
||||
|
||||
//#define _LBMQI(level, ii,ij,ik, is, lunused) ( ((is)*mLevel[level].lOffsz) + (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
|
||||
#define _LBMQI(level, ii,ij,ik, is, lunused) ( (mLevel[level].lOffsy*(ik)) + (mLevel[level].lOffsx*(ij)) + (ii) )
|
||||
#define _QCELL(level,xx,yy,zz,set,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx),(yy),(zz),(set), l)*dTotalNum +(l)])
|
||||
#define _QCELL_NB(level,xx,yy,zz,set, dir,l) (mLevel[level].mprsCells[(set)][ LBMQI((level),(xx)+this->dfVecX[dir],(yy)+this->dfVecY[dir],(zz)+this->dfVecZ[dir],set, l)*dTotalNum +(l)])
|
||||
@ -835,12 +875,17 @@ inline LbmFloat LbmFsgrSolver::getCollideEq(int l, LbmFloat rho, LbmFloat ux, L
|
||||
+ rho - (3.0/2.0*(ux*ux + uy*uy + uz*uz))
|
||||
+ 3.0 *tmp
|
||||
+ 9.0/2.0 *(tmp*tmp) )
|
||||
);
|
||||
); // */
|
||||
};
|
||||
|
||||
/*****************************************************************************/
|
||||
/* init a given cell with flag, density, mass and equilibrium dist. funcs */
|
||||
|
||||
void LbmFsgrSolver::forceChangeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
|
||||
// also overwrite persisting flags
|
||||
// function is useful for tracking accesses...
|
||||
RFLAG(level,xx,yy,zz,set) = newflag;
|
||||
}
|
||||
void LbmFsgrSolver::changeFlag(int level, int xx,int yy,int zz,int set,CellFlagType newflag) {
|
||||
CellFlagType pers = RFLAG(level,xx,yy,zz,set) & CFPersistMask;
|
||||
RFLAG(level,xx,yy,zz,set) = newflag | pers;
|
||||
@ -857,7 +902,6 @@ LbmFsgrSolver::initEmptyCell(int level, int i,int j,int k, CellFlagType flag, Lb
|
||||
RAC(ecel, dMass) = mass;
|
||||
RAC(ecel, dFfrac) = mass/rho;
|
||||
RAC(ecel, dFlux) = FLUX_INIT;
|
||||
//RFLAG(level, i,j,k, workSet)= flag;
|
||||
changeFlag(level, i,j,k, workSet, flag);
|
||||
|
||||
workSet ^= 1;
|
||||
@ -875,7 +919,6 @@ LbmFsgrSolver::initVelocityCell(int level, int i,int j,int k, CellFlagType flag,
|
||||
RAC(ecel, dMass) = mass;
|
||||
RAC(ecel, dFfrac) = mass/rho;
|
||||
RAC(ecel, dFlux) = FLUX_INIT;
|
||||
//RFLAG(level, i,j,k, workSet) = flag;
|
||||
changeFlag(level, i,j,k, workSet, flag);
|
||||
|
||||
workSet ^= 1;
|
||||
|
||||
@ -16,6 +16,11 @@
|
||||
// all boundary types at once
|
||||
#define FGI_ALLBOUNDS ( FGI_BNDNO | FGI_BNDFREE | FGI_BNDPART | FGI_MBNDINFLOW | FGI_MBNDOUTFLOW )
|
||||
|
||||
// helper for 2d init
|
||||
#define SWAPYZ(vec) { \
|
||||
const LbmFloat tmp = (vec)[2]; \
|
||||
(vec)[2] = (vec)[1]; (vec)[1] = tmp; }
|
||||
|
||||
|
||||
/*****************************************************************************/
|
||||
//! common variables
|
||||
@ -302,15 +307,14 @@ LbmFsgrSolver::LbmFsgrSolver() :
|
||||
mpPreviewSurface(NULL),
|
||||
mTimeAdap(true), mForceTimeStepReduce(false),
|
||||
mFVHeight(0.0), mFVArea(1.0), mUpdateFVHeight(false),
|
||||
mInitSurfaceSmoothing(0),
|
||||
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),
|
||||
mDropMode(1), mDropSize(0.15), mDropSpeed(0.0),
|
||||
mObjectSpeeds(), mObjectPartslips(), mObjectMassMovnd(),
|
||||
mMOIVertices(), mMOIVerticesOld(),
|
||||
mMOIVertices(), mMOIVerticesOld(), mMOINormals(),
|
||||
mIsoWeightMethod(1),
|
||||
mMaxRefine(1),
|
||||
mDfScaleUp(-1.0), mDfScaleDown(-1.0),
|
||||
@ -319,6 +323,7 @@ LbmFsgrSolver::LbmFsgrSolver() :
|
||||
mLastOmega(1e10), mLastGravity(1e10),
|
||||
mNumInvIfTotal(0), mNumFsgrChanges(0),
|
||||
mDisableStandingFluidInit(0),
|
||||
mInit2dYZ(false),
|
||||
mForceTadapRefine(-1), mCutoff(-1)
|
||||
{
|
||||
// not much to do here...
|
||||
@ -456,12 +461,21 @@ void LbmFsgrSolver::parseAttrList()
|
||||
this->mSmoothSurface = this->mpAttrs->readFloat("smoothsurface", this->mSmoothSurface, "SimulationLbm","mSmoothSurface", false );
|
||||
this->mSmoothNormals = this->mpAttrs->readFloat("smoothnormals", this->mSmoothNormals, "SimulationLbm","mSmoothNormals", false );
|
||||
|
||||
mFsSurfGenSetting = this->mpAttrs->readInt("fssurfgen", mFsSurfGenSetting, "SimulationLbm","mFsSurfGenSetting", false );
|
||||
if(mFsSurfGenSetting==-1) {
|
||||
// all on
|
||||
mFsSurfGenSetting =
|
||||
fssgNormal | fssgNoNorth | fssgNoSouth | fssgNoEast |
|
||||
fssgNoWest | fssgNoTop | fssgNoBottom | fssgNoObs ;
|
||||
}
|
||||
|
||||
// refinement
|
||||
mMaxRefine = this->mRefinementDesired;
|
||||
mMaxRefine = this->mpAttrs->readInt("maxrefine", mMaxRefine ,"LbmFsgrSolver", "mMaxRefine", false);
|
||||
if(mMaxRefine<0) mMaxRefine=0;
|
||||
if(mMaxRefine>FSGR_MAXNOOFLEVELS) mMaxRefine=FSGR_MAXNOOFLEVELS-1;
|
||||
mDisableStandingFluidInit = this->mpAttrs->readInt("disable_stfluidinit", mDisableStandingFluidInit,"LbmFsgrSolver", "mDisableStandingFluidInit", false);
|
||||
mInit2dYZ = this->mpAttrs->readBool("init2dyz", mInit2dYZ,"LbmFsgrSolver", "mInit2dYZ", false);
|
||||
mForceTadapRefine = this->mpAttrs->readInt("forcetadaprefine", mForceTadapRefine,"LbmFsgrSolver", "mForceTadapRefine", false);
|
||||
|
||||
// demo mode settings
|
||||
@ -489,7 +503,7 @@ void LbmFsgrSolver::parseAttrList()
|
||||
mUseTestdata = 0;
|
||||
if(mFarFieldSize>=2.) mUseTestdata=1; // equiv. to test solver check
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS!=1
|
||||
if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
|
||||
//if(mUseTestdata) { mMaxRefine=0; } // force fsgr off
|
||||
|
||||
if(mMaxRefine==0) mInitialCsmago=0.02;
|
||||
mInitialCsmago = this->mpAttrs->readFloat("csmago", mInitialCsmago, "SimulationLbm","mInitialCsmago", false );
|
||||
@ -513,6 +527,7 @@ void LbmFsgrSolver::initLevelOmegas()
|
||||
this->mOmega = this->mpParam->calculateOmega(mSimulationTime);
|
||||
this->mGravity = vec2L( this->mpParam->calculateGravity(mSimulationTime) );
|
||||
this->mSurfaceTension = 0.; //this->mpParam->calculateSurfaceTension(); // unused
|
||||
if(mInit2dYZ) { SWAPYZ(mGravity); }
|
||||
|
||||
// check if last init was ok
|
||||
LbmFloat gravDelta = norm(this->mGravity-mLastGravity);
|
||||
@ -831,16 +846,33 @@ bool LbmFsgrSolver::initializeSolverMemory()
|
||||
isoend[2] = se;
|
||||
twodOff = 2;
|
||||
}
|
||||
int isosx = this->mSizex+2;
|
||||
int isosy = this->mSizey+2;
|
||||
int isosz = this->mSizez+2+twodOff;
|
||||
|
||||
// MPT
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
|
||||
int xfac = 2;
|
||||
isosx *= xfac;
|
||||
isoend[0] = isostart[0] + (isoend[0]-isostart[0])*(LbmFloat)xfac;
|
||||
}
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
|
||||
//errMsg(" SETISO ", " "<<isostart<<" - "<<isoend<<" "<<(((isoend[0]-isostart[0]) / (LbmFloat)(this->mSizex+1.0))*0.5)<<" "<<(LbmFloat)(this->mSizex+1.0)<<" " );
|
||||
this->mpIso->setStart( vec2G(isostart) );
|
||||
this->mpIso->setEnd( vec2G(isoend) );
|
||||
mpIso->setStart( vec2G(isostart) );
|
||||
mpIso->setEnd( vec2G(isoend) );
|
||||
LbmVec isodist = isoend-isostart;
|
||||
this->mpIso->initializeIsosurface( this->mSizex+2, this->mSizey+2, this->mSizez+2+twodOff, vec2G(isodist) );
|
||||
for(int ak=0;ak<this->mSizez+2+twodOff;ak++)
|
||||
for(int aj=0;aj<this->mSizey+2;aj++)
|
||||
for(int ai=0;ai<this->mSizex+2;ai++) { *this->mpIso->getData(ai,aj,ak) = 0.0; }
|
||||
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) = RFLAG(lev,i,j,k,0) = 0; // reset for changeFlag usage
|
||||
@ -968,6 +1000,45 @@ bool LbmFsgrSolver::initializeSolverGrids() {
|
||||
}
|
||||
}
|
||||
// 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; }
|
||||
|
||||
@ -1083,13 +1154,14 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
|
||||
|
||||
#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; \
|
||||
const int k=0; \
|
||||
|
||||
#else // LBMDIM -> 3
|
||||
#define POS2GRID_CHECK(vec,n) \
|
||||
@ -1117,7 +1189,25 @@ bool LbmFsgrSolver::initializeSolverPostinit() {
|
||||
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, this->dfInv[l] )<<" partv="<<partv<<",id="<<(int)(mnbf>>24)<<" newval="<<newval ); / part slip debug */ \
|
||||
/* m[l] = (RAC(ccel, this->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
|
||||
@ -1128,6 +1218,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
// 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 + this->mpParam->getTimestep();
|
||||
@ -1180,6 +1271,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
|
||||
otype = ntype = CFInvalid;
|
||||
switch(obj->getGeoInitType()) {
|
||||
/*
|
||||
case FGI_BNDPART:
|
||||
case FGI_BNDFREE:
|
||||
if(!staticInit) {
|
||||
@ -1191,16 +1283,15 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
}
|
||||
break;
|
||||
// off */
|
||||
/*
|
||||
case FGI_BNDPART: rhomass = BND_FILL;
|
||||
otype = ntype = CFBnd|CFBndPartslip;
|
||||
otype = ntype = CFBnd|CFBndPartslip|(OId<<24);
|
||||
break;
|
||||
case FGI_BNDFREE: rhomass = BND_FILL;
|
||||
otype = ntype = CFBnd|CFBndFreeslip;
|
||||
otype = ntype = CFBnd|CFBndFreeslip|(OId<<24);
|
||||
break;
|
||||
// off */
|
||||
case FGI_BNDNO: rhomass = BND_FILL;
|
||||
otype = ntype = CFBnd|CFBndNoslip;
|
||||
otype = ntype = CFBnd|CFBndNoslip|(OId<<24);
|
||||
break;
|
||||
case FGI_FLUID:
|
||||
otype = ntype = CFFluid;
|
||||
@ -1223,26 +1314,32 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
mObjectSpeeds[OId] = vec2L(this->mpParam->calculateLattVelocityFromRw( vec2P( (*this->mpGiObjects)[OId]->getInitialVelocity(mSimulationTime) )));
|
||||
debMsgStd("LbmFsgrSolver::initMovingObstacles",DM_MSG,"id"<<OId<<" "<<obj->getName()<<" inivel set to "<< mObjectSpeeds[OId]<<", unscaled:"<< (*this->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, mLevel[mMaxRefine].nodeSize, this->mvGeoStart, this->mvGeoEnd);
|
||||
obj->initMovingPointsAnim(sourceTime,mMOIVerticesOld, targetTime, mMOIVertices, pNormals, mLevel[mMaxRefine].nodeSize, this->mvGeoStart, this->mvGeoEnd);
|
||||
monTrafo += mMOIVerticesOld.size();
|
||||
obj->applyTransformation(sourceTime, &mMOIVerticesOld,NULL, 0, mMOIVerticesOld.size(), false );
|
||||
obj->applyTransformation(sourceTime, &mMOIVerticesOld,pNormals, 0, mMOIVerticesOld.size(), false );
|
||||
monTrafo += mMOIVertices.size();
|
||||
obj->applyTransformation(targetTime, &mMOIVertices,NULL, 0, mMOIVertices.size(), false );
|
||||
obj->applyTransformation(targetTime, &mMOIVertices,NULL /* no old normals needed */, 0, mMOIVertices.size(), false );
|
||||
} else {
|
||||
// only do transform update
|
||||
obj->getMovingPoints(mMOIVertices);
|
||||
obj->getMovingPoints(mMOIVertices,pNormals);
|
||||
mMOIVerticesOld = mMOIVertices;
|
||||
// WARNING - assumes mSimulationTime is global!?
|
||||
obj->applyTransformation(targetTime, &mMOIVertices,NULL, 0, mMOIVertices.size(), false );
|
||||
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, 0, mMOIVerticesOld.size(), false );
|
||||
obj->applyTransformation(sourceTime, &mMOIVerticesOld, NULL /* no old normals needed */, 0, mMOIVerticesOld.size(), false );
|
||||
monTrafo += mMOIVerticesOld.size();
|
||||
}
|
||||
|
||||
@ -1263,33 +1360,22 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
LbmFloat massCheck = 0.;
|
||||
int massReinits=0;
|
||||
bool fillCells = (mObjectMassMovnd[OId]<=-1.);
|
||||
LbmFloat massCScale; //, massCScalePos, massCScaleNeg;
|
||||
//if(mInitialMass>0.) {
|
||||
//massCScale = (mInitialMass-mObjectMassMovnd[OId])/mInitialMass;
|
||||
massCScale = 1.-(mObjectMassMovnd[OId]/(LbmFloat)(mMOIVertices.size()/10) );
|
||||
//massCScalePos = MIN(massCScale,1.);
|
||||
//massCScaleNeg = MAX(massCScale,1.);
|
||||
//} else {
|
||||
//massCScale = massCScalePos = massCScaleNeg = 1.;
|
||||
//}
|
||||
|
||||
// first pass - set new obs. cells
|
||||
if(active) {
|
||||
for(size_t n=0; n<mMOIVertices.size(); n++) {
|
||||
//errMsg("AAABB","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
|
||||
//errMsg("initMovingObstacles_Debug","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK);
|
||||
POS2GRID_CHECK(mMOIVertices,n);
|
||||
//if(i==30 && j==14) { errMsg("AAABB","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK<<" "); }
|
||||
//{ 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);
|
||||
}
|
||||
FORDF0 { massCheck -= QCELL(level, i,j,k, workSet, l); }
|
||||
massReinits++;
|
||||
}
|
||||
if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
|
||||
else if(RFLAG(level, i,j,k, workSet)&(CFInter)) {
|
||||
massCheck -= QCELL(level, i,j,k, workSet, dMass);
|
||||
massReinits++;
|
||||
}
|
||||
@ -1314,60 +1400,94 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
const LbmFloat ux = this->dfDvecX[l]*objvel[0];
|
||||
const LbmFloat uy = this->dfDvecY[l]*objvel[1];
|
||||
const LbmFloat uz = this->dfDvecZ[l]*objvel[2];
|
||||
LbmFloat factor = 2.0*this->dfLength[l]* 3.0 *(ux+uy+uz); // rhoTest, dont multiply by density...
|
||||
|
||||
/*LbmFloat *rhodstCell = RACPNT(level, i+dfVecX[l],j+dfVecY[l],k+dfVecZ[l],workSet);
|
||||
LbmFloat rhodst = RAC(rhodstCell,dC);
|
||||
for(int ll=1; ll<19; ll++) { rhodst += RAC(rhodstCell,ll); }
|
||||
rhodst = 1.; // DEBUG TEST! */
|
||||
|
||||
LbmFloat factor = 2. * this->dfLength[l] * 3.0 * (ux+uy+uz); //
|
||||
/* FSTEST
|
||||
*/
|
||||
if(ntype&(CFBndFreeslip|CFBndPartslip)) {
|
||||
// missing, diag mass checks...
|
||||
//if(l<=LBMDIM*2) massCheck += factor;
|
||||
|
||||
//if(l<=LBMDIM*2) factor *= 1.4142;
|
||||
//if(l>LBMDIM*2) factor = 0.;
|
||||
//else
|
||||
//factor = 0.; // TODO, optimize
|
||||
factor *= 2.0; // TODO, optimize
|
||||
} else {
|
||||
factor *= 1.2; // TODO, optimize
|
||||
}
|
||||
RAC(dstCell,l) = factor;
|
||||
massCheck += RAC(dstCell,l);
|
||||
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; }
|
||||
}
|
||||
//errMsg("AAABB","OId"<<OId<<" n"<<n<<" -> "<<PRINT_IJK" objvel"<<objvel<<" ul"<<PRINT_VEC(ux,uy,uz) );
|
||||
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)) {
|
||||
//? unused ntlVec3Gfx objvel= (mMOIVertices[n]-mMOIVerticesOld[n]);
|
||||
// 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,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;
|
||||
//LbmFloat avgrho=0.0, avgux=0.0, avguy=0.0, avguz=0.0;
|
||||
if(nbored&CFFluid) {
|
||||
//if(nbored&(CFFluid|CFInter)) {
|
||||
settype = CFInter|CFNoInterpolSrc;
|
||||
if(fillCells) rhomass = 1.0;
|
||||
else rhomass = 0.;
|
||||
rhomass = 1.5;
|
||||
if(!fillCells) rhomass = 0.;
|
||||
//settype = CFFluid|CFNoInterpolSrc; rhomass=1.; // test
|
||||
//rhomass = 1.01; // DEBUGT
|
||||
|
||||
//interpolateCellValues(level,i,j,k, workSet, avgrho,avgux,avguy,avguz);
|
||||
//LbmVec speed(avgux,avguy,avguz);
|
||||
//initVelocityCell(level, i,j,k, settype, avgrho, rhomass, speed );
|
||||
OBJVEL_CALC;
|
||||
initVelocityCell(level, i,j,k, settype, 1., rhomass, objvel );
|
||||
if(!(nbored&CFEmpty)) { settype=CFFluid|CFNoInterpolSrc; rhomass=1.; }
|
||||
//settype=CFFluid; // rhomass = 1.; objvel = LbmVec(0.); // DEBUG TEST
|
||||
|
||||
// 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);
|
||||
//objvel = LbmVec(avgux,avguy,avguz);
|
||||
//avgrho=1.;
|
||||
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) );
|
||||
|
||||
// old - fixed init
|
||||
//initVelocityCell(level, i,j,k, settype, 1., rhomass, objvel );
|
||||
massCheck += rhomass;
|
||||
}
|
||||
/*else if((nbored&CFInter)&&(fillCells)) {
|
||||
settype = CFInter|CFNoInterpolSrc; rhomass = 1.0;
|
||||
_interpolateCellValues(level,i,j,k, workSet, avgrho,avgux,avguy,avguz);
|
||||
} // */
|
||||
else {
|
||||
} else {
|
||||
settype = CFEmpty; rhomass = 0.0;
|
||||
initEmptyCell(level, i,j,k, settype, 1., rhomass );
|
||||
}
|
||||
//settype = CFBnd|CFBndNoslip; rhomass = 0.0;
|
||||
//avgux=avguy=avguz=0.0; avgrho=1.0;
|
||||
monFluids++;
|
||||
massReinits++;
|
||||
} // flag & simtime
|
||||
@ -1376,8 +1496,8 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
|
||||
// only compute mass transfer when init is done
|
||||
if(this->mInitDone) {
|
||||
errMsg("initMov","Massd test "<<obj->getName()<<" dccd massCheck="<<massCheck<<" massReinits"<<massReinits<<
|
||||
" fillCells"<<fillCells<<" massmovbnd:"<<mObjectMassMovnd[OId]<<" massCScale"<<massCScale<<" inim:"<<mInitialMass );
|
||||
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
|
||||
@ -1437,9 +1557,10 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
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, set2Flag);
|
||||
forceChangeFlag(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);
|
||||
forceChangeFlag(level, i,j,k, workSet,
|
||||
(RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
|
||||
}
|
||||
}
|
||||
} // second static inflow pass
|
||||
@ -1453,7 +1574,7 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
// 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)) {
|
||||
changeFlag(level, i,j,k, workSet, CFMbndOutflow); }
|
||||
forceChangeFlag(level, i,j,k, workSet, CFMbndOutflow); }
|
||||
continue;
|
||||
}
|
||||
monFluids++;
|
||||
@ -1481,9 +1602,10 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
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, set2Flag);
|
||||
forceChangeFlag(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);
|
||||
forceChangeFlag(level, i,j,k, workSet,
|
||||
(RFLAG(level, i,j,k, workSet)&CFNoPersistMask)|set2Flag);
|
||||
}
|
||||
}
|
||||
} // second static outflow pass
|
||||
@ -1511,6 +1633,15 @@ void LbmFsgrSolver::initMovingObstacles(bool staticInit) {
|
||||
}
|
||||
|
||||
|
||||
// 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) */
|
||||
/*****************************************************************************/
|
||||
@ -1580,6 +1711,7 @@ bool LbmFsgrSolver::initGeometryFlags() {
|
||||
if(LBMDIM==2) {
|
||||
dvec[2] = 0.0;
|
||||
iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0, 0.0, (this->mvGeoEnd[2]-this->mvGeoStart[2])*0.5 ))+(dvec*0.5);
|
||||
//if(mInit2dYZ) { SWAPYZ(mGravity); for(int lev=0; lev<=mMaxRefine; lev++){ SWAPYZ( mLevel[lev].gravity ); } }
|
||||
} else {
|
||||
iniPos =(this->mvGeoStart + ntlVec3Gfx( 0.0 ))+(dvec*0.5);
|
||||
}
|
||||
@ -1587,16 +1719,13 @@ bool LbmFsgrSolver::initGeometryFlags() {
|
||||
|
||||
// first init boundary conditions
|
||||
// invalid cells are set to empty afterwards
|
||||
#define GETPOS(i,j,k) \
|
||||
ntlVec3Gfx( iniPos[0]+ dvec[0]*(gfxReal)(i), \
|
||||
iniPos[1]+ dvec[1]*(gfxReal)(j), \
|
||||
iniPos[2]+ dvec[2]*(gfxReal)(k) )
|
||||
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;
|
||||
|
||||
const bool inside = this->geoInitCheckPointInside( GETPOS(i,j,k) , FGI_ALLBOUNDS, OId, distance);
|
||||
GETPOS(i,j,k);
|
||||
const bool inside = this->geoInitCheckPointInside( ggpos, FGI_ALLBOUNDS, OId, distance);
|
||||
if(inside) {
|
||||
pObj = (*this->mpGiObjects)[OId];
|
||||
switch( pObj->getGeoInitType() ){
|
||||
@ -1620,6 +1749,9 @@ bool LbmFsgrSolver::initGeometryFlags() {
|
||||
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;
|
||||
@ -1662,9 +1794,8 @@ bool LbmFsgrSolver::initGeometryFlags() {
|
||||
if(!(RFLAG(level, i,j,k, mLevel[level].setCurr)==CFEmpty)) continue;
|
||||
ntype = CFInvalid;
|
||||
int inits = 0;
|
||||
//if((i==1) && (j==31) && (k==48)) globGeoInitDebug=1;
|
||||
//else globGeoInitDebug=0;
|
||||
const bool inside = this->geoInitCheckPointInside( GETPOS(i,j,k) , FGI_FLUID, OId, distance);
|
||||
GETPOS(i,j,k);
|
||||
const bool inside = this->geoInitCheckPointInside( ggpos, FGI_FLUID, OId, distance);
|
||||
if(inside) {
|
||||
ntype = CFFluid;
|
||||
}
|
||||
@ -1723,23 +1854,16 @@ void LbmFsgrSolver::initFreeSurfaces() {
|
||||
// set interface cells
|
||||
FSGR_FORIJK1(mMaxRefine) {
|
||||
if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, mLevel[mMaxRefine].setCurr), CFFluid )) {
|
||||
//int initInter = 0; // check for neighboring empty cells
|
||||
FORDF1 {
|
||||
int ni=i+this->dfVecX[l], nj=j+this->dfVecY[l], nk=k+this->dfVecZ[l];
|
||||
if( TESTFLAG( 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) );
|
||||
initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
|
||||
// unnecessary? initEmptyCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill);
|
||||
initVelocityCell(mMaxRefine, ni,nj,nk, CFInter, arho, interfaceFill, LbmVec(aux,auy,auz) );
|
||||
//initEmptyCell(level, i,j,k, settype, avgrho, rhomass, speed ); */
|
||||
//initInter = 1;
|
||||
}
|
||||
}
|
||||
/*if(initInter) {
|
||||
QCELL(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr, dMass) = interfaceFill;
|
||||
RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setCurr) = RFLAG(mMaxRefine,i,j,k,mLevel[mMaxRefine].setOther) = CFInter;
|
||||
} // */
|
||||
}
|
||||
}
|
||||
|
||||
@ -1950,10 +2074,6 @@ void LbmFsgrSolver::initStandingFluidGradient() {
|
||||
if(haveStandingFluid) {
|
||||
int rhoworkSet = mLevel[lev].setCurr;
|
||||
myTime_t timestart = getTime(); // FIXME use user time here?
|
||||
LbmFloat lcsmqo;
|
||||
#if OPT3D==1
|
||||
LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
|
||||
#endif // OPT3D==true
|
||||
|
||||
GRAVLOOP {
|
||||
int i = gravIndex[0], j = gravIndex[1], k = gravIndex[2];
|
||||
@ -1980,87 +2100,21 @@ void LbmFsgrSolver::initStandingFluidGradient() {
|
||||
}
|
||||
|
||||
} // 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);
|
||||
|
||||
#if ELBEEM_PLUGIN!=1 && LBMDIM==3
|
||||
/*int lowj = 0;
|
||||
for(int k=1;k<mLevel[lev].lSizez-1;++k) {
|
||||
for(int i=1;i<mLevel[lev].lSizex-1;++i) {
|
||||
LbmFloat rho = 1.0+ (fluidHeight) * (mLevel[lev].gravity[maxGravComp])* (-3.0/1.0)*(mLevel[lev].omega);
|
||||
RFLAG(lev, i,lowj,k, rhoworkSet^1) =
|
||||
RFLAG(lev, i,lowj,k, rhoworkSet) = CFFluid;
|
||||
FORDF0 { QCELL(lev, i,lowj,k, rhoworkSet, l) = this->dfEquil[l]*rho; }
|
||||
QCELL(lev, i,lowj,k, rhoworkSet, dMass) = rho;
|
||||
} } // */
|
||||
#endif
|
||||
|
||||
int preinitSteps = (haveStandingFluid* ((mLevel[lev].lSizey+mLevel[lev].lSizez+mLevel[lev].lSizex)/3) );
|
||||
preinitSteps = (haveStandingFluid>>2); // not much use...?
|
||||
//preinitSteps = 4; // DEBUG!!!!
|
||||
//this->mInitDone = 1; // GRAVTEST
|
||||
//preinitSteps = 0;
|
||||
debMsgNnl("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
|
||||
debMsgStd("Standing fluid preinit", DM_MSG, "Performing "<<preinitSteps<<" prerelaxations ",10);
|
||||
for(int s=0; s<preinitSteps; s++) {
|
||||
int workSet = SRCS(lev); //mLevel[lev].setCurr;
|
||||
int otherSet = TSET(lev); //mLevel[lev].setOther;
|
||||
debMsgDirect(".");
|
||||
if(debugStandingPreinit) debMsgStd("Standing fluid preinit", DM_MSG, "s="<<s<<" curset="<<workSet<<" srcs"<<SRCS(lev), 10);
|
||||
LbmFloat *ccel;
|
||||
LbmFloat *tcel;
|
||||
LbmFloat m[LBM_DFNUM];
|
||||
|
||||
// grav loop not necessary here
|
||||
#define NBFLAG(l) (nbflag[(l)])
|
||||
LbmFloat rho, ux,uy,uz, usqr;
|
||||
int kstart=getForZMinBnd(), kend=getForZMaxBnd(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=0;j<mLevel[lev].lSizey-0;++j) {
|
||||
for(int i=0;i<mLevel[lev].lSizex-0;++i) {
|
||||
const CellFlagType currFlag = RFLAG(lev, i,j,k,workSet);
|
||||
if( (currFlag & (CFEmpty|CFBnd)) ) continue;
|
||||
ccel = RACPNT(lev, i,j,k,workSet);
|
||||
tcel = RACPNT(lev, i,j,k,otherSet);
|
||||
|
||||
if( (currFlag & (CFInter)) ) {
|
||||
// copy all values
|
||||
for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
|
||||
continue;
|
||||
}
|
||||
|
||||
if( (currFlag & CFNoBndFluid)) {
|
||||
OPTIMIZED_STREAMCOLLIDE;
|
||||
} else {
|
||||
FORDF1 {
|
||||
nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
|
||||
}
|
||||
DEFAULT_STREAM;
|
||||
//ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2];
|
||||
DEFAULT_COLLIDEG(mLevel[lev].gravity);
|
||||
}
|
||||
for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
|
||||
} } } // GRAVLOOP
|
||||
|
||||
mLevel[lev].setOther = mLevel[lev].setCurr;
|
||||
mLevel[lev].setCurr ^= 1;
|
||||
// in solver main cpp
|
||||
standingFluidPreinit();
|
||||
}
|
||||
//this->mInitDone = 0; // GRAVTEST
|
||||
// */
|
||||
|
||||
myTime_t timeend = getTime();
|
||||
debMsgDirect(" done, "<<getTimeString(timeend-timestart)<<" \n");
|
||||
debMsgStd("Standing fluid preinit", DM_MSG, " done, "<<getTimeString(timeend-timestart), 9);
|
||||
#undef NBFLAG
|
||||
}
|
||||
}
|
||||
@ -2150,17 +2204,15 @@ LbmFsgrSolver::interpolateCellValues(
|
||||
LbmFloat avgrho = 0.0;
|
||||
LbmFloat avgux = 0.0, avguy = 0.0, avguz = 0.0;
|
||||
LbmFloat cellcnt = 0.0;
|
||||
//LbmFloat avgnbdf[LBM_DFNUM];
|
||||
//FORDF0M { avgnbdf[m]= 0.0; }
|
||||
|
||||
for(int nbl=1; nbl< this->cDfNum ; ++nbl) {
|
||||
if( (RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFFluid) ||
|
||||
((!(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFNoInterpolSrc) ) &&
|
||||
(RFLAG_NB(level,ei,ej,ek,workSet,nbl) & CFInter) )) {
|
||||
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< this->cDfNum ; ++rl) {
|
||||
LbmFloat nbdf = QCELL_NB(level,ei,ej,ek, workSet,nbl, rl);
|
||||
//avgnbdf[rl] += nbdf;
|
||||
avgux += (this->dfDvecX[rl]*nbdf);
|
||||
avguy += (this->dfDvecY[rl]*nbdf);
|
||||
avguz += (this->dfDvecZ[rl]*nbdf);
|
||||
@ -2171,20 +2223,18 @@ LbmFsgrSolver::interpolateCellValues(
|
||||
|
||||
if(cellcnt<=0.0) {
|
||||
// no nbs? just use eq.
|
||||
//FORDF0 { QCELL(level,ei,ej,ek, workSet, l) = this->dfEquil[l]; }
|
||||
avgrho = 1.0;
|
||||
avgux = avguy = avguz = 0.0;
|
||||
//TTT mNumProblems++;
|
||||
#if ELBEEM_PLUGIN!=1
|
||||
//this->mPanic=1;
|
||||
// this can happen for worst case moving obj scenarios...
|
||||
errMsg("LbmFsgrSolver::interpolateCellValues","Cellcnt<=0.0 at "<<PRINT_VEC(ei,ej,ek)); //,SIMWORLD_GENERICERROR);
|
||||
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;
|
||||
//FORDF0M { avgnbdf[m] /= cellcnt; } // CHECK FIXME test?
|
||||
}
|
||||
|
||||
retrho = avgrho;
|
||||
|
||||
@ -285,7 +285,11 @@ void LbmSolverInterface::initGeoTree() {
|
||||
if(mpGiTree != NULL) delete mpGiTree;
|
||||
char treeFlag = (1<<(this->mLbmInitId+4));
|
||||
mpGiTree = new ntlTree(
|
||||
# if FSGR_STRICT_DEBUG!=1
|
||||
15, 8, // TREEwarning - fixed values for depth & maxtriangles here...
|
||||
# else // FSGR_STRICT_DEBUG==1
|
||||
10, 20, // reduced/slower debugging values
|
||||
# endif // FSGR_STRICT_DEBUG==1
|
||||
scene, treeFlag );
|
||||
}
|
||||
|
||||
|
||||
@ -128,6 +128,7 @@ typedef int BubbleId;
|
||||
|
||||
// above 24 is used to encode in/outflow object type
|
||||
#define CFPersistMask (0xFF000000 | CFMbndInflow | CFMbndOutflow)
|
||||
#define CFNoPersistMask (~CFPersistMask)
|
||||
|
||||
|
||||
// nk
|
||||
|
||||
@ -10,15 +10,60 @@
|
||||
#include "solver_class.h"
|
||||
#include "solver_relax.h"
|
||||
#include "particletracer.h"
|
||||
#include "loop_tools.h"
|
||||
|
||||
/*****************************************************************************/
|
||||
/*! perform a single LBM step */
|
||||
/*****************************************************************************/
|
||||
|
||||
double globdfcnt;
|
||||
double globdfavg[19];
|
||||
double globdfmax[19];
|
||||
double globdfmin[19];
|
||||
|
||||
void LbmFsgrSolver::step() {
|
||||
initLevelOmegas();
|
||||
stepMain();
|
||||
|
||||
// intlbm test
|
||||
if(0) {
|
||||
if(this->mStepCnt<5) {
|
||||
// init
|
||||
globdfcnt=0.;
|
||||
FORDF0{
|
||||
globdfavg[l] = 0.;
|
||||
globdfmax[l] = -1000.; //this->dfEquil[l];
|
||||
globdfmin[l] = 1000.; // this->dfEquil[l];
|
||||
}
|
||||
} else {
|
||||
|
||||
int workSet = mLevel[mMaxRefine].setCurr;
|
||||
for(int k= getForZMinBnd(); k< getForZMaxBnd(mMaxRefine); ++k) {
|
||||
for(int j=0;j<mLevel[mMaxRefine].lSizey-0;j++) {
|
||||
for(int i=0;i<mLevel[mMaxRefine].lSizex-0;i++) {
|
||||
//float val = 0.;
|
||||
if(RFLAG(mMaxRefine, i,j,k, workSet) & CFFluid) {
|
||||
//val = QCELL(mMaxRefine,i,j,k, workSet,dFfrac);
|
||||
FORDF0{
|
||||
const double df = (double)QCELL(mMaxRefine,i,j,k, workSet,l);
|
||||
globdfavg[l] += df;
|
||||
if(df>globdfmax[l]) globdfmax[l] = df;
|
||||
//if(df<=0.01) { errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l); }
|
||||
if(df<globdfmin[l]) globdfmin[l] = df;
|
||||
//errMsg("GLOBDFERR"," at "<<PRINT_IJK<<" "<<l<<" "<<df<<" "<<globdfmin[l]);
|
||||
}
|
||||
globdfcnt+=1.;
|
||||
}
|
||||
}
|
||||
}
|
||||
}
|
||||
if(this->mStepCnt%10==0) {
|
||||
FORDF0{ errMsg("GLOBDF","l="<<l<<" avg="<<(globdfavg[l]/globdfcnt)<<" max="<<globdfmax[l]<<" min="<<globdfmin[l]<<" "<<globdfcnt); }
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
// intlbm test */
|
||||
}
|
||||
|
||||
void LbmFsgrSolver::stepMain()
|
||||
@ -189,6 +234,7 @@ void LbmFsgrSolver::stepMain()
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
if((mUseTestdata)&&(this->mInitDone)) { handleTestdata(); }
|
||||
testXchng();
|
||||
#endif
|
||||
|
||||
// advance positions with current grid
|
||||
@ -213,8 +259,13 @@ void LbmFsgrSolver::stepMain()
|
||||
|
||||
// output total step time
|
||||
timeend = getTime();
|
||||
double mdelta = (lastMass-mCurrentMass);
|
||||
if(ABS(mdelta)<1e-12) mdelta=0.;
|
||||
debMsgStd("LbmFsgrSolver::stepMain",DM_MSG,"step:"<<this->mStepCnt
|
||||
<<": dccd="<< mCurrentMass<<",d"<<(lastMass-mCurrentMass)<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
|
||||
<<": dccd="<< mCurrentMass
|
||||
<<",d"<<mdelta
|
||||
<<",ds"<<(mCurrentMass-mObjectMassMovnd[1])
|
||||
<<"/"<<mCurrentVolume<<"(fix="<<this->mFixMass<<",ini="<<mInitialMass<<"), "
|
||||
<<" totst:"<<getTimeString(timeend-timestart), 3);
|
||||
// nicer output
|
||||
debMsgDirect(std::endl);
|
||||
@ -248,7 +299,7 @@ void LbmFsgrSolver::fineAdvance()
|
||||
// time adaptivity
|
||||
this->mpParam->setSimulationMaxSpeed( sqrt(mMaxVlen / 1.5) );
|
||||
//if(mStartSymm) { checkSymmetry("step2"); } // DEBUG
|
||||
if(!this->mSilent){ errMsg("fineAdvance"," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps) ); }
|
||||
if(!this->mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," stepped from "<<mLevel[mMaxRefine].setCurr<<" to "<<mLevel[mMaxRefine].setOther<<" step"<< (mLevel[mMaxRefine].lsteps), 3 ); }
|
||||
|
||||
// update other set
|
||||
mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
|
||||
@ -257,62 +308,117 @@ void LbmFsgrSolver::fineAdvance()
|
||||
|
||||
// flag init... (work on current set, to simplify flag checks)
|
||||
reinitFlags( mLevel[mMaxRefine].setCurr );
|
||||
if(!this->mSilent){ errMsg("fineAdvance"," flags reinit on set "<< mLevel[mMaxRefine].setCurr ); }
|
||||
if(!this->mSilent){ debMsgStd("fineAdvance",DM_NOTIFY," flags reinit on set "<< mLevel[mMaxRefine].setCurr, 3 ); }
|
||||
|
||||
// DEBUG VEL CHECK
|
||||
if(0) {
|
||||
int lev = mMaxRefine;
|
||||
int workSet = mLevel[lev].setCurr;
|
||||
int mi=0,mj=0,mk=0;
|
||||
LbmFloat compRho=0.;
|
||||
LbmFloat compUx=0.;
|
||||
LbmFloat compUy=0.;
|
||||
LbmFloat compUz=0.;
|
||||
LbmFloat maxUlen=0.;
|
||||
LbmVec maxU(0.);
|
||||
LbmFloat maxRho=-100.;
|
||||
int ri=0,rj=0,rk=0;
|
||||
|
||||
FSGR_FORIJK1(lev) {
|
||||
if( (RFLAG(lev,i,j,k, workSet) & (CFFluid| CFInter| CFGrFromCoarse| CFGrFromFine| CFGrNorm)) ) {
|
||||
compRho=QCELL(lev, i,j,k,workSet, dC);
|
||||
compUx = compUy = compUz = 0.0;
|
||||
for(int l=1; l<this->cDfNum; l++) {
|
||||
LbmFloat df = QCELL(lev, i,j,k,workSet, l);
|
||||
compRho += df;
|
||||
compUx += (this->dfDvecX[l]*df);
|
||||
compUy += (this->dfDvecY[l]*df);
|
||||
compUz += (this->dfDvecZ[l]*df);
|
||||
}
|
||||
LbmVec u(compUx,compUy,compUz);
|
||||
LbmFloat nu = norm(u);
|
||||
if(nu>maxUlen) {
|
||||
maxU = u;
|
||||
maxUlen = nu;
|
||||
mi=i; mj=j; mk=k;
|
||||
}
|
||||
if(compRho>maxRho) {
|
||||
maxRho=compRho;
|
||||
ri=i; rj=j; rk=k;
|
||||
}
|
||||
} else {
|
||||
continue;
|
||||
}
|
||||
}
|
||||
|
||||
errMsg("MAXVELCHECK"," at "<<PRINT_VEC(mi,mj,mk)<<" norm:"<<maxUlen<<" u:"<<maxU);
|
||||
errMsg("MAXRHOCHECK"," at "<<PRINT_VEC(ri,rj,rk)<<" rho:"<<maxRho);
|
||||
printLbmCell(lev, 30,36,23, -1);
|
||||
} // DEBUG VEL CHECK
|
||||
|
||||
}
|
||||
|
||||
|
||||
/*****************************************************************************/
|
||||
//! fine step function
|
||||
/*****************************************************************************/
|
||||
|
||||
// fine step defines
|
||||
|
||||
// access to own dfs during step (may be changed to local array)
|
||||
#define MYDF(l) RAC(ccel, l)
|
||||
|
||||
// drop model definitions
|
||||
#define RWVEL_THRESH 1.5
|
||||
#define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
|
||||
|
||||
#if LBMDIM==3
|
||||
// normal
|
||||
#define SLOWDOWNREGION (this->mSizez/4)
|
||||
#else // LBMDIM==2
|
||||
// off
|
||||
#define SLOWDOWNREGION 10
|
||||
#endif // LBMDIM==2
|
||||
#define P_LCSMQO 0.01
|
||||
|
||||
/*****************************************************************************/
|
||||
//! fine step function
|
||||
/*****************************************************************************/
|
||||
void
|
||||
LbmFsgrSolver::mainLoop(int lev)
|
||||
{
|
||||
// loops over _only inner_ cells -----------------------------------------------------------------------------------
|
||||
LbmFloat calcCurrentMass = 0.0;
|
||||
LbmFloat calcCurrentVolume = 0.0;
|
||||
int calcCellsFilled = this->mNumFilledCells;
|
||||
int calcCellsEmptied = this->mNumEmptiedCells;
|
||||
int calcNumUsedCells = this->mNumUsedCells;
|
||||
|
||||
// slow surf regions smooth (if below)
|
||||
const LbmFloat smoothStrength = 0.0; //0.01;
|
||||
const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
|
||||
const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
|
||||
|
||||
const int cutMin = 1;
|
||||
const int cutConst = mCutoff+2;
|
||||
|
||||
|
||||
# if LBM_INCLUDE_TESTSOLVERS==1
|
||||
// 3d region off... quit
|
||||
if((mUseTestdata)&&(mpTest->mDebugvalue1>0.0)) { return; }
|
||||
#endif // ELBEEM_PLUGIN!=1
|
||||
//printLbmCell(lev, 6,6,16, mLevel[lev].setCurr ); // DEBUG
|
||||
# endif // ELBEEM_PLUGIN!=1
|
||||
|
||||
// main loop region
|
||||
const bool doReduce = true;
|
||||
const int gridLoopBound=1;
|
||||
GRID_REGION_INIT();
|
||||
#if PARALLEL==1
|
||||
#include "paraloop.h"
|
||||
#include "paraloopstart.h"
|
||||
GRID_REGION_START();
|
||||
#else // PARALLEL==1
|
||||
{ // main loop region
|
||||
int kstart=getForZMin1(), kend=getForZMax1(mMaxRefine);
|
||||
//{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIMcDimension<<" mlsz:"<< mLevel[mMaxRefine].lSizez<<" zmax1:"<<getForZMax1(mMaxRefine) ); } // DEBUG
|
||||
#define PERFORM_USQRMAXCHECK USQRMAXCHECK(usqr,ux,uy,uz, mMaxVlen, mMxvx,mMxvy,mMxvz);
|
||||
#define LIST_EMPTY(x) mListEmpty.push_back( x );
|
||||
#define LIST_FULL(x) mListFull.push_back( x );
|
||||
GRID_REGION_START();
|
||||
#endif // PARALLEL==1
|
||||
|
||||
// local to loop
|
||||
// local to main
|
||||
CellFlagType nbflag[LBM_DFNUM];
|
||||
LbmFloat *ccel = NULL;
|
||||
LbmFloat *tcel = NULL;
|
||||
int oldFlag;
|
||||
int newFlag;
|
||||
int nbored;
|
||||
int oldFlag, newFlag, nbored;
|
||||
LbmFloat m[LBM_DFNUM];
|
||||
LbmFloat rho, ux, uy, uz, tmp, usqr;
|
||||
LbmFloat mass, change, lcsmqo=0.0;
|
||||
usqr = tmp = 0.0;
|
||||
#if OPT3D==1
|
||||
LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
|
||||
#endif // OPT3D==true
|
||||
|
||||
// smago vars
|
||||
LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
|
||||
|
||||
// ifempty cell conversion flags
|
||||
bool iffilled, ifemptied;
|
||||
@ -320,72 +426,22 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
int recons[LBM_DFNUM]; // reconstruct this DF?
|
||||
int numRecons; // how many are reconstructed?
|
||||
|
||||
// slow surf regions smooth (if below)
|
||||
const LbmFloat smoothStrength = 0.0; //0.01;
|
||||
const LbmFloat sssUsqrLimit = 1.5 * 0.03*0.03;
|
||||
const LbmFloat sssUsqrLimitInv = 1.0/sssUsqrLimit;
|
||||
|
||||
CellFlagType *pFlagSrc;
|
||||
CellFlagType *pFlagDst;
|
||||
pFlagSrc = &RFLAG(lev, 0,1, kstart,SRCS(lev)); // omp
|
||||
pFlagDst = &RFLAG(lev, 0,1, kstart,TSET(lev)); // omp
|
||||
ccel = RACPNT(lev, 0,1, kstart ,SRCS(lev)); // omp
|
||||
tcel = RACPNT(lev, 0,1, kstart ,TSET(lev)); // omp
|
||||
//CellFlagType *pFlagTar = NULL;
|
||||
int pFlagTarOff;
|
||||
if(mLevel[lev].setOther==1) pFlagTarOff = mLevel[lev].lOffsz;
|
||||
else pFlagTarOff = -mLevel[lev].lOffsz;
|
||||
#define ADVANCE_POINTERS(p) \
|
||||
ccel += (QCELLSTEP*(p)); \
|
||||
tcel += (QCELLSTEP*(p)); \
|
||||
pFlagSrc+= (p); \
|
||||
pFlagDst+= (p); \
|
||||
i+= (p);
|
||||
LbmFloat mass=0., change=0., lcsmqo=0.;
|
||||
rho= ux= uy= uz= usqr= tmp= 0.;
|
||||
lcsmqadd = lcsmomega = 0.;
|
||||
FORDF0{ lcsmeq[l] = 0.; }
|
||||
|
||||
// ---
|
||||
// now stream etc.
|
||||
|
||||
//{ errMsg("LbmFsgrSolver::mainLoop","Err MAINADVANCE0 ks:"<< kstart<<" ke:"<<kend<<" dim:"<<LBMDIM<<" mlsz:"<<mLevel[mMaxRefine].lSizez ); } // DEBUG
|
||||
|
||||
// use //template functions for 2D/3D
|
||||
#if COMPRESSGRIDS==0
|
||||
for(int k=kstart;k<kend;++k) {
|
||||
for(int j=1;j<mLevel[lev].lSizey-1;++j) {
|
||||
for(int i=0;i<mLevel[lev].lSizex-2; ) {
|
||||
#else // COMPRESSGRIDS==0
|
||||
int kdir = 1; // COMPRT ON
|
||||
if(mLevel[mMaxRefine].setCurr==1) {
|
||||
kdir = -1;
|
||||
int temp = kend;
|
||||
kend = kstart-1;
|
||||
kstart = temp-1;
|
||||
} // COMPRT
|
||||
|
||||
#if PARALLEL==0
|
||||
//const int id = 0, Nthrds = 1;
|
||||
const int jstart = 0;
|
||||
const int jend = mLevel[mMaxRefine].lSizey;
|
||||
//#endif // PARALLEL==1
|
||||
#else // PARALLEL==1
|
||||
PARA_INITIALIZE();
|
||||
errMsg("LbmFsgrSolver::mainLoop","id="<<id<<" js="<<jstart<<" je="<<jend<<" jdir="<<(1) ); // debug
|
||||
#endif // PARALLEL==1
|
||||
GRID_LOOP_START();
|
||||
// loop start
|
||||
// stream from current set to other, then collide and store
|
||||
//errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
|
||||
|
||||
for(int k=kstart;k!=kend;k+=kdir) {
|
||||
|
||||
//errMsg("LbmFsgrSolver::mainLoop","k="<<k<<" ks="<<kstart<<" ke="<<kend<<" kdir="<<kdir<<" x*y="<<mLevel[mMaxRefine].lSizex*mLevel[mMaxRefine].lSizey*dTotalNum ); // debug
|
||||
pFlagSrc = &RFLAG(lev, 0, jstart, k, SRCS(lev)); // omp test // COMPRT ON
|
||||
pFlagDst = &RFLAG(lev, 0, jstart, k, TSET(lev)); // omp test
|
||||
ccel = RACPNT(lev, 0, jstart, k, SRCS(lev)); // omp test
|
||||
tcel = RACPNT(lev, 0, jstart, k, TSET(lev)); // omp test // COMPRT ON
|
||||
|
||||
//for(int j=1;j<mLevel[lev].lSizey-1;++j) {
|
||||
for(int j=jstart;j!=jend;++j) {
|
||||
for(int i=0;i<mLevel[lev].lSizex-2; ) {
|
||||
#endif // COMPRESSGRIDS==0
|
||||
|
||||
ADVANCE_POINTERS(1);
|
||||
#if FSGR_STRICT_DEBUG==1
|
||||
# if FSGR_STRICT_DEBUG==1
|
||||
// safety pointer check
|
||||
rho = ux = uy = uz = tmp = usqr = -100.0; // DEBUG
|
||||
if( (&RFLAG(lev, i,j,k,mLevel[lev].setCurr) != pFlagSrc) ||
|
||||
(&RFLAG(lev, i,j,k,mLevel[lev].setOther) != pFlagDst) ) {
|
||||
@ -405,13 +461,11 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
);
|
||||
CAUSE_PANIC;
|
||||
}
|
||||
#endif
|
||||
# endif
|
||||
oldFlag = *pFlagSrc;
|
||||
//DEBUG if(LBMDIM==2) { errMsg("LbmFsgrSolver::mainLoop","Err flagp "<<PRINT_IJK<<"="<< RFLAG(lev, i,j,k,mLevel[lev].setCurr)<<" "); }
|
||||
// stream from current set to other, then collide and store
|
||||
|
||||
// old INTCFCOARSETEST==1
|
||||
if( (oldFlag & (CFGrFromCoarse)) ) { // interpolateFineFromCoarse test!
|
||||
if( (oldFlag & (CFGrFromCoarse)) ) {
|
||||
if(( this->mStepCnt & (1<<(mMaxRefine-lev)) ) ==1) {
|
||||
FORDF0 { RAC(tcel,l) = RAC(ccel,l); }
|
||||
} else {
|
||||
@ -419,7 +473,7 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
calcNumUsedCells++;
|
||||
}
|
||||
continue; // interpolateFineFromCoarse test!
|
||||
} // interpolateFineFromCoarse test! old INTCFCOARSETEST==1
|
||||
} // interpolateFineFromCoarse test!
|
||||
|
||||
if(oldFlag & (CFMbndInflow)) {
|
||||
// fluid & if are ok, fill if later on
|
||||
@ -434,9 +488,10 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
RAC(tcel, dMass) = RAC(tcel, dFfrac) = iniRho;
|
||||
RAC(tcel, dFlux) = FLUX_INIT;
|
||||
changeFlag(lev, i,j,k, TSET(lev), CFInter);
|
||||
calcCurrentMass += iniRho; calcCurrentVolume += 1.0; calcNumUsedCells++;
|
||||
calcCurrentMass += iniRho;
|
||||
calcCurrentVolume += 1.0;
|
||||
calcNumUsedCells++;
|
||||
mInitialMass += iniRho;
|
||||
//errMsg("INFLOW_DEBUG","ini at "<<PRINT_IJK<<" v="<<vel<<" inirho="<<iniRho);
|
||||
// dont treat cell until next step
|
||||
continue;
|
||||
}
|
||||
@ -456,11 +511,6 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
// same as ifemptied for if below
|
||||
LbmPoint oemptyp; oemptyp.flag = 0;
|
||||
oemptyp.x = i; oemptyp.y = j; oemptyp.z = k;
|
||||
#if PARALLEL==1
|
||||
//calcListEmpty[id].push_back( oemptyp );
|
||||
#else // PARALLEL==1
|
||||
//mListEmpty.push_back( oemptyp );
|
||||
#endif // PARALLEL==1
|
||||
LIST_EMPTY(oemptyp);
|
||||
calcCellsEmptied++;
|
||||
continue;
|
||||
@ -597,7 +647,7 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
// for fluid cells - just the f_i difference from streaming to empty cells ----
|
||||
numRecons = 0;
|
||||
bool onlyBndnb = ((!(oldFlag&CFNoBndFluid))&&(oldFlag&CFNoNbFluid)&&(nbored&CFBndNoslip));
|
||||
//onlyBndnb = false; // DEBUG test off
|
||||
//onlyBndnb = false; // DEBUG test off
|
||||
|
||||
FORDF1 { // dfl loop
|
||||
recons[l] = 0;
|
||||
@ -690,13 +740,13 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
if(nbflag[dN] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
|
||||
if(nbflag[dS] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsx*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
|
||||
ny = 0.5* (nv2-nv1);
|
||||
#if LBMDIM==3
|
||||
# if LBMDIM==3
|
||||
if(nbflag[dT] &(CFFluid|CFInter)){ nv1 = RAC((ccel+(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv1 = 0.0;
|
||||
if(nbflag[dB] &(CFFluid|CFInter)){ nv2 = RAC((ccel-(mLevel[lev].lOffsy*QCELLSTEP)),dFfrac); } else nv2 = 0.0;
|
||||
nz = 0.5* (nv2-nv1);
|
||||
#else // LBMDIM==3
|
||||
# else // LBMDIM==3
|
||||
nz = 0.0;
|
||||
#endif // LBMDIM==3
|
||||
# endif // LBMDIM==3
|
||||
|
||||
if( (ABS(nx)+ABS(ny)+ABS(nz)) > LBM_EPSILON) {
|
||||
// normal ok and usable...
|
||||
@ -712,16 +762,27 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
// calculate macroscopic cell values
|
||||
LbmFloat oldUx, oldUy, oldUz;
|
||||
LbmFloat oldRho; // OLD rho = ccel->rho;
|
||||
#if OPT3D==0
|
||||
oldRho=RAC(ccel,0);
|
||||
oldUx = oldUy = oldUz = 0.0;
|
||||
for(int l=1; l<this->cDfNum; l++) {
|
||||
oldRho += RAC(ccel,l);
|
||||
oldUx += (this->dfDvecX[l]*RAC(ccel,l));
|
||||
oldUy += (this->dfDvecY[l]*RAC(ccel,l));
|
||||
oldUz += (this->dfDvecZ[l]*RAC(ccel,l));
|
||||
# define REFERENCE_PRESSURE 1.0 // always atmosphere...
|
||||
# if OPT3D==0
|
||||
oldRho=RAC(ccel,0);
|
||||
oldUx = oldUy = oldUz = 0.0;
|
||||
for(int l=1; l<this->cDfNum; l++) {
|
||||
oldRho += RAC(ccel,l);
|
||||
oldUx += (this->dfDvecX[l]*RAC(ccel,l));
|
||||
oldUy += (this->dfDvecY[l]*RAC(ccel,l));
|
||||
oldUz += (this->dfDvecZ[l]*RAC(ccel,l));
|
||||
}
|
||||
// reconstruct dist funcs from empty cells
|
||||
FORDF1 {
|
||||
if(recons[ l ]) {
|
||||
m[ this->dfInv[l] ] =
|
||||
this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) +
|
||||
this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz)
|
||||
- MYDF( l );
|
||||
}
|
||||
#else // OPT3D==0
|
||||
}
|
||||
usqr = 1.5 * (oldUx*oldUx + oldUy*oldUy + oldUz*oldUz); // needed later on
|
||||
# else // OPT3D==0
|
||||
oldRho = + RAC(ccel,dC) + RAC(ccel,dN )
|
||||
+ RAC(ccel,dS ) + RAC(ccel,dE )
|
||||
+ RAC(ccel,dW ) + RAC(ccel,dT )
|
||||
@ -733,39 +794,25 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
+ RAC(ccel,dEB) + RAC(ccel,dWT)
|
||||
+ RAC(ccel,dWB);
|
||||
|
||||
oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
|
||||
oldUx = + RAC(ccel,dE) - RAC(ccel,dW)
|
||||
+ RAC(ccel,dNE) - RAC(ccel,dNW)
|
||||
+ RAC(ccel,dSE) - RAC(ccel,dSW)
|
||||
+ RAC(ccel,dET) + RAC(ccel,dEB)
|
||||
- RAC(ccel,dWT) - RAC(ccel,dWB);
|
||||
|
||||
oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
|
||||
oldUy = + RAC(ccel,dN) - RAC(ccel,dS)
|
||||
+ RAC(ccel,dNE) + RAC(ccel,dNW)
|
||||
- RAC(ccel,dSE) - RAC(ccel,dSW)
|
||||
+ RAC(ccel,dNT) + RAC(ccel,dNB)
|
||||
- RAC(ccel,dST) - RAC(ccel,dSB);
|
||||
|
||||
oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
|
||||
oldUz = + RAC(ccel,dT) - RAC(ccel,dB)
|
||||
+ RAC(ccel,dNT) - RAC(ccel,dNB)
|
||||
+ RAC(ccel,dST) - RAC(ccel,dSB)
|
||||
+ RAC(ccel,dET) - RAC(ccel,dEB)
|
||||
+ RAC(ccel,dWT) - RAC(ccel,dWB);
|
||||
#endif
|
||||
|
||||
// now reconstruction
|
||||
#define REFERENCE_PRESSURE 1.0 // always atmosphere...
|
||||
#if OPT3D==0
|
||||
// NOW - construct dist funcs from empty cells
|
||||
FORDF1 {
|
||||
if(recons[ l ]) {
|
||||
m[ this->dfInv[l] ] =
|
||||
this->getCollideEq(l, REFERENCE_PRESSURE, oldUx,oldUy,oldUz) +
|
||||
this->getCollideEq(this->dfInv[l], REFERENCE_PRESSURE, oldUx,oldUy,oldUz)
|
||||
- MYDF( l );
|
||||
}
|
||||
}
|
||||
usqr = 1.5 * (oldUx*oldUx + oldUy*oldUy + oldUz*oldUz); // needed later on
|
||||
#else
|
||||
ux=oldUx, uy=oldUy, uz=oldUz; // no local vars, only for usqr
|
||||
rho = REFERENCE_PRESSURE;
|
||||
usqr = 1.5 * (ux*ux + uy*uy + uz*uz); // needed later on
|
||||
@ -787,7 +834,7 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
if(recons[dEB]) { m[dWT] = EQEB + EQWT - MYDF(dEB); }
|
||||
if(recons[dWT]) { m[dEB] = EQWT + EQEB - MYDF(dWT); }
|
||||
if(recons[dWB]) { m[dET] = EQWB + EQET - MYDF(dWB); }
|
||||
#endif
|
||||
# endif
|
||||
|
||||
|
||||
// inflow bc handling
|
||||
@ -835,9 +882,6 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
&& (this->mPartGenProb>0.0)) {
|
||||
bool doAdd = true;
|
||||
bool bndOk=true;
|
||||
//if( (i<1+mCutoff)||(i>this->mSizex-1-mCutoff)||
|
||||
//(j<1+mCutoff)||(j>this->mSizey-1-mCutoff)||
|
||||
//(k<1+mCutoff)||(k>this->mSizez-1-mCutoff) ) { bndOk=false; }
|
||||
if( (i<cutMin)||(i>this->mSizex-cutMin)||
|
||||
(j<cutMin)||(j>this->mSizey-cutMin)||
|
||||
(k<cutMin)||(k>this->mSizez-cutMin) ) { bndOk=false; }
|
||||
@ -853,27 +897,14 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
LbmFloat prob = (rand()/(RAND_MAX+1.0));
|
||||
LbmFloat basethresh = this->mPartGenProb*lcsmqo*rl;
|
||||
|
||||
// reduce probability in outer region
|
||||
//if( (i<cutConst)||(i>this->mSizex-cutConst) ){ prob *= 0.5; }
|
||||
//if( (j<cutConst)||(j>this->mSizey-cutConst) ){ prob *= 0.5; }
|
||||
//if( (k<cutConst)||(k>this->mSizez-cutConst) ){ prob *= 0.5; }
|
||||
// NEW TEST 040627
|
||||
//if( (i<cutConst)||(i>this->mSizex-cutConst) ){ doAdd=false; }
|
||||
//if( (j<cutConst)||(j>this->mSizey-cutConst) ){ doAdd=false; }
|
||||
//if( (k<cutConst)||(k>this->mSizez-cutConst) ){ doAdd=false; }
|
||||
// reduce probability in outer region?
|
||||
const int pibord = mLevel[mMaxRefine].lSizex/2-cutConst;
|
||||
const int pjbord = mLevel[mMaxRefine].lSizey/2-cutConst;
|
||||
LbmFloat pifac = 1.-(LbmFloat)(ABS(i-pibord)) / (LbmFloat)(pibord);
|
||||
LbmFloat pjfac = 1.-(LbmFloat)(ABS(j-pjbord)) / (LbmFloat)(pjbord);
|
||||
if(pifac<0.) pifac=0.;
|
||||
if(pjfac<0.) pjfac=0.;
|
||||
//const LbmFloat pkfac = 1.-(LbmFloat)(ABS(k-mLevel[mMaxRefine].lSizez/2))/(LbmFloat)(mLevel[mMaxRefine].lSizez/2);
|
||||
//errMsg("PROBTTT"," at "<<PRINT_IJK<<" prob"<<prob<<" pifac"<<pifac<<" pjfac"<<pjfac<<" "<<(basethresh*rl*pifac*pjfac));
|
||||
//prob *= pifac*pjfac; //*pkfac;
|
||||
|
||||
//#define RWVEL_THRESH 1.0
|
||||
#define RWVEL_THRESH 1.5
|
||||
#define RWVEL_WINDTHRESH (RWVEL_THRESH*0.5)
|
||||
//if( (prob< (basethresh*rl)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
|
||||
if( (prob< (basethresh*rl*pifac*pjfac)) && (lcsmqo>0.0095) && (rl>RWVEL_THRESH) ) {
|
||||
// add
|
||||
@ -881,24 +912,12 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
doAdd = false; // dont...
|
||||
}
|
||||
|
||||
//#define SLOWDOWNREGION (2*mCutoff)
|
||||
#if LBMDIM==3
|
||||
// normal
|
||||
#define SLOWDOWNREGION (this->mSizez/4)
|
||||
#else // LBMDIM==2
|
||||
// off
|
||||
#define SLOWDOWNREGION 10
|
||||
#endif // LBMDIM==2
|
||||
#define P_LCSMQO 0.01
|
||||
|
||||
// "wind" disturbance
|
||||
// use realworld relative velocity here instead?
|
||||
if( (doAdd &&
|
||||
((rl>RWVEL_WINDTHRESH) && (lcsmqo<P_LCSMQO)) )// normal checks
|
||||
||(k>this->mSizez-SLOWDOWNREGION) ) {
|
||||
LbmFloat nuz = uz;
|
||||
//? LbmFloat jdf; // = 0.05 * (rand()/(RAND_MAX+1.0));
|
||||
//? if(rl>RWVEL_WINDTHRESH) jdf *= (rl-RWVEL_WINDTHRESH);
|
||||
if(k>this->mSizez-SLOWDOWNREGION) {
|
||||
// special case
|
||||
LbmFloat zfac = (LbmFloat)( k-(this->mSizez-SLOWDOWNREGION) );
|
||||
@ -917,7 +936,6 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
const LbmFloat add = this->dfLength[l]*(-ux*this->dfDvecX[l]-uy*this->dfDvecY[l]-nuz*this->dfDvecZ[l])*jdf;
|
||||
RAC(tcel,l) += add; }
|
||||
}
|
||||
//errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<PRINT_VEC(ux,uy,nuz)<<" rwv"<<PRINT_VEC(rux,ruy,ruz) );
|
||||
//errMsg("TOPDOWNCORR"," jdf:"<<jdf<<" rl"<<rl<<" vel "<<norm(LbmVec(ux,uy,nuz))<<" rwv"<<norm(LbmVec(rux,ruy,ruz)) );
|
||||
} // wind disturbance
|
||||
|
||||
@ -972,17 +990,13 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
mpParticles->getLast()->setType(type);
|
||||
//if((s%3)==2) mpParticles->getLast()->setType(PART_FLOAT);
|
||||
mpParticles->getLast()->setSize(size);
|
||||
//errMsg("NEWPART"," at "<<PRINT_IJK<<" u="<<PRINT_VEC(ux,uy,uz) <<" RWu="<<PRINT_VEC(rux,ruy,ruz)<<" add"<<doAdd<<" pvel="<<pv );
|
||||
//errMsg("NEWPART"," at "<<PRINT_IJK<<" u="<<norm(LbmVec(ux,uy,uz)) <<" RWu="<<norm(LbmVec(rux,ruy,ruz))<<" add"<<doAdd<<" pvel="<<norm(pv) );
|
||||
//if(!bndOk){ mass -= val2fac*size*0.02; } // FORCEDISSOLVE
|
||||
//const LbmFloat val2fac = 1.0; //? TODO N test? /(this->mPartGenProb); // FIXME remove factor!
|
||||
//mass -= val2fac*size*0.0015; // NTEST!
|
||||
//mass -= val2fac*size*0.001; // NTEST!
|
||||
//mass -= size*0.001; // NTEST!
|
||||
mass -= size*0.0020; // NTEST!
|
||||
#if LBMDIM==2
|
||||
# if LBMDIM==2
|
||||
mpParticles->getLast()->setVel(pv[0],pv[1],0.0);
|
||||
mpParticles->getLast()->setPos(ntlVec3Gfx(srci,srcj,0.5));
|
||||
#endif // LBMDIM==2
|
||||
# endif // LBMDIM==2
|
||||
//errMsg("PIT","NEW pit"<<mpParticles->getLast()->getId()<<" pos:"<< mpParticles->getLast()->getPos()<<" status:"<<convertFlags2String(mpParticles->getLast()->getFlags())<<" vel:"<< mpParticles->getLast()->getVel() );
|
||||
} // multiple parts
|
||||
} // doAdd
|
||||
@ -1004,8 +1018,8 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
}
|
||||
|
||||
// looks much nicer... LISTTRICK
|
||||
#if FSGR_LISTTRICK==1
|
||||
if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
|
||||
# if FSGR_LISTTRICK==1
|
||||
//if((oldFlag&CFNoNbEmpty)&&(newFlag&CFNoNbEmpty)) { TEST_IF_CHECK; }
|
||||
if(newFlag&CFNoBndFluid) { // test NEW TEST
|
||||
if(!iffilled) {
|
||||
// remove cells independent from amount of change...
|
||||
@ -1023,7 +1037,7 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
}
|
||||
}
|
||||
} // nobndfluid only */
|
||||
#endif
|
||||
# endif
|
||||
//iffilled = ifemptied = 0; // DEBUG!!!!!!!!!!!!!!!
|
||||
|
||||
|
||||
@ -1032,11 +1046,6 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
LbmPoint filledp; filledp.flag=0;
|
||||
if(!(newFlag&CFNoBndFluid)) filledp.flag |= 1; // NEWSURFT
|
||||
filledp.x = i; filledp.y = j; filledp.z = k;
|
||||
#if PARALLEL==1
|
||||
//calcListFull[id].push_back( filledp );
|
||||
#else // PARALLEL==1
|
||||
//mListFull.push_back( filledp );
|
||||
#endif // PARALLEL==1
|
||||
LIST_FULL(filledp);
|
||||
//this->mNumFilledCells++; // DEBUG
|
||||
calcCellsFilled++;
|
||||
@ -1045,11 +1054,6 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
LbmPoint emptyp; emptyp.flag=0;
|
||||
if(!(newFlag&CFNoBndFluid)) emptyp.flag |= 1; // NEWSURFT
|
||||
emptyp.x = i; emptyp.y = j; emptyp.z = k;
|
||||
#if PARALLEL==1
|
||||
//calcListEmpty[id].push_back( emptyp );
|
||||
#else // PARALLEL==1
|
||||
//mListEmpty.push_back( emptyp );
|
||||
#endif // PARALLEL==1
|
||||
LIST_EMPTY(emptyp);
|
||||
//this->mNumEmptiedCells++; // DEBUG
|
||||
calcCellsEmptied++;
|
||||
@ -1087,40 +1091,135 @@ LbmFsgrSolver::mainLoop(int lev)
|
||||
calcCurrentVolume += RAC(tcel,dFfrac);
|
||||
|
||||
// interface cell handling done...
|
||||
} // i
|
||||
int i=0; //dummy
|
||||
ADVANCE_POINTERS(2);
|
||||
} // j
|
||||
|
||||
#if COMPRESSGRIDS==1
|
||||
#if PARALLEL==1
|
||||
//frintf(stderr," (id=%d k=%d) ",id,k);
|
||||
# pragma omp barrier
|
||||
//#if PARALLEL==1
|
||||
//#include "paraloopend.h"
|
||||
//GRID_REGION_END();
|
||||
//#else // PARALLEL==1
|
||||
//GRID_LOOPREG_END();
|
||||
//GRID_REGION_END();
|
||||
//#endif // PARALLEL==1
|
||||
#if PARALLEL!=1
|
||||
GRID_LOOPREG_END();
|
||||
#else // PARALLEL==1
|
||||
#include "paraloopend.h" // = GRID_LOOPREG_END();
|
||||
#endif // PARALLEL==1
|
||||
#else // COMPRESSGRIDS==1
|
||||
int i=0; //dummy
|
||||
ADVANCE_POINTERS(mLevel[lev].lSizex*2);
|
||||
#endif // COMPRESSGRIDS==1
|
||||
} // all cell loop k,j,i
|
||||
|
||||
} // main loop region
|
||||
|
||||
// write vars from parallel computations to class
|
||||
//errMsg("DFINI"," maxr l"<<mMaxRefine<<" cm="<<calcCurrentMass<<" cv="<<calcCurrentVolume );
|
||||
// write vars from computations to class
|
||||
mLevel[lev].lmass = calcCurrentMass;
|
||||
mLevel[lev].lvolume = calcCurrentVolume;
|
||||
this->mNumFilledCells = calcCellsFilled;
|
||||
this->mNumEmptiedCells = calcCellsEmptied;
|
||||
this->mNumUsedCells = calcNumUsedCells;
|
||||
#if PARALLEL==1
|
||||
PARA_FINISH();
|
||||
#endif // PARALLEL==1
|
||||
|
||||
// check other vars...?
|
||||
}
|
||||
|
||||
|
||||
|
||||
void
|
||||
LbmFsgrSolver::preinitGrids()
|
||||
{
|
||||
const int lev = mMaxRefine;
|
||||
const bool doReduce = false;
|
||||
const int gridLoopBound=0;
|
||||
|
||||
// touch both grids
|
||||
for(int s=0; s<2; s++) {
|
||||
|
||||
GRID_REGION_INIT();
|
||||
#if PARALLEL==1
|
||||
#include "paraloopstart.h"
|
||||
#endif // PARALLEL==1
|
||||
GRID_REGION_START();
|
||||
GRID_LOOP_START();
|
||||
FORDF0{ RAC(ccel,l) = 0.; }
|
||||
*pFlagSrc =0;
|
||||
*pFlagDst =0;
|
||||
//errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
|
||||
#if PARALLEL!=1
|
||||
GRID_LOOPREG_END();
|
||||
#else // PARALLEL==1
|
||||
#include "paraloopend.h" // = GRID_LOOPREG_END();
|
||||
#endif // PARALLEL==1
|
||||
//GRID_REGION_END();
|
||||
// TEST! */
|
||||
|
||||
/* dummy remove warnings */
|
||||
calcCurrentMass = calcCurrentVolume = 0.;
|
||||
calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
|
||||
|
||||
// change grid
|
||||
mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
|
||||
mLevel[mMaxRefine].setCurr ^= 1;
|
||||
}
|
||||
}
|
||||
|
||||
void
|
||||
LbmFsgrSolver::standingFluidPreinit()
|
||||
{
|
||||
const int lev = mMaxRefine;
|
||||
const bool doReduce = false;
|
||||
const int gridLoopBound=1;
|
||||
|
||||
GRID_REGION_INIT();
|
||||
#if PARALLEL==1
|
||||
#include "paraloopstart.h"
|
||||
#endif // PARALLEL==1
|
||||
GRID_REGION_START();
|
||||
|
||||
LbmFloat rho, ux,uy,uz, usqr;
|
||||
CellFlagType nbflag[LBM_DFNUM];
|
||||
LbmFloat m[LBM_DFNUM];
|
||||
LbmFloat lcsmqo;
|
||||
# if OPT3D==1
|
||||
LbmFloat lcsmqadd, lcsmeq[LBM_DFNUM], lcsmomega;
|
||||
CellFlagType nbored=0;
|
||||
# endif // OPT3D==true
|
||||
|
||||
GRID_LOOP_START();
|
||||
//FORDF0{ RAC(ccel,l) = 0.; }
|
||||
//*pFlagSrc =0;
|
||||
//*pFlagDst =0;
|
||||
//errMsg("l1"," at "<<PRINT_IJK<<" id"<<id);
|
||||
const CellFlagType currFlag = *pFlagSrc; //RFLAG(lev, i,j,k,workSet);
|
||||
if( (currFlag & (CFEmpty|CFBnd)) ) continue;
|
||||
//ccel = RACPNT(lev, i,j,k,workSet);
|
||||
//tcel = RACPNT(lev, i,j,k,otherSet);
|
||||
|
||||
if( (currFlag & (CFInter)) ) {
|
||||
// copy all values
|
||||
for(int l=0; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
|
||||
continue;
|
||||
}
|
||||
|
||||
if( (currFlag & CFNoBndFluid)) {
|
||||
OPTIMIZED_STREAMCOLLIDE;
|
||||
} else {
|
||||
FORDF1 {
|
||||
nbflag[l] = RFLAG_NB(lev, i,j,k, SRCS(lev),l);
|
||||
}
|
||||
DEFAULT_STREAM;
|
||||
//ux = [0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2];
|
||||
DEFAULT_COLLIDEG(mLevel[lev].gravity);
|
||||
}
|
||||
for(int l=LBM_DFNUM; l<dTotalNum;l++) { RAC(tcel,l) = RAC(ccel,l); }
|
||||
#if PARALLEL!=1
|
||||
GRID_LOOPREG_END();
|
||||
#else // PARALLEL==1
|
||||
#include "paraloopend.h" // = GRID_LOOPREG_END();
|
||||
#endif // PARALLEL==1
|
||||
//GRID_REGION_END();
|
||||
// TEST! */
|
||||
|
||||
/* dummy remove warnings */
|
||||
calcCurrentMass = calcCurrentVolume = 0.;
|
||||
calcCellsFilled = calcCellsEmptied = calcNumUsedCells = 0;
|
||||
|
||||
// change grid
|
||||
mLevel[mMaxRefine].setOther = mLevel[mMaxRefine].setCurr;
|
||||
mLevel[mMaxRefine].setCurr ^= 1;
|
||||
}
|
||||
|
||||
|
||||
/******************************************************************************
|
||||
* work on lists from updateCellMass to reinit cell flags
|
||||
*****************************************************************************/
|
||||
|
||||
File diff suppressed because it is too large
Load Diff
@ -11,16 +11,34 @@
|
||||
#include "solver_relax.h"
|
||||
#include "particletracer.h"
|
||||
|
||||
// MPT
|
||||
#include "ntl_world.h"
|
||||
#include "simulation_object.h"
|
||||
|
||||
|
||||
/******************************************************************************
|
||||
* helper functions
|
||||
*****************************************************************************/
|
||||
|
||||
// try to enhance surface?
|
||||
#define SURFACE_ENH 2
|
||||
|
||||
//! for raytracing
|
||||
void LbmFsgrSolver::prepareVisualization( void ) {
|
||||
int lev = mMaxRefine;
|
||||
int workSet = mLevel[lev].setCurr;
|
||||
|
||||
int mainGravDir=0;
|
||||
LbmFloat mainGravLen = 0.;
|
||||
FORDF1{
|
||||
LbmFloat thisGravLen = dot(LbmVec(dfVecX[l],dfVecY[l],dfVecZ[l]), getNormalized(mLevel[mMaxRefine].gravity) );
|
||||
if(thisGravLen>mainGravLen) {
|
||||
mainGravLen = thisGravLen;
|
||||
mainGravDir = l;
|
||||
}
|
||||
//errMsg("GRAVT","l"<<l<<" dfv"<<LbmVec(dfVecX[l],dfVecY[l],dfVecZ[l])<<" g"<< mLevel[mMaxRefine].gravity<< " mgl"<<mainGravLen<<" mgd"<<mainGravDir);
|
||||
}
|
||||
|
||||
//make same prepareVisualization and getIsoSurface...
|
||||
#if LBMDIM==2
|
||||
// 2d, place in the middle of isofield slice (k=2)
|
||||
@ -45,8 +63,17 @@ void LbmFsgrSolver::prepareVisualization( void ) {
|
||||
}
|
||||
#endif // LBMDIM==2
|
||||
|
||||
// MPT, ignore
|
||||
//if( strstr( this->getName().c_str(), "mpfluid2" ) != NULL) {
|
||||
//errMsg("MPT TESTX","skipping mpfluid2");
|
||||
//return;
|
||||
//}
|
||||
if( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) {
|
||||
mpIso->resetAll(0.);
|
||||
}
|
||||
|
||||
LbmFloat minval = this->mIsoValue*1.1; // / mIsoWeight[13];
|
||||
|
||||
LbmFloat minval = this->mIsoValue*1.05; // / mIsoWeight[13];
|
||||
// add up...
|
||||
float val = 0.0;
|
||||
for(int k= getForZMin1(); k< getForZMax1(lev); ++k)
|
||||
@ -54,60 +81,67 @@ void LbmFsgrSolver::prepareVisualization( void ) {
|
||||
for(int i=1;i<mLevel[lev].lSizex-1;i++) {
|
||||
const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
|
||||
//if(cflag&(CFBnd|CFEmpty)) {
|
||||
if(cflag&(CFBnd)) {
|
||||
|
||||
#if SURFACE_ENH==0
|
||||
|
||||
// no enhancements...
|
||||
if( (cflag&(CFFluid|CFUnused)) ) {
|
||||
val = 1.;
|
||||
} else if( (cflag&CFInter) ) {
|
||||
val = (QCELL(lev, i,j,k,workSet, dFfrac));
|
||||
} else {
|
||||
continue;
|
||||
}
|
||||
|
||||
} else if( (cflag&CFEmpty) ) {
|
||||
//continue; // OFF DEBUG
|
||||
int noslipbnd = 0;
|
||||
int intercnt = 0;
|
||||
CellFlagType nbored;
|
||||
FORDF1 {
|
||||
const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
|
||||
if((l==6)&&(nbflag&CFBnd)){ noslipbnd=1; }
|
||||
if(nbflag&CFInter){ intercnt++; }
|
||||
nbored |= nbflag;
|
||||
}
|
||||
//? val = (QCELL(lev, i,j,k,workSet, dFfrac));
|
||||
if((noslipbnd)&&(intercnt>6)) {
|
||||
//if(val<minval) val = minval;
|
||||
//*this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] );
|
||||
*this->mpIso->lbmGetData(i,j,ZKOFF) += minval;
|
||||
} else if((noslipbnd)&&(intercnt>0)) {
|
||||
*this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.95;
|
||||
} else {
|
||||
}
|
||||
#else // SURFACE_ENH!=1
|
||||
if(cflag&CFBnd) {
|
||||
// treated in second loop
|
||||
continue;
|
||||
|
||||
//} else if( (cflag&CFInter) ) {
|
||||
|
||||
} else if( (cflag&CFFluid) && (cflag&CFNoBndFluid) ) {
|
||||
} else if(cflag&CFUnused) {
|
||||
val = 1.;
|
||||
} else if( (cflag&CFFluid) && (cflag&CFNoBndFluid)) {
|
||||
// optimized fluid
|
||||
val = 1.;
|
||||
|
||||
} else if( (cflag&(CFInter|CFFluid)) ) {
|
||||
} else if( (cflag&(CFEmpty|CFInter|CFFluid)) ) {
|
||||
int noslipbnd = 0;
|
||||
int intercnt = 0;
|
||||
FORDF1 {
|
||||
const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
|
||||
if((l==6)&&(nbflag&CFBnd)){ noslipbnd=1; l=100; }
|
||||
if(nbflag&CFInter){ intercnt++; }
|
||||
|
||||
if(l!=mainGravDir) continue; // only check bnd along main grav. dir
|
||||
//if((nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; }
|
||||
if((nbflag&CFBnd)){ noslipbnd=1; }
|
||||
}
|
||||
// no empty nb interface cells are treated as full
|
||||
if(cflag&(CFNoNbEmpty|CFFluid)) {
|
||||
|
||||
if(cflag&CFEmpty) {
|
||||
// special empty treatment
|
||||
if((noslipbnd)&&(intercnt>6)) {
|
||||
*this->mpIso->lbmGetData(i,j,ZKOFF) += minval;
|
||||
} else if((noslipbnd)&&(intercnt>0)) {
|
||||
// necessary?
|
||||
*this->mpIso->lbmGetData(i,j,ZKOFF) += this->mIsoValue*0.9;
|
||||
} else {
|
||||
// nothing to do...
|
||||
}
|
||||
|
||||
continue;
|
||||
} else if(cflag&(CFNoNbEmpty|CFFluid)) {
|
||||
// no empty nb interface cells are treated as full
|
||||
val=1.0;
|
||||
} else {
|
||||
val = (QCELL(lev, i,j,k,workSet, dFfrac));
|
||||
}
|
||||
val = (QCELL(lev, i,j,k,workSet, dFfrac));
|
||||
|
||||
if(noslipbnd) {
|
||||
//errMsg("NEWVAL", PRINT_IJK<<" val"<<val <<" f"<< convertCellFlagType2String(cflag)<<" "<<noslipbnd); //" nbfl"<<convertCellFlagType2String(nbored) );
|
||||
if(val<minval) val = minval;
|
||||
*this->mpIso->lbmGetData(i,j,ZKOFF) += minval-( val * mIsoWeight[13] );
|
||||
}
|
||||
// */
|
||||
#endif // SURFACE_ENH>0
|
||||
|
||||
} else { // unused?
|
||||
} else { // all others, unused?
|
||||
continue;
|
||||
// old fluid val = 1.0;
|
||||
} // */
|
||||
}
|
||||
|
||||
*this->mpIso->lbmGetData( i-1 , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[0] );
|
||||
*this->mpIso->lbmGetData( i , j-1 ,ZKOFF-ZKD1) += ( val * mIsoWeight[1] );
|
||||
@ -148,6 +182,105 @@ void LbmFsgrSolver::prepareVisualization( void ) {
|
||||
*this->mpIso->lbmGetData( i+1 , j+1 ,ZKOFF+ZKD1) += ( val * mIsoWeight[26] );
|
||||
}
|
||||
|
||||
// TEST!?
|
||||
#if SURFACE_ENH>=2
|
||||
|
||||
if(mFsSurfGenSetting&fssgNoObs) {
|
||||
for(int k= getForZMin1(); k< getForZMax1(lev); ++k)
|
||||
for(int j=1;j<mLevel[lev].lSizey-1;j++)
|
||||
for(int i=1;i<mLevel[lev].lSizex-1;i++) {
|
||||
const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
|
||||
if(cflag&(CFBnd)) {
|
||||
CellFlagType nbored=0;
|
||||
LbmFloat avgfill=0.,avgfcnt=0.;
|
||||
FORDF1 {
|
||||
const int ni = i+this->dfVecX[l];
|
||||
const int nj = j+this->dfVecY[l];
|
||||
const int nk = ZKOFF+this->dfVecZ[l];
|
||||
|
||||
const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet);
|
||||
nbored |= nbflag;
|
||||
if(nbflag&CFInter) {
|
||||
avgfill += QCELL(lev, ni,nj,nk, workSet,dFfrac); avgfcnt += 1.;
|
||||
} else if(nbflag&CFFluid) {
|
||||
avgfill += 1.; avgfcnt += 1.;
|
||||
} else if(nbflag&CFEmpty) {
|
||||
avgfill += 0.; avgfcnt += 1.;
|
||||
}
|
||||
|
||||
//if( (ni<0) || (nj<0) || (nk<0)
|
||||
//|| (ni>=mLevel[mMaxRefine].lSizex)
|
||||
//|| (nj>=mLevel[mMaxRefine].lSizey)
|
||||
//|| (nk>=mLevel[mMaxRefine].lSizez) ) continue;
|
||||
}
|
||||
|
||||
if(nbored&CFInter) {
|
||||
if(avgfcnt>0.) avgfill/=avgfcnt;
|
||||
*this->mpIso->lbmGetData(i,j,ZKOFF) = avgfill; continue;
|
||||
}
|
||||
else if(nbored&CFFluid) {
|
||||
*this->mpIso->lbmGetData(i,j,ZKOFF) = 1.; continue;
|
||||
}
|
||||
|
||||
}
|
||||
}
|
||||
|
||||
// move surface towards inner "row" of obstacle
|
||||
// cells if necessary (all obs cells without fluid/inter
|
||||
// nbs (=iso==0) next to obstacles...)
|
||||
for(int k= getForZMin1(); k< getForZMax1(lev); ++k)
|
||||
for(int j=1;j<mLevel[lev].lSizey-1;j++)
|
||||
for(int i=1;i<mLevel[lev].lSizex-1;i++) {
|
||||
const CellFlagType cflag = RFLAG(lev, i,j,k,workSet);
|
||||
if( (cflag&(CFBnd)) && (*this->mpIso->lbmGetData(i,j,ZKOFF)==0.)) {
|
||||
int bndnbcnt=0;
|
||||
FORDF1 {
|
||||
const int ni = i+this->dfVecX[l];
|
||||
const int nj = j+this->dfVecY[l];
|
||||
const int nk = ZKOFF+this->dfVecZ[l];
|
||||
const CellFlagType nbflag = RFLAG(lev, ni,nj,nk, workSet);
|
||||
if(nbflag&CFBnd) bndnbcnt++;
|
||||
}
|
||||
if(bndnbcnt>0) *this->mpIso->lbmGetData(i,j,ZKOFF)=this->mIsoValue*0.95;
|
||||
}
|
||||
}
|
||||
}
|
||||
// */
|
||||
|
||||
if(mFsSurfGenSetting&fssgNoNorth)
|
||||
for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k)
|
||||
for(int j=0;j<mLevel[lev].lSizey-0;j++) {
|
||||
*this->mpIso->lbmGetData(0, j,ZKOFF) = *this->mpIso->lbmGetData(1, j,ZKOFF);
|
||||
}
|
||||
if(mFsSurfGenSetting&fssgNoEast)
|
||||
for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k)
|
||||
for(int i=0;i<mLevel[lev].lSizex-0;i++) {
|
||||
*this->mpIso->lbmGetData(i,0, ZKOFF) = *this->mpIso->lbmGetData(i,1, ZKOFF);
|
||||
}
|
||||
if(mFsSurfGenSetting&fssgNoSouth)
|
||||
for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k)
|
||||
for(int j=0;j<mLevel[lev].lSizey-0;j++) {
|
||||
*this->mpIso->lbmGetData(mLevel[lev].lSizex-1,j,ZKOFF) = *this->mpIso->lbmGetData(mLevel[lev].lSizex-2,j,ZKOFF);
|
||||
}
|
||||
if(mFsSurfGenSetting&fssgNoWest)
|
||||
for(int k= getForZMinBnd(); k< getForZMaxBnd(lev); ++k)
|
||||
for(int i=0;i<mLevel[lev].lSizex-0;i++) {
|
||||
*this->mpIso->lbmGetData(i,mLevel[lev].lSizey-1,ZKOFF) = *this->mpIso->lbmGetData(i,mLevel[lev].lSizey-2,ZKOFF);
|
||||
}
|
||||
if(LBMDIM>2) {
|
||||
if(mFsSurfGenSetting&fssgNoBottom)
|
||||
for(int j=0;j<mLevel[lev].lSizey-0;j++)
|
||||
for(int i=0;i<mLevel[lev].lSizex-0;i++) {
|
||||
*this->mpIso->lbmGetData(i,j,0 ) = *this->mpIso->lbmGetData(i,j,1 );
|
||||
}
|
||||
if(mFsSurfGenSetting&fssgNoTop)
|
||||
for(int j=0;j<mLevel[lev].lSizey-0;j++)
|
||||
for(int i=0;i<mLevel[lev].lSizex-0;i++) {
|
||||
*this->mpIso->lbmGetData(i,j,mLevel[lev].lSizez-1) = *this->mpIso->lbmGetData(i,j,mLevel[lev].lSizez-2);
|
||||
}
|
||||
}
|
||||
#endif // SURFACE_ENH>=2
|
||||
|
||||
|
||||
// update preview, remove 2d?
|
||||
if((this->mOutputSurfacePreview)&&(LBMDIM==3)) {
|
||||
@ -215,7 +348,47 @@ void LbmFsgrSolver::prepareVisualization( void ) {
|
||||
}
|
||||
}
|
||||
|
||||
// correction
|
||||
// MPT
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
if( ( strstr( this->getName().c_str(), "mpfluid1" ) != NULL) && (mInitDone) ) {
|
||||
errMsg("MPT TESTX","special mpfluid1");
|
||||
vector<SimulationObject*> *sims = mpGlob->getSims();
|
||||
for(int simi=0; simi<(int)(sims->size()); simi++) {
|
||||
SimulationObject *sim = (*sims)[simi];
|
||||
if( strstr( sim->getName().c_str(), "mpfluid2" ) != NULL) {
|
||||
LbmFsgrSolver *fsgr = (LbmFsgrSolver *)(sim->getSolver());
|
||||
int msyLev = mMaxRefine;
|
||||
int simLev = fsgr->mMaxRefine;
|
||||
|
||||
IsoSurface* simIso = fsgr->mpIso;
|
||||
int idst = mLevel[msyLev].lSizex-2 -1; // start at boundary
|
||||
int isrc = 0 +2;
|
||||
if((simIso)&&(fsgr->mInitDone)) {
|
||||
fsgr->prepareVisualization();
|
||||
|
||||
for(int k=0;k<=mLevel[simLev].lSizez-1;++k) {
|
||||
for(int j=0;j<=mLevel[simLev].lSizey-1;++j) {
|
||||
for(int i=isrc; i<mLevel[simLev].lSizex-1; i++) {
|
||||
//LbmFloat *cceldst = &QCELL( msyLev ,idst+i,j,k, msySet,0);
|
||||
//LbmFloat *ccelsrc = &SIM_QCELL(fsgr, simLev ,isrc+i,j,k, simSet,0);
|
||||
//for(int l=0; l<dTotalNum; l++) {
|
||||
//RAC(cceldst,l) = RAC(ccelsrc,l);
|
||||
//}
|
||||
// both flag sets, make sure curr/other are same!?
|
||||
//RFLAG(msyLev ,idst+i,j,k, 0) = SIM_RFLAG(fsgr, simLev ,isrc+i,j,k, 0);
|
||||
//RFLAG(msyLev ,idst+i,j,k, 1) = SIM_RFLAG(fsgr, simLev ,isrc+i,j,k, 1);
|
||||
errMsg("TESTX"," "<<PRINT_VEC(idst+i,j,k)<<" < "<<PRINT_VEC(i,j,k) );
|
||||
*mpIso->lbmGetData(idst+i,j,k) = *simIso->lbmGetData(i,j,k);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
} // simIso
|
||||
}
|
||||
}
|
||||
}
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
|
||||
return;
|
||||
}
|
||||
|
||||
@ -315,7 +488,7 @@ int LbmFsgrSolver::initParticles() {
|
||||
//if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid) )
|
||||
//&& TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFNoNbFluid )
|
||||
//if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter|CFMbndInflow) ) {
|
||||
if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid) ) {
|
||||
if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFNoBndFluid|CFUnused) ) {
|
||||
bool cellOk = true;
|
||||
//? FORDF1 { if(!(RFLAG_NB(mMaxRefine,i,j,k,workSet, l) & CFFluid)) cellOk = false; }
|
||||
if(!cellOk) continue;
|
||||
@ -426,10 +599,22 @@ int LbmFsgrSolver::initParticles() {
|
||||
/* errMsg("PIT","U pit"<<(p)->getId()<<" pos:"<< (p)->getPos()<<" status:"<<convertFlags2String((p)->getFlags())<<" to "<< (newtype) ); */ \
|
||||
p->setType(newtype);
|
||||
|
||||
// tracer defines
|
||||
#define TRACE_JITTER 0.025
|
||||
#define TRACE_RAND (rand()/(RAND_MAX+1.0))*TRACE_JITTER-(TRACE_JITTER*0.5)
|
||||
#define FFGET_NORM(var,dl) \
|
||||
if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFInter)){ (var) = QCELL_NB(lev,i,j,k,workSet,dl,dFfrac); } \
|
||||
else if(RFLAG_NB(lev,i,j,k,workSet, dl) &(CFFluid|CFUnused)){ (var) = 1.; } else (var) = 0.0;
|
||||
|
||||
// float jitter
|
||||
#define FLOAT_JITTER_BND (FLOAT_JITTER*2.0)
|
||||
#define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5))
|
||||
|
||||
|
||||
void LbmFsgrSolver::advanceParticles() {
|
||||
int workSet = mLevel[mMaxRefine].setCurr;
|
||||
LbmFloat vx=0.0,vy=0.0,vz=0.0;
|
||||
LbmFloat rho, df[27]; //feq[27];
|
||||
//LbmFloat df[27]; //feq[27];
|
||||
#define DEL_PART { \
|
||||
/*errMsg("PIT","DEL AT "<< __LINE__<<" type:"<<p->getType()<<" "); */ \
|
||||
p->setActive( false ); \
|
||||
@ -453,10 +638,6 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
const bool useff = (mFarFieldSize>2.); // if(mpTest->mDebugvalue1>0.0){
|
||||
|
||||
// TODO use timestep size
|
||||
//bool isIn,isOut,isInZ;
|
||||
//const int cutval = 1+mCutoff/2; // TODO FIXME add half border!
|
||||
//const int cutval = mCutoff/2; // TODO FIXME add half border!
|
||||
//const int cutval = 0; // TODO FIXME add half border!
|
||||
const int cutval = mCutoff; // use full border!?
|
||||
int actCnt=0;
|
||||
if(this->mStepCnt%50==49) { mpParticles->cleanup(); }
|
||||
@ -518,27 +699,39 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
(p->getType()==PART_TRACER) ) {
|
||||
|
||||
// no interpol
|
||||
rho = vx = vy = vz = 0.0;
|
||||
vx = vy = vz = 0.0;
|
||||
if(p->getStatus()&PART_IN) { // IN
|
||||
if(k>=cutval) {
|
||||
if(k>this->mSizez-1-cutval) DEL_PART;
|
||||
FORDF0{
|
||||
LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
|
||||
df[l] = cdf;
|
||||
rho += cdf;
|
||||
vx += (this->dfDvecX[l]*cdf);
|
||||
vy += (this->dfDvecY[l]*cdf);
|
||||
vz += (this->dfDvecZ[l]*cdf);
|
||||
}
|
||||
|
||||
// remove gravity influence
|
||||
const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les
|
||||
vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5;
|
||||
vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5;
|
||||
vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5;
|
||||
|
||||
if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid) ) {
|
||||
if( RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFUnused) ) {
|
||||
// still ok
|
||||
int partLev = mMaxRefine;
|
||||
int si=i, sj=j, sk=k;
|
||||
while(partLev>0 && RFLAG(partLev, si,sj,sk, workSet)&(CFUnused)) {
|
||||
partLev--;
|
||||
si/=2;
|
||||
sj/=2;
|
||||
sk/=2;
|
||||
}
|
||||
//LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
|
||||
LbmFloat *ccel = RACPNT(partLev, si,sj,sk, workSet);
|
||||
FORDF1{
|
||||
//LbmFloat cdf = QCELL(mMaxRefine, i,j,k, workSet, l);
|
||||
LbmFloat cdf = RAC(ccel, l);
|
||||
// TODO update below
|
||||
//df[l] = cdf;
|
||||
vx += (this->dfDvecX[l]*cdf);
|
||||
vy += (this->dfDvecY[l]*cdf);
|
||||
vz += (this->dfDvecZ[l]*cdf);
|
||||
}
|
||||
|
||||
// remove gravity influence
|
||||
const LbmFloat lesomega = mLevel[mMaxRefine].omega; // no les
|
||||
vx -= mLevel[mMaxRefine].gravity[0] * lesomega*0.5;
|
||||
vy -= mLevel[mMaxRefine].gravity[1] * lesomega*0.5;
|
||||
vz -= mLevel[mMaxRefine].gravity[2] * lesomega*0.5;
|
||||
|
||||
} else { // OUT
|
||||
// out of bounds, deactivate...
|
||||
// FIXME make fsgr treatment
|
||||
@ -548,12 +741,12 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
// below 3d region, just rise
|
||||
}
|
||||
} else { // OUT
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
# if LBM_INCLUDE_TESTSOLVERS==1
|
||||
if(useff) { mpTest->handleParticle(p, i,j,k); }
|
||||
else DEL_PART;
|
||||
#else // LBM_INCLUDE_TESTSOLVERS==1
|
||||
# else // LBM_INCLUDE_TESTSOLVERS==1
|
||||
DEL_PART;
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
# endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
// TODO use x,y vel...?
|
||||
}
|
||||
|
||||
@ -577,15 +770,11 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
//LbmVec change = (fb+fd) *timestep / (pvolume*rhoAir) *(timestep/cellsize);
|
||||
//actCnt++; // should be after active test
|
||||
if(actCnt<0) {
|
||||
errMsg("\nPIT","BTEST1 vol="<<pvolume<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel);
|
||||
errMsg("PIT","BTEST1 vol="<<pvolume<<" radius="<<radius<<" vn="<<velRelNorm<<" velPart="<<velPart<<" velRel"<<velRel);
|
||||
errMsg("PIT","BTEST2 cellsize="<<cellsize<<" timestep="<<timestep<<" viscW="<<viscWater<<" ss/mb="<<(timestep/(pvolume*rhoAir)));
|
||||
errMsg("PIT","BTEST2 grav="<<rwgrav<<" " );
|
||||
errMsg("PIT","BTEST2 change="<<(change)<<" fb="<<(fb)<<" fd="<<(fd)<<" ");
|
||||
errMsg("PIT","BTEST2 change="<<norm(change)<<" fb="<<norm(fb)<<" fd="<<norm(fd)<<" ");
|
||||
#if LOOPTEST==1
|
||||
errMsg("PIT","BTEST2 n="<<n<<" "); // LOOPTEST! DEBUG
|
||||
#endif // LOOPTEST==1
|
||||
errMsg("PIT","\n");
|
||||
}
|
||||
|
||||
//v += change;
|
||||
@ -601,25 +790,15 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
if(fflag&(CFFluid|CFInter) ) { p->setInFluid(true);
|
||||
} else { p->setInFluid(false); }
|
||||
|
||||
//if(( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) {
|
||||
if( (( fflag&CFFluid ) && ( fflag&CFNoBndFluid )) ||
|
||||
(( fflag&CFInter ) && (!(fflag&CFNoNbFluid)) ) ) {
|
||||
// only real fluid
|
||||
#define TRACE_JITTER 0.025
|
||||
#define TRACE_RAND (rand()/(RAND_MAX+1.0))*TRACE_JITTER-(TRACE_JITTER*0.5)
|
||||
#define __TRACE_RAND 0.
|
||||
#if LBMDIM==3
|
||||
# if LBMDIM==3
|
||||
p->advance( TRACE_RAND,TRACE_RAND,TRACE_RAND);
|
||||
#else
|
||||
# else
|
||||
p->advance( TRACE_RAND,TRACE_RAND, 0.);
|
||||
#endif
|
||||
# endif
|
||||
|
||||
//} // test!?
|
||||
//if( fflag&(CFNoBndFluid) ) {
|
||||
// ok
|
||||
//} else if( fflag&(CFEmpty|CFNoNbFluid) ) {
|
||||
//} else if( fflag&(CFEmpty) ) {
|
||||
//v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity);
|
||||
} else {
|
||||
// move inwards along normal, make sure normal is valid first
|
||||
const int lev = mMaxRefine;
|
||||
@ -629,25 +808,21 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
if(i>=this->mSizex-1) { nx = 1.; nonorm = true; }
|
||||
if(j<=0) { ny = -1.; nonorm = true; }
|
||||
if(j>=this->mSizey-1) { ny = 1.; nonorm = true; }
|
||||
#if LBMDIM==3
|
||||
# if LBMDIM==3
|
||||
if(k<=0) { nz = -1.; nonorm = true; }
|
||||
if(k>=this->mSizez-1) { nz = 1.; nonorm = true; }
|
||||
#endif // LBMDIM==3
|
||||
//mynbfac = QCELL_NB(lev, i,j,k,SRCS(lev),l, dFlux) / QCELL(lev, i,j,k,SRCS(lev), dFlux);
|
||||
# endif // LBMDIM==3
|
||||
if(!nonorm) {
|
||||
if(RFLAG_NB(lev,i,j,k,workSet, dE) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dE,dFfrac); } else nv1 = 0.0;
|
||||
if(RFLAG_NB(lev,i,j,k,workSet, dW) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dW,dFfrac); } else nv2 = 0.0;
|
||||
FFGET_NORM(nv1,dE); FFGET_NORM(nv2,dW);
|
||||
nx = 0.5* (nv2-nv1);
|
||||
if(RFLAG_NB(lev,i,j,k,workSet, dN) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dN,dFfrac); } else nv1 = 0.0;
|
||||
if(RFLAG_NB(lev,i,j,k,workSet, dS) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dS,dFfrac); } else nv2 = 0.0;
|
||||
FFGET_NORM(nv1,dN); FFGET_NORM(nv2,dS);
|
||||
ny = 0.5* (nv2-nv1);
|
||||
#if LBMDIM==3
|
||||
if(RFLAG_NB(lev,i,j,k,workSet, dT) &(CFFluid|CFInter)){ nv1 = QCELL_NB(lev,i,j,k,workSet,dT,dFfrac); } else nv1 = 0.0;
|
||||
if(RFLAG_NB(lev,i,j,k,workSet, dB) &(CFFluid|CFInter)){ nv2 = QCELL_NB(lev,i,j,k,workSet,dB,dFfrac); } else nv2 = 0.0;
|
||||
# if LBMDIM==3
|
||||
FFGET_NORM(nv1,dT); FFGET_NORM(nv2,dB);
|
||||
nz = 0.5* (nv2-nv1);
|
||||
#else //LBMDIM==3
|
||||
nz = 0.0;
|
||||
#endif //LBMDIM==3
|
||||
# else // LBMDIM==3
|
||||
nz = 0.;
|
||||
# endif // LBMDIM==3
|
||||
} else {
|
||||
v = p->getVel() + vec2G(mLevel[mMaxRefine].gravity);
|
||||
}
|
||||
@ -708,7 +883,7 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
if( RFLAG(mMaxRefine, i,j,k, workSet)& (CFEmpty|CFInter|CFBnd)) {
|
||||
// still ok
|
||||
// shipt3 } else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFInter) ){
|
||||
} else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid) ){
|
||||
} else if( RFLAG(mMaxRefine, i,j,k, workSet) & (CFFluid|CFUnused) ){
|
||||
// FIXME make fsgr treatment
|
||||
P_CHANGETYPE(p, PART_FLOAT ); continue;
|
||||
// jitter in cell to prevent stacking when hitting a steep surface
|
||||
@ -721,12 +896,12 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
}
|
||||
}
|
||||
} else { // OUT
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
# if LBM_INCLUDE_TESTSOLVERS==1
|
||||
if(useff) { mpTest->handleParticle(p, i,j,k); }
|
||||
else{ DEL_PART; }
|
||||
#else // LBM_INCLUDE_TESTSOLVERS==1
|
||||
# else // LBM_INCLUDE_TESTSOLVERS==1
|
||||
{ DEL_PART; }
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
# endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
}
|
||||
|
||||
} // air particle
|
||||
@ -741,7 +916,7 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
|
||||
if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFInter ) ) {
|
||||
// still ok
|
||||
} else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
|
||||
} else if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid|CFUnused) ) ) {
|
||||
//errMsg("PIT","NEWBUB pit "<< (*pit).getPos()<<" status:"<<convertFlags2String((*pit).getFlags()) );
|
||||
|
||||
//P_CHANGETYPE(p, PART_BUBBLE ); continue;
|
||||
@ -767,7 +942,7 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
// test - delte on boundary!?
|
||||
//if( (i<=cutval)||(i>=this->mSizex-1-cutval)|| (j<=cutval)||(j>=this->mSizey-1-cutval)) { DEL_PART; } // DEBUG TEST
|
||||
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
# if LBM_INCLUDE_TESTSOLVERS==1
|
||||
/*LbmFloat prob = 1.0;
|
||||
// vanishing
|
||||
prob = (rand()/(RAND_MAX+1.0));
|
||||
@ -783,30 +958,26 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
//? if(k>this->mSizez*3/4) { if(prob<3.0*mLevel[mMaxRefine].timestep*0.1) DEL_PART;}
|
||||
// */
|
||||
|
||||
#else // LBM_INCLUDE_TESTSOLVERS==1
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
# endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
|
||||
if(p->getStatus()&PART_IN) { // IN
|
||||
//if((k<cutval)||(k>this->mSizez-1-cutval)) DEL_PART;
|
||||
if(k<cutval) DEL_PART;
|
||||
if(k>this->mSizez-1-cutval) { P_CHANGETYPE(p, PART_DROP ); continue; } // create new drop
|
||||
//ntlVec3Gfx v = getVelocityAt(i,j,k);
|
||||
rho = vx = vy = vz = 0.0;
|
||||
vx = vy = vz = 0.0;
|
||||
|
||||
// define from particletracer.h
|
||||
#if MOVE_FLOATS==1
|
||||
const int DEPTH_AVG=3; // only average interface vels
|
||||
int ccnt=0;
|
||||
for(int kk=0;kk<DEPTH_AVG;kk+=1) {
|
||||
//for(int kk=1;kk<DEPTH_AVG;kk+=1) {
|
||||
if((k-kk)<1) continue;
|
||||
//if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFFluid|CFInter)) {} else continue;
|
||||
if(RFLAG(mMaxRefine, i,j,k, workSet)&(CFInter)) {} else continue;
|
||||
ccnt++;
|
||||
FORDF0{
|
||||
FORDF1{
|
||||
LbmFloat cdf = QCELL(mMaxRefine, i,j,k-kk, workSet, l);
|
||||
df[l] = cdf;
|
||||
//rho += cdf;
|
||||
//df[l] = cdf;
|
||||
vx += (this->dfDvecX[l]*cdf);
|
||||
vy += (this->dfDvecY[l]*cdf);
|
||||
vz += (this->dfDvecZ[l]*cdf);
|
||||
@ -824,7 +995,7 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
vy += (rand()/(RAND_MAX+1.0))*(FLOAT_JITTER*0.2)-(FLOAT_JITTER*0.2*0.5);
|
||||
|
||||
//bool delfloat = false;
|
||||
if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), CFFluid ) ) {
|
||||
if( TESTFLAG( RFLAG(mMaxRefine, i,j,k, workSet), (CFFluid|CFUnused) ) ) {
|
||||
// in fluid cell
|
||||
/*if((1) && (k<this->mSizez-3) &&
|
||||
(
|
||||
@ -888,16 +1059,10 @@ void LbmFsgrSolver::advanceParticles() {
|
||||
// use half butoff border 1/8
|
||||
int maxdw = (int)(mLevel[mMaxRefine].lSizex*0.125*0.5);
|
||||
if(maxdw<3) maxdw=3;
|
||||
#define FLOAT_JITTER_BND (FLOAT_JITTER*2.0)
|
||||
#define FLOAT_JITTBNDRAND(x) ((rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND*(1.-(x/(LbmFloat)maxdw))-(FLOAT_JITTER_BND*(1.-(x)/(LbmFloat)maxdw)*0.5))
|
||||
//if(ABS(i-( cutval))<maxdw) { p->advance( (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5) ,0.,0.); }
|
||||
//if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance( (rand()/(RAND_MAX+1.0))*FLOAT_JITTER_BND-(FLOAT_JITTER_BND*0.5) ,0.,0.); }
|
||||
if((j>=0)&&(j<=this->mSizey-1)) {
|
||||
if(ABS(i-( cutval))<maxdw) { p->advance( FLOAT_JITTBNDRAND( ABS(i-( cutval))), 0.,0.); }
|
||||
if(ABS(i-(this->mSizex-1-cutval))<maxdw) { p->advance( FLOAT_JITTBNDRAND( ABS(i-(this->mSizex-1-cutval))), 0.,0.); }
|
||||
}
|
||||
//#undef FLOAT_JITTER_BND
|
||||
#undef FLOAT_JITTBNDRAND
|
||||
//if( (i<cutval)||(i>this->mSizex-1-cutval)|| //(j<cutval)||(j>this->mSizey-1-cutval)
|
||||
}
|
||||
} // PART_FLOAT
|
||||
@ -1306,11 +1471,15 @@ int LbmFsgrSolver::debLBMGI(int level, int ii,int ij,int ik, int is) {
|
||||
if(level < 0){ errMsg("LbmStrict::debLBMGI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(level > mMaxRefine){ errMsg("LbmStrict::debLBMGI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
|
||||
if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if((ii==-1)&&(ij==0)) {
|
||||
// special case for main loop, ok
|
||||
} else {
|
||||
if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
}
|
||||
if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
@ -1338,11 +1507,15 @@ int LbmFsgrSolver::debLBMQI(int level, int ii,int ij,int ik, int is, int l) {
|
||||
if(level < 0){ errMsg("LbmStrict::debLBMQI"," invLev- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(level > mMaxRefine){ errMsg("LbmStrict::debLBMQI"," invLev+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
|
||||
if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if((ii==-1)&&(ij==0)) {
|
||||
// special case for main loop, ok
|
||||
} else {
|
||||
if(ii<0){ errMsg("LbmStrict"," invX- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij<0){ errMsg("LbmStrict"," invY- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
}
|
||||
if(ik<0){ errMsg("LbmStrict"," invZ- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ii>mLevel[level].lSizex-1){ errMsg("LbmStrict"," invX+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ij>mLevel[level].lSizey-1){ errMsg("LbmStrict"," invY+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(ik>mLevel[level].lSizez-1){ errMsg("LbmStrict"," invZ+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(is<0){ errMsg("LbmStrict"," invS- l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
if(is>1){ errMsg("LbmStrict"," invS+ l"<<level<<"|"<<ii<<","<<ij<<","<<ik<<" s"<<is); STRICT_EXIT; }
|
||||
@ -1561,7 +1734,10 @@ void LbmFsgrSolver::debugDisplayNode(int dispset, CellIdentifierInterface* cell
|
||||
void LbmFsgrSolver::lbmDebugDisplay(int dispset) {
|
||||
// DEBUG always display testdata
|
||||
#if LBM_INCLUDE_TESTSOLVERS==1
|
||||
if(mUseTestdata){ mpTest->testDebugDisplay(dispset); }
|
||||
if(mUseTestdata){
|
||||
cpDebugDisplay(dispset);
|
||||
mpTest->testDebugDisplay(dispset);
|
||||
}
|
||||
#endif // LBM_INCLUDE_TESTSOLVERS==1
|
||||
if(dispset<=FLUIDDISPNothing) return;
|
||||
//if(!dispset->on) return;
|
||||
|
||||
Loading…
Reference in New Issue
Block a user