From 87dd4319a0a5ff3b1a7ca0be647aa50369486f81 Mon Sep 17 00:00:00 2001
From: "LEGOLAS\\lenz" <lenz@irmb.tu-bs.de>
Date: Thu, 5 Dec 2019 13:40:12 +0100
Subject: [PATCH] implements non-equilibrium initialization
---
src/Core/Input/ConfigData/ConfigData.h | 2 +
src/Core/Input/ConfigData/ConfigDataImp.cpp | 16 +
src/Core/Input/ConfigData/ConfigDataImp.h | 5 +
.../ConfigFileReader/ConfigFileReader.cpp | 3 +
src/VirtualFluids_GPU/LBM/LB.h | 2 +-
src/VirtualFluids_GPU/Parameter/Parameter.cpp | 13 +
src/VirtualFluids_GPU/Parameter/Parameter.h | 2 +
.../InitCompSP27/InitCompSP27.cu | 47 ++-
.../InitCompSP27/InitCompSP27_Device.cu | 298 ++++++++++++++++++
.../InitCompSP27/InitCompSP27_Device.cuh | 14 +
targets/apps/LBM/TGV_3D/TGV_3D.cpp | 2 +
11 files changed, 391 insertions(+), 13 deletions(-)
diff --git a/src/Core/Input/ConfigData/ConfigData.h b/src/Core/Input/ConfigData/ConfigData.h
index 66d1c4800..3ea009cb5 100644
--- a/src/Core/Input/ConfigData/ConfigData.h
+++ b/src/Core/Input/ConfigData/ConfigData.h
@@ -26,6 +26,7 @@ public:
virtual bool getStreetVelocityFile() = 0;
virtual bool getUseMeasurePoints() = 0;
virtual bool getUseWale() = 0;
+ virtual bool getUseInitNeq() = 0;
virtual bool getSimulatePorousMedia() = 0;
virtual uint getD3Qxx() = 0;
virtual uint getTEnd() = 0;
@@ -107,6 +108,7 @@ public:
virtual bool isStreetVelocityFileInConfigFile() = 0;
virtual bool isUseMeasurePointsInConfigFile() = 0;
virtual bool isUseWaleInConfigFile() = 0;
+ virtual bool isUseInitNeqInConfigFile() = 0;
virtual bool isSimulatePorousMediaInConfigFile() = 0;
virtual bool isD3QxxInConfigFile() = 0;
virtual bool isTEndInConfigFile() = 0;
diff --git a/src/Core/Input/ConfigData/ConfigDataImp.cpp b/src/Core/Input/ConfigData/ConfigDataImp.cpp
index 4221c9155..0b6b88e91 100644
--- a/src/Core/Input/ConfigData/ConfigDataImp.cpp
+++ b/src/Core/Input/ConfigData/ConfigDataImp.cpp
@@ -174,6 +174,11 @@ bool ConfigDataImp::getUseWale()
return this->useWale;
}
+bool ConfigDataImp::getUseInitNeq()
+{
+ return this->useInitNeq;
+}
+
bool ConfigDataImp::getSimulatePorousMedia()
{
return this->simulatePorousMedia;
@@ -576,6 +581,12 @@ void ConfigDataImp::setUseWale(bool useWale)
this->isUseWale = true;
}
+void ConfigDataImp::setUseInitNeq(bool useInitNeq)
+{
+ this->useInitNeq = useInitNeq;
+ this->isUseInitNeq = true;
+}
+
void ConfigDataImp::setSimulatePorousMedia(bool simulatePorousMedia)
{
this->simulatePorousMedia = simulatePorousMedia;
@@ -987,6 +998,11 @@ bool ConfigDataImp::isUseWaleInConfigFile()
return this->isUseWale;
}
+bool ConfigDataImp::isUseInitNeqInConfigFile()
+{
+ return this->isUseInitNeq;
+}
+
bool ConfigDataImp::isSimulatePorousMediaInConfigFile()
{
return this->isSimulatePorousMedia;
diff --git a/src/Core/Input/ConfigData/ConfigDataImp.h b/src/Core/Input/ConfigData/ConfigDataImp.h
index b2c6c5158..99f62fc9b 100644
--- a/src/Core/Input/ConfigData/ConfigDataImp.h
+++ b/src/Core/Input/ConfigData/ConfigDataImp.h
@@ -29,6 +29,7 @@ public:
bool getStreetVelocityFile();
bool getUseMeasurePoints();
bool getUseWale();
+ bool getUseInitNeq();
bool getSimulatePorousMedia();
uint getD3Qxx();
uint getTEnd();
@@ -110,6 +111,7 @@ public:
void setStreetVelocityFile(bool streetVelocityFile);
void setUseMeasurePoints(bool useMeasurePoints);
void setUseWale(bool useWale);
+ void setUseInitNeq(bool useInitNeq);
void setSimulatePorousMedia(bool simulatePorousMedia);
void setD3Qxx(uint d3Qxx);
void setTEnd(uint tEnd);
@@ -192,6 +194,7 @@ public:
bool isStreetVelocityFileInConfigFile();
bool isUseMeasurePointsInConfigFile();
bool isUseWaleInConfigFile();
+ bool isUseInitNeqInConfigFile();
bool isSimulatePorousMediaInConfigFile();
bool isD3QxxInConfigFile();
bool isTEndInConfigFile();
@@ -276,6 +279,7 @@ private:
bool streetVelocityFile;
bool useMeasurePoints;
bool useWale;
+ bool useInitNeq;
bool simulatePorousMedia;
uint d3Qxx;
uint tEnd;
@@ -358,6 +362,7 @@ private:
bool isStreetVelocityFile;
bool isUseMeasurePoints;
bool isUseWale;
+ bool isUseInitNeq;
bool isSimulatePorousMedia;
bool isD3Qxx;
bool isTEnd;
diff --git a/src/Core/Input/ConfigFileReader/ConfigFileReader.cpp b/src/Core/Input/ConfigFileReader/ConfigFileReader.cpp
index d5863d508..75958b286 100644
--- a/src/Core/Input/ConfigFileReader/ConfigFileReader.cpp
+++ b/src/Core/Input/ConfigFileReader/ConfigFileReader.cpp
@@ -85,6 +85,9 @@ VF_PUBLIC std::shared_ptr<ConfigData> ConfigFileReader::readConfigFile(const std
if (input->getValue("UseWale") != "")
data->setUseWale(StringUtil::toBool(input->getValue("UseWale")));
+ if (input->getValue("UseInitNeq") != "")
+ data->setUseInitNeq(StringUtil::toBool(input->getValue("UseInitNeq")));
+
if (input->getValue("SimulatePorousMedia") != "")
data->setSimulatePorousMedia(StringUtil::toBool(input->getValue("SimulatePorousMedia")));
diff --git a/src/VirtualFluids_GPU/LBM/LB.h b/src/VirtualFluids_GPU/LBM/LB.h
index b4af44419..0aa989457 100644
--- a/src/VirtualFluids_GPU/LBM/LB.h
+++ b/src/VirtualFluids_GPU/LBM/LB.h
@@ -96,7 +96,7 @@ typedef struct InitCond{
real vis, vis_ratio;
real u0, u0_ratio;
real delta_rho, delta_press;
- bool printFiles, readGeo, doRestart, doCheckPoint, isGeo, isProp, isCp, calcMedian, GeometryValues, isConc, is2ndOrderMoments, is3rdOrderMoments, isHighOrderMoments, isWale, isMeasurePoints;
+ bool printFiles, readGeo, doRestart, doCheckPoint, isGeo, isProp, isCp, calcMedian, GeometryValues, isConc, is2ndOrderMoments, is3rdOrderMoments, isHighOrderMoments, isWale, isMeasurePoints, isInitNeq;
bool isGeoNormal, isInflowNormal, isOutflowNormal;
bool simulatePorousMedia;
bool streetVelocityFile;
diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/VirtualFluids_GPU/Parameter/Parameter.cpp
index b24ef75cf..dc25ced87 100644
--- a/src/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -116,6 +116,11 @@ Parameter::Parameter(SPtr<ConfigData> configData, Communicator* comm)
else
this->setUseWale(false);
//////////////////////////////////////////////////////////////////////////
+ if (configData->isUseInitNeqInConfigFile())
+ this->setUseInitNeq(configData->getUseInitNeq());
+ else
+ this->setUseInitNeq(false);
+ //////////////////////////////////////////////////////////////////////////
if (configData->isSimulatePorousMediaInConfigFile())
this->setSimulatePorousMedia(configData->getSimulatePorousMedia());
else
@@ -3617,6 +3622,10 @@ void Parameter::setUseWale(bool useWale)
{
ic.isWale = useWale;
}
+void Parameter::setUseInitNeq(bool useInitNeq)
+{
+ ic.isInitNeq = useInitNeq;
+}
void Parameter::setSimulatePorousMedia(bool simulatePorousMedia)
{
ic.simulatePorousMedia = simulatePorousMedia;
@@ -4942,6 +4951,10 @@ bool Parameter::getUseWale()
{
return ic.isWale;
}
+bool Parameter::getUseInitNeq()
+{
+ return ic.isInitNeq;
+}
bool Parameter::getSimulatePorousMedia()
{
return ic.simulatePorousMedia;
diff --git a/src/VirtualFluids_GPU/Parameter/Parameter.h b/src/VirtualFluids_GPU/Parameter/Parameter.h
index 6ac497489..8089d3fbb 100644
--- a/src/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/VirtualFluids_GPU/Parameter/Parameter.h
@@ -709,6 +709,7 @@ public:
void setStreetVelocityFile(bool streetVelocityFile);
void setUseMeasurePoints(bool useMeasurePoints);
void setUseWale(bool useWale);
+ void setUseInitNeq(bool useInitNeq);
void setSimulatePorousMedia(bool simulatePorousMedia);
void setclockCycleForMP(real clockCycleForMP);
void setDevices(std::vector<uint> devices);
@@ -952,6 +953,7 @@ public:
bool isStreetVelocityFile();
bool getUseMeasurePoints();
bool getUseWale();
+ bool getUseInitNeq();
bool getSimulatePorousMedia();
double getMemsizeGPU();
//1D domain decomposition
diff --git a/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27.cu b/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27.cu
index 69df08643..c02b4ffc8 100644
--- a/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27.cu
+++ b/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27.cu
@@ -28,18 +28,41 @@ void InitCompSP27::init(int level)
dim3 grid(Grid1, Grid2);
dim3 threads(numberOfThreads, 1, 1);
- LB_Init_Comp_SP_27 << < grid, threads >> >( para->getParD(level)->neighborX_SP,
- para->getParD(level)->neighborY_SP,
- para->getParD(level)->neighborZ_SP,
- para->getParD(level)->geoSP,
- para->getParD(level)->rho_SP,
- para->getParD(level)->vx_SP,
- para->getParD(level)->vy_SP,
- para->getParD(level)->vz_SP,
- para->getParD(level)->size_Mat_SP,
- para->getParD(level)->d0SP.f[0],
- para->getParD(level)->evenOrOdd);
- getLastCudaError("LBInitSP27 execution failed");
+ if( ! para->getUseInitNeq() )
+ {
+ LB_Init_Comp_SP_27 <<< grid, threads >>> (para->getParD(level)->neighborX_SP,
+ para->getParD(level)->neighborY_SP,
+ para->getParD(level)->neighborZ_SP,
+ para->getParD(level)->geoSP,
+ para->getParD(level)->rho_SP,
+ para->getParD(level)->vx_SP,
+ para->getParD(level)->vy_SP,
+ para->getParD(level)->vz_SP,
+ para->getParD(level)->size_Mat_SP,
+ para->getParD(level)->d0SP.f[0],
+ para->getParD(level)->evenOrOdd);
+ getLastCudaError("LBInitSP27 execution failed");
+ }
+ else
+ {
+ LB_Init_Comp_Neq_SP_27 <<< grid, threads >>> (para->getParD(level)->neighborX_SP,
+ para->getParD(level)->neighborY_SP,
+ para->getParD(level)->neighborZ_SP,
+ para->getParD(level)->neighborWSB_SP,
+ para->getParD(level)->geoSP,
+ para->getParD(level)->rho_SP,
+ para->getParD(level)->vx_SP,
+ para->getParD(level)->vy_SP,
+ para->getParD(level)->vz_SP,
+ para->getParD(level)->size_Mat_SP,
+ para->getParD(level)->d0SP.f[0],
+ para->getParD(level)->omega,
+ para->getParD(level)->evenOrOdd);
+ cudaDeviceSynchronize();
+ getLastCudaError("LBInitNeqSP27 execution failed");
+ }
+
+
}
diff --git a/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cu b/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cu
index 90ab9c5bd..6ecedde3d 100644
--- a/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cu
+++ b/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cu
@@ -3,6 +3,8 @@
#include "Core/RealConstants.h"
#include "math.h"
+#include <stdio.h>
+
extern "C" __global__ void LB_Init_Comp_SP_27(unsigned int* neighborX,
unsigned int* neighborY,
unsigned int* neighborZ,
@@ -162,4 +164,300 @@ extern "C" __global__ void LB_Init_Comp_SP_27(unsigned int* neighborX,
(D.f[dirTNW ])[ktnw ] = c1o216*(drho+ (c1o1+drho) * (c3o1*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq));
}
}
+}
+
+
+
+
+
+
+
+
+
+
+////////////////////////////////////////////////////////////////////////////////
+extern "C" __global__ void LB_Init_Comp_Neq_SP_27( unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int* neighborWSB,
+ unsigned int* geoD,
+ real* rho,
+ real* ux,
+ real* uy,
+ real* uz,
+ unsigned int size_Mat,
+ real* DD,
+ real omega,
+ 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<size_Mat)
+ {
+ ////////////////////////////////////////////////////////////////////////////////
+ unsigned int BC;
+ BC = geoD[k];
+
+ if( BC != GEO_SOLID && BC != GEO_VOID)
+ {
+ 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[dirBSW ] = &DD[dirTNE *size_Mat];
+ D.f[dirBNE ] = &DD[dirTSW *size_Mat];
+ D.f[dirBNW ] = &DD[dirTSE *size_Mat];
+ D.f[dirBSE ] = &DD[dirTNW *size_Mat];
+ D.f[dirTSW ] = &DD[dirBNE *size_Mat];
+ D.f[dirTNE ] = &DD[dirBSW *size_Mat];
+ D.f[dirTNW ] = &DD[dirBSE *size_Mat];
+ D.f[dirTSE ] = &DD[dirBNW *size_Mat];
+ }
+ //////////////////////////////////////////////////////////////////////////
+ real drho = rho[k];//0.0f;//
+ real vx1 = ux[k]; //0.0f;//
+ real vx2 = uy[k]; //0.0f;//
+ real vx3 = uz[k]; //0.0f;//
+ //////////////////////////////////////////////////////////////////////////
+ //index
+ //////////////////////////////////////////////////////////////////////////
+ unsigned int kzero= k;
+ unsigned int ke = k;
+ unsigned int kw = neighborX[k];
+ unsigned int kn = k;
+ unsigned int ks = neighborY[k];
+ unsigned int kt = k;
+ unsigned int kb = neighborZ[k];
+ unsigned int ksw = neighborY[kw];
+ unsigned int kne = k;
+ unsigned int kse = ks;
+ unsigned int knw = kw;
+ unsigned int kbw = neighborZ[kw];
+ unsigned int kte = k;
+ unsigned int kbe = kb;
+ unsigned int ktw = kw;
+ unsigned int kbs = neighborZ[ks];
+ unsigned int ktn = k;
+ 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 = k;
+ unsigned int kbsw = neighborZ[ksw];
+ //////////////////////////////////////////////////////////////////////////////
+ //neighbor index
+ uint kPx = neighborX[k];
+ uint kPy = neighborY[k];
+ uint kPz = neighborZ[k];
+ uint kMxyz = neighborWSB[k];
+ uint kMx = neighborZ[neighborY[kMxyz]];
+ uint kMy = neighborZ[neighborX[kMxyz]];
+ uint kMz = neighborY[neighborX[kMxyz]];
+ //////////////////////////////////////////////////////////////////////////
+ //getVeloX//
+ real vx1NeighborPx = ux[kPx];
+ real vx1NeighborMx = ux[kMx];
+ real vx1NeighborPy = ux[kPy];
+ real vx1NeighborMy = ux[kMy];
+ real vx1NeighborPz = ux[kPz];
+ real vx1NeighborMz = ux[kMz];
+ //getVeloY//
+ real vx2NeighborPx = uy[kPx];
+ real vx2NeighborMx = uy[kMx];
+ real vx2NeighborPy = uy[kPy];
+ real vx2NeighborMy = uy[kMy];
+ real vx2NeighborPz = uy[kPz];
+ real vx2NeighborMz = uy[kMz];
+ //getVeloZ//
+ real vx3NeighborPx = uz[kPx];
+ real vx3NeighborMx = uz[kMx];
+ real vx3NeighborPy = uz[kPy];
+ real vx3NeighborMy = uz[kMy];
+ real vx3NeighborPz = uz[kPz];
+ real vx3NeighborMz = uz[kMz];
+ //////////////////////////////////////////////////////////////////////////
+
+ real dvx1dx = (vx1NeighborPx - vx1NeighborMx) / c2o1;
+ real dvx1dy = (vx1NeighborPy - vx1NeighborMy) / c2o1;
+ real dvx1dz = (vx1NeighborPz - vx1NeighborMz) / c2o1;
+
+ real dvx2dx = (vx2NeighborPx - vx2NeighborMx) / c2o1;
+ real dvx2dy = (vx2NeighborPy - vx2NeighborMy) / c2o1;
+ real dvx2dz = (vx2NeighborPz - vx2NeighborMz) / c2o1;
+
+ real dvx3dx = (vx3NeighborPx - vx3NeighborMx) / c2o1;
+ real dvx3dy = (vx3NeighborPy - vx3NeighborMy) / c2o1;
+ real dvx3dz = (vx3NeighborPz - vx3NeighborMz) / c2o1;
+
+ //////////////////////////////////////////////////////////////////////////
+
+ // the following code is copy and pasted from VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
+ // i.e. Konstantins code
+
+ real ax = dvx1dx;
+ real ay = dvx1dy;
+ real az = dvx1dz;
+
+ real bx = dvx2dx;
+ real by = dvx2dy;
+ real bz = dvx2dz;
+
+ real cx = dvx3dx;
+ real cy = dvx3dy;
+ real cz = dvx3dz;
+
+ real eps_new = c1o1;
+ real op = c1o1;
+ real o = omega;
+
+ real f_E = eps_new *((5.*ax*o + 5.*by*o + 5.*cz*o - 8.*ax*op + 4.*by*op + 4.*cz*op)/(54.*o*op));
+
+ real f_N = f_E + eps_new *((2.*(ax - by))/(9.*o));
+ real f_T = f_E + eps_new *((2.*(ax - cz))/(9.*o));
+ real f_NE = eps_new *(-(5.*cz*o + 3.*(ay + bx)*op - 2.*cz*op + ax*(5.*o + op) + by*(5.*o + op))/(54.*o*op));
+ real f_SE = f_NE + eps_new *(( ay + bx )/(9.*o));
+ real f_TE = eps_new *(-(5.*cz*o + by*(5.*o - 2.*op) + 3.*(az + cx)*op + cz*op + ax*(5.*o + op))/(54.*o*op));
+ real f_BE = f_TE + eps_new *(( az + cx )/(9.*o));
+ real f_TN = eps_new *(-(5.*ax*o + 5.*by*o + 5.*cz*o - 2.*ax*op + by*op + 3.*bz*op + 3.*cy*op + cz*op)/(54.*o*op));
+ real f_BN = f_TN + eps_new *(( bz + cy )/(9.*o));
+ real f_ZERO = eps_new *((5.*(ax + by + cz))/(9.*op));
+ real f_TNE = eps_new *(-(ay + az + bx + bz + cx + cy)/(72.*o));
+ real f_TSW = - f_TNE - eps_new *((ay + bx)/(36.*o));
+ real f_TSE = - f_TNE - eps_new *((az + cx)/(36.*o));
+ real f_TNW = - f_TNE - eps_new *((bz + cy)/(36.*o));
+
+ //////////////////////////////////////////////////////////////////////////
+ real cu_sq=c3o2*(vx1*vx1+vx2*vx2+vx3*vx3);
+
+ (D.f[dirZERO])[kzero] = c8o27* (drho-cu_sq*(c1o1+drho));
+ (D.f[dirE ])[ke ] = c2o27* (drho+ (c1o1+drho) * (c3o1*( vx1 )+c9o2*( vx1 )*( vx1 )-cu_sq));
+ (D.f[dirW ])[kw ] = c2o27* (drho+ (c1o1+drho) * (c3o1*(-vx1 )+c9o2*(-vx1 )*(-vx1 )-cu_sq));
+ (D.f[dirN ])[kn ] = c2o27* (drho+ (c1o1+drho) * (c3o1*( vx2 )+c9o2*( vx2 )*( vx2 )-cu_sq));
+ (D.f[dirS ])[ks ] = c2o27* (drho+ (c1o1+drho) * (c3o1*( -vx2 )+c9o2*( -vx2 )*( -vx2 )-cu_sq));
+ (D.f[dirT ])[kt ] = c2o27* (drho+ (c1o1+drho) * (c3o1*( vx3)+c9o2*( vx3)*( vx3)-cu_sq));
+ (D.f[dirB ])[kb ] = c2o27* (drho+ (c1o1+drho) * (c3o1*( -vx3)+c9o2*( -vx3)*( -vx3)-cu_sq));
+ (D.f[dirNE ])[kne ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( vx1+vx2 )+c9o2*( vx1+vx2 )*( vx1+vx2 )-cu_sq));
+ (D.f[dirSW ])[ksw ] = c1o54* (drho+ (c1o1+drho) * (c3o1*(-vx1-vx2 )+c9o2*(-vx1-vx2 )*(-vx1-vx2 )-cu_sq));
+ (D.f[dirSE ])[kse ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( vx1-vx2 )+c9o2*( vx1-vx2 )*( vx1-vx2 )-cu_sq));
+ (D.f[dirNW ])[knw ] = c1o54* (drho+ (c1o1+drho) * (c3o1*(-vx1+vx2 )+c9o2*(-vx1+vx2 )*(-vx1+vx2 )-cu_sq));
+ (D.f[dirTE ])[kte ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( vx1 +vx3)+c9o2*( vx1 +vx3)*( vx1 +vx3)-cu_sq));
+ (D.f[dirBW ])[kbw ] = c1o54* (drho+ (c1o1+drho) * (c3o1*(-vx1 -vx3)+c9o2*(-vx1 -vx3)*(-vx1 -vx3)-cu_sq));
+ (D.f[dirBE ])[kbe ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( vx1 -vx3)+c9o2*( vx1 -vx3)*( vx1 -vx3)-cu_sq));
+ (D.f[dirTW ])[ktw ] = c1o54* (drho+ (c1o1+drho) * (c3o1*(-vx1 +vx3)+c9o2*(-vx1 +vx3)*(-vx1 +vx3)-cu_sq));
+ (D.f[dirTN ])[ktn ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( vx2+vx3)+c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq));
+ (D.f[dirBS ])[kbs ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( -vx2-vx3)+c9o2*( -vx2-vx3)*( -vx2-vx3)-cu_sq));
+ (D.f[dirBN ])[kbn ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( vx2-vx3)+c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq));
+ (D.f[dirTS ])[kts ] = c1o54* (drho+ (c1o1+drho) * (c3o1*( -vx2+vx3)+c9o2*( -vx2+vx3)*( -vx2+vx3)-cu_sq));
+ (D.f[dirTNE ])[ktne ] = c1o216*(drho+ (c1o1+drho) * (c3o1*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq));
+ (D.f[dirBSW ])[kbsw ] = c1o216*(drho+ (c1o1+drho) * (c3o1*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq));
+ (D.f[dirBNE ])[kbne ] = c1o216*(drho+ (c1o1+drho) * (c3o1*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq));
+ (D.f[dirTSW ])[ktsw ] = c1o216*(drho+ (c1o1+drho) * (c3o1*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq));
+ (D.f[dirTSE ])[ktse ] = c1o216*(drho+ (c1o1+drho) * (c3o1*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq));
+ (D.f[dirBNW ])[kbnw ] = c1o216*(drho+ (c1o1+drho) * (c3o1*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq));
+ (D.f[dirBSE ])[kbse ] = c1o216*(drho+ (c1o1+drho) * (c3o1*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq));
+ (D.f[dirTNW ])[ktnw ] = c1o216*(drho+ (c1o1+drho) * (c3o1*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq));
+
+ //////////////////////////////////////////////////////////////////////////
+
+ (D.f[dirZERO])[kzero] += (c1o1+drho) * f_ZERO;
+ (D.f[dirE ])[ke ] += (c1o1+drho) * f_E ;
+ (D.f[dirW ])[kw ] += (c1o1+drho) * f_E ;
+ (D.f[dirN ])[kn ] += (c1o1+drho) * f_N ;
+ (D.f[dirS ])[ks ] += (c1o1+drho) * f_N ;
+ (D.f[dirT ])[kt ] += (c1o1+drho) * f_T ;
+ (D.f[dirB ])[kb ] += (c1o1+drho) * f_T ;
+ (D.f[dirNE ])[kne ] += (c1o1+drho) * f_NE ;
+ (D.f[dirSW ])[ksw ] += (c1o1+drho) * f_NE ;
+ (D.f[dirSE ])[kse ] += (c1o1+drho) * f_SE ;
+ (D.f[dirNW ])[knw ] += (c1o1+drho) * f_SE ;
+ (D.f[dirTE ])[kte ] += (c1o1+drho) * f_TE ;
+ (D.f[dirBW ])[kbw ] += (c1o1+drho) * f_TE ;
+ (D.f[dirBE ])[kbe ] += (c1o1+drho) * f_BE ;
+ (D.f[dirTW ])[ktw ] += (c1o1+drho) * f_BE ;
+ (D.f[dirTN ])[ktn ] += (c1o1+drho) * f_TN ;
+ (D.f[dirBS ])[kbs ] += (c1o1+drho) * f_TN ;
+ (D.f[dirBN ])[kbn ] += (c1o1+drho) * f_BN ;
+ (D.f[dirTS ])[kts ] += (c1o1+drho) * f_BN ;
+ (D.f[dirTNE ])[ktne ] += (c1o1+drho) * f_TNE ;
+ (D.f[dirBSW ])[kbsw ] += (c1o1+drho) * f_TNE ;
+ (D.f[dirBNE ])[kbne ] += (c1o1+drho) * f_TSW ;
+ (D.f[dirTSW ])[ktsw ] += (c1o1+drho) * f_TSW ;
+ (D.f[dirTSE ])[ktse ] += (c1o1+drho) * f_TSE ;
+ (D.f[dirBNW ])[kbnw ] += (c1o1+drho) * f_TSE ;
+ (D.f[dirBSE ])[kbse ] += (c1o1+drho) * f_TNW ;
+ (D.f[dirTNW ])[ktnw ] += (c1o1+drho) * f_TNW ;
+
+ //////////////////////////////////////////////////////////////////////////
+ }
+ else
+ {
+ //////////////////////////////////////////////////////////////////////////
+ Distributions27 D;
+ D.f[dirZERO] = &DD[dirZERO*size_Mat];
+ //////////////////////////////////////////////////////////////////////////
+ (D.f[dirZERO])[k] = c96o1;
+ //////////////////////////////////////////////////////////////////////////
+ }
+ }
}
\ No newline at end of file
diff --git a/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cuh b/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cuh
index c6c784186..dd9dbd7d0 100644
--- a/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cuh
+++ b/src/VirtualFluids_GPU/PreProcessor/PreProcessorStrategy/InitCompSP27/InitCompSP27_Device.cuh
@@ -16,4 +16,18 @@ extern "C" __global__ void LB_Init_Comp_SP_27(unsigned int* neighborX,
real* DD,
bool EvenOrOdd);
+extern "C" __global__ void LB_Init_Comp_Neq_SP_27(unsigned int* neighborX,
+ unsigned int* neighborY,
+ unsigned int* neighborZ,
+ unsigned int* neighborWSB,
+ unsigned int* geoD,
+ real* rho,
+ real* ux,
+ real* uy,
+ real* uz,
+ unsigned int size_Mat,
+ real* DD,
+ real omega,
+ bool EvenOrOdd);
+
#endif
\ No newline at end of file
diff --git a/targets/apps/LBM/TGV_3D/TGV_3D.cpp b/targets/apps/LBM/TGV_3D/TGV_3D.cpp
index fac2cc777..4b138f72d 100644
--- a/targets/apps/LBM/TGV_3D/TGV_3D.cpp
+++ b/targets/apps/LBM/TGV_3D/TGV_3D.cpp
@@ -245,6 +245,8 @@ void multipleLevel(const std::string& configPath)
if( useWale )
para->setUseWale( true );
+ //para->setUseInitNeq( true );
+
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
--
GitLab