diff options
author | Nils Thuerey <nils@thuerey.de> | 2006-08-22 14:45:19 +0400 |
---|---|---|
committer | Nils Thuerey <nils@thuerey.de> | 2006-08-22 14:45:19 +0400 |
commit | b6257305c929464ef50902ab7c98a7c37dd17fd7 (patch) | |
tree | 8bb33edeb0726760795f0ed9dcf27b4af6f36b92 /intern/elbeem | |
parent | d6d2d6f063994e1481ce498c5e245d8470ef4b6d (diff) |
elbeem update:
- slightly improved and optimized
handling of moving obstacles
- cleanup of debug messages and minor fixes
Diffstat (limited to 'intern/elbeem')
-rw-r--r-- | intern/elbeem/intern/elbeem.h | 1 | ||||
-rw-r--r-- | intern/elbeem/intern/isosurface.cpp | 10 | ||||
-rw-r--r-- | intern/elbeem/intern/isosurface.h | 3 | ||||
-rw-r--r-- | intern/elbeem/intern/loop_tools.h | 105 | ||||
-rw-r--r-- | intern/elbeem/intern/ntl_geometryobject.cpp | 143 | ||||
-rw-r--r-- | intern/elbeem/intern/ntl_geometryobject.h | 5 | ||||
-rw-r--r-- | intern/elbeem/intern/simulation_object.cpp | 2 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_adap.cpp | 44 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_class.h | 89 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_init.cpp | 378 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_interface.cpp | 4 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_interface.h | 1 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_main.cpp | 507 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_relax.h | 1077 | ||||
-rw-r--r-- | intern/elbeem/intern/solver_util.cpp | 422 |
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; |