From af4593612b409481fc2773e0b868016b0585fc85 Mon Sep 17 00:00:00 2001
From: "LEGOLAS\\lenz" <lenz@irmb.tu-bs.de>
Date: Fri, 14 Sep 2018 17:39:53 +0200
Subject: [PATCH] implements thin wall
---
src/GridGenerator/grid/Grid.h | 6 +-
.../grid/GridBuilder/GridBuilder.h | 4 +-
.../grid/GridBuilder/LevelGridBuilder.cpp | 6 +-
.../grid/GridBuilder/LevelGridBuilder.h | 7 +-
.../grid/GridBuilder/MultipleGridBuilder.cpp | 9 +-
.../grid/GridBuilder/MultipleGridBuilder.h | 2 +-
src/GridGenerator/grid/GridImp.cu | 159 +++-
src/GridGenerator/grid/GridImp.h | 23 +-
src/GridGenerator/grid/GridMocks.h | 10 +-
.../GridCpuStrategy/GridCpuStrategy.cpp | 8 +-
.../GridGpuStrategy/GridGpuStrategy.cpp | 47 +-
.../GridReaderGenerator/GridGenerator.cpp | 3 +
.../GPU/CudaMemoryManager.cpp | 20 +
src/VirtualFluids_GPU/GPU/CudaMemoryManager.h | 4 +
src/VirtualFluids_GPU/GPU/GPU_Interface.h | 20 +
src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh | 34 +
src/VirtualFluids_GPU/GPU/LBMKernel.cu | 71 ++
src/VirtualFluids_GPU/GPU/VelocityBCs27.cu | 833 ++++++++++++++++++
src/VirtualFluids_GPU/LBM/Simulation.cpp | 32 +-
targets/apps/HULC/main.cpp | 10 +-
20 files changed, 1223 insertions(+), 85 deletions(-)
diff --git a/src/GridGenerator/grid/Grid.h b/src/GridGenerator/grid/Grid.h
index 86ea77c97..19615b95b 100644
--- a/src/GridGenerator/grid/Grid.h
+++ b/src/GridGenerator/grid/Grid.h
@@ -53,6 +53,7 @@ public:
HOST virtual int *getNeighborsX() const = 0;
HOST virtual int *getNeighborsY() const = 0;
HOST virtual int *getNeighborsZ() const = 0;
+ HOST virtual int *getNeighborsNegative() const = 0;
HOST virtual uint* getCF_coarse() const = 0;
HOST virtual uint* getCF_fine() const = 0;
@@ -67,7 +68,7 @@ public:
HOST virtual int getStartDirection() const = 0;
HOST virtual int getEndDirection() const = 0;
- HOST virtual void getNodeValues(real *xCoords, real *yCoords, real *zCoords, unsigned int *neighborX, unsigned int *neighborY, unsigned int *neighborZ, unsigned int *geo) const = 0;
+ HOST virtual void getNodeValues(real *xCoords, real *yCoords, real *zCoords, uint *neighborX, uint *neighborY, uint *neighborZ, uint *neighborNegative, uint *geo) const = 0;
HOST virtual SPtr<GridStrategy> getGridStrategy() const = 0;
HOSTDEVICE virtual void transIndexToCoords(uint index, real &x, real &y, real &z) const = 0;
@@ -78,6 +79,9 @@ public:
HOST virtual void setOddStart( bool xOddStart, bool yOddStart, bool zOddStart ) = 0;
HOST virtual void findGridInterface(SPtr<Grid> grid, LbmOrGks lbmOrGks) = 0;
+
+ HOST virtual void enableFindSolidBoundaryNodes() = 0;
+ HOST virtual void enableComputeQs() = 0;
HOST virtual void mesh(TriangularMesh& geometry) = 0;
HOST virtual void mesh(Object* object) = 0;
diff --git a/src/GridGenerator/grid/GridBuilder/GridBuilder.h b/src/GridGenerator/grid/GridBuilder/GridBuilder.h
index fef10e893..02453626f 100644
--- a/src/GridGenerator/grid/GridBuilder/GridBuilder.h
+++ b/src/GridGenerator/grid/GridBuilder/GridBuilder.h
@@ -47,7 +47,9 @@ public:
virtual SPtr<Grid> getGrid(uint level) = 0;
virtual unsigned int getNumberOfNodes(unsigned int level) const = 0;
- virtual void getNodeValues(real *xCoords, real *yCoords, real *zCoords, unsigned int *nx, unsigned int *ny, unsigned int *nz, unsigned int *geo, const int level) const = 0;
+ virtual void getNodeValues(real *xCoords, real *yCoords, real *zCoords,
+ uint *neighborX, uint *neighborY, uint *neighborZ, uint *neighborNegative,
+ uint *geo, const int level) const = 0;
virtual void getDimensions(int &nx, int &ny, int &nz, const int level) const = 0;
virtual uint getNumberOfNodesCF(int level) = 0;
virtual uint getNumberOfNodesFC(int level) = 0;
diff --git a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
index 02d3374f3..42c598d77 100644
--- a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
+++ b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.cpp
@@ -248,9 +248,11 @@ void LevelGridBuilder::getDimensions(int &nx, int &ny, int &nz, const int level)
nz = grids[level]->getNumberOfNodesZ();
}
-void LevelGridBuilder::getNodeValues(real *xCoords, real *yCoords, real *zCoords, unsigned int *neighborX, unsigned int *neighborY, unsigned int *neighborZ, unsigned int *geo, const int level) const
+void LevelGridBuilder::getNodeValues(real *xCoords, real *yCoords, real *zCoords,
+ uint *neighborX, uint *neighborY, uint *neighborZ, uint *neighborNegative,
+ uint *geo, const int level) const
{
- grids[level]->getNodeValues(xCoords, yCoords, zCoords, neighborX, neighborY, neighborZ, geo);
+ grids[level]->getNodeValues(xCoords, yCoords, zCoords, neighborX, neighborY, neighborZ, neighborNegative, geo);
}
diff --git a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
index e3dfb810c..eacc8b3a8 100644
--- a/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
+++ b/src/GridGenerator/grid/GridBuilder/LevelGridBuilder.h
@@ -52,11 +52,12 @@ public:
VF_PUBLIC virtual std::shared_ptr<Grid> getGrid(int level, int box);
- VF_PUBLIC virtual unsigned int getNumberOfNodes(unsigned int level) const;;
+ VF_PUBLIC virtual unsigned int getNumberOfNodes(unsigned int level) const;
- VF_PUBLIC virtual void getNodeValues(real *xCoords, real *yCoords, real *zCoords, unsigned int *nx,
- unsigned int *ny, unsigned int *nz, unsigned int *geo, const int level) const;
+ VF_PUBLIC virtual void getNodeValues(real *xCoords, real *yCoords, real *zCoords,
+ uint *neighborX, uint *neighborY, uint *neighborZ, uint *neighborNegative,
+ uint *geo, const int level) const override;
VF_PUBLIC virtual void getDimensions(int &nx, int &ny, int &nz, const int level) const;
diff --git a/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp b/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
index c64e76a77..501305e93 100644
--- a/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
+++ b/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.cpp
@@ -405,7 +405,7 @@ std::vector<SPtr<Grid> > MultipleGridBuilder::getGrids() const
return this->grids;
}
-void MultipleGridBuilder::buildGrids( LbmOrGks lbmOrGks )
+void MultipleGridBuilder::buildGrids( LbmOrGks lbmOrGks, bool enableThinWalls )
{
//////////////////////////////////////////////////////////////////////////
@@ -444,6 +444,13 @@ void MultipleGridBuilder::buildGrids( LbmOrGks lbmOrGks )
*logging::out << logging::Logger::INFO_INTERMEDIATE << "Start with Q Computation\n";
grids[grids.size() - 1]->mesh(solidObject);
+
+ if( enableThinWalls ){
+ grids[grids.size() - 1]->enableFindSolidBoundaryNodes();
+ grids[grids.size() - 1]->findQs(solidObject);
+ grids[grids.size() - 1]->enableComputeQs();
+ }
+
grids[grids.size() - 1]->findQs(solidObject);
//for (size_t i = 0; i < grids.size(); i++){
diff --git a/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h b/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h
index 148bd1df0..5d50e3a29 100644
--- a/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h
+++ b/src/GridGenerator/grid/GridBuilder/MultipleGridBuilder.h
@@ -42,7 +42,7 @@ public:
VF_PUBLIC real getEndZ(uint level) const;
VF_PUBLIC std::vector<SPtr<Grid> > getGrids() const;
- VF_PUBLIC void buildGrids(LbmOrGks lbmOrGks);
+ VF_PUBLIC void buildGrids(LbmOrGks lbmOrGks, bool enableThinWalls = false);
VF_PUBLIC void setNumberOfLayers( uint numberOfLayersFine, uint numberOfLayersBetweenLevels );
diff --git a/src/GridGenerator/grid/GridImp.cu b/src/GridGenerator/grid/GridImp.cu
index 3b4d619f1..c2f3c150f 100644
--- a/src/GridGenerator/grid/GridImp.cu
+++ b/src/GridGenerator/grid/GridImp.cu
@@ -49,7 +49,8 @@ HOST GridImp::GridImp(Object* object, real startX, real startY, real startZ, rea
level(level),
gridInterface(nullptr),
innerRegionFromFinerGrid(false),
- numberOfLayers(0)
+ numberOfLayers(0),
+ qComputationStage(qComputationStageType::ComputeQs)
{
initalNumberOfNodesAndSize();
}
@@ -802,9 +803,10 @@ HOSTDEVICE void GridImp::setNeighborIndices(uint index)
real x, y, z;
this->transIndexToCoords(index, x, y, z);
- neighborIndexX[index] = -1;
- neighborIndexY[index] = -1;
- neighborIndexZ[index] = -1;
+ this->neighborIndexX[index] = -1;
+ this->neighborIndexY[index] = -1;
+ this->neighborIndexZ[index] = -1;
+ this->neighborIndexNegative[index] = -1;
if (this->field.isStopper(index) || this->field.is(index, STOPPER_OUT_OF_GRID_BOUNDARY))
{
@@ -822,9 +824,13 @@ HOSTDEVICE void GridImp::setNeighborIndices(uint index)
const int neighborY = getSparseIndex(x, neighborYCoord, z);
const int neighborZ = getSparseIndex(x, y, neighborZCoord);
- this->neighborIndexX[index] = neighborX;
- this->neighborIndexY[index] = neighborY;
- this->neighborIndexZ[index] = neighborZ;
+ this->getNegativeNeighborCoords(neighborXCoord, neighborYCoord, neighborZCoord, x, y, z);
+ const int neighborNegative = getSparseIndex(neighborXCoord, neighborYCoord, neighborZCoord);
+
+ this->neighborIndexX[index] = neighborX;
+ this->neighborIndexY[index] = neighborY;
+ this->neighborIndexZ[index] = neighborZ;
+ this->neighborIndexNegative[index] = neighborNegative;
}
HOSTDEVICE void GridImp::setStopperNeighborCoords(uint index)
@@ -840,6 +846,14 @@ HOSTDEVICE void GridImp::setStopperNeighborCoords(uint index)
if (vf::Math::lessEqual(z + delta, endZ) && !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x, y, z + delta)))
neighborIndexZ[index] = getSparseIndex(x, y, z + delta);
+
+ if (vf::Math::greaterEqual(x - delta, endX) &&
+ vf::Math::greaterEqual(y - delta, endY) &&
+ vf::Math::greaterEqual(z - delta, endZ) &&
+ !this->field.isInvalidOutOfGrid(this->transCoordToIndex(x - delta, y - delta, z - delta)))
+ {
+ neighborIndexNegative[index] = getSparseIndex(x - delta, y - delta, z - delta);
+ }
}
HOSTDEVICE void GridImp::getNeighborCoords(real &neighborX, real &neighborY, real &neighborZ, real x, real y, real z) const
@@ -867,6 +881,31 @@ HOSTDEVICE real GridImp::getNeighborCoord(bool periodicity, real startCoord, rea
return coords[direction] + delta;
}
+HOSTDEVICE void GridImp::getNegativeNeighborCoords(real &neighborX, real &neighborY, real &neighborZ, real x, real y, real z) const
+{
+ real coords[3] = { x, y, z };
+ neighborX = getNegativeNeighborCoord(periodicityX, startX, coords, 0);
+ neighborY = getNegativeNeighborCoord(periodicityY, startY, coords, 1);
+ neighborZ = getNegativeNeighborCoord(periodicityZ, startZ, coords, 2);
+}
+
+HOSTDEVICE real GridImp::getNegativeNeighborCoord(bool periodicity, real startCoord, real coords[3], int direction) const
+{
+ if (periodicity)
+ {
+ real neighborCoords[3] = {coords[0], coords[1] , coords[2] };
+ neighborCoords[direction] = neighborCoords[direction] - delta;
+ const int neighborIndex = this->transCoordToIndex(neighborCoords[0], neighborCoords[1], neighborCoords[2]);
+
+ if(!field.isStopperOutOfGrid(neighborIndex) && !field.is(neighborIndex, STOPPER_OUT_OF_GRID_BOUNDARY) )
+ return coords[direction] - delta;
+
+ return getLastFluidNode(coords, direction, startCoord);
+ }
+
+ return coords[direction] - delta;
+}
+
HOSTDEVICE real GridImp::getLastFluidNode(real coords[3], int direction, real startCoord) const
{
@@ -998,7 +1037,10 @@ HOST void GridImp::findQs(TriangularMesh &triangularMesh)
{
const clock_t begin = clock();
- gridStrategy->allocateQs(shared_from_this());
+ if( this->qComputationStage == qComputationStageType::ComputeQs ){
+ gridStrategy->allocateQs(shared_from_this());
+ }
+
gridStrategy->findQs(shared_from_this(), triangularMesh);
const clock_t end = clock();
@@ -1026,11 +1068,23 @@ HOSTDEVICE void GridImp::findQs(Triangle &triangle)
const Vertex point(x, y, z);
- if(this->field.is(index, BC_SOLID)) //TODO: not working with thin walls
+ if( this->qComputationStage == qComputationStageType::ComputeQs ){
+ if(this->field.is(index, BC_SOLID)) //TODO: not working with thin walls
+ {
+ calculateQs(index, point, triangle);
+ }
+ }
+ else if( this->qComputationStage == qComputationStageType::FindSolidBoundaryNodes )
{
- //calculateQs(point, triangle);
- calculateQs(index, point, triangle);
- }
+ if( this->field.is(index, BC_SOLID) || this->field.is(index, STOPPER_SOLID ) ) continue;
+
+ if( checkIfAtLeastOneValidQ(index, point, triangle) )
+ {
+ // similar as in void GridImp::findBoundarySolidNode(uint index)
+ this->field.setFieldEntry( index, BC_SOLID );
+ this->qIndices[index] = this->numberOfSolidBoundaryNodes++;
+ }
+ }
}
}
}
@@ -1045,7 +1099,7 @@ HOSTDEVICE void GridImp::setDebugPoint(uint index, int pointValue)
field.setFieldEntry(index, pointValue);
}
-HOSTDEVICE void GridImp::calculateQs(const Vertex &point, const Triangle &triangle) const
+HOSTDEVICE void GridImp::calculateQs(const Vertex &point, const Triangle &triangle) const // NOT USED !!!!
{
Vertex pointOnTriangle, direction;
//VertexInteger solid_node;
@@ -1095,12 +1149,14 @@ HOSTDEVICE void GridImp::calculateQs(const uint index, const Vertex &point, cons
point.y + direction.y * this->delta,
point.z + direction.z * this->delta);
if (neighborIndex == INVALID_INDEX) continue;
- if (!this->field.is(neighborIndex, STOPPER_SOLID)) continue;
- if (this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] < -0.5)
- {
- this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] = 1.0;
- }
+ if ( this->field.is(neighborIndex, STOPPER_SOLID) )
+ {
+ if (this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] < -0.5)
+ {
+ this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] = 1.0;
+ }
+ }
error = triangle.getTriangleIntersection(point, direction, pointOnTriangle, subdistance);
@@ -1108,25 +1164,46 @@ HOSTDEVICE void GridImp::calculateQs(const uint index, const Vertex &point, cons
if (error == 0 && vf::Math::lessEqual(subdistance, 1.0) && vf::Math::greaterEqual(subdistance, 0.0))
{
- if (subdistance < this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]])
+ if ( -0.5 > this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] ||
+ subdistance < this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] )
{
this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] = subdistance;
//printf("%d %f \n", this->qIndices[index], subdistance);
}
}
+ }
+}
+
+HOSTDEVICE bool GridImp::checkIfAtLeastOneValidQ(const uint index, const Vertex & point, const Triangle & triangle) const
+{
+ Vertex pointOnTriangle, direction;
+ //VertexInteger solid_node;
+ real subdistance;
+ int error;
+ for (int i = distribution.dir_start; i <= distribution.dir_end; i++)
+ {
+#if defined(__CUDA_ARCH__)
+ direction = Vertex(DIRECTIONS[i][0], DIRECTIONS[i][1], DIRECTIONS[i][2]);
+#else
+ direction = Vertex(real(distribution.dirs[i * DIMENSION + 0]), real(distribution.dirs[i * DIMENSION + 1]),
+ real(distribution.dirs[i * DIMENSION + 2]));
+#endif
- //if (error == 0 && vf::Math::lessEqual(subdistance, 1.0) && vf::Math::greaterEqual(subdistance, 0.0))
- //{
- // this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] = subdistance;
- //}
- //else
- //{
- // if (this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] < -0.5)
- // {
- // this->qValues[i*this->numberOfSolidBoundaryNodes + this->qIndices[index]] = 1.0;
- // }
- //}
+ uint neighborIndex = this->transCoordToIndex(point.x + direction.x * this->delta,
+ point.y + direction.y * this->delta,
+ point.z + direction.z * this->delta);
+ if (neighborIndex == INVALID_INDEX) continue;
+
+ error = triangle.getTriangleIntersection(point, direction, pointOnTriangle, subdistance);
+
+ subdistance /= this->delta;
+
+ if (error == 0 && vf::Math::lessEqual(subdistance, 1.0) && vf::Math::greaterEqual(subdistance, 0.0))
+ {
+ return true;
+ }
}
+ return false;
}
@@ -1338,6 +1415,11 @@ int* GridImp::getNeighborsZ() const
return this->neighborIndexZ;
}
+int* GridImp::getNeighborsNegative() const
+{
+ return this->neighborIndexNegative;
+}
+
uint GridImp::getNumberOfNodesCF() const
{
@@ -1400,7 +1482,7 @@ void GridImp::getGridInterface(uint* gridInterfaceList, const uint* oldGridInter
#define GEOFLUID 19
#define GEOSOLID 16
-HOST void GridImp::getNodeValues(real *xCoords, real *yCoords, real *zCoords, uint *neighborX, uint *neighborY, uint *neighborZ, uint *geo) const
+HOST void GridImp::getNodeValues(real *xCoords, real *yCoords, real *zCoords, uint *neighborX, uint *neighborY, uint *neighborZ, uint *neighborNegative, uint *geo) const
{
xCoords[0] = 0;
yCoords[0] = 0;
@@ -1419,9 +1501,10 @@ HOST void GridImp::getNodeValues(real *xCoords, real *yCoords, real *zCoords, ui
real x, y, z;
this->transIndexToCoords(i, x, y, z);
- const uint neighborXIndex = uint(this->neighborIndexX[i] + 1);
- const uint neighborYIndex = uint(this->neighborIndexY[i] + 1);
- const uint neighborZIndex = uint(this->neighborIndexZ[i] + 1);
+ const uint neighborXIndex = uint(this->neighborIndexX[i] + 1);
+ const uint neighborYIndex = uint(this->neighborIndexY[i] + 1);
+ const uint neighborZIndex = uint(this->neighborIndexZ[i] + 1);
+ const uint neighborNegativeIndex = uint(this->neighborIndexNegative[i] + 1);
const char type2 = this->field.getFieldEntry(i);
@@ -1431,9 +1514,11 @@ HOST void GridImp::getNodeValues(real *xCoords, real *yCoords, real *zCoords, ui
yCoords[nodeNumber + 1] = y;
zCoords[nodeNumber + 1] = z;
- neighborX[nodeNumber + 1] = neighborXIndex;
- neighborY[nodeNumber + 1] = neighborYIndex;
- neighborZ[nodeNumber + 1] = neighborZIndex;
+ neighborX [nodeNumber + 1] = neighborXIndex;
+ neighborY [nodeNumber + 1] = neighborYIndex;
+ neighborZ [nodeNumber + 1] = neighborZIndex;
+ neighborNegative[nodeNumber + 1] = neighborNegativeIndex;
+
geo[nodeNumber + 1] = type;
nodeNumber++;
}
diff --git a/src/GridGenerator/grid/GridImp.h b/src/GridGenerator/grid/GridImp.h
index 4b0f7bd7f..690449c17 100644
--- a/src/GridGenerator/grid/GridImp.h
+++ b/src/GridGenerator/grid/GridImp.h
@@ -68,7 +68,7 @@ private:
Object* object;
GridInterface* gridInterface;
- int *neighborIndexX, *neighborIndexY, *neighborIndexZ;
+ int *neighborIndexX, *neighborIndexY, *neighborIndexZ, *neighborIndexNegative;
int *sparseIndices;
uint *qIndices;
@@ -177,7 +177,7 @@ public:
HOSTDEVICE uint getNumberOfNodesX() const override;
HOSTDEVICE uint getNumberOfNodesY() const override;
HOSTDEVICE uint getNumberOfNodesZ() const override;
- HOST void getNodeValues(real *xCoords, real *yCoords, real *zCoords, uint *neighborX, uint *neighborY, uint *neighborZ, uint *geo) const override;
+ HOST void getNodeValues(real *xCoords, real *yCoords, real *zCoords, uint *neighborX, uint *neighborY, uint *neighborZ, uint *neighborNegative, uint *geo) const override;
HOSTDEVICE uint getNumberOfNodesCF() const override;
HOSTDEVICE uint getNumberOfNodesFC() const override;
@@ -187,6 +187,7 @@ public:
int* getNeighborsX() const override;
int* getNeighborsY() const override;
int* getNeighborsZ() const override;
+ int* getNeighborsNegative() const override;
HOST uint* getCF_coarse() const override;
HOST uint* getCF_fine() const override;
@@ -214,6 +215,8 @@ private:
HOSTDEVICE void setStopperNeighborCoords(uint index);
HOSTDEVICE void getNeighborCoords(real &neighborX, real &neighborY, real &neighborZ, real x, real y, real z) const;
HOSTDEVICE real getNeighborCoord(bool periodicity, real endCoord, real coords[3], int direction) const;
+ HOSTDEVICE void getNegativeNeighborCoords(real &neighborX, real &neighborY, real &neighborZ, real x, real y, real z) const;
+ HOSTDEVICE real getNegativeNeighborCoord(bool periodicity, real endCoord, real coords[3], int direction) const;
HOSTDEVICE int getSparseIndex(const real &expectedX, const real &expectedY, const real &expectedZ) const;
@@ -233,10 +236,22 @@ public:
HOST void findQs(Object* object) override;
HOST void findQs(TriangularMesh &triangularMesh);
HOSTDEVICE void findQs(Triangle &triangle);
+private:
+
+ enum class qComputationStageType{
+ FindSolidBoundaryNodes,
+ ComputeQs
+ } qComputationStage;
+
+public:
+ HOST void enableFindSolidBoundaryNodes(){ qComputationStage = qComputationStageType::FindSolidBoundaryNodes; }
+ HOST void enableComputeQs(){ qComputationStage = qComputationStageType::ComputeQs; }
+
private:
HOSTDEVICE void setDebugPoint(uint index, int pointValue);
- HOSTDEVICE void calculateQs(const Vertex &point, const Triangle &actualTriangle) const;
- HOSTDEVICE void calculateQs(const uint index, const Vertex &point, const Triangle &actualTriangle) const;
+ HOSTDEVICE void calculateQs(const Vertex &point, const Triangle &triangle) const;
+ HOSTDEVICE void calculateQs(const uint index, const Vertex &point, const Triangle &triangle) const;
+ HOSTDEVICE bool checkIfAtLeastOneValidQ(const uint index, const Vertex &point, const Triangle &triangle) const;
diff --git a/src/GridGenerator/grid/GridMocks.h b/src/GridGenerator/grid/GridMocks.h
index 38e40c4e9..35d3ff3b0 100644
--- a/src/GridGenerator/grid/GridMocks.h
+++ b/src/GridGenerator/grid/GridMocks.h
@@ -52,19 +52,25 @@ public:
virtual int* getDirection() const override { return nullptr; }
virtual int getStartDirection() const override { return 0; }
virtual int getEndDirection() const override { return 0; }
- virtual void getNodeValues(real* xCoords, real* yCoords, real* zCoords, unsigned* neighborX, unsigned* neighborY,
- unsigned* neighborZ, unsigned* geo) const override {}
+ virtual void getNodeValues(real* xCoords, real* yCoords, real* zCoords,
+ uint* neighborX, uint* neighborY, uint* neighborZ, uint* neighborNegative,
+ uint* geo) const override {}
virtual SPtr<GridStrategy> getGridStrategy() const override { return nullptr; }
virtual void transIndexToCoords(uint index, real& x, real& y, real& z) const override {}
virtual void setPeriodicity(bool periodicityX, bool periodicityY, bool periodicityZ) override {}
virtual void freeMemory() override {}
virtual void findGridInterface(SPtr<Grid> grid, LbmOrGks lbmOrGks) override {}
+
+ virtual void enableFindSolidBoundaryNodes() override {}
+ virtual void enableComputeQs() override {}
+
virtual void mesh(TriangularMesh& geometry) override {}
virtual uint transCoordToIndex(const real& x, const real& y, const real& z) const override { return 0; }
virtual int* getNeighborsX() const override { return nullptr; }
virtual int* getNeighborsY() const override { return nullptr; }
virtual int* getNeighborsZ() const override { return nullptr; }
+ virtual int* getNeighborsNegative() const override { return nullptr; }
virtual void inital(const SPtr<Grid> fineGrid, uint numberOfLayers) override {};
virtual void setOddStart( bool xOddStart, bool yOddStart, bool zOddStart ) {};
virtual bool nodeInCellIs(Cell& cell, char type) const override { return false; }
diff --git a/src/GridGenerator/grid/GridStrategy/GridCpuStrategy/GridCpuStrategy.cpp b/src/GridGenerator/grid/GridStrategy/GridCpuStrategy/GridCpuStrategy.cpp
index a223d8f23..5a2a14bc0 100644
--- a/src/GridGenerator/grid/GridStrategy/GridCpuStrategy/GridCpuStrategy.cpp
+++ b/src/GridGenerator/grid/GridStrategy/GridCpuStrategy/GridCpuStrategy.cpp
@@ -19,9 +19,10 @@
void GridCpuStrategy::allocateGridMemory(SPtr<GridImp> grid)
{
- grid->neighborIndexX = new int[grid->size];
- grid->neighborIndexY = new int[grid->size];
- grid->neighborIndexZ = new int[grid->size];
+ grid->neighborIndexX = new int[grid->size];
+ grid->neighborIndexY = new int[grid->size];
+ grid->neighborIndexZ = new int[grid->size];
+ grid->neighborIndexNegative = new int[grid->size];
grid->sparseIndices = new int[grid->size];
@@ -245,6 +246,7 @@ void GridCpuStrategy::freeMemory(SPtr<GridImp> grid)
delete[] grid->neighborIndexX;
delete[] grid->neighborIndexY;
delete[] grid->neighborIndexZ;
+ delete[] grid->neighborIndexNegative;
delete[] grid->sparseIndices;
delete[] grid->qIndices;
delete[] grid->qValues;
diff --git a/src/GridGenerator/grid/GridStrategy/GridGpuStrategy/GridGpuStrategy.cpp b/src/GridGenerator/grid/GridStrategy/GridGpuStrategy/GridGpuStrategy.cpp
index dc8395b5a..b6aa03934 100644
--- a/src/GridGenerator/grid/GridStrategy/GridGpuStrategy/GridGpuStrategy.cpp
+++ b/src/GridGenerator/grid/GridStrategy/GridGpuStrategy/GridGpuStrategy.cpp
@@ -190,6 +190,7 @@ void GridGpuStrategy::freeMemory(SPtr<GridImp> grid)
delete[] grid->neighborIndexX;
delete[] grid->neighborIndexY;
delete[] grid->neighborIndexZ;
+ delete[] grid->neighborIndexNegative;
delete[] grid->sparseIndices;
delete[] grid->distribution.f;
@@ -264,13 +265,18 @@ void GridGpuStrategy::allocDistribution(SPtr<GridImp> grid)
void GridGpuStrategy::allocNeighborsIndices(SPtr<GridImp> grid)
{
int size_in_bytes_neighbors = grid->size * sizeof(int);
- int *neighborIndexX, *neighborIndexY, *neighborIndexZ;
- CudaSafeCall(cudaMalloc(&neighborIndexX, size_in_bytes_neighbors));
- CudaSafeCall(cudaMalloc(&neighborIndexY, size_in_bytes_neighbors));
- CudaSafeCall(cudaMalloc(&neighborIndexZ, size_in_bytes_neighbors));
- grid->neighborIndexX = neighborIndexX;
- grid->neighborIndexY = neighborIndexY;
- grid->neighborIndexZ = neighborIndexZ;
+ int *neighborIndexX, *neighborIndexY, *neighborIndexZ, *neighborIndexNegative;;
+
+ CudaSafeCall(cudaMalloc(&neighborIndexX, size_in_bytes_neighbors));
+ CudaSafeCall(cudaMalloc(&neighborIndexY, size_in_bytes_neighbors));
+ CudaSafeCall(cudaMalloc(&neighborIndexZ, size_in_bytes_neighbors));
+ CudaSafeCall(cudaMalloc(&neighborIndexNegative, size_in_bytes_neighbors));
+
+ grid->neighborIndexX = neighborIndexX;
+ grid->neighborIndexY = neighborIndexY;
+ grid->neighborIndexZ = neighborIndexZ;
+ grid->neighborIndexNegative = neighborIndexNegative;
+
CudaCheckError();
}
@@ -356,21 +362,26 @@ void GridGpuStrategy::copyAndFreeDistributiondFromGPU(SPtr<GridImp> grid)
void GridGpuStrategy::copyAndFreeNeighborsToCPU(SPtr<GridImp> grid)
{
int size_in_bytes_neighbors = grid->size * sizeof(int);
- int *neighborIndexX_h, *neighborIndexY_h, *neighborIndexZ_h;
- neighborIndexX_h = new int[grid->size];
- neighborIndexY_h = new int[grid->size];
- neighborIndexZ_h = new int[grid->size];
-
- CudaSafeCall(cudaMemcpy(neighborIndexX_h, grid->neighborIndexX, size_in_bytes_neighbors, cudaMemcpyDeviceToHost));
- CudaSafeCall(cudaMemcpy(neighborIndexY_h, grid->neighborIndexY, size_in_bytes_neighbors, cudaMemcpyDeviceToHost));
- CudaSafeCall(cudaMemcpy(neighborIndexZ_h, grid->neighborIndexZ, size_in_bytes_neighbors, cudaMemcpyDeviceToHost));
+ int *neighborIndexX_h, *neighborIndexY_h, *neighborIndexZ_h, *neighborIndexNegative_h;
+ neighborIndexX_h = new int[grid->size];
+ neighborIndexY_h = new int[grid->size];
+ neighborIndexZ_h = new int[grid->size];
+ neighborIndexNegative_h = new int[grid->size];
+
+ CudaSafeCall(cudaMemcpy(neighborIndexX_h, grid->neighborIndexX, size_in_bytes_neighbors, cudaMemcpyDeviceToHost));
+ CudaSafeCall(cudaMemcpy(neighborIndexY_h, grid->neighborIndexY, size_in_bytes_neighbors, cudaMemcpyDeviceToHost));
+ CudaSafeCall(cudaMemcpy(neighborIndexZ_h, grid->neighborIndexZ, size_in_bytes_neighbors, cudaMemcpyDeviceToHost));
+ CudaSafeCall(cudaMemcpy(neighborIndexNegative_h, grid->neighborIndexNegative, size_in_bytes_neighbors, cudaMemcpyDeviceToHost));
+
CudaSafeCall(cudaFree(grid->neighborIndexX));
CudaSafeCall(cudaFree(grid->neighborIndexY));
CudaSafeCall(cudaFree(grid->neighborIndexZ));
+ CudaSafeCall(cudaFree(grid->neighborIndexNegative));
- grid->neighborIndexX = neighborIndexX_h;
- grid->neighborIndexY = neighborIndexY_h;
- grid->neighborIndexZ = neighborIndexZ_h;
+ grid->neighborIndexX = neighborIndexX_h;
+ grid->neighborIndexY = neighborIndexY_h;
+ grid->neighborIndexZ = neighborIndexZ_h;
+ grid->neighborIndexNegative = neighborIndexNegative_h;
CudaCheckError();
}
diff --git a/src/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp b/src/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
index 36fcb3170..43c5b7e86 100644
--- a/src/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
+++ b/src/VirtualFluids_GPU/DataStructureInitializer/GridReaderGenerator/GridGenerator.cpp
@@ -60,6 +60,7 @@ void GridGenerator::allocArrays_CoordNeighborGeo()
cudaMemoryManager->cudaAllocCoord(level);
cudaMemoryManager->cudaAllocSP(level);
+ cudaMemoryManager->cudaAllocNeighborWSB(level);
builder->getNodeValues(
para->getParH(level)->coordX_SP,
@@ -68,11 +69,13 @@ void GridGenerator::allocArrays_CoordNeighborGeo()
para->getParH(level)->neighborX_SP,
para->getParH(level)->neighborY_SP,
para->getParH(level)->neighborZ_SP,
+ para->getParH(level)->neighborWSB_SP,
para->getParH(level)->geoSP,
level);
setInitalNodeValues(numberOfNodesPerLevel, level);
+ cudaMemoryManager->cudaCopyNeighborWSB(level);
cudaMemoryManager->cudaCopySP(level);
cudaMemoryManager->cudaCopyCoord(level);
diff --git a/src/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp b/src/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
index 208db9ccb..80de21b60 100644
--- a/src/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
+++ b/src/VirtualFluids_GPU/GPU/CudaMemoryManager.cpp
@@ -95,6 +95,26 @@ void CudaMemoryManager::cudaFreeSP(int lev)
checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->neighborY_SP));
checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->neighborZ_SP));
}
+//negative neighbor (WSB)
+void CudaMemoryManager::cudaAllocNeighborWSB(int lev)
+{
+ //Host
+ checkCudaErrors( cudaMallocHost((void**) &(parameter->getParH(lev)->neighborWSB_SP ), parameter->getParH(lev)->mem_size_int_SP ));
+ //Device
+ checkCudaErrors( cudaMalloc((void**) &(parameter->getParD(lev)->neighborWSB_SP ), parameter->getParD(lev)->mem_size_int_SP ));
+ //////////////////////////////////////////////////////////////////////////
+ double tmp = (double)parameter->getParH(lev)->mem_size_int_SP;
+ parameter->setMemsizeGPU(tmp, false);
+}
+void CudaMemoryManager::cudaCopyNeighborWSB(int lev)
+{
+ //copy host to device
+ checkCudaErrors( cudaMemcpy(parameter->getParD(lev)->neighborWSB_SP, parameter->getParH(lev)->neighborWSB_SP, parameter->getParH(lev)->mem_size_int_SP , cudaMemcpyHostToDevice));
+}
+void CudaMemoryManager::cudaFreeNeighborWSB(int lev)
+{
+ checkCudaErrors( cudaFreeHost(parameter->getParH(lev)->neighborWSB_SP));
+}
//Velo
void CudaMemoryManager::cudaAllocVeloBC(int lev)
{
diff --git a/src/VirtualFluids_GPU/GPU/CudaMemoryManager.h b/src/VirtualFluids_GPU/GPU/CudaMemoryManager.h
index 09af493b3..263fb37b8 100644
--- a/src/VirtualFluids_GPU/GPU/CudaMemoryManager.h
+++ b/src/VirtualFluids_GPU/GPU/CudaMemoryManager.h
@@ -28,6 +28,10 @@ public:
void cudaCopySP(int lev);
void cudaFreeSP(int lev);
+ void cudaAllocNeighborWSB(int lev);
+ void cudaCopyNeighborWSB(int lev);
+ void cudaFreeNeighborWSB(int lev);
+
void cudaAllocVeloBC(int lev);
void cudaCopyVeloBC(int lev);
void cudaFreeVeloBC(int lev);
diff --git a/src/VirtualFluids_GPU/GPU/GPU_Interface.h b/src/VirtualFluids_GPU/GPU/GPU_Interface.h
index 245f1d2e2..d5cbff5c3 100644
--- a/src/VirtualFluids_GPU/GPU/GPU_Interface.h
+++ b/src/VirtualFluids_GPU/GPU/GPU_Interface.h
@@ -1047,6 +1047,26 @@ extern "C" void QVelDevComp27(unsigned int numberOfThreads,
unsigned int size_Mat,
bool evenOrOdd);
+extern "C" void QVelDevCompThinWalls27(unsigned int numberOfThreads,
+ int nx,
+ int ny,
+ real* vx,
+ real* vy,
+ real* vz,
+ real* DD,
+ int* k_Q,
+ real* QQ,
+ unsigned int sizeQ,
+ unsigned int kQ,
+ real om1,
+ unsigned int* geom,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int* neighborWSB,
+ unsigned int size_Mat,
+ bool evenOrOdd);
+
extern "C" void QVelDevCompZeroPress27(unsigned int numberOfThreads,
int nx,
int ny,
diff --git a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
index 080afdb02..1d8201296 100644
--- a/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
+++ b/src/VirtualFluids_GPU/GPU/GPU_Kernels.cuh
@@ -965,6 +965,40 @@ extern "C" __global__ void QVelDeviceComp27(int inx,
unsigned int size_Mat,
bool evenOrOdd);
+
+extern "C" __global__ void QVelDeviceCompThinWallsPartOne27(int inx,
+ int iny,
+ real* vx,
+ real* vy,
+ real* vz,
+ real* DD,
+ int* k_Q,
+ real* QQ,
+ unsigned int sizeQ,
+ int kQ,
+ real om1,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int size_Mat,
+ bool evenOrOdd);
+
+extern "C" __global__ void QVelDeviceCompThinWallsPartTwo27(int inx,
+ int iny,
+ real* DD,
+ int* k_Q,
+ real* QQ,
+ unsigned int sizeQ,
+ int kQ,
+ real om1,
+ unsigned int* geom,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int* neighborWSB,
+ unsigned int size_Mat,
+ bool evenOrOdd);
+
extern "C" __global__ void QVelDeviceCompZeroPress27( int inx,
int iny,
real* vx,
diff --git a/src/VirtualFluids_GPU/GPU/LBMKernel.cu b/src/VirtualFluids_GPU/GPU/LBMKernel.cu
index 6c1a9638e..02de7b7da 100644
--- a/src/VirtualFluids_GPU/GPU/LBMKernel.cu
+++ b/src/VirtualFluids_GPU/GPU/LBMKernel.cu
@@ -3920,6 +3920,77 @@ extern "C" void QVelDevComp27(unsigned int numberOfThreads,
getLastCudaError("QVelDeviceComp27 execution failed");
}
//////////////////////////////////////////////////////////////////////////
+extern "C" void QVelDevCompThinWalls27(unsigned int numberOfThreads,
+ int nx,
+ int ny,
+ real* vx,
+ real* vy,
+ real* vz,
+ real* DD,
+ int* k_Q,
+ real* QQ,
+ unsigned int sizeQ,
+ unsigned int kQ,
+ real om1,
+ unsigned int* geom,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int* neighborWSB,
+ unsigned int size_Mat,
+ bool evenOrOdd)
+{
+ int Grid = (kQ / numberOfThreads)+1;
+ int Grid1, Grid2;
+ if (Grid>512)
+ {
+ Grid1 = 512;
+ Grid2 = (Grid/Grid1)+1;
+ }
+ else
+ {
+ Grid1 = 1;
+ Grid2 = Grid;
+ }
+ dim3 gridQ(Grid1, Grid2);
+ dim3 threads(numberOfThreads, 1, 1 );
+
+ QVelDeviceCompThinWallsPartOne27<<< gridQ, threads >>> (nx,
+ ny,
+ vx,
+ vy,
+ vz,
+ DD,
+ k_Q,
+ QQ,
+ sizeQ,
+ kQ,
+ om1,
+ neighborX,
+ neighborY,
+ neighborZ,
+ size_Mat,
+ evenOrOdd);
+ getLastCudaError("QVelDeviceCompThinWallsPartOne27 execution failed");
+
+ QVelDeviceCompThinWallsPartTwo27<<< gridQ, threads >>> (nx,
+ ny,
+ DD,
+ k_Q,
+ QQ,
+ sizeQ,
+ kQ,
+ om1,
+ geom,
+ neighborX,
+ neighborY,
+ neighborZ,
+ neighborWSB,
+ size_Mat,
+ evenOrOdd);
+ getLastCudaError("QVelDeviceCompThinWallsPartTwo27 execution failed");
+}
+//////////////////////////////////////////////////////////////////////////
extern "C" void QVelDevCompZeroPress27( unsigned int numberOfThreads,
int nx,
int ny,
diff --git a/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu b/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu
index 809c65ada..16ec07ee2 100644
--- a/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu
+++ b/src/VirtualFluids_GPU/GPU/VelocityBCs27.cu
@@ -1,6 +1,839 @@
/* Device code */
#include "LBM/D3Q27.h"
#include "GPU/constant.h"
+
+//////////////////////////////////////////////////////////////////////////////
+extern "C" __global__ void QVelDeviceCompThinWallsPartOne27(int inx,
+ int iny,
+ real* vx,
+ real* vy,
+ real* vz,
+ real* DD,
+ int* k_Q,
+ real* QQ,
+ unsigned int sizeQ,
+ int kQ,
+ real om1,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int size_Mat,
+ bool evenOrOdd)
+{
+ Distributions27 D;
+ if (evenOrOdd==true)
+ {
+ D.f[dirE ] = &DD[dirE *size_Mat];
+ D.f[dirW ] = &DD[dirW *size_Mat];
+ D.f[dirN ] = &DD[dirN *size_Mat];
+ D.f[dirS ] = &DD[dirS *size_Mat];
+ D.f[dirT ] = &DD[dirT *size_Mat];
+ D.f[dirB ] = &DD[dirB *size_Mat];
+ D.f[dirNE ] = &DD[dirNE *size_Mat];
+ D.f[dirSW ] = &DD[dirSW *size_Mat];
+ D.f[dirSE ] = &DD[dirSE *size_Mat];
+ D.f[dirNW ] = &DD[dirNW *size_Mat];
+ D.f[dirTE ] = &DD[dirTE *size_Mat];
+ D.f[dirBW ] = &DD[dirBW *size_Mat];
+ D.f[dirBE ] = &DD[dirBE *size_Mat];
+ D.f[dirTW ] = &DD[dirTW *size_Mat];
+ D.f[dirTN ] = &DD[dirTN *size_Mat];
+ D.f[dirBS ] = &DD[dirBS *size_Mat];
+ D.f[dirBN ] = &DD[dirBN *size_Mat];
+ D.f[dirTS ] = &DD[dirTS *size_Mat];
+ D.f[dirZERO] = &DD[dirZERO*size_Mat];
+ D.f[dirTNE ] = &DD[dirTNE *size_Mat];
+ D.f[dirTSW ] = &DD[dirTSW *size_Mat];
+ D.f[dirTSE ] = &DD[dirTSE *size_Mat];
+ D.f[dirTNW ] = &DD[dirTNW *size_Mat];
+ D.f[dirBNE ] = &DD[dirBNE *size_Mat];
+ D.f[dirBSW ] = &DD[dirBSW *size_Mat];
+ D.f[dirBSE ] = &DD[dirBSE *size_Mat];
+ D.f[dirBNW ] = &DD[dirBNW *size_Mat];
+ }
+ else
+ {
+ D.f[dirW ] = &DD[dirE *size_Mat];
+ D.f[dirE ] = &DD[dirW *size_Mat];
+ D.f[dirS ] = &DD[dirN *size_Mat];
+ D.f[dirN ] = &DD[dirS *size_Mat];
+ D.f[dirB ] = &DD[dirT *size_Mat];
+ D.f[dirT ] = &DD[dirB *size_Mat];
+ D.f[dirSW ] = &DD[dirNE *size_Mat];
+ D.f[dirNE ] = &DD[dirSW *size_Mat];
+ D.f[dirNW ] = &DD[dirSE *size_Mat];
+ D.f[dirSE ] = &DD[dirNW *size_Mat];
+ D.f[dirBW ] = &DD[dirTE *size_Mat];
+ D.f[dirTE ] = &DD[dirBW *size_Mat];
+ D.f[dirTW ] = &DD[dirBE *size_Mat];
+ D.f[dirBE ] = &DD[dirTW *size_Mat];
+ D.f[dirBS ] = &DD[dirTN *size_Mat];
+ D.f[dirTN ] = &DD[dirBS *size_Mat];
+ D.f[dirTS ] = &DD[dirBN *size_Mat];
+ D.f[dirBN ] = &DD[dirTS *size_Mat];
+ D.f[dirZERO] = &DD[dirZERO*size_Mat];
+ D.f[dirTNE ] = &DD[dirBSW *size_Mat];
+ D.f[dirTSW ] = &DD[dirBNE *size_Mat];
+ D.f[dirTSE ] = &DD[dirBNW *size_Mat];
+ D.f[dirTNW ] = &DD[dirBSE *size_Mat];
+ D.f[dirBNE ] = &DD[dirTSW *size_Mat];
+ D.f[dirBSW ] = &DD[dirTNE *size_Mat];
+ D.f[dirBSE ] = &DD[dirTNW *size_Mat];
+ D.f[dirBNW ] = &DD[dirTSE *size_Mat];
+ }
+ ////////////////////////////////////////////////////////////////////////////////
+ const unsigned x = threadIdx.x; // Globaler x-Index
+ const unsigned y = blockIdx.x; // Globaler y-Index
+ const unsigned z = blockIdx.y; // Globaler z-Index
+
+ const unsigned nx = blockDim.x;
+ const unsigned ny = gridDim.x;
+
+ const unsigned k = nx*(ny*z + y) + x;
+ //////////////////////////////////////////////////////////////////////////
+
+ if(k<kQ)
+ {
+ ////////////////////////////////////////////////////////////////////////////////
+ real VeloX = vx[k];
+ real VeloY = vy[k];
+ real VeloZ = vz[k]; //(16.0*(u0*2.0)*bbx*bby*(grid_nx-bbx)*(grid_ny-bby))/(grid_nx*grid_nx*grid_ny*grid_ny)
+ ////////////////////////////////////////////////////////////////////////////////
+ real *q_dirE, *q_dirW, *q_dirN, *q_dirS, *q_dirT, *q_dirB,
+ *q_dirNE, *q_dirSW, *q_dirSE, *q_dirNW, *q_dirTE, *q_dirBW,
+ *q_dirBE, *q_dirTW, *q_dirTN, *q_dirBS, *q_dirBN, *q_dirTS,
+ *q_dirTNE, *q_dirTSW, *q_dirTSE, *q_dirTNW, *q_dirBNE, *q_dirBSW,
+ *q_dirBSE, *q_dirBNW;
+ q_dirE = &QQ[dirE *sizeQ];
+ q_dirW = &QQ[dirW *sizeQ];
+ q_dirN = &QQ[dirN *sizeQ];
+ q_dirS = &QQ[dirS *sizeQ];
+ q_dirT = &QQ[dirT *sizeQ];
+ q_dirB = &QQ[dirB *sizeQ];
+ q_dirNE = &QQ[dirNE *sizeQ];
+ q_dirSW = &QQ[dirSW *sizeQ];
+ q_dirSE = &QQ[dirSE *sizeQ];
+ q_dirNW = &QQ[dirNW *sizeQ];
+ q_dirTE = &QQ[dirTE *sizeQ];
+ q_dirBW = &QQ[dirBW *sizeQ];
+ q_dirBE = &QQ[dirBE *sizeQ];
+ q_dirTW = &QQ[dirTW *sizeQ];
+ q_dirTN = &QQ[dirTN *sizeQ];
+ q_dirBS = &QQ[dirBS *sizeQ];
+ q_dirBN = &QQ[dirBN *sizeQ];
+ q_dirTS = &QQ[dirTS *sizeQ];
+ q_dirTNE = &QQ[dirTNE *sizeQ];
+ q_dirTSW = &QQ[dirTSW *sizeQ];
+ q_dirTSE = &QQ[dirTSE *sizeQ];
+ q_dirTNW = &QQ[dirTNW *sizeQ];
+ q_dirBNE = &QQ[dirBNE *sizeQ];
+ q_dirBSW = &QQ[dirBSW *sizeQ];
+ q_dirBSE = &QQ[dirBSE *sizeQ];
+ q_dirBNW = &QQ[dirBNW *sizeQ];
+ ////////////////////////////////////////////////////////////////////////////////
+ //index
+ unsigned int KQK = k_Q[k];
+ unsigned int kzero= KQK;
+ unsigned int ke = KQK;
+ unsigned int kw = neighborX[KQK];
+ unsigned int kn = KQK;
+ unsigned int ks = neighborY[KQK];
+ unsigned int kt = KQK;
+ unsigned int kb = neighborZ[KQK];
+ unsigned int ksw = neighborY[kw];
+ unsigned int kne = KQK;
+ unsigned int kse = ks;
+ unsigned int knw = kw;
+ unsigned int kbw = neighborZ[kw];
+ unsigned int kte = KQK;
+ unsigned int kbe = kb;
+ unsigned int ktw = kw;
+ unsigned int kbs = neighborZ[ks];
+ unsigned int ktn = KQK;
+ unsigned int kbn = kb;
+ unsigned int kts = ks;
+ unsigned int ktse = ks;
+ unsigned int kbnw = kbw;
+ unsigned int ktnw = kw;
+ unsigned int kbse = kbs;
+ unsigned int ktsw = ksw;
+ unsigned int kbne = kb;
+ unsigned int ktne = KQK;
+ unsigned int kbsw = neighborZ[ksw];
+ ////////////////////////////////////////////////////////////////////////////////
+ real f_E, f_W, f_N, f_S, f_T, f_B, f_NE, f_SW, f_SE, f_NW, f_TE, f_BW, f_BE,
+ f_TW, f_TN, f_BS, f_BN, f_TS, f_TNE, f_TSW, f_TSE, f_TNW, f_BNE, f_BSW, f_BSE, f_BNW;
+
+ f_W = (D.f[dirE ])[ke ];
+ f_E = (D.f[dirW ])[kw ];
+ f_S = (D.f[dirN ])[kn ];
+ f_N = (D.f[dirS ])[ks ];
+ f_B = (D.f[dirT ])[kt ];
+ f_T = (D.f[dirB ])[kb ];
+ f_SW = (D.f[dirNE ])[kne ];
+ f_NE = (D.f[dirSW ])[ksw ];
+ f_NW = (D.f[dirSE ])[kse ];
+ f_SE = (D.f[dirNW ])[knw ];
+ f_BW = (D.f[dirTE ])[kte ];
+ f_TE = (D.f[dirBW ])[kbw ];
+ f_TW = (D.f[dirBE ])[kbe ];
+ f_BE = (D.f[dirTW ])[ktw ];
+ f_BS = (D.f[dirTN ])[ktn ];
+ f_TN = (D.f[dirBS ])[kbs ];
+ f_TS = (D.f[dirBN ])[kbn ];
+ f_BN = (D.f[dirTS ])[kts ];
+ f_BSW = (D.f[dirTNE ])[ktne ];
+ f_BNE = (D.f[dirTSW ])[ktsw ];
+ f_BNW = (D.f[dirTSE ])[ktse ];
+ f_BSE = (D.f[dirTNW ])[ktnw ];
+ f_TSW = (D.f[dirBNE ])[kbne ];
+ f_TNE = (D.f[dirBSW ])[kbsw ];
+ f_TNW = (D.f[dirBSE ])[kbse ];
+ f_TSE = (D.f[dirBNW ])[kbnw ];
+ ////////////////////////////////////////////////////////////////////////////////
+ real vx1, vx2, vx3, drho, feq, q;
+ drho = f_TSE + f_TNW + f_TNE + f_TSW + f_BSE + f_BNW + f_BNE + f_BSW +
+ f_BN + f_TS + f_TN + f_BS + f_BE + f_TW + f_TE + f_BW + f_SE + f_NW + f_NE + f_SW +
+ f_T + f_B + f_N + f_S + f_E + f_W + ((D.f[dirZERO])[kzero]);
+
+ vx1 = (((f_TSE - f_BNW) - (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
+ ((f_BE - f_TW) + (f_TE - f_BW)) + ((f_SE - f_NW) + (f_NE - f_SW)) +
+ (f_E - f_W)) / (one + drho);
+
+
+ vx2 = ((-(f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) - (f_TSW - f_BNE)) +
+ ((f_BN - f_TS) + (f_TN - f_BS)) + (-(f_SE - f_NW) + (f_NE - f_SW)) +
+ (f_N - f_S)) / (one + drho);
+
+ vx3 = (((f_TSE - f_BNW) + (f_TNW - f_BSE)) + ((f_TNE - f_BSW) + (f_TSW - f_BNE)) +
+ (-(f_BN - f_TS) + (f_TN - f_BS)) + ((f_TE - f_BW) - (f_BE - f_TW)) +
+ (f_T - f_B)) / (one + drho);
+
+ real cu_sq=c3o2*(vx1*vx1+vx2*vx2+vx3*vx3) * (one + drho);
+
+ //////////////////////////////////////////////////////////////////////////
+ //if (evenOrOdd==false)
+ //{
+ // D.f[dirE ] = &DD[dirE *size_Mat];
+ // D.f[dirW ] = &DD[dirW *size_Mat];
+ // D.f[dirN ] = &DD[dirN *size_Mat];
+ // D.f[dirS ] = &DD[dirS *size_Mat];
+ // D.f[dirT ] = &DD[dirT *size_Mat];
+ // D.f[dirB ] = &DD[dirB *size_Mat];
+ // D.f[dirNE ] = &DD[dirNE *size_Mat];
+ // D.f[dirSW ] = &DD[dirSW *size_Mat];
+ // D.f[dirSE ] = &DD[dirSE *size_Mat];
+ // D.f[dirNW ] = &DD[dirNW *size_Mat];
+ // D.f[dirTE ] = &DD[dirTE *size_Mat];
+ // D.f[dirBW ] = &DD[dirBW *size_Mat];
+ // D.f[dirBE ] = &DD[dirBE *size_Mat];
+ // D.f[dirTW ] = &DD[dirTW *size_Mat];
+ // D.f[dirTN ] = &DD[dirTN *size_Mat];
+ // D.f[dirBS ] = &DD[dirBS *size_Mat];
+ // D.f[dirBN ] = &DD[dirBN *size_Mat];
+ // D.f[dirTS ] = &DD[dirTS *size_Mat];
+ // D.f[dirZERO] = &DD[dirZERO*size_Mat];
+ // D.f[dirTNE ] = &DD[dirTNE *size_Mat];
+ // D.f[dirTSW ] = &DD[dirTSW *size_Mat];
+ // D.f[dirTSE ] = &DD[dirTSE *size_Mat];
+ // D.f[dirTNW ] = &DD[dirTNW *size_Mat];
+ // D.f[dirBNE ] = &DD[dirBNE *size_Mat];
+ // D.f[dirBSW ] = &DD[dirBSW *size_Mat];
+ // D.f[dirBSE ] = &DD[dirBSE *size_Mat];
+ // D.f[dirBNW ] = &DD[dirBNW *size_Mat];
+ //}
+ //else
+ //{
+ // D.f[dirW ] = &DD[dirE *size_Mat];
+ // D.f[dirE ] = &DD[dirW *size_Mat];
+ // D.f[dirS ] = &DD[dirN *size_Mat];
+ // D.f[dirN ] = &DD[dirS *size_Mat];
+ // D.f[dirB ] = &DD[dirT *size_Mat];
+ // D.f[dirT ] = &DD[dirB *size_Mat];
+ // D.f[dirSW ] = &DD[dirNE *size_Mat];
+ // D.f[dirNE ] = &DD[dirSW *size_Mat];
+ // D.f[dirNW ] = &DD[dirSE *size_Mat];
+ // D.f[dirSE ] = &DD[dirNW *size_Mat];
+ // D.f[dirBW ] = &DD[dirTE *size_Mat];
+ // D.f[dirTE ] = &DD[dirBW *size_Mat];
+ // D.f[dirTW ] = &DD[dirBE *size_Mat];
+ // D.f[dirBE ] = &DD[dirTW *size_Mat];
+ // D.f[dirBS ] = &DD[dirTN *size_Mat];
+ // D.f[dirTN ] = &DD[dirBS *size_Mat];
+ // D.f[dirTS ] = &DD[dirBN *size_Mat];
+ // D.f[dirBN ] = &DD[dirTS *size_Mat];
+ // D.f[dirZERO] = &DD[dirZERO*size_Mat];
+ // D.f[dirTNE ] = &DD[dirBSW *size_Mat];
+ // D.f[dirTSW ] = &DD[dirBNE *size_Mat];
+ // D.f[dirTSE ] = &DD[dirBNW *size_Mat];
+ // D.f[dirTNW ] = &DD[dirBSE *size_Mat];
+ // D.f[dirBNE ] = &DD[dirTSW *size_Mat];
+ // D.f[dirBSW ] = &DD[dirTNE *size_Mat];
+ // D.f[dirBSE ] = &DD[dirTNW *size_Mat];
+ // D.f[dirBNW ] = &DD[dirTSE *size_Mat];
+ //}
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //Test
+ //(D.f[dirZERO])[k]=c1o10;
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ //ToDo anders Klammern
+
+ q = q_dirE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c2over27* (drho/*+three*( vx1 )*/+c9over2*( vx1 )*( vx1 ) * (one + drho)-cu_sq);
+ (D.f[dirW])[kw]=(one-q)/(one+q)*(f_E-f_W+(f_E+f_W-two*feq*om1)/(one-om1))*c1o2+(q*(f_E+f_W)-six*c2over27*( VeloX ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
+ //(D.f[dirW])[kw]=zero;
+ }
+
+ q = q_dirW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c2over27* (drho/*+three*(-vx1 )*/+c9over2*(-vx1 )*(-vx1 ) * (one + drho)-cu_sq);
+ (D.f[dirE])[ke]=(one-q)/(one+q)*(f_W-f_E+(f_W+f_E-two*feq*om1)/(one-om1))*c1o2+(q*(f_W+f_E)-six*c2over27*(-VeloX ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
+ //(D.f[dirE])[ke]=zero;
+ }
+
+ q = q_dirN[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c2over27* (drho/*+three*( vx2 )*/+c9over2*( vx2 )*( vx2 ) * (one + drho)-cu_sq);
+ (D.f[dirS])[ks]=(one-q)/(one+q)*(f_N-f_S+(f_N+f_S-two*feq*om1)/(one-om1))*c1o2+(q*(f_N+f_S)-six*c2over27*( VeloY ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
+ //(D.f[dirS])[ks]=zero;
+ }
+
+ q = q_dirS[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c2over27* (drho/*+three*( -vx2 )*/+c9over2*( -vx2 )*( -vx2 ) * (one + drho)-cu_sq);
+ (D.f[dirN])[kn]=(one-q)/(one+q)*(f_S-f_N+(f_S+f_N-two*feq*om1)/(one-om1))*c1o2+(q*(f_S+f_N)-six*c2over27*(-VeloY ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
+ //(D.f[dirN])[kn]=zero;
+ }
+
+ q = q_dirT[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c2over27* (drho/*+three*( vx3)*/+c9over2*( vx3)*( vx3) * (one + drho)-cu_sq);
+ (D.f[dirB])[kb]=(one-q)/(one+q)*(f_T-f_B+(f_T+f_B-two*feq*om1)/(one-om1))*c1o2+(q*(f_T+f_B)-six*c2over27*( VeloZ )/* * (one + drho)*/)/(one+q);// - c2over27 * drho;
+ //(D.f[dirB])[kb]=one;
+ }
+
+ q = q_dirB[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c2over27* (drho/*+three*( -vx3)*/+c9over2*( -vx3)*( -vx3) * (one + drho)-cu_sq);
+ (D.f[dirT])[kt]=(one-q)/(one+q)*(f_B-f_T+(f_B+f_T-two*feq*om1)/(one-om1))*c1o2+(q*(f_B+f_T)-six*c2over27*(-VeloZ ) /** (one + drho)*/)/(one+q);// - c2over27 * drho;
+ //(D.f[dirT])[kt]=zero;
+ }
+
+ q = q_dirNE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( vx1+vx2 )*/+c9over2*( vx1+vx2 )*( vx1+vx2 ) * (one + drho)-cu_sq);
+ (D.f[dirSW])[ksw]=(one-q)/(one+q)*(f_NE-f_SW+(f_NE+f_SW-two*feq*om1)/(one-om1))*c1o2+(q*(f_NE+f_SW)-six*c1over54*(VeloX+VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirSW])[ksw]=zero;
+ }
+
+ q = q_dirSW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*(-vx1-vx2 )*/+c9over2*(-vx1-vx2 )*(-vx1-vx2 ) * (one + drho)-cu_sq);
+ (D.f[dirNE])[kne]=(one-q)/(one+q)*(f_SW-f_NE+(f_SW+f_NE-two*feq*om1)/(one-om1))*c1o2+(q*(f_SW+f_NE)-six*c1over54*(-VeloX-VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirNE])[kne]=zero;
+ }
+
+ q = q_dirSE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( vx1-vx2 )*/+c9over2*( vx1-vx2 )*( vx1-vx2 ) * (one + drho)-cu_sq);
+ (D.f[dirNW])[knw]=(one-q)/(one+q)*(f_SE-f_NW+(f_SE+f_NW-two*feq*om1)/(one-om1))*c1o2+(q*(f_SE+f_NW)-six*c1over54*( VeloX-VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirNW])[knw]=zero;
+ }
+
+ q = q_dirNW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*(-vx1+vx2 )*/+c9over2*(-vx1+vx2 )*(-vx1+vx2 ) * (one + drho)-cu_sq);
+ (D.f[dirSE])[kse]=(one-q)/(one+q)*(f_NW-f_SE+(f_NW+f_SE-two*feq*om1)/(one-om1))*c1o2+(q*(f_NW+f_SE)-six*c1over54*(-VeloX+VeloY) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirSE])[kse]=zero;
+ }
+
+ q = q_dirTE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( vx1 +vx3)*/+c9over2*( vx1 +vx3)*( vx1 +vx3) * (one + drho)-cu_sq);
+ (D.f[dirBW])[kbw]=(one-q)/(one+q)*(f_TE-f_BW+(f_TE+f_BW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TE+f_BW)-six*c1over54*( VeloX+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirBW])[kbw]=zero;
+ }
+
+ q = q_dirBW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*(-vx1 -vx3)*/+c9over2*(-vx1 -vx3)*(-vx1 -vx3) * (one + drho)-cu_sq);
+ (D.f[dirTE])[kte]=(one-q)/(one+q)*(f_BW-f_TE+(f_BW+f_TE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BW+f_TE)-six*c1over54*(-VeloX-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirTE])[kte]=zero;
+ }
+
+ q = q_dirBE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( vx1 -vx3)*/+c9over2*( vx1 -vx3)*( vx1 -vx3) * (one + drho)-cu_sq);
+ (D.f[dirTW])[ktw]=(one-q)/(one+q)*(f_BE-f_TW+(f_BE+f_TW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BE+f_TW)-six*c1over54*( VeloX-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirTW])[ktw]=zero;
+ }
+
+ q = q_dirTW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*(-vx1 +vx3)*/+c9over2*(-vx1 +vx3)*(-vx1 +vx3) * (one + drho)-cu_sq);
+ (D.f[dirBE])[kbe]=(one-q)/(one+q)*(f_TW-f_BE+(f_TW+f_BE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TW+f_BE)-six*c1over54*(-VeloX+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirBE])[kbe]=zero;
+ }
+
+ q = q_dirTN[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( vx2+vx3)*/+c9over2*( vx2+vx3)*( vx2+vx3) * (one + drho)-cu_sq);
+ (D.f[dirBS])[kbs]=(one-q)/(one+q)*(f_TN-f_BS+(f_TN+f_BS-two*feq*om1)/(one-om1))*c1o2+(q*(f_TN+f_BS)-six*c1over54*( VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirBS])[kbs]=zero;
+ }
+
+ q = q_dirBS[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( -vx2-vx3)*/+c9over2*( -vx2-vx3)*( -vx2-vx3) * (one + drho)-cu_sq);
+ (D.f[dirTN])[ktn]=(one-q)/(one+q)*(f_BS-f_TN+(f_BS+f_TN-two*feq*om1)/(one-om1))*c1o2+(q*(f_BS+f_TN)-six*c1over54*( -VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirTN])[ktn]=zero;
+ }
+
+ q = q_dirBN[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( vx2-vx3)*/+c9over2*( vx2-vx3)*( vx2-vx3) * (one + drho)-cu_sq);
+ (D.f[dirTS])[kts]=(one-q)/(one+q)*(f_BN-f_TS+(f_BN+f_TS-two*feq*om1)/(one-om1))*c1o2+(q*(f_BN+f_TS)-six*c1over54*( VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirTS])[kts]=zero;
+ }
+
+ q = q_dirTS[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over54* (drho/*+three*( -vx2+vx3)*/+c9over2*( -vx2+vx3)*( -vx2+vx3) * (one + drho)-cu_sq);
+ (D.f[dirBN])[kbn]=(one-q)/(one+q)*(f_TS-f_BN+(f_TS+f_BN-two*feq*om1)/(one-om1))*c1o2+(q*(f_TS+f_BN)-six*c1over54*( -VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over54 * drho;
+ //(D.f[dirBN])[kbn]=zero;
+ }
+
+ q = q_dirTNE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*( vx1+vx2+vx3)*/+c9over2*( vx1+vx2+vx3)*( vx1+vx2+vx3) * (one + drho)-cu_sq);
+ (D.f[dirBSW])[kbsw]=(one-q)/(one+q)*(f_TNE-f_BSW+(f_TNE+f_BSW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TNE+f_BSW)-six*c1over216*( VeloX+VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirBSW])[kbsw]=zero;
+ }
+
+ q = q_dirBSW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*(-vx1-vx2-vx3)*/+c9over2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3) * (one + drho)-cu_sq);
+ (D.f[dirTNE])[ktne]=(one-q)/(one+q)*(f_BSW-f_TNE+(f_BSW+f_TNE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BSW+f_TNE)-six*c1over216*(-VeloX-VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirTNE])[ktne]=zero;
+ }
+
+ q = q_dirBNE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*( vx1+vx2-vx3)*/+c9over2*( vx1+vx2-vx3)*( vx1+vx2-vx3) * (one + drho)-cu_sq);
+ (D.f[dirTSW])[ktsw]=(one-q)/(one+q)*(f_BNE-f_TSW+(f_BNE+f_TSW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BNE+f_TSW)-six*c1over216*( VeloX+VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirTSW])[ktsw]=zero;
+ }
+
+ q = q_dirTSW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*(-vx1-vx2+vx3)*/+c9over2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3) * (one + drho)-cu_sq);
+ (D.f[dirBNE])[kbne]=(one-q)/(one+q)*(f_TSW-f_BNE+(f_TSW+f_BNE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TSW+f_BNE)-six*c1over216*(-VeloX-VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirBNE])[kbne]=zero;
+ }
+
+ q = q_dirTSE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*( vx1-vx2+vx3)*/+c9over2*( vx1-vx2+vx3)*( vx1-vx2+vx3) * (one + drho)-cu_sq);
+ (D.f[dirBNW])[kbnw]=(one-q)/(one+q)*(f_TSE-f_BNW+(f_TSE+f_BNW-two*feq*om1)/(one-om1))*c1o2+(q*(f_TSE+f_BNW)-six*c1over216*( VeloX-VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirBNW])[kbnw]=zero;
+ }
+
+ q = q_dirBNW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*(-vx1+vx2-vx3)*/+c9over2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3) * (one + drho)-cu_sq);
+ (D.f[dirTSE])[ktse]=(one-q)/(one+q)*(f_BNW-f_TSE+(f_BNW+f_TSE-two*feq*om1)/(one-om1))*c1o2+(q*(f_BNW+f_TSE)-six*c1over216*(-VeloX+VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirTSE])[ktse]=zero;
+ }
+
+ q = q_dirBSE[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*( vx1-vx2-vx3)*/+c9over2*( vx1-vx2-vx3)*( vx1-vx2-vx3) * (one + drho)-cu_sq);
+ (D.f[dirTNW])[ktnw]=(one-q)/(one+q)*(f_BSE-f_TNW+(f_BSE+f_TNW-two*feq*om1)/(one-om1))*c1o2+(q*(f_BSE+f_TNW)-six*c1over216*( VeloX-VeloY-VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirTNW])[ktnw]=zero;
+ }
+
+ q = q_dirTNW[k];
+ if (q>=zero && q<=one)
+ {
+ feq=c1over216*(drho/*+three*(-vx1+vx2+vx3)*/+c9over2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3) * (one + drho)-cu_sq);
+ (D.f[dirBSE])[kbse]=(one-q)/(one+q)*(f_TNW-f_BSE+(f_TNW+f_BSE-two*feq*om1)/(one-om1))*c1o2+(q*(f_TNW+f_BSE)-six*c1over216*(-VeloX+VeloY+VeloZ) /** (one + drho)*/)/(one+q);// - c1over216 * drho;
+ //(D.f[dirBSE])[kbse]=zero;
+ }
+ }
+}
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+//////////////////////////////////////////////////////////////////////////////
+extern "C" __global__ void QVelDeviceCompThinWallsPartTwo27(int inx,
+ int iny,
+ real* DD,
+ int* k_Q,
+ real* QQ,
+ unsigned int sizeQ,
+ int kQ,
+ real om1,
+ unsigned int* geom,
+ unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int* neighborWSB,
+ unsigned int size_Mat,
+ bool evenOrOdd)
+{
+ ////////////////////////////////////////////////////////////////////////////////
+ const unsigned x = threadIdx.x; // Globaler x-Index
+ const unsigned y = blockIdx.x; // Globaler y-Index
+ const unsigned z = blockIdx.y; // Globaler z-Index
+
+ const unsigned nx = blockDim.x;
+ const unsigned ny = gridDim.x;
+
+ const unsigned k = nx*(ny*z + y) + x;
+ //////////////////////////////////////////////////////////////////////////
+
+ if(k<kQ)
+ {
+ ////////////////////////////////////////////////////////////////////////////////
+ real *q_dirE, *q_dirW, *q_dirN, *q_dirS, *q_dirT, *q_dirB,
+ *q_dirNE, *q_dirSW, *q_dirSE, *q_dirNW, *q_dirTE, *q_dirBW,
+ *q_dirBE, *q_dirTW, *q_dirTN, *q_dirBS, *q_dirBN, *q_dirTS,
+ *q_dirTNE, *q_dirTSW, *q_dirTSE, *q_dirTNW, *q_dirBNE, *q_dirBSW,
+ *q_dirBSE, *q_dirBNW;
+ q_dirE = &QQ[dirE *sizeQ];
+ q_dirW = &QQ[dirW *sizeQ];
+ q_dirN = &QQ[dirN *sizeQ];
+ q_dirS = &QQ[dirS *sizeQ];
+ q_dirT = &QQ[dirT *sizeQ];
+ q_dirB = &QQ[dirB *sizeQ];
+ q_dirNE = &QQ[dirNE *sizeQ];
+ q_dirSW = &QQ[dirSW *sizeQ];
+ q_dirSE = &QQ[dirSE *sizeQ];
+ q_dirNW = &QQ[dirNW *sizeQ];
+ q_dirTE = &QQ[dirTE *sizeQ];
+ q_dirBW = &QQ[dirBW *sizeQ];
+ q_dirBE = &QQ[dirBE *sizeQ];
+ q_dirTW = &QQ[dirTW *sizeQ];
+ q_dirTN = &QQ[dirTN *sizeQ];
+ q_dirBS = &QQ[dirBS *sizeQ];
+ q_dirBN = &QQ[dirBN *sizeQ];
+ q_dirTS = &QQ[dirTS *sizeQ];
+ q_dirTNE = &QQ[dirTNE *sizeQ];
+ q_dirTSW = &QQ[dirTSW *sizeQ];
+ q_dirTSE = &QQ[dirTSE *sizeQ];
+ q_dirTNW = &QQ[dirTNW *sizeQ];
+ q_dirBNE = &QQ[dirBNE *sizeQ];
+ q_dirBSW = &QQ[dirBSW *sizeQ];
+ q_dirBSE = &QQ[dirBSE *sizeQ];
+ q_dirBNW = &QQ[dirBNW *sizeQ];
+ ////////////////////////////////////////////////////////////////////////////////
+ //index
+ unsigned int KQK = k_Q[k];
+ unsigned int kzero= KQK;
+ unsigned int ke = KQK;
+ unsigned int kw = neighborX[KQK];
+ unsigned int kn = KQK;
+ unsigned int ks = neighborY[KQK];
+ unsigned int kt = KQK;
+ unsigned int kb = neighborZ[KQK];
+ unsigned int ksw = neighborY[kw];
+ unsigned int kne = KQK;
+ unsigned int kse = ks;
+ unsigned int knw = kw;
+ unsigned int kbw = neighborZ[kw];
+ unsigned int kte = KQK;
+ unsigned int kbe = kb;
+ unsigned int ktw = kw;
+ unsigned int kbs = neighborZ[ks];
+ unsigned int ktn = KQK;
+ unsigned int kbn = kb;
+ unsigned int kts = ks;
+ unsigned int ktse = ks;
+ unsigned int kbnw = kbw;
+ unsigned int ktnw = kw;
+ unsigned int kbse = kbs;
+ unsigned int ktsw = ksw;
+ unsigned int kbne = kb;
+ unsigned int ktne = KQK;
+ unsigned int kbsw = neighborZ[ksw];
+ ////////////////////////////////////////////////////////////////////////////////
+ //anti ET intermediate steps
+ unsigned int kmmm = neighborWSB[KQK]; // -1 -1 -1
+ unsigned int k0mm = neighborX[kmmm]; // 0 -1 -1
+ unsigned int km0m = neighborY[kmmm]; // -1 0 -1
+ unsigned int kmm0 = neighborZ[kmmm]; // -1 -1 0
+ unsigned int k0m0 = neighborX[kmm0]; // 0 -1 0
+ unsigned int km00 = neighborY[kmm0]; // -1 0 0
+ /////////////////////////////////////////////////
+ //final indices for anti ET
+ unsigned int kpmm = neighborX[k0mm]; // 1 -1 -1
+ unsigned int kmpm = neighborY[km0m]; // -1 1 -1
+ unsigned int kmmp = neighborZ[kmm0]; // -1 -1 1
+ unsigned int kmp0 = neighborY[km00]; // -1 1 0
+ unsigned int km0p = neighborZ[km00]; // -1 0 1
+ unsigned int k0mp = neighborZ[k0m0]; // 0 -1 1
+ ////////////////////////////////////////////////////////////////////////////////
+ Distributions27 D, DN;
+ if (evenOrOdd == true)
+ {
+ D.f[dirE] = &DD[dirE *size_Mat];
+ D.f[dirW] = &DD[dirW *size_Mat];
+ D.f[dirN] = &DD[dirN *size_Mat];
+ D.f[dirS] = &DD[dirS *size_Mat];
+ D.f[dirT] = &DD[dirT *size_Mat];
+ D.f[dirB] = &DD[dirB *size_Mat];
+ D.f[dirNE] = &DD[dirNE *size_Mat];
+ D.f[dirSW] = &DD[dirSW *size_Mat];
+ D.f[dirSE] = &DD[dirSE *size_Mat];
+ D.f[dirNW] = &DD[dirNW *size_Mat];
+ D.f[dirTE] = &DD[dirTE *size_Mat];
+ D.f[dirBW] = &DD[dirBW *size_Mat];
+ D.f[dirBE] = &DD[dirBE *size_Mat];
+ D.f[dirTW] = &DD[dirTW *size_Mat];
+ D.f[dirTN] = &DD[dirTN *size_Mat];
+ D.f[dirBS] = &DD[dirBS *size_Mat];
+ D.f[dirBN] = &DD[dirBN *size_Mat];
+ D.f[dirTS] = &DD[dirTS *size_Mat];
+ D.f[dirZERO] = &DD[dirZERO*size_Mat];
+ D.f[dirTNE] = &DD[dirTNE *size_Mat];
+ D.f[dirTSW] = &DD[dirTSW *size_Mat];
+ D.f[dirTSE] = &DD[dirTSE *size_Mat];
+ D.f[dirTNW] = &DD[dirTNW *size_Mat];
+ D.f[dirBNE] = &DD[dirBNE *size_Mat];
+ D.f[dirBSW] = &DD[dirBSW *size_Mat];
+ D.f[dirBSE] = &DD[dirBSE *size_Mat];
+ D.f[dirBNW] = &DD[dirBNW *size_Mat];
+ }
+ else
+ {
+ D.f[dirW] = &DD[dirE *size_Mat];
+ D.f[dirE] = &DD[dirW *size_Mat];
+ D.f[dirS] = &DD[dirN *size_Mat];
+ D.f[dirN] = &DD[dirS *size_Mat];
+ D.f[dirB] = &DD[dirT *size_Mat];
+ D.f[dirT] = &DD[dirB *size_Mat];
+ D.f[dirSW] = &DD[dirNE *size_Mat];
+ D.f[dirNE] = &DD[dirSW *size_Mat];
+ D.f[dirNW] = &DD[dirSE *size_Mat];
+ D.f[dirSE] = &DD[dirNW *size_Mat];
+ D.f[dirBW] = &DD[dirTE *size_Mat];
+ D.f[dirTE] = &DD[dirBW *size_Mat];
+ D.f[dirTW] = &DD[dirBE *size_Mat];
+ D.f[dirBE] = &DD[dirTW *size_Mat];
+ D.f[dirBS] = &DD[dirTN *size_Mat];
+ D.f[dirTN] = &DD[dirBS *size_Mat];
+ D.f[dirTS] = &DD[dirBN *size_Mat];
+ D.f[dirBN] = &DD[dirTS *size_Mat];
+ D.f[dirZERO] = &DD[dirZERO*size_Mat];
+ D.f[dirTNE] = &DD[dirBSW *size_Mat];
+ D.f[dirTSW] = &DD[dirBNE *size_Mat];
+ D.f[dirTSE] = &DD[dirBNW *size_Mat];
+ D.f[dirTNW] = &DD[dirBSE *size_Mat];
+ D.f[dirBNE] = &DD[dirTSW *size_Mat];
+ D.f[dirBSW] = &DD[dirTNE *size_Mat];
+ D.f[dirBSE] = &DD[dirTNW *size_Mat];
+ D.f[dirBNW] = &DD[dirTSE *size_Mat];
+ }
+ if (evenOrOdd==false)
+ {
+ DN.f[dirE ] = &DD[dirE *size_Mat];
+ DN.f[dirW ] = &DD[dirW *size_Mat];
+ DN.f[dirN ] = &DD[dirN *size_Mat];
+ DN.f[dirS ] = &DD[dirS *size_Mat];
+ DN.f[dirT ] = &DD[dirT *size_Mat];
+ DN.f[dirB ] = &DD[dirB *size_Mat];
+ DN.f[dirNE ] = &DD[dirNE *size_Mat];
+ DN.f[dirSW ] = &DD[dirSW *size_Mat];
+ DN.f[dirSE ] = &DD[dirSE *size_Mat];
+ DN.f[dirNW ] = &DD[dirNW *size_Mat];
+ DN.f[dirTE ] = &DD[dirTE *size_Mat];
+ DN.f[dirBW ] = &DD[dirBW *size_Mat];
+ DN.f[dirBE ] = &DD[dirBE *size_Mat];
+ DN.f[dirTW ] = &DD[dirTW *size_Mat];
+ DN.f[dirTN ] = &DD[dirTN *size_Mat];
+ DN.f[dirBS ] = &DD[dirBS *size_Mat];
+ DN.f[dirBN ] = &DD[dirBN *size_Mat];
+ DN.f[dirTS ] = &DD[dirTS *size_Mat];
+ DN.f[dirZERO] = &DD[dirZERO*size_Mat];
+ DN.f[dirTNE ] = &DD[dirTNE *size_Mat];
+ DN.f[dirTSW ] = &DD[dirTSW *size_Mat];
+ DN.f[dirTSE ] = &DD[dirTSE *size_Mat];
+ DN.f[dirTNW ] = &DD[dirTNW *size_Mat];
+ DN.f[dirBNE ] = &DD[dirBNE *size_Mat];
+ DN.f[dirBSW ] = &DD[dirBSW *size_Mat];
+ DN.f[dirBSE ] = &DD[dirBSE *size_Mat];
+ DN.f[dirBNW ] = &DD[dirBNW *size_Mat];
+ }
+ else
+ {
+ DN.f[dirW ] = &DD[dirE *size_Mat];
+ DN.f[dirE ] = &DD[dirW *size_Mat];
+ DN.f[dirS ] = &DD[dirN *size_Mat];
+ DN.f[dirN ] = &DD[dirS *size_Mat];
+ DN.f[dirB ] = &DD[dirT *size_Mat];
+ DN.f[dirT ] = &DD[dirB *size_Mat];
+ DN.f[dirSW ] = &DD[dirNE *size_Mat];
+ DN.f[dirNE ] = &DD[dirSW *size_Mat];
+ DN.f[dirNW ] = &DD[dirSE *size_Mat];
+ DN.f[dirSE ] = &DD[dirNW *size_Mat];
+ DN.f[dirBW ] = &DD[dirTE *size_Mat];
+ DN.f[dirTE ] = &DD[dirBW *size_Mat];
+ DN.f[dirTW ] = &DD[dirBE *size_Mat];
+ DN.f[dirBE ] = &DD[dirTW *size_Mat];
+ DN.f[dirBS ] = &DD[dirTN *size_Mat];
+ DN.f[dirTN ] = &DD[dirBS *size_Mat];
+ DN.f[dirTS ] = &DD[dirBN *size_Mat];
+ DN.f[dirBN ] = &DD[dirTS *size_Mat];
+ DN.f[dirZERO] = &DD[dirZERO*size_Mat];
+ DN.f[dirTNE ] = &DD[dirBSW *size_Mat];
+ DN.f[dirTSW ] = &DD[dirBNE *size_Mat];
+ DN.f[dirTSE ] = &DD[dirBNW *size_Mat];
+ DN.f[dirTNW ] = &DD[dirBSE *size_Mat];
+ DN.f[dirBNE ] = &DD[dirTSW *size_Mat];
+ DN.f[dirBSW ] = &DD[dirTNE *size_Mat];
+ DN.f[dirBSE ] = &DD[dirTNW *size_Mat];
+ DN.f[dirBNW ] = &DD[dirTSE *size_Mat];
+ }
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //directions allways exchange
+ //(-1 -1 -1) (-1 0 0) ( 0 -1 0) ( 0 0 -1) (-1 -1 0) (-1 0 -1) ( 0 -1 -1) ( 1 1 -1) ( 1 -1 1) (-1 1 1) ( 1 -1 0) ( 1 0 -1) ( 0 1 -1)
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //directions exchange if solid neighbor
+ //( 1 1 1) ( 1 0 0) ( 0 1 0) ( 0 0 1) ( 1 1 0) ( 1 0 1) ( 0 1 1) (-1 -1 1) (-1 1 -1) ( 1 -1 -1) (-1 1 0) (-1 0 1) ( 0 -1 1)
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ real q, tmp;
+ q = q_dirE[k]; if (q>=zero && q<=one){ if (geom[kw ] < GEO_FLUID){tmp = (DN.f[dirW ])[kw ]; (DN.f[dirW ])[kw ]=(D.f[dirW ])[kw ]; (D.f[dirW ])[kw ]=tmp;}}
+ q = q_dirW[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirE ])[ke ]; (DN.f[dirE ])[ke ]=(D.f[dirE ])[ke ]; (D.f[dirE ])[ke ]=tmp;}}
+ q = q_dirN[k]; if (q>=zero && q<=one){ if (geom[ks ] < GEO_FLUID){tmp = (DN.f[dirS ])[ks ]; (DN.f[dirS ])[ks ]=(D.f[dirS ])[ks ]; (D.f[dirS ])[ks ]=tmp;}}
+ q = q_dirS[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirN ])[kn ]; (DN.f[dirN ])[kn ]=(D.f[dirN ])[kn ]; (D.f[dirN ])[kn ]=tmp;}}
+ q = q_dirT[k]; if (q>=zero && q<=one){ if (geom[kb ] < GEO_FLUID){tmp = (DN.f[dirB ])[kb ]; (DN.f[dirB ])[kb ]=(D.f[dirB ])[kb ]; (D.f[dirB ])[kb ]=tmp;}}
+ q = q_dirB[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirT ])[kt ]; (DN.f[dirT ])[kt ]=(D.f[dirT ])[kt ]; (D.f[dirT ])[kt ]=tmp;}}
+ q = q_dirNE[k]; if (q>=zero && q<=one){ if (geom[ksw ] < GEO_FLUID){tmp = (DN.f[dirSW ])[ksw ]; (DN.f[dirSW ])[ksw ]=(D.f[dirSW ])[ksw ]; (D.f[dirSW ])[ksw ]=tmp;}}
+ q = q_dirSW[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirNE ])[kne ]; (DN.f[dirNE ])[kne ]=(D.f[dirNE ])[kne ]; (D.f[dirNE ])[kne ]=tmp;}}
+ q = q_dirSE[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirNW ])[knw ]; (DN.f[dirNW ])[knw ]=(D.f[dirNW ])[knw ]; (D.f[dirNW ])[knw ]=tmp;}}
+ q = q_dirNW[k]; if (q>=zero && q<=one){ if (geom[kmp0] < GEO_FLUID){tmp = (DN.f[dirSE ])[kse ]; (DN.f[dirSE ])[kse ]=(D.f[dirSE ])[kse ]; (D.f[dirSE ])[kse ]=tmp;}}
+ q = q_dirTE[k]; if (q>=zero && q<=one){ if (geom[kbw ] < GEO_FLUID){tmp = (DN.f[dirBW ])[kbw ]; (DN.f[dirBW ])[kbw ]=(D.f[dirBW ])[kbw ]; (D.f[dirBW ])[kbw ]=tmp;}}
+ q = q_dirBW[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirTE ])[kte ]; (DN.f[dirTE ])[kte ]=(D.f[dirTE ])[kte ]; (D.f[dirTE ])[kte ]=tmp;}}
+ q = q_dirBE[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirTW ])[ktw ]; (DN.f[dirTW ])[ktw ]=(D.f[dirTW ])[ktw ]; (D.f[dirTW ])[ktw ]=tmp;}}
+ q = q_dirTW[k]; if (q>=zero && q<=one){ if (geom[km0p] < GEO_FLUID){tmp = (DN.f[dirBE ])[kbe ]; (DN.f[dirBE ])[kbe ]=(D.f[dirBE ])[kbe ]; (D.f[dirBE ])[kbe ]=tmp;}}
+ q = q_dirTN[k]; if (q>=zero && q<=one){ if (geom[kbs ] < GEO_FLUID){tmp = (DN.f[dirBS ])[kbs ]; (DN.f[dirBS ])[kbs ]=(D.f[dirBS ])[kbs ]; (D.f[dirBS ])[kbs ]=tmp;}}
+ q = q_dirBS[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirTN ])[ktn ]; (DN.f[dirTN ])[ktn ]=(D.f[dirTN ])[ktn ]; (D.f[dirTN ])[ktn ]=tmp;}}
+ q = q_dirBN[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirTS ])[kts ]; (DN.f[dirTS ])[kts ]=(D.f[dirTS ])[kts ]; (D.f[dirTS ])[kts ]=tmp;}}
+ q = q_dirTS[k]; if (q>=zero && q<=one){ if (geom[k0mp] < GEO_FLUID){tmp = (DN.f[dirBN ])[kbn ]; (DN.f[dirBN ])[kbn ]=(D.f[dirBN ])[kbn ]; (D.f[dirBN ])[kbn ]=tmp;}}
+ q = q_dirTNE[k]; if (q>=zero && q<=one){ if (geom[kbsw] < GEO_FLUID){tmp = (DN.f[dirBSW])[kbsw]; (DN.f[dirBSW])[kbsw]=(D.f[dirBSW])[kbsw]; (D.f[dirBSW])[kbsw]=tmp;}}
+ q = q_dirBSW[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirTNE])[ktne]; (DN.f[dirTNE])[ktne]=(D.f[dirTNE])[ktne]; (D.f[dirTNE])[ktne]=tmp;}}
+ q = q_dirBNE[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirTSW])[ktsw]; (DN.f[dirTSW])[ktsw]=(D.f[dirTSW])[ktsw]; (D.f[dirTSW])[ktsw]=tmp;}}
+ q = q_dirTSW[k]; if (q>=zero && q<=one){ if (geom[kmmp] < GEO_FLUID){tmp = (DN.f[dirBNE])[kbne]; (DN.f[dirBNE])[kbne]=(D.f[dirBNE])[kbne]; (D.f[dirBNE])[kbne]=tmp;}}
+ q = q_dirTSE[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirBNW])[kbnw]; (DN.f[dirBNW])[kbnw]=(D.f[dirBNW])[kbnw]; (D.f[dirBNW])[kbnw]=tmp;}}
+ q = q_dirBNW[k]; if (q>=zero && q<=one){ if (geom[kmpm] < GEO_FLUID){tmp = (DN.f[dirTSE])[ktse]; (DN.f[dirTSE])[ktse]=(D.f[dirTSE])[ktse]; (D.f[dirTSE])[ktse]=tmp;}}
+ q = q_dirBSE[k]; if (q>=zero && q<=one){ if (geom[kpmm] < GEO_FLUID){tmp = (DN.f[dirTNW])[ktnw]; (DN.f[dirTNW])[ktnw]=(D.f[dirTNW])[ktnw]; (D.f[dirTNW])[ktnw]=tmp;}}
+ q = q_dirTNW[k]; if (q>=zero && q<=one){ {tmp = (DN.f[dirBSE])[kbse]; (DN.f[dirBSE])[kbse]=(D.f[dirBSE])[kbse]; (D.f[dirBSE])[kbse]=tmp;}}
+ }
+}
+////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
//////////////////////////////////////////////////////////////////////////////
extern "C" __global__ void QVelDeviceCompPlusSlip27(int inx,
diff --git a/src/VirtualFluids_GPU/LBM/Simulation.cpp b/src/VirtualFluids_GPU/LBM/Simulation.cpp
index 6057ed193..5ecfbbaae 100644
--- a/src/VirtualFluids_GPU/LBM/Simulation.cpp
+++ b/src/VirtualFluids_GPU/LBM/Simulation.cpp
@@ -1256,13 +1256,30 @@ else
// ////////////////////////////////////////////////////////////////////////////
if ((para->getParD(0)->QGeom.kQ) > 0 && (para->getIsGeometryValues()))
{
- QVelDevCompZeroPress27( para->getParD(0)->numberofthreads, para->getParD(0)->nx, para->getParD(0)->ny,
- para->getParD(0)->QGeom.Vx, para->getParD(0)->QGeom.Vy, para->getParD(0)->QGeom.Vz,
- para->getParD(0)->d0SP.f[0], para->getParD(0)->QGeom.k, para->getParD(0)->QGeom.q27[0],
- para->getParD(0)->QGeom.kQ, para->getParD(0)->QGeom.kQ, para->getParD(0)->omega,
- para->getParD(0)->neighborX_SP, para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
- para->getParD(0)->size_Mat_SP, para->getParD(0)->evenOrOdd);
- getLastCudaError("QVelDevCompZeroPress27 execution failed");
+ //////////////////////////////////////////////////////////////////////////
+ //QVelDevCompZeroPress27( para->getParD(0)->numberofthreads, para->getParD(0)->nx, para->getParD(0)->ny,
+ // para->getParD(0)->QGeom.Vx, para->getParD(0)->QGeom.Vy, para->getParD(0)->QGeom.Vz,
+ // para->getParD(0)->d0SP.f[0], para->getParD(0)->QGeom.k, para->getParD(0)->QGeom.q27[0],
+ // para->getParD(0)->QGeom.kQ, para->getParD(0)->QGeom.kQ, para->getParD(0)->omega,
+ // para->getParD(0)->neighborX_SP, para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+ // para->getParD(0)->size_Mat_SP, para->getParD(0)->evenOrOdd);
+ //getLastCudaError("QVelDevCompZeroPress27 execution failed");
+ //////////////////////////////////////////////////////////////////////////
+ // QVelDevComp27(para->getParD(0)->numberofthreads, para->getParD(0)->nx, para->getParD(0)->ny,
+ //para->getParD(0)->QGeom.Vx, para->getParD(0)->QGeom.Vy, para->getParD(0)->QGeom.Vz,
+ //para->getParD(0)->d0SP.f[0], para->getParD(0)->QGeom.k, para->getParD(0)->QGeom.q27[0],
+ //para->getParD(0)->QGeom.kQ, para->getParD(0)->QGeom.kQ, para->getParD(0)->omega,
+ //para->getParD(0)->neighborX_SP, para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP,
+ //para->getParD(0)->size_Mat_SP, para->getParD(0)->evenOrOdd);
+ // getLastCudaError("QVelDevComp27 execution failed");
+ //////////////////////////////////////////////////////////////////////////
+ QVelDevCompThinWalls27(para->getParD(0)->numberofthreads, para->getParD(0)->nx, para->getParD(0)->ny,
+ para->getParD(0)->QGeom.Vx, para->getParD(0)->QGeom.Vy, para->getParD(0)->QGeom.Vz,
+ para->getParD(0)->d0SP.f[0], para->getParD(0)->QGeom.k, para->getParD(0)->QGeom.q27[0],
+ para->getParD(0)->QGeom.kQ, para->getParD(0)->QGeom.kQ, para->getParD(0)->omega, para->getParD(0)->geoSP,
+ para->getParD(0)->neighborX_SP, para->getParD(0)->neighborY_SP, para->getParD(0)->neighborZ_SP, para->getParD(0)->neighborWSB_SP,
+ para->getParD(0)->size_Mat_SP, para->getParD(0)->evenOrOdd);
+ getLastCudaError("QVelDevComp27 execution failed");
//////////////////////////////////////////////////////////////////////////
//QVelDevCompPlusSlip27( para->getParD(0)->numberofthreads, para->getParD(0)->nx, para->getParD(0)->ny,
// para->getParD(0)->QGeom.Vx, para->getParD(0)->QGeom.Vy, para->getParD(0)->QGeom.Vz,
@@ -2681,6 +2698,7 @@ else
//para->cudaFreeFull(lev);
para->cudaFreeCoord(lev);
para->cudaFreeSP(lev);
+ para->cudaFreeNeighborWSB(lev);
if (para->getCalcMedian())
{
para->cudaFreeMedianSP(lev);
diff --git a/targets/apps/HULC/main.cpp b/targets/apps/HULC/main.cpp
index 219ca56f2..9485ed252 100644
--- a/targets/apps/HULC/main.cpp
+++ b/targets/apps/HULC/main.cpp
@@ -341,12 +341,12 @@ void multipleLevel(const std::string& configPath)
// Testing Thin Wall
//////////////////////////////////////////////////////////////////////////
- real dx = 0.05;
+ real dx = 0.0125;
real vx = 0.05;
- TriangularMesh* triangularMesh = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/ThinWallTest.stl");
+ TriangularMesh* triangularMesh = TriangularMesh::make("C:/Users/lenz/Desktop/Work/gridGenerator/stl/ThinWallTest2.stl");
- gridBuilder->addCoarseGrid(-2, -2, -2+0.5*dx, 3, 2, 2+0.5*dx, dx);
+ gridBuilder->addCoarseGrid(-2, -1.5, -0.5, 4, 1.5, 0.5, dx);
gridBuilder->addGeometry( triangularMesh );
@@ -373,7 +373,7 @@ void multipleLevel(const std::string& configPath)
gridBuilder->setPeriodicBoundaryCondition(false, false, false);
- gridBuilder->buildGrids(LBM); // buildGrids() has to be called before setting the BCs!!!!
+ gridBuilder->buildGrids(LBM, true); // buildGrids() has to be called before setting the BCs!!!!
///////////////////////////////////////////////////////////////////////////
@@ -413,7 +413,7 @@ void multipleLevel(const std::string& configPath)
gridBuilder->writeGridsToVtk("C:/Users/lenz/Desktop/Work/gridGenerator/out/Test_");
//gridBuilder->writeGridsToVtk("M:/TestGridGeneration/results/CylinderTest_");
- //gridBuilder->writeArrows ("C:/Users/lenz/Desktop/Work/gridGenerator/Sphere_Arrow");
+ gridBuilder->writeArrows ("C:/Users/lenz/Desktop/Work/gridGenerator/Test_Arrow");
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
--
GitLab