Welcome to mirror list, hosted at ThFree Co, Russian Federation.

git.blender.org/blender.git - Unnamed repository; edit this file 'description' to name the repository.
summaryrefslogtreecommitdiff
path: root/intern
diff options
context:
space:
mode:
authorNils Thuerey <nils@thuerey.de>2006-08-22 14:45:19 +0400
committerNils Thuerey <nils@thuerey.de>2006-08-22 14:45:19 +0400
commitb6257305c929464ef50902ab7c98a7c37dd17fd7 (patch)
tree8bb33edeb0726760795f0ed9dcf27b4af6f36b92 /intern
parentd6d2d6f063994e1481ce498c5e245d8470ef4b6d (diff)
elbeem update:
- slightly improved and optimized handling of moving obstacles - cleanup of debug messages and minor fixes
Diffstat (limited to 'intern')
-rw-r--r--intern/elbeem/intern/elbeem.h1
-rw-r--r--intern/elbeem/intern/isosurface.cpp10
-rw-r--r--intern/elbeem/intern/isosurface.h3
-rw-r--r--intern/elbeem/intern/loop_tools.h105
-rw-r--r--intern/elbeem/intern/ntl_geometryobject.cpp143
-rw-r--r--intern/elbeem/intern/ntl_geometryobject.h5
-rw-r--r--intern/elbeem/intern/simulation_object.cpp2
-rw-r--r--intern/elbeem/intern/solver_adap.cpp44
-rw-r--r--intern/elbeem/intern/solver_class.h89
-rw-r--r--intern/elbeem/intern/solver_init.cpp378
-rw-r--r--intern/elbeem/intern/solver_interface.cpp4
-rw-r--r--intern/elbeem/intern/solver_interface.h1
-rw-r--r--intern/elbeem/intern/solver_main.cpp507
-rw-r--r--intern/elbeem/intern/solver_relax.h1077
-rw-r--r--intern/elbeem/intern/solver_util.cpp422
15 files changed, 1624 insertions, 1167 deletions
diff --git a/intern/elbeem/intern/elbeem.h b/intern/elbeem/intern/elbeem.h
index 8eb8a5433b1..ea21e807a6f 100644
--- a/intern/elbeem/intern/elbeem.h
+++ b/intern/elbeem/intern/elbeem.h
@@ -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];
diff --git a/intern/elbeem/intern/isosurface.cpp b/intern/elbeem/intern/isosurface.cpp
index 373a6f74761..283fd17eb94 100644
--- a/intern/elbeem/intern/isosurface.cpp
+++ b/intern/elbeem/intern/isosurface.cpp
@@ -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;
diff --git a/intern/elbeem/intern/isosurface.h b/intern/elbeem/intern/isosurface.h
index d1ff1529069..5b5198171f3 100644
--- a/intern/elbeem/intern/isosurface.h
+++ b/intern/elbeem/intern/isosurface.h
@@ -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 );
diff --git a/intern/elbeem/intern/loop_tools.h b/intern/elbeem/intern/loop_tools.h
new file mode 100644
index 00000000000..927263b1204
--- /dev/null
+++ b/intern/elbeem/intern/loop_tools.h
@@ -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; ) {
+
diff --git a/intern/elbeem/intern/ntl_geometryobject.cpp b/intern/elbeem/intern/ntl_geometryobject.cpp
index 90c3b3437ea..d815fb55287 100644
--- a/intern/elbeem/intern/ntl_geometryobject.cpp
+++ b/intern/elbeem/intern/ntl_geometryobject.cpp
@@ -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;
}
}
diff --git a/intern/elbeem/intern/ntl_geometryobject.h b/intern/elbeem/intern/ntl_geometryobject.h
index 5fe994682e3..a9a8109ba1b 100644
--- a/intern/elbeem/intern/ntl_geometryobject.h
+++ b/intern/elbeem/intern/ntl_geometryobject.h
@@ -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? */
diff --git a/intern/elbeem/intern/simulation_object.cpp b/intern/elbeem/intern/simulation_object.cpp
index e34030180a5..3f954f1a4c7 100644
--- a/intern/elbeem/intern/simulation_object.cpp
+++ b/intern/elbeem/intern/simulation_object.cpp
@@ -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;
}
diff --git a/intern/elbeem/intern/solver_adap.cpp b/intern/elbeem/intern/solver_adap.cpp
index aade9b30545..99892081fca 100644
--- a/intern/elbeem/intern/solver_adap.cpp
+++ b/intern/elbeem/intern/solver_adap.cpp
@@ -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);
}
diff --git a/intern/elbeem/intern/solver_class.h b/intern/elbeem/intern/solver_class.h
index 46462ec52de..578e1079cd7 100644
--- a/intern/elbeem/intern/solver_class.h
+++ b/intern/elbeem/intern/solver_class.h
@@ -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
@@ -298,8 +303,10 @@ class LbmFsgrSolver :
void printLbmCell(int level, int i, int j, int k,int set);
// 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 );
+
+ 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
- public: // former LbmModelLBGK functions
+ // 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;
diff --git a/intern/elbeem/intern/solver_init.cpp b/intern/elbeem/intern/solver_init.cpp
index 5fb3f84743a..1a19e1deabf 100644
--- a/intern/elbeem/intern/solver_init.cpp
+++ b/intern/elbeem/intern/solver_init.cpp
@@ -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;
diff --git a/intern/elbeem/intern/solver_interface.cpp b/intern/elbeem/intern/solver_interface.cpp
index 7b3f71e438c..e9db2a10db7 100644
--- a/intern/elbeem/intern/solver_interface.cpp
+++ b/intern/elbeem/intern/solver_interface.cpp
@@ -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 );
}
diff --git a/intern/elbeem/intern/solver_interface.h b/intern/elbeem/intern/solver_interface.h
index 4cbdd2945e5..b95ce0358cb 100644
--- a/intern/elbeem/intern/solver_interface.h
+++ b/intern/elbeem/intern/solver_interface.h
@@ -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
diff --git a/intern/elbeem/intern/solver_main.cpp b/intern/elbeem/intern/solver_main.cpp
index 8fa3b8d8a65..964192420e4 100644
--- a/intern/elbeem/intern/solver_main.cpp
+++ b/intern/elbeem/intern/solver_main.cpp
@@ -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
- 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
+ GRID_LOOP_START();
+ // loop start
+ // stream from current set to other, then collide and store
+ //errMsg("l2"," at "<<PRINT_IJK<<" id"<<id);
- 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));
- }
-#else // OPT3D==0
+# 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 );
+ }
+ }
+ 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,38 +1091,133 @@ 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;
+}
+
+
+
+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
- PARA_FINISH();
+#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! */
- // check other vars...?
+ /* 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;
+}
/******************************************************************************
diff --git a/intern/elbeem/intern/solver_relax.h b/intern/elbeem/intern/solver_relax.h
index 1a1e2c01172..2c99a853c62 100644
--- a/intern/elbeem/intern/solver_relax.h
+++ b/intern/elbeem/intern/solver_relax.h
@@ -15,17 +15,51 @@
#define CAUSE_PANIC { this->mPanic=1; } /*set flag*/
#endif // FSGR_STRICT_DEBUG==1
+#if LBM_INCLUDE_TESTSOLVERS!=1
+#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
+ ux += (grav)[0]; \
+ uy += (grav)[1]; \
+ uz += (grav)[2];
+#define TEST_IF_CHECK
+#else // LBM_INCLUDE_TESTSOLVERS!=1
+// defined in test.h
+
+#define NEWDIRVELMOTEST 0
+#if NEWDIRVELMOTEST==1
// off for non testing
+#undef PRECOLLIDE_MODS
#define PRECOLLIDE_MODS(rho,ux,uy,uz, grav) \
ux += (grav)[0]; \
uy += (grav)[1]; \
- uz += (grav)[2]; \
-
-#define TEST_IF_CHECK
+ uz += (grav)[2]; \
+ { \
+ int lev = mMaxRefine, nomb=0; \
+ LbmFloat bcnt = 0.,nux=0.,nuy=0.,nuz=0.; \
+ for(int l=1; l<this->cDfNum; l++) { \
+ if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { \
+ if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBndMoving) { \
+ nux += QCELL_NB(lev, i,j,k,SRCS(lev),l, dMass); \
+ nuy += QCELL_NB(lev, i,j,k,SRCS(lev),l, dFfrac); \
+ bcnt += 1.; \
+ } else { \
+ nomb++; \
+ } \
+ } \
+ } \
+ if((bcnt>0.)&&(nomb==0)) { \
+ ux = nux/bcnt; \
+ uy = nuy/bcnt; \
+ uz = nuz/bcnt; \
+ } \
+ }
+#else // NEWDIRVELMOTEST==1
+// off for non testing
+#endif // NEWDIRVELMOTEST==1
+#endif // LBM_INCLUDE_TESTSOLVERS!=1
/******************************************************************************
@@ -180,65 +214,114 @@
// target set
#define TSET(l) mLevel[(l)].setOther
+// handle mov. obj
+#if FSGR_STRICT_DEBUG==1
+
+#define LBMDS_ADDMOV(linv,l) \
+ \
+ if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \
+ \
+ LbmFloat dte=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)-(mSimulationTime+this->mpParam->getTimestep()); \
+ if( ABS(dte)< 1e-15 ) { \
+ m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
+ } else { \
+ const int sdx = i+this->dfVecX[linv], sdy = j+this->dfVecY[linv], sdz = k+this->dfVecZ[linv]; \
+ \
+ errMsg("INVALID_MOV_OBJ_TIME"," at "<<PRINT_IJK<<" from l"<<l<<" "<<PRINT_VEC(sdx,sdy,sdz)<<" t="<<(mSimulationTime+this->mpParam->getTimestep())<<" ct="<<QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)<<" dte="<<dte); \
+ debugMarkCell(lev,sdx,sdy,sdz); \
+ } \
+ } \
+
+
+
+#else // FSGR_STRICT_DEBUG==1
+
+#define LBMDS_ADDMOV(linv,l) \
+ \
+ \
+ if((nbflag[linv]&CFBndMoving)&&(!(nbflag[l]&CFBnd))){ \
+ \
+ m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
+ } \
+
+
+
+#endif // !FSGR_STRICT_DEBUG==1
+
// treatment of freeslip reflection
// used both for OPT and nonOPT
-#define DEFAULT_STREAM_FREESLIP(l,invl,mnbf) \
- /*const int inv_l = this->dfInv[l];*/ \
- int nb1 = 0, nb2 = 0; /* is neighbor in this direction an obstacle? */\
- LbmFloat newval = 0.0; /* new value for m[l], differs for free/part slip */\
- const int dx = this->dfVecX[invl], dy = this->dfVecY[invl], dz = this->dfVecZ[invl]; \
- if(dz==0) { \
- nb1 = !(RFLAG(lev, i, j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
- nb2 = !(RFLAG(lev, i+dx,j, k, SRCS(lev))&(CFFluid|CFInter)); \
- if((nb1)&&(!nb2)) { \
- /* x reflection */\
- newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
- } else \
- if((!nb1)&&(nb2)) { \
- /* y reflection */\
- newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
- } else { \
- /* normal no slip in all other cases */\
- newval = QCELL(lev, i,j,k,SRCS(lev), invl); \
- } \
- } else /* z=0 */\
- if(dy==0) { \
- nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
- nb2 = !(RFLAG(lev, i+dx,j,k, SRCS(lev))&(CFFluid|CFInter)); \
- if((nb1)&&(!nb2)) { \
- /* x reflection */\
- newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
- } else \
- if((!nb1)&&(nb2)) { \
- /* z reflection */\
- newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
- } else { \
- /* normal no slip in all other cases */\
- newval = ( QCELL(lev, i,j,k,SRCS(lev), invl) ); \
- } \
- /* end y=0 */ \
- } else { \
- /* x=0 */\
- nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
- nb2 = !(RFLAG(lev, i,j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
- if((nb1)&&(!nb2)) { \
- /* y reflection */\
- newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
- } else \
- if((!nb1)&&(nb2)) { \
- /* z reflection */\
- newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
- } else { \
- /* normal no slip in all other cases */\
- newval = ( QCELL(lev, i,j,k,SRCS(lev), invl) ); \
- } \
- } \
- if(mnbf & CFBndPartslip) { /* part slip interpolation */ \
- const LbmFloat partv = mObjectPartslips[(int)(mnbf>>24)]; \
- m[l] = RAC(ccel, this->dfInv[l] ) * partv + newval * (1.0-partv); /* part slip */ \
- } else {\
- m[l] = newval; /* normal free slip*/\
- }\
+#define DEFAULT_STREAM_FREESLIP(l,invl,mnbf) \
+ \
+ int nb1 = 0, nb2 = 0; \
+ LbmFloat newval = 0.0; \
+ const int dx = this->dfVecX[invl], dy = this->dfVecY[invl], dz = this->dfVecZ[invl]; \
+ \
+ \
+ \
+ const LbmFloat movadd = ( \
+ ((nbflag[invl]&CFBndMoving)&&(!(nbflag[l]&CFBnd))) ? \
+ (QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l)) : 0.); \
+ \
+ if(dz==0) { \
+ nb1 = !(RFLAG(lev, i, j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
+ nb2 = !(RFLAG(lev, i+dx,j, k, SRCS(lev))&(CFFluid|CFInter)); \
+ if((nb1)&&(!nb2)) { \
+ \
+ newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
+ } else \
+ if((!nb1)&&(nb2)) { \
+ \
+ newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
+ } else { \
+ \
+ newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \
+ } \
+ } else \
+ if(dy==0) { \
+ nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
+ nb2 = !(RFLAG(lev, i+dx,j,k, SRCS(lev))&(CFFluid|CFInter)); \
+ if((nb1)&&(!nb2)) { \
+ \
+ newval = QCELL(lev, i+dx,j,k,SRCS(lev), this->dfRefX[l]); \
+ } else \
+ if((!nb1)&&(nb2)) { \
+ \
+ newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
+ } else { \
+ \
+ newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \
+ } \
+ \
+ } else \
+ \
+ { \
+ \
+ nb1 = !(RFLAG(lev, i,j,k+dz, SRCS(lev))&(CFFluid|CFInter)); \
+ nb2 = !(RFLAG(lev, i,j+dy,k, SRCS(lev))&(CFFluid|CFInter)); \
+ if((nb1)&&(!nb2)) { \
+ \
+ newval = QCELL(lev, i,j+dy,k,SRCS(lev), this->dfRefY[l]); \
+ } else \
+ if((!nb1)&&(nb2)) { \
+ \
+ newval = QCELL(lev, i,j,k+dz,SRCS(lev), this->dfRefZ[l]); \
+ } else { \
+ \
+ newval = RAC(ccel, this->dfInv[l] ) +movadd /* */; \
+ } \
+ } \
+ \
+ if(mnbf & CFBndPartslip) { \
+ const LbmFloat partv = mObjectPartslips[(int)(mnbf>>24)]; \
+ \
+ m[l] = (RAC(ccel, this->dfInv[l] ) +movadd /* d *(1./1.) */ ) * partv + newval * (1.0-partv); \
+ } else { \
+ m[l] = newval; \
+ } \
+ \
+
+
+
// complete default stream&collide, 2d/3d
/* read distribution funtions of adjacent cells = sweep step */
@@ -248,16 +331,18 @@
#define MARKCELLCHECK \
debugMarkCell(lev,i,j,k); CAUSE_PANIC;
#define STREAMCHECK(id,ni,nj,nk,nl) \
- if((m[nl] < -1.0) || (m[nl]>1.0)) {\
+ if((!(m[nl] > -1.0) && (m[nl]<1.0)) ) {\
errMsg("STREAMCHECK","ID"<<id<<" Invalid streamed DF nl"<<nl<<" value:"<<m[nl]<<" at "<<PRINT_IJK<<" from "<<PRINT_VEC(ni,nj,nk)<<" nl"<<(nl)<<\
" nfc"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setCurr)<<" nfo"<< RFLAG(lev, ni,nj,nk, mLevel[lev].setOther) ); \
/*FORDF0{ errMsg("STREAMCHECK"," at "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); } */ \
MARKCELLCHECK; \
+ m[nl] = dfEquil[nl]; /* REPAIR */ \
}
#define COLLCHECK \
if( (rho>2.0) || (rho<-1.0) || (ABS(ux)>1.0) || (ABS(uy)>1.0) |(ABS(uz)>1.0) ) {\
errMsg("COLLCHECK","Invalid collision values r:"<<rho<<" u:"PRINT_VEC(ux,uy,uz)<<" at? "<<PRINT_IJK ); \
/*FORDF0{ errMsg("COLLCHECK"," at? "<<PRINT_IJK<<" df "<<l<<"="<<m[l] ); }*/ \
+ rho=ux=uy=uz= 0.; /* REPAIR */ \
MARKCELLCHECK; \
}
#else
@@ -266,455 +351,444 @@
#endif
// careful ux,uy,uz need to be inited before!
+#define DEFAULT_STREAM \
+ m[dC] = RAC(ccel,dC); \
+ STREAMCHECK(1, i,j,k, dC); \
+ FORDF1 { \
+ CellFlagType nbf = nbflag[ this->dfInv[l] ]; \
+ if(nbf & CFBnd) { \
+ if(nbf & CFBndNoslip) { \
+ \
+ m[l] = RAC(ccel, this->dfInv[l] ); \
+ LBMDS_ADDMOV(this->dfInv[l],l); \
+ STREAMCHECK(2, i,j,k, l); \
+ } else if(nbf & (CFBndFreeslip|CFBndPartslip)) { \
+ \
+ if(l<=LBMDIM*2) { \
+ m[l] = RAC(ccel, this->dfInv[l] ); STREAMCHECK(3, i,j,k, l); \
+ LBMDS_ADDMOV(this->dfInv[l],l); \
+ } else { \
+ const int inv_l = this->dfInv[l]; \
+ DEFAULT_STREAM_FREESLIP(l,inv_l,nbf); \
+ } \
+ \
+ } \
+ else { \
+ errMsg("LbmFsgrSolver","Invalid Bnd type at "<<PRINT_IJK<<" f"<<convertCellFlagType2String(nbf)<<",nbdir"<<this->dfInv[l] ); \
+ } \
+ } else { \
+ m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
+ STREAMCHECK(4, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
+ } \
+ } \
+
-#define DEFAULT_STREAM \
- m[dC] = RAC(ccel,dC); \
- STREAMCHECK(1, i,j,k, dC); \
- FORDF1 { \
- CellFlagType nbf = nbflag[ this->dfInv[l] ];\
- if(nbf & CFBnd) { \
- if(nbf & CFBndNoslip) { \
- /* no slip, default */ \
- m[l] = RAC(ccel, this->dfInv[l] ); /* noslip */ \
- STREAMCHECK(2, i,j,k, l); \
- } else if(nbf & (CFBndFreeslip|CFBndPartslip)) { \
- /* free slip */ \
- if(l<=LBMDIM*2) { \
- m[l] = RAC(ccel, this->dfInv[l] ); STREAMCHECK(3, i,j,k, l); /* noslip for <dim*2 */ \
- } else { \
- const int inv_l = this->dfInv[l]; \
- DEFAULT_STREAM_FREESLIP(l,inv_l,nbf); \
- } /* l>2*dim free slip */ \
- \
- } /* type reflect */\
- else {\
- errMsg("LbmFsgrSolver","Invalid Bnd type at "<<PRINT_IJK<<" f"<<convertCellFlagType2String(nbf)<<",nbdir"<<this->dfInv[l] ); \
- } \
- if(nbf&CFBndMoving) m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); /* obs. speed*/ \
- } else { \
- m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); \
- STREAMCHECK(4, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
- } \
- }
// careful ux,uy,uz need to be inited before!
-#define DEFAULT_COLLIDEG(grav) \
- this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), grav, mLevel[lev].lcsmago, &mDebugOmegaRet, &lcsmqo ); \
- CSMOMEGA_STATS(lev,mDebugOmegaRet); \
- FORDF0 { RAC(tcel,l) = m[l]; } \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- COLLCHECK;
-#define OPTIMIZED_STREAMCOLLIDE \
- m[0] = RAC(ccel,0); \
- FORDF1 { /* df0 is set later on... */ \
- /* FIXME CHECK INV ? */\
- if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC; \
- } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
- STREAMCHECK(8, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
- } \
- rho=m[0]; \
- /* ux = mLevel[lev].gravity[0]; uy = [1]; uz = mLevel[lev].gravity[2]; */ \
- this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo ); \
- CSMOMEGA_STATS(lev,mDebugOmegaRet); \
- FORDF0 { RAC(tcel,l) = m[l]; } \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- COLLCHECK;
+#define DEFAULT_COLLIDEG(grav) \
+ this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), grav, mLevel[lev].lcsmago, &mDebugOmegaRet, &lcsmqo ); \
+ CSMOMEGA_STATS(lev,mDebugOmegaRet); \
+ FORDF0 { RAC(tcel,l) = m[l]; } \
+ usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+ COLLCHECK; \
-#else // 3D, opt OPT3D==true
-// handle mov. obj
-#if FSGR_STRICT_DEBUG==1
-#define LBMDS_ADDMOV(linv,l) if(nbflag[linv]&CFBndMoving){ \
- if(QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)== (mSimulationTime+this->mpParam->getTimestep()) ) { \
- m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); /* obs. speed*/ \
- } else { \
- CAUSE_PANIC; errMsg("INVALID_MOV_OBJ_TIME"," at "<<PRINT_IJK<<" from l"<<l<<" t="<<mSimulationTime<<" ct="<<QCELL_NBINV(lev, i, j, k, SRCS(lev), l,dFlux)); \
- } \
- }
-#else // FSGR_STRICT_DEBUG==1
-#define LBMDS_ADDMOV(linv,l) if(nbflag[linv]&CFBndMoving){ m[l]+=QCELL_NBINV(lev, i, j, k, SRCS(lev), l,l); } /* obs. speed*/
-#endif // FSGR_STRICT_DEBUG==1
+#define OPTIMIZED_STREAMCOLLIDE \
+ m[0] = RAC(ccel,0); \
+ FORDF1 { \
+ \
+ if(RFLAG_NBINV(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC; \
+ } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
+ STREAMCHECK(8, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
+ } \
+ rho=m[0]; \
+ DEFAULT_COLLIDEG(mLevel[lev].gravity) \
+
+
+
+#define OPTIMIZED_STREAMCOLLIDE___UNUSED \
+ \
+ this->collideArrays(lev, i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo ); \
+ CSMOMEGA_STATS(lev,mDebugOmegaRet); \
+ FORDF0 { RAC(tcel,l) = m[l]; } \
+ usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+ COLLCHECK; \
+
+
+
+#else // 3D, opt OPT3D==true
+
// default stream opt3d add moving bc val
-#define DEFAULT_STREAM \
- m[dC] = RAC(ccel,dC); \
- /* explicit streaming */ \
- if((!nbored & CFBnd)) { \
- /* no boundary near?, no real speed diff.? */ \
- m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
- m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
- m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
- m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
- m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
- m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
- } else { \
- /* explicit streaming, normal velocity always zero for obstacles */ \
- if(nbflag[dS ]&CFBnd) { m[dN ] = RAC(ccel,dS ); LBMDS_ADDMOV(dS ,dN ); } else { m[dN ] = CSRC_N ; } \
- if(nbflag[dN ]&CFBnd) { m[dS ] = RAC(ccel,dN ); LBMDS_ADDMOV(dN ,dS ); } else { m[dS ] = CSRC_S ; } \
- if(nbflag[dW ]&CFBnd) { m[dE ] = RAC(ccel,dW ); LBMDS_ADDMOV(dW ,dE ); } else { m[dE ] = CSRC_E ; } \
- if(nbflag[dE ]&CFBnd) { m[dW ] = RAC(ccel,dE ); LBMDS_ADDMOV(dE ,dW ); } else { m[dW ] = CSRC_W ; } \
- if(nbflag[dB ]&CFBnd) { m[dT ] = RAC(ccel,dB ); LBMDS_ADDMOV(dB ,dT ); } else { m[dT ] = CSRC_T ; } \
- if(nbflag[dT ]&CFBnd) { m[dB ] = RAC(ccel,dT ); LBMDS_ADDMOV(dT ,dB ); } else { m[dB ] = CSRC_B ; } \
- \
- /* also treat free slip here */ \
- if(nbflag[dSW]&CFBnd) { if(nbflag[dSW]&CFBndNoslip){ m[dNE] = RAC(ccel,dSW); LBMDS_ADDMOV(dSW,dNE); }else{ DEFAULT_STREAM_FREESLIP(dNE,dSW,nbflag[dSW]);} LBMDS_ADDMOV(dSW,dNE); } else { m[dNE] = CSRC_NE; } \
- if(nbflag[dSE]&CFBnd) { if(nbflag[dSE]&CFBndNoslip){ m[dNW] = RAC(ccel,dSE); LBMDS_ADDMOV(dSE,dNW); }else{ DEFAULT_STREAM_FREESLIP(dNW,dSE,nbflag[dSE]);} LBMDS_ADDMOV(dSE,dNW); } else { m[dNW] = CSRC_NW; } \
- if(nbflag[dNW]&CFBnd) { if(nbflag[dNW]&CFBndNoslip){ m[dSE] = RAC(ccel,dNW); LBMDS_ADDMOV(dNW,dSE); }else{ DEFAULT_STREAM_FREESLIP(dSE,dNW,nbflag[dNW]);} LBMDS_ADDMOV(dNW,dSE); } else { m[dSE] = CSRC_SE; } \
- if(nbflag[dNE]&CFBnd) { if(nbflag[dNE]&CFBndNoslip){ m[dSW] = RAC(ccel,dNE); LBMDS_ADDMOV(dNE,dSW); }else{ DEFAULT_STREAM_FREESLIP(dSW,dNE,nbflag[dNE]);} LBMDS_ADDMOV(dNE,dSW); } else { m[dSW] = CSRC_SW; } \
- if(nbflag[dSB]&CFBnd) { if(nbflag[dSB]&CFBndNoslip){ m[dNT] = RAC(ccel,dSB); LBMDS_ADDMOV(dSB,dNT); }else{ DEFAULT_STREAM_FREESLIP(dNT,dSB,nbflag[dSB]);} LBMDS_ADDMOV(dSB,dNT); } else { m[dNT] = CSRC_NT; } \
- if(nbflag[dST]&CFBnd) { if(nbflag[dST]&CFBndNoslip){ m[dNB] = RAC(ccel,dST); LBMDS_ADDMOV(dST,dNB); }else{ DEFAULT_STREAM_FREESLIP(dNB,dST,nbflag[dST]);} LBMDS_ADDMOV(dST,dNB); } else { m[dNB] = CSRC_NB; } \
- if(nbflag[dNB]&CFBnd) { if(nbflag[dNB]&CFBndNoslip){ m[dST] = RAC(ccel,dNB); LBMDS_ADDMOV(dNB,dST); }else{ DEFAULT_STREAM_FREESLIP(dST,dNB,nbflag[dNB]);} LBMDS_ADDMOV(dNB,dST); } else { m[dST] = CSRC_ST; } \
- if(nbflag[dNT]&CFBnd) { if(nbflag[dNT]&CFBndNoslip){ m[dSB] = RAC(ccel,dNT); LBMDS_ADDMOV(dNT,dSB); }else{ DEFAULT_STREAM_FREESLIP(dSB,dNT,nbflag[dNT]);} LBMDS_ADDMOV(dNT,dSB); } else { m[dSB] = CSRC_SB; } \
- if(nbflag[dWB]&CFBnd) { if(nbflag[dWB]&CFBndNoslip){ m[dET] = RAC(ccel,dWB); LBMDS_ADDMOV(dWB,dET); }else{ DEFAULT_STREAM_FREESLIP(dET,dWB,nbflag[dWB]);} LBMDS_ADDMOV(dWB,dET); } else { m[dET] = CSRC_ET; } \
- if(nbflag[dWT]&CFBnd) { if(nbflag[dWT]&CFBndNoslip){ m[dEB] = RAC(ccel,dWT); LBMDS_ADDMOV(dWT,dEB); }else{ DEFAULT_STREAM_FREESLIP(dEB,dWT,nbflag[dWT]);} LBMDS_ADDMOV(dWT,dEB); } else { m[dEB] = CSRC_EB; } \
- if(nbflag[dEB]&CFBnd) { if(nbflag[dEB]&CFBndNoslip){ m[dWT] = RAC(ccel,dEB); LBMDS_ADDMOV(dEB,dWT); }else{ DEFAULT_STREAM_FREESLIP(dWT,dEB,nbflag[dEB]);} LBMDS_ADDMOV(dEB,dWT); } else { m[dWT] = CSRC_WT; } \
- if(nbflag[dET]&CFBnd) { if(nbflag[dET]&CFBndNoslip){ m[dWB] = RAC(ccel,dET); LBMDS_ADDMOV(dET,dWB); }else{ DEFAULT_STREAM_FREESLIP(dWB,dET,nbflag[dET]);} LBMDS_ADDMOV(dET,dWB); } else { m[dWB] = CSRC_WB; } \
- }
-
-
-
-#define COLL_CALCULATE_DFEQ(dstarray) \
- dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \
- dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \
- dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \
- dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \
- dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \
- dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB;
-#define COLL_CALCULATE_NONEQTENSOR(csolev, srcArray ) \
- lcsmqadd = (srcArray##NE - lcsmeq[ dNE ]); \
- lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \
- lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \
- lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
- lcsmqo = (lcsmqadd* lcsmqadd); \
- lcsmqadd = (srcArray##ET - lcsmeq[ dET ]); \
- lcsmqadd -= (srcArray##EB - lcsmeq[ dEB ]); \
- lcsmqadd -= (srcArray##WT - lcsmeq[ dWT ]); \
- lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \
- lcsmqo += (lcsmqadd* lcsmqadd); \
- lcsmqadd = (srcArray##NT - lcsmeq[ dNT ]); \
- lcsmqadd -= (srcArray##NB - lcsmeq[ dNB ]); \
- lcsmqadd -= (srcArray##ST - lcsmeq[ dST ]); \
- lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \
- lcsmqo += (lcsmqadd* lcsmqadd); \
- lcsmqo *= 2.0; \
- lcsmqadd = (srcArray##E - lcsmeq[ dE ]); \
- lcsmqadd += (srcArray##W - lcsmeq[ dW ]); \
- lcsmqadd += (srcArray##NE - lcsmeq[ dNE ]); \
- lcsmqadd += (srcArray##NW - lcsmeq[ dNW ]); \
- lcsmqadd += (srcArray##SE - lcsmeq[ dSE ]); \
- lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
- lcsmqadd += (srcArray##ET - lcsmeq[ dET ]); \
- lcsmqadd += (srcArray##EB - lcsmeq[ dEB ]); \
- lcsmqadd += (srcArray##WT - lcsmeq[ dWT ]); \
- lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \
- lcsmqo += (lcsmqadd* lcsmqadd); \
- lcsmqadd = (srcArray##N - lcsmeq[ dN ]); \
- lcsmqadd += (srcArray##S - lcsmeq[ dS ]); \
- lcsmqadd += (srcArray##NE - lcsmeq[ dNE ]); \
- lcsmqadd += (srcArray##NW - lcsmeq[ dNW ]); \
- lcsmqadd += (srcArray##SE - lcsmeq[ dSE ]); \
- lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
- lcsmqadd += (srcArray##NT - lcsmeq[ dNT ]); \
- lcsmqadd += (srcArray##NB - lcsmeq[ dNB ]); \
- lcsmqadd += (srcArray##ST - lcsmeq[ dST ]); \
- lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \
- lcsmqo += (lcsmqadd* lcsmqadd); \
- lcsmqadd = (srcArray##T - lcsmeq[ dT ]); \
- lcsmqadd += (srcArray##B - lcsmeq[ dB ]); \
- lcsmqadd += (srcArray##NT - lcsmeq[ dNT ]); \
- lcsmqadd += (srcArray##NB - lcsmeq[ dNB ]); \
- lcsmqadd += (srcArray##ST - lcsmeq[ dST ]); \
- lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \
- lcsmqadd += (srcArray##ET - lcsmeq[ dET ]); \
- lcsmqadd += (srcArray##EB - lcsmeq[ dEB ]); \
- lcsmqadd += (srcArray##WT - lcsmeq[ dWT ]); \
- lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \
- lcsmqo += (lcsmqadd* lcsmqadd); \
- lcsmqo = sqrt(lcsmqo); /* FIXME check effect of sqrt*/ \
+#define DEFAULT_STREAM \
+ m[dC] = RAC(ccel,dC); \
+ \
+ if((!nbored & CFBnd)) { \
+ \
+ m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
+ m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
+ m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
+ m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
+ m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
+ m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
+ } else { \
+ \
+ if(nbflag[dS ]&CFBnd) { m[dN ] = RAC(ccel,dS ); LBMDS_ADDMOV(dS ,dN ); } else { m[dN ] = CSRC_N ; } \
+ if(nbflag[dN ]&CFBnd) { m[dS ] = RAC(ccel,dN ); LBMDS_ADDMOV(dN ,dS ); } else { m[dS ] = CSRC_S ; } \
+ if(nbflag[dW ]&CFBnd) { m[dE ] = RAC(ccel,dW ); LBMDS_ADDMOV(dW ,dE ); } else { m[dE ] = CSRC_E ; } \
+ if(nbflag[dE ]&CFBnd) { m[dW ] = RAC(ccel,dE ); LBMDS_ADDMOV(dE ,dW ); } else { m[dW ] = CSRC_W ; } \
+ if(nbflag[dB ]&CFBnd) { m[dT ] = RAC(ccel,dB ); LBMDS_ADDMOV(dB ,dT ); } else { m[dT ] = CSRC_T ; } \
+ if(nbflag[dT ]&CFBnd) { m[dB ] = RAC(ccel,dT ); LBMDS_ADDMOV(dT ,dB ); } else { m[dB ] = CSRC_B ; } \
+ \
+ \
+ if(nbflag[dSW]&CFBnd) { if(nbflag[dSW]&CFBndNoslip){ m[dNE] = RAC(ccel,dSW); LBMDS_ADDMOV(dSW,dNE); }else{ DEFAULT_STREAM_FREESLIP(dNE,dSW,nbflag[dSW]);} } else { m[dNE] = CSRC_NE; } \
+ if(nbflag[dSE]&CFBnd) { if(nbflag[dSE]&CFBndNoslip){ m[dNW] = RAC(ccel,dSE); LBMDS_ADDMOV(dSE,dNW); }else{ DEFAULT_STREAM_FREESLIP(dNW,dSE,nbflag[dSE]);} } else { m[dNW] = CSRC_NW; } \
+ if(nbflag[dNW]&CFBnd) { if(nbflag[dNW]&CFBndNoslip){ m[dSE] = RAC(ccel,dNW); LBMDS_ADDMOV(dNW,dSE); }else{ DEFAULT_STREAM_FREESLIP(dSE,dNW,nbflag[dNW]);} } else { m[dSE] = CSRC_SE; } \
+ if(nbflag[dNE]&CFBnd) { if(nbflag[dNE]&CFBndNoslip){ m[dSW] = RAC(ccel,dNE); LBMDS_ADDMOV(dNE,dSW); }else{ DEFAULT_STREAM_FREESLIP(dSW,dNE,nbflag[dNE]);} } else { m[dSW] = CSRC_SW; } \
+ if(nbflag[dSB]&CFBnd) { if(nbflag[dSB]&CFBndNoslip){ m[dNT] = RAC(ccel,dSB); LBMDS_ADDMOV(dSB,dNT); }else{ DEFAULT_STREAM_FREESLIP(dNT,dSB,nbflag[dSB]);} } else { m[dNT] = CSRC_NT; } \
+ if(nbflag[dST]&CFBnd) { if(nbflag[dST]&CFBndNoslip){ m[dNB] = RAC(ccel,dST); LBMDS_ADDMOV(dST,dNB); }else{ DEFAULT_STREAM_FREESLIP(dNB,dST,nbflag[dST]);} } else { m[dNB] = CSRC_NB; } \
+ if(nbflag[dNB]&CFBnd) { if(nbflag[dNB]&CFBndNoslip){ m[dST] = RAC(ccel,dNB); LBMDS_ADDMOV(dNB,dST); }else{ DEFAULT_STREAM_FREESLIP(dST,dNB,nbflag[dNB]);} } else { m[dST] = CSRC_ST; } \
+ if(nbflag[dNT]&CFBnd) { if(nbflag[dNT]&CFBndNoslip){ m[dSB] = RAC(ccel,dNT); LBMDS_ADDMOV(dNT,dSB); }else{ DEFAULT_STREAM_FREESLIP(dSB,dNT,nbflag[dNT]);} } else { m[dSB] = CSRC_SB; } \
+ if(nbflag[dWB]&CFBnd) { if(nbflag[dWB]&CFBndNoslip){ m[dET] = RAC(ccel,dWB); LBMDS_ADDMOV(dWB,dET); }else{ DEFAULT_STREAM_FREESLIP(dET,dWB,nbflag[dWB]);} } else { m[dET] = CSRC_ET; } \
+ if(nbflag[dWT]&CFBnd) { if(nbflag[dWT]&CFBndNoslip){ m[dEB] = RAC(ccel,dWT); LBMDS_ADDMOV(dWT,dEB); }else{ DEFAULT_STREAM_FREESLIP(dEB,dWT,nbflag[dWT]);} } else { m[dEB] = CSRC_EB; } \
+ if(nbflag[dEB]&CFBnd) { if(nbflag[dEB]&CFBndNoslip){ m[dWT] = RAC(ccel,dEB); LBMDS_ADDMOV(dEB,dWT); }else{ DEFAULT_STREAM_FREESLIP(dWT,dEB,nbflag[dEB]);} } else { m[dWT] = CSRC_WT; } \
+ if(nbflag[dET]&CFBnd) { if(nbflag[dET]&CFBndNoslip){ m[dWB] = RAC(ccel,dET); LBMDS_ADDMOV(dET,dWB); }else{ DEFAULT_STREAM_FREESLIP(dWB,dET,nbflag[dET]);} } else { m[dWB] = CSRC_WB; } \
+ } \
+
+
+
+
+
+#define COLL_CALCULATE_DFEQ(dstarray) \
+ dstarray[dN ] = EQN ; dstarray[dS ] = EQS ; \
+ dstarray[dE ] = EQE ; dstarray[dW ] = EQW ; \
+ dstarray[dT ] = EQT ; dstarray[dB ] = EQB ; \
+ dstarray[dNE] = EQNE; dstarray[dNW] = EQNW; dstarray[dSE] = EQSE; dstarray[dSW] = EQSW; \
+ dstarray[dNT] = EQNT; dstarray[dNB] = EQNB; dstarray[dST] = EQST; dstarray[dSB] = EQSB; \
+ dstarray[dET] = EQET; dstarray[dEB] = EQEB; dstarray[dWT] = EQWT; dstarray[dWB] = EQWB; \
+
+
+
+#define COLL_CALCULATE_NONEQTENSOR(csolev, srcArray ) \
+ lcsmqadd = (srcArray##NE - lcsmeq[ dNE ]); \
+ lcsmqadd -= (srcArray##NW - lcsmeq[ dNW ]); \
+ lcsmqadd -= (srcArray##SE - lcsmeq[ dSE ]); \
+ lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
+ lcsmqo = (lcsmqadd* lcsmqadd); \
+ lcsmqadd = (srcArray##ET - lcsmeq[ dET ]); \
+ lcsmqadd -= (srcArray##EB - lcsmeq[ dEB ]); \
+ lcsmqadd -= (srcArray##WT - lcsmeq[ dWT ]); \
+ lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \
+ lcsmqo += (lcsmqadd* lcsmqadd); \
+ lcsmqadd = (srcArray##NT - lcsmeq[ dNT ]); \
+ lcsmqadd -= (srcArray##NB - lcsmeq[ dNB ]); \
+ lcsmqadd -= (srcArray##ST - lcsmeq[ dST ]); \
+ lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \
+ lcsmqo += (lcsmqadd* lcsmqadd); \
+ lcsmqo *= 2.0; \
+ lcsmqadd = (srcArray##E - lcsmeq[ dE ]); \
+ lcsmqadd += (srcArray##W - lcsmeq[ dW ]); \
+ lcsmqadd += (srcArray##NE - lcsmeq[ dNE ]); \
+ lcsmqadd += (srcArray##NW - lcsmeq[ dNW ]); \
+ lcsmqadd += (srcArray##SE - lcsmeq[ dSE ]); \
+ lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
+ lcsmqadd += (srcArray##ET - lcsmeq[ dET ]); \
+ lcsmqadd += (srcArray##EB - lcsmeq[ dEB ]); \
+ lcsmqadd += (srcArray##WT - lcsmeq[ dWT ]); \
+ lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \
+ lcsmqo += (lcsmqadd* lcsmqadd); \
+ lcsmqadd = (srcArray##N - lcsmeq[ dN ]); \
+ lcsmqadd += (srcArray##S - lcsmeq[ dS ]); \
+ lcsmqadd += (srcArray##NE - lcsmeq[ dNE ]); \
+ lcsmqadd += (srcArray##NW - lcsmeq[ dNW ]); \
+ lcsmqadd += (srcArray##SE - lcsmeq[ dSE ]); \
+ lcsmqadd += (srcArray##SW - lcsmeq[ dSW ]); \
+ lcsmqadd += (srcArray##NT - lcsmeq[ dNT ]); \
+ lcsmqadd += (srcArray##NB - lcsmeq[ dNB ]); \
+ lcsmqadd += (srcArray##ST - lcsmeq[ dST ]); \
+ lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \
+ lcsmqo += (lcsmqadd* lcsmqadd); \
+ lcsmqadd = (srcArray##T - lcsmeq[ dT ]); \
+ lcsmqadd += (srcArray##B - lcsmeq[ dB ]); \
+ lcsmqadd += (srcArray##NT - lcsmeq[ dNT ]); \
+ lcsmqadd += (srcArray##NB - lcsmeq[ dNB ]); \
+ lcsmqadd += (srcArray##ST - lcsmeq[ dST ]); \
+ lcsmqadd += (srcArray##SB - lcsmeq[ dSB ]); \
+ lcsmqadd += (srcArray##ET - lcsmeq[ dET ]); \
+ lcsmqadd += (srcArray##EB - lcsmeq[ dEB ]); \
+ lcsmqadd += (srcArray##WT - lcsmeq[ dWT ]); \
+ lcsmqadd += (srcArray##WB - lcsmeq[ dWB ]); \
+ lcsmqo += (lcsmqadd* lcsmqadd); \
+ lcsmqo = sqrt(lcsmqo); \
+
+
// COLL_CALCULATE_CSMOMEGAVAL(csolev, lcsmomega);
// careful - need lcsmqo
-#define COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega ) \
- dstomega = 1.0/\
- ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*(\
- -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \
- / (6.0*mLevel[(csolev)].lcsmago_sqr)) \
- ) +0.5 );
-
-#define DEFAULT_COLLIDE_LES(grav) \
- rho = + MSRC_C + MSRC_N \
- + MSRC_S + MSRC_E \
- + MSRC_W + MSRC_T \
- + MSRC_B + MSRC_NE \
- + MSRC_NW + MSRC_SE \
- + MSRC_SW + MSRC_NT \
- + MSRC_NB + MSRC_ST \
- + MSRC_SB + MSRC_ET \
- + MSRC_EB + MSRC_WT \
- + MSRC_WB; \
- \
- ux = MSRC_E - MSRC_W \
- + MSRC_NE - MSRC_NW \
- + MSRC_SE - MSRC_SW \
- + MSRC_ET + MSRC_EB \
- - MSRC_WT - MSRC_WB ; \
- \
- uy = MSRC_N - MSRC_S \
- + MSRC_NE + MSRC_NW \
- - MSRC_SE - MSRC_SW \
- + MSRC_NT + MSRC_NB \
- - MSRC_ST - MSRC_SB ; \
- \
- uz = MSRC_T - MSRC_B \
- + MSRC_NT - MSRC_NB \
- + MSRC_ST - MSRC_SB \
- + MSRC_ET - MSRC_EB \
- + MSRC_WT - MSRC_WB ; \
- PRECOLLIDE_MODS(rho,ux,uy,uz, grav); \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- COLL_CALCULATE_DFEQ(lcsmeq); \
- COLL_CALCULATE_NONEQTENSOR(lev, MSRC_); \
- COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
- CSMOMEGA_STATS(lev,lcsmomega); \
- \
- RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C + lcsmomega*EQC ; \
- \
- RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N + lcsmomega*lcsmeq[ dN ]; \
- RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S + lcsmomega*lcsmeq[ dS ]; \
- RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E + lcsmomega*lcsmeq[ dE ]; \
- RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W + lcsmomega*lcsmeq[ dW ]; \
- RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T + lcsmomega*lcsmeq[ dT ]; \
- RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B + lcsmomega*lcsmeq[ dB ]; \
- \
- RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
- RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
- RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
- RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
- RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
- RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
- RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
- RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
- RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
- RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
- RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
- RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB];
-
-#define DEFAULT_COLLIDE_NOLES(grav) \
- rho = + MSRC_C + MSRC_N \
- + MSRC_S + MSRC_E \
- + MSRC_W + MSRC_T \
- + MSRC_B + MSRC_NE \
- + MSRC_NW + MSRC_SE \
- + MSRC_SW + MSRC_NT \
- + MSRC_NB + MSRC_ST \
- + MSRC_SB + MSRC_ET \
- + MSRC_EB + MSRC_WT \
- + MSRC_WB; \
- \
- ux = MSRC_E - MSRC_W \
- + MSRC_NE - MSRC_NW \
- + MSRC_SE - MSRC_SW \
- + MSRC_ET + MSRC_EB \
- - MSRC_WT - MSRC_WB ; \
- \
- uy = MSRC_N - MSRC_S \
- + MSRC_NE + MSRC_NW \
- - MSRC_SE - MSRC_SW \
- + MSRC_NT + MSRC_NB \
- - MSRC_ST - MSRC_SB ; \
- \
- uz = MSRC_T - MSRC_B \
- + MSRC_NT - MSRC_NB \
- + MSRC_ST - MSRC_SB \
- + MSRC_ET - MSRC_EB \
- + MSRC_WT - MSRC_WB ; \
- PRECOLLIDE_MODS(rho, ux,uy,uz, grav); \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- \
- RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C + OMEGA(lev)*EQC ; \
- \
- RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N + OMEGA(lev)*EQN ; \
- RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S + OMEGA(lev)*EQS ; \
- RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E + OMEGA(lev)*EQE ; \
- RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W + OMEGA(lev)*EQW ; \
- RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T + OMEGA(lev)*EQT ; \
- RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B + OMEGA(lev)*EQB ; \
- \
- RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \
- RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \
- RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \
- RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \
- RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \
- RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \
- RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \
- RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \
- RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \
- RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \
- RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \
- RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB;
-
-
-
-#define OPTIMIZED_STREAMCOLLIDE_LES \
- /* only surrounded by fluid cells...!, so safe streaming here... */ \
- m[dC ] = CSRC_C ; \
- m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
- m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
- m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
- m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
- m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
- m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
- \
- rho = MSRC_C + MSRC_N + MSRC_S + MSRC_E + MSRC_W + MSRC_T \
- + MSRC_B + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \
- + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \
- ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \
- + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB; \
- uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \
- + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB; \
- uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \
- + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB; \
- PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- COLL_CALCULATE_DFEQ(lcsmeq); \
- COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \
- COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
- CSMOMEGA_STATS(lev,lcsmomega); \
- \
- RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C + lcsmomega*EQC ; \
- RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N + lcsmomega*lcsmeq[ dN ]; \
- RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S + lcsmomega*lcsmeq[ dS ]; \
- RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E + lcsmomega*lcsmeq[ dE ]; \
- RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W + lcsmomega*lcsmeq[ dW ]; \
- RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T + lcsmomega*lcsmeq[ dT ]; \
- RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B + lcsmomega*lcsmeq[ dB ]; \
- \
- RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
- RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
- RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
- RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
- \
- RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
- RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
- RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
- RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
- \
- RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
- RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
- RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
- RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \
-
-#define OPTIMIZED_STREAMCOLLIDE_UNUSED \
- /* only surrounded by fluid cells...!, so safe streaming here... */ \
- rho = CSRC_C + CSRC_N + CSRC_S + CSRC_E + CSRC_W + CSRC_T \
- + CSRC_B + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
- + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
- ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
- + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
- uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
- + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
- uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
- + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
- PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- COLL_CALCULATE_DFEQ(lcsmeq); \
- COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \
- COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
- \
- RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C + lcsmomega*EQC ; \
- RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N + lcsmomega*lcsmeq[ dN ]; \
- RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S + lcsmomega*lcsmeq[ dS ]; \
- RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E + lcsmomega*lcsmeq[ dE ]; \
- RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W + lcsmomega*lcsmeq[ dW ]; \
- RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T + lcsmomega*lcsmeq[ dT ]; \
- RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B + lcsmomega*lcsmeq[ dB ]; \
- \
- RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE]; \
- RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW]; \
- RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE]; \
- RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \
- \
- RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT]; \
- RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB]; \
- RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST]; \
- RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \
- \
- RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET]; \
- RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \
- RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT]; \
- RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB]; \
-
-#define OPTIMIZED_STREAMCOLLIDE_NOLES \
- /* only surrounded by fluid cells...!, so safe streaming here... */ \
- rho = CSRC_C + CSRC_N + CSRC_S + CSRC_E + CSRC_W + CSRC_T \
- + CSRC_B + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
- + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
- ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
- + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
- uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
- + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
- uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
- + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
- PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C + OMEGA(lev)*EQC ; \
- RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N + OMEGA(lev)*EQN ; \
- RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S + OMEGA(lev)*EQS ; \
- RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E + OMEGA(lev)*EQE ; \
- RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W + OMEGA(lev)*EQW ; \
- RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T + OMEGA(lev)*EQT ; \
- RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B + OMEGA(lev)*EQB ; \
- \
- RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE; \
- RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW; \
- RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE; \
- RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \
- \
- RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT; \
- RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB; \
- RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST; \
- RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \
- \
- RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET; \
- RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \
- RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT; \
- RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB; \
-
-
-// debug version1
-#define STREAMCHECK(ni,nj,nk,nl)
-#define COLLCHECK
-#define OPTIMIZED_STREAMCOLLIDE_DEBUG \
- m[0] = RAC(ccel,0); \
- FORDF1 { /* df0 is set later on... */ \
- if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl"); CAUSE_PANIC;\
- } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
- STREAMCHECK(9, i+this->dfVecX[this->dfInv[l]], j+this->dfVecY[this->dfInv[l]],k+this->dfVecZ[this->dfInv[l]], l); \
- } \
- /* rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; */ \
- this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].gravity, mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo ); \
- CSMOMEGA_STATS(lev,mDebugOmegaRet); \
- FORDF0 { RAC(tcel,l) = m[l]; } \
- usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
- COLLCHECK;
-
-
-
-// more debugging
-/*DEBUG \
- m[0] = RAC(ccel,0); \
- FORDF1 { \
- if(RFLAG_NB(lev, i,j,k,SRCS(lev),l)&CFBnd) { errMsg("???", "bnd-err-nobndfl");CAUSE_PANIC; \
- } else { m[l] = QCELL_NBINV(lev, i, j, k, SRCS(lev), l, l); } \
- } \
-errMsg("T","QSDM at %d,%d,%d lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lcsmomega ); \
- rho=m[0]; ux = mLevel[lev].gravity[0]; uy = mLevel[lev].gravity[1]; uz = mLevel[lev].gravity[2]; \
- this->collideArrays(i,j,k, m, rho,ux,uy,uz, OMEGA(lev), mLevel[lev].lcsmago , &mDebugOmegaRet, &lcsmqo ); \
- CSMOMEGA_STATS(lev,mDebugOmegaRet); \
- */
+#define COLL_CALCULATE_CSMOMEGAVAL(csolev, dstomega ) \
+ dstomega = 1.0/ \
+ ( 3.0*( mLevel[(csolev)].lcnu+mLevel[(csolev)].lcsmago_sqr*( \
+ -mLevel[(csolev)].lcnu + sqrt( mLevel[(csolev)].lcnu*mLevel[(csolev)].lcnu + 18.0*mLevel[(csolev)].lcsmago_sqr* lcsmqo ) \
+ / (6.0*mLevel[(csolev)].lcsmago_sqr)) \
+ ) +0.5 ); \
+
+
+
+#define DEFAULT_COLLIDE_LES(grav) \
+ rho = + MSRC_C + MSRC_N \
+ + MSRC_S + MSRC_E \
+ + MSRC_W + MSRC_T \
+ + MSRC_B + MSRC_NE \
+ + MSRC_NW + MSRC_SE \
+ + MSRC_SW + MSRC_NT \
+ + MSRC_NB + MSRC_ST \
+ + MSRC_SB + MSRC_ET \
+ + MSRC_EB + MSRC_WT \
+ + MSRC_WB; \
+ \
+ ux = MSRC_E - MSRC_W \
+ + MSRC_NE - MSRC_NW \
+ + MSRC_SE - MSRC_SW \
+ + MSRC_ET + MSRC_EB \
+ - MSRC_WT - MSRC_WB ; \
+ \
+ uy = MSRC_N - MSRC_S \
+ + MSRC_NE + MSRC_NW \
+ - MSRC_SE - MSRC_SW \
+ + MSRC_NT + MSRC_NB \
+ - MSRC_ST - MSRC_SB ; \
+ \
+ uz = MSRC_T - MSRC_B \
+ + MSRC_NT - MSRC_NB \
+ + MSRC_ST - MSRC_SB \
+ + MSRC_ET - MSRC_EB \
+ + MSRC_WT - MSRC_WB ; \
+ PRECOLLIDE_MODS(rho,ux,uy,uz, grav); \
+ usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+ COLL_CALCULATE_DFEQ(lcsmeq); \
+ COLL_CALCULATE_NONEQTENSOR(lev, MSRC_); \
+ COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
+ CSMOMEGA_STATS(lev,lcsmomega); \
+ \
+ RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C + lcsmomega*EQC ; \
+ \
+ RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N + lcsmomega*lcsmeq[ dN ]; \
+ RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S + lcsmomega*lcsmeq[ dS ]; \
+ RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E + lcsmomega*lcsmeq[ dE ]; \
+ RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W + lcsmomega*lcsmeq[ dW ]; \
+ RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T + lcsmomega*lcsmeq[ dT ]; \
+ RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B + lcsmomega*lcsmeq[ dB ]; \
+ \
+ RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
+ RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
+ RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
+ RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
+ RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
+ RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
+ RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
+ RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
+ RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
+ RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
+ RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
+ RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \
+
+
+
+#define DEFAULT_COLLIDE_NOLES(grav) \
+ rho = + MSRC_C + MSRC_N \
+ + MSRC_S + MSRC_E \
+ + MSRC_W + MSRC_T \
+ + MSRC_B + MSRC_NE \
+ + MSRC_NW + MSRC_SE \
+ + MSRC_SW + MSRC_NT \
+ + MSRC_NB + MSRC_ST \
+ + MSRC_SB + MSRC_ET \
+ + MSRC_EB + MSRC_WT \
+ + MSRC_WB; \
+ \
+ ux = MSRC_E - MSRC_W \
+ + MSRC_NE - MSRC_NW \
+ + MSRC_SE - MSRC_SW \
+ + MSRC_ET + MSRC_EB \
+ - MSRC_WT - MSRC_WB ; \
+ \
+ uy = MSRC_N - MSRC_S \
+ + MSRC_NE + MSRC_NW \
+ - MSRC_SE - MSRC_SW \
+ + MSRC_NT + MSRC_NB \
+ - MSRC_ST - MSRC_SB ; \
+ \
+ uz = MSRC_T - MSRC_B \
+ + MSRC_NT - MSRC_NB \
+ + MSRC_ST - MSRC_SB \
+ + MSRC_ET - MSRC_EB \
+ + MSRC_WT - MSRC_WB ; \
+ PRECOLLIDE_MODS(rho, ux,uy,uz, grav); \
+ usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+ \
+ RAC(tcel,dC ) = (1.0-OMEGA(lev))*MSRC_C + OMEGA(lev)*EQC ; \
+ \
+ RAC(tcel,dN ) = (1.0-OMEGA(lev))*MSRC_N + OMEGA(lev)*EQN ; \
+ RAC(tcel,dS ) = (1.0-OMEGA(lev))*MSRC_S + OMEGA(lev)*EQS ; \
+ RAC(tcel,dE ) = (1.0-OMEGA(lev))*MSRC_E + OMEGA(lev)*EQE ; \
+ RAC(tcel,dW ) = (1.0-OMEGA(lev))*MSRC_W + OMEGA(lev)*EQW ; \
+ RAC(tcel,dT ) = (1.0-OMEGA(lev))*MSRC_T + OMEGA(lev)*EQT ; \
+ RAC(tcel,dB ) = (1.0-OMEGA(lev))*MSRC_B + OMEGA(lev)*EQB ; \
+ \
+ RAC(tcel,dNE) = (1.0-OMEGA(lev))*MSRC_NE + OMEGA(lev)*EQNE; \
+ RAC(tcel,dNW) = (1.0-OMEGA(lev))*MSRC_NW + OMEGA(lev)*EQNW; \
+ RAC(tcel,dSE) = (1.0-OMEGA(lev))*MSRC_SE + OMEGA(lev)*EQSE; \
+ RAC(tcel,dSW) = (1.0-OMEGA(lev))*MSRC_SW + OMEGA(lev)*EQSW; \
+ RAC(tcel,dNT) = (1.0-OMEGA(lev))*MSRC_NT + OMEGA(lev)*EQNT; \
+ RAC(tcel,dNB) = (1.0-OMEGA(lev))*MSRC_NB + OMEGA(lev)*EQNB; \
+ RAC(tcel,dST) = (1.0-OMEGA(lev))*MSRC_ST + OMEGA(lev)*EQST; \
+ RAC(tcel,dSB) = (1.0-OMEGA(lev))*MSRC_SB + OMEGA(lev)*EQSB; \
+ RAC(tcel,dET) = (1.0-OMEGA(lev))*MSRC_ET + OMEGA(lev)*EQET; \
+ RAC(tcel,dEB) = (1.0-OMEGA(lev))*MSRC_EB + OMEGA(lev)*EQEB; \
+ RAC(tcel,dWT) = (1.0-OMEGA(lev))*MSRC_WT + OMEGA(lev)*EQWT; \
+ RAC(tcel,dWB) = (1.0-OMEGA(lev))*MSRC_WB + OMEGA(lev)*EQWB; \
+
+
+
+
+
+#define OPTIMIZED_STREAMCOLLIDE_LES \
+ \
+ m[dC ] = CSRC_C ; \
+ m[dN ] = CSRC_N ; m[dS ] = CSRC_S ; \
+ m[dE ] = CSRC_E ; m[dW ] = CSRC_W ; \
+ m[dT ] = CSRC_T ; m[dB ] = CSRC_B ; \
+ m[dNE] = CSRC_NE; m[dNW] = CSRC_NW; m[dSE] = CSRC_SE; m[dSW] = CSRC_SW; \
+ m[dNT] = CSRC_NT; m[dNB] = CSRC_NB; m[dST] = CSRC_ST; m[dSB] = CSRC_SB; \
+ m[dET] = CSRC_ET; m[dEB] = CSRC_EB; m[dWT] = CSRC_WT; m[dWB] = CSRC_WB; \
+ \
+ rho = MSRC_C + MSRC_N + MSRC_S + MSRC_E + MSRC_W + MSRC_T \
+ + MSRC_B + MSRC_NE + MSRC_NW + MSRC_SE + MSRC_SW + MSRC_NT \
+ + MSRC_NB + MSRC_ST + MSRC_SB + MSRC_ET + MSRC_EB + MSRC_WT + MSRC_WB; \
+ ux = MSRC_E - MSRC_W + MSRC_NE - MSRC_NW + MSRC_SE - MSRC_SW \
+ + MSRC_ET + MSRC_EB - MSRC_WT - MSRC_WB; \
+ uy = MSRC_N - MSRC_S + MSRC_NE + MSRC_NW - MSRC_SE - MSRC_SW \
+ + MSRC_NT + MSRC_NB - MSRC_ST - MSRC_SB; \
+ uz = MSRC_T - MSRC_B + MSRC_NT - MSRC_NB + MSRC_ST - MSRC_SB \
+ + MSRC_ET - MSRC_EB + MSRC_WT - MSRC_WB; \
+ PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
+ usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+ COLL_CALCULATE_DFEQ(lcsmeq); \
+ COLL_CALCULATE_NONEQTENSOR(lev, MSRC_) \
+ COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
+ CSMOMEGA_STATS(lev,lcsmomega); \
+ \
+ RAC(tcel,dC ) = (1.0-lcsmomega)*MSRC_C + lcsmomega*EQC ; \
+ RAC(tcel,dN ) = (1.0-lcsmomega)*MSRC_N + lcsmomega*lcsmeq[ dN ]; \
+ RAC(tcel,dS ) = (1.0-lcsmomega)*MSRC_S + lcsmomega*lcsmeq[ dS ]; \
+ RAC(tcel,dE ) = (1.0-lcsmomega)*MSRC_E + lcsmomega*lcsmeq[ dE ]; \
+ RAC(tcel,dW ) = (1.0-lcsmomega)*MSRC_W + lcsmomega*lcsmeq[ dW ]; \
+ RAC(tcel,dT ) = (1.0-lcsmomega)*MSRC_T + lcsmomega*lcsmeq[ dT ]; \
+ RAC(tcel,dB ) = (1.0-lcsmomega)*MSRC_B + lcsmomega*lcsmeq[ dB ]; \
+ \
+ RAC(tcel,dNE) = (1.0-lcsmomega)*MSRC_NE + lcsmomega*lcsmeq[ dNE]; \
+ RAC(tcel,dNW) = (1.0-lcsmomega)*MSRC_NW + lcsmomega*lcsmeq[ dNW]; \
+ RAC(tcel,dSE) = (1.0-lcsmomega)*MSRC_SE + lcsmomega*lcsmeq[ dSE]; \
+ RAC(tcel,dSW) = (1.0-lcsmomega)*MSRC_SW + lcsmomega*lcsmeq[ dSW]; \
+ \
+ RAC(tcel,dNT) = (1.0-lcsmomega)*MSRC_NT + lcsmomega*lcsmeq[ dNT]; \
+ RAC(tcel,dNB) = (1.0-lcsmomega)*MSRC_NB + lcsmomega*lcsmeq[ dNB]; \
+ RAC(tcel,dST) = (1.0-lcsmomega)*MSRC_ST + lcsmomega*lcsmeq[ dST]; \
+ RAC(tcel,dSB) = (1.0-lcsmomega)*MSRC_SB + lcsmomega*lcsmeq[ dSB]; \
+ \
+ RAC(tcel,dET) = (1.0-lcsmomega)*MSRC_ET + lcsmomega*lcsmeq[ dET]; \
+ RAC(tcel,dEB) = (1.0-lcsmomega)*MSRC_EB + lcsmomega*lcsmeq[ dEB]; \
+ RAC(tcel,dWT) = (1.0-lcsmomega)*MSRC_WT + lcsmomega*lcsmeq[ dWT]; \
+ RAC(tcel,dWB) = (1.0-lcsmomega)*MSRC_WB + lcsmomega*lcsmeq[ dWB]; \
+
+
+
+#define OPTIMIZED_STREAMCOLLIDE_UNUSED \
+ \
+ rho = CSRC_C + CSRC_N + CSRC_S + CSRC_E + CSRC_W + CSRC_T \
+ + CSRC_B + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
+ + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
+ ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
+ + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
+ uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
+ + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
+ uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
+ + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
+ PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
+ usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+ COLL_CALCULATE_DFEQ(lcsmeq); \
+ COLL_CALCULATE_NONEQTENSOR(lev, CSRC_) \
+ COLL_CALCULATE_CSMOMEGAVAL(lev, lcsmomega); \
+ \
+ RAC(tcel,dC ) = (1.0-lcsmomega)*CSRC_C + lcsmomega*EQC ; \
+ RAC(tcel,dN ) = (1.0-lcsmomega)*CSRC_N + lcsmomega*lcsmeq[ dN ]; \
+ RAC(tcel,dS ) = (1.0-lcsmomega)*CSRC_S + lcsmomega*lcsmeq[ dS ]; \
+ RAC(tcel,dE ) = (1.0-lcsmomega)*CSRC_E + lcsmomega*lcsmeq[ dE ]; \
+ RAC(tcel,dW ) = (1.0-lcsmomega)*CSRC_W + lcsmomega*lcsmeq[ dW ]; \
+ RAC(tcel,dT ) = (1.0-lcsmomega)*CSRC_T + lcsmomega*lcsmeq[ dT ]; \
+ RAC(tcel,dB ) = (1.0-lcsmomega)*CSRC_B + lcsmomega*lcsmeq[ dB ]; \
+ \
+ RAC(tcel,dNE) = (1.0-lcsmomega)*CSRC_NE + lcsmomega*lcsmeq[ dNE]; \
+ RAC(tcel,dNW) = (1.0-lcsmomega)*CSRC_NW + lcsmomega*lcsmeq[ dNW]; \
+ RAC(tcel,dSE) = (1.0-lcsmomega)*CSRC_SE + lcsmomega*lcsmeq[ dSE]; \
+ RAC(tcel,dSW) = (1.0-lcsmomega)*CSRC_SW + lcsmomega*lcsmeq[ dSW]; \
+ \
+ RAC(tcel,dNT) = (1.0-lcsmomega)*CSRC_NT + lcsmomega*lcsmeq[ dNT]; \
+ RAC(tcel,dNB) = (1.0-lcsmomega)*CSRC_NB + lcsmomega*lcsmeq[ dNB]; \
+ RAC(tcel,dST) = (1.0-lcsmomega)*CSRC_ST + lcsmomega*lcsmeq[ dST]; \
+ RAC(tcel,dSB) = (1.0-lcsmomega)*CSRC_SB + lcsmomega*lcsmeq[ dSB]; \
+ \
+ RAC(tcel,dET) = (1.0-lcsmomega)*CSRC_ET + lcsmomega*lcsmeq[ dET]; \
+ RAC(tcel,dEB) = (1.0-lcsmomega)*CSRC_EB + lcsmomega*lcsmeq[ dEB]; \
+ RAC(tcel,dWT) = (1.0-lcsmomega)*CSRC_WT + lcsmomega*lcsmeq[ dWT]; \
+ RAC(tcel,dWB) = (1.0-lcsmomega)*CSRC_WB + lcsmomega*lcsmeq[ dWB]; \
+
+
+
+#define OPTIMIZED_STREAMCOLLIDE_NOLES \
+ \
+ rho = CSRC_C + CSRC_N + CSRC_S + CSRC_E + CSRC_W + CSRC_T \
+ + CSRC_B + CSRC_NE + CSRC_NW + CSRC_SE + CSRC_SW + CSRC_NT \
+ + CSRC_NB + CSRC_ST + CSRC_SB + CSRC_ET + CSRC_EB + CSRC_WT + CSRC_WB; \
+ ux = CSRC_E - CSRC_W + CSRC_NE - CSRC_NW + CSRC_SE - CSRC_SW \
+ + CSRC_ET + CSRC_EB - CSRC_WT - CSRC_WB; \
+ uy = CSRC_N - CSRC_S + CSRC_NE + CSRC_NW - CSRC_SE - CSRC_SW \
+ + CSRC_NT + CSRC_NB - CSRC_ST - CSRC_SB; \
+ uz = CSRC_T - CSRC_B + CSRC_NT - CSRC_NB + CSRC_ST - CSRC_SB \
+ + CSRC_ET - CSRC_EB + CSRC_WT - CSRC_WB; \
+ PRECOLLIDE_MODS(rho, ux,uy,uz, mLevel[lev].gravity); \
+ usqr = 1.5 * (ux*ux + uy*uy + uz*uz); \
+ RAC(tcel,dC ) = (1.0-OMEGA(lev))*CSRC_C + OMEGA(lev)*EQC ; \
+ RAC(tcel,dN ) = (1.0-OMEGA(lev))*CSRC_N + OMEGA(lev)*EQN ; \
+ RAC(tcel,dS ) = (1.0-OMEGA(lev))*CSRC_S + OMEGA(lev)*EQS ; \
+ RAC(tcel,dE ) = (1.0-OMEGA(lev))*CSRC_E + OMEGA(lev)*EQE ; \
+ RAC(tcel,dW ) = (1.0-OMEGA(lev))*CSRC_W + OMEGA(lev)*EQW ; \
+ RAC(tcel,dT ) = (1.0-OMEGA(lev))*CSRC_T + OMEGA(lev)*EQT ; \
+ RAC(tcel,dB ) = (1.0-OMEGA(lev))*CSRC_B + OMEGA(lev)*EQB ; \
+ \
+ RAC(tcel,dNE) = (1.0-OMEGA(lev))*CSRC_NE + OMEGA(lev)*EQNE; \
+ RAC(tcel,dNW) = (1.0-OMEGA(lev))*CSRC_NW + OMEGA(lev)*EQNW; \
+ RAC(tcel,dSE) = (1.0-OMEGA(lev))*CSRC_SE + OMEGA(lev)*EQSE; \
+ RAC(tcel,dSW) = (1.0-OMEGA(lev))*CSRC_SW + OMEGA(lev)*EQSW; \
+ \
+ RAC(tcel,dNT) = (1.0-OMEGA(lev))*CSRC_NT + OMEGA(lev)*EQNT; \
+ RAC(tcel,dNB) = (1.0-OMEGA(lev))*CSRC_NB + OMEGA(lev)*EQNB; \
+ RAC(tcel,dST) = (1.0-OMEGA(lev))*CSRC_ST + OMEGA(lev)*EQST; \
+ RAC(tcel,dSB) = (1.0-OMEGA(lev))*CSRC_SB + OMEGA(lev)*EQSB; \
+ \
+ RAC(tcel,dET) = (1.0-OMEGA(lev))*CSRC_ET + OMEGA(lev)*EQET; \
+ RAC(tcel,dEB) = (1.0-OMEGA(lev))*CSRC_EB + OMEGA(lev)*EQEB; \
+ RAC(tcel,dWT) = (1.0-OMEGA(lev))*CSRC_WT + OMEGA(lev)*EQWT; \
+ RAC(tcel,dWB) = (1.0-OMEGA(lev))*CSRC_WB + OMEGA(lev)*EQWB; \
+
+
+
+
+
+// LES switching for OPT3D
#if USE_LES==1
#define DEFAULT_COLLIDEG(grav) DEFAULT_COLLIDE_LES(grav)
#define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_LES
@@ -723,7 +797,7 @@ errMsg("T","QSDM at %d,%d,%d lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lc
#define OPTIMIZED_STREAMCOLLIDE OPTIMIZED_STREAMCOLLIDE_NOLES
#endif
-#endif
+#endif // 3D, opt OPT3D==true
#define USQRMAXCHECK(Cusqr,Cux,Cuy,Cuz, CmMaxVlen,CmMxvx,CmMxvy,CmMxvz) \
if(Cusqr>CmMaxVlen) { \
@@ -779,10 +853,6 @@ errMsg("T","QSDM at %d,%d,%d lcsmqo=%25.15f, lcsmomega=%f \n", i,j,k, lcsmqo,lc
}
// end ADD_INT_DFSCHECK
- //if( !(RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) & CFUnused) ){
- //errMsg("INTFLAGUNU", PRINT_VEC(i,j,k)<<" child at "<<PRINT_VEC((ix),(iy),(iz)) );
- //if(iy==15) errMsg("IFFC", PRINT_VEC(i,j,k)<<" child interpolated at "<<PRINT_VEC((ix),(iy),(iz)) );
- //if(((ix)>10)&&(iy>5)&&(iz>5)) { debugMarkCell(lev+1, (ix),(iy),(iz) ); }
#define INTUNUTCHECK(ix,iy,iz) \
if( (RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) != (CFFluid|CFGrFromCoarse)) ){\
errMsg("INTFLAGUNU_CHECK", PRINT_VEC(i,j,k)<<" child not unused at l"<<(lev+1)<<" "<<PRINT_VEC((ix),(iy),(iz))<<" flag: "<< RFLAG(lev+1, (ix),(iy),(iz), mLevel[lev+1].setCurr) ); \
@@ -1078,7 +1148,7 @@ inline LbmFloat LbmFsgrSolver::getLesOmega(LbmFloat omega, LbmFloat csmago, LbmF
// "normal" collision
inline void LbmFsgrSolver::collideArrays(
- int i, int j, int k, // position - more for debugging
+ int lev, int i, int j, int k, // position - more for debugging
LbmFloat df[],
LbmFloat &outrho, // out only!
// velocity modifiers (returns actual velocity!)
@@ -1104,6 +1174,7 @@ inline void LbmFsgrSolver::collideArrays(
uz += (this->dfDvecZ[l]*df[l]);
}
+
PRECOLLIDE_MODS(rho,ux,uy,uz, gravity);
for(l=0; l<this->cDfNum; l++) {
feq[l] = getCollideEq(l,rho,ux,uy,uz);
diff --git a/intern/elbeem/intern/solver_util.cpp b/intern/elbeem/intern/solver_util.cpp
index a38307a3b02..6b7f6f2ed9c 100644
--- a/intern/elbeem/intern/solver_util.cpp
+++ b/intern/elbeem/intern/solver_util.cpp
@@ -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
+#else // SURFACE_ENH!=1
+ if(cflag&CFBnd) {
+ // treated in second loop
+ continue;
+ } else if(cflag&CFUnused) {
+ val = 1.;
+ } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid)) {
+ // optimized fluid
+ val = 1.;
+ } else if( (cflag&(CFEmpty|CFInter|CFFluid)) ) {
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 {
- }
- continue;
- //} else if( (cflag&CFInter) ) {
+ if(l!=mainGravDir) continue; // only check bnd along main grav. dir
+ //if((nbflag&CFBnd)&&(nbflag&CFBndNoslip)){ noslipbnd=1; }
+ if((nbflag&CFBnd)){ noslipbnd=1; }
+ }
- } else if( (cflag&CFFluid) && (cflag&CFNoBndFluid) ) {
- // optimized fluid
- val = 1.;
+ 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...
+ }
- } else if( (cflag&(CFInter|CFFluid)) ) {
- int noslipbnd = 0;
- FORDF1 {
- const CellFlagType nbflag = RFLAG_NB(lev, i,j,k, workSet,l);
- if((l==6)&&(nbflag&CFBnd)){ noslipbnd=1; l=100; }
- }
- // no empty nb interface cells are treated as full
- if(cflag&(CFNoNbEmpty|CFFluid)) {
+ 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;
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-23 22:35:25 +0300