diff --git a/.clang-format b/.clang-format
index 1067e3c24a255eab69b5149f5d40d4b9ce3a596b..cca4fe23c1d1eb12fe71b5d097d7e5ec4146839b 100644
--- a/.clang-format
+++ b/.clang-format
@@ -16,7 +16,7 @@ AllowShortBlocksOnASingleLine: Never
AllowShortCaseLabelsOnASingleLine: false
AllowShortFunctionsOnASingleLine: None
AllowShortLambdasOnASingleLine: All
-AllowShortIfStatementsOnASingleLine: WithoutElse
+AllowShortIfStatementsOnASingleLine: false
AllowShortLoopsOnASingleLine: false
AlwaysBreakAfterDefinitionReturnType: None
AlwaysBreakAfterReturnType: None
@@ -99,7 +99,7 @@ PenaltyBreakString: 1000
PenaltyBreakTemplateDeclaration: 10
PenaltyExcessCharacter: 1000000
PenaltyReturnTypeOnItsOwnLine: 60
-PointerAlignment: Right
+PointerAlignment: Left
ReflowComments: true
SortIncludes: true
SortUsingDeclarations: true
@@ -129,5 +129,6 @@ StatementMacros:
TabWidth: 4
UseCRLF: false
UseTab: Never
+SpaceBeforeCpp11BracedList: true
...
diff --git a/CMakeLists.txt b/CMakeLists.txt
index 34edea623975db57b6b3a23b0a1b49721367d67b..ec5f5dbedbdb2f6b7d38fc39b6b864f3e33559e2 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -1,5 +1,5 @@
#################################################################################
-# _ ___ __ __________ _ __
+# _ ___ __ __________ _ __
# | | / (_)____/ /___ ______ _/ / ____/ /_ __(_)___/ /____
# | | / / / ___/ __/ / / / __ `/ / /_ / / / / / / __ / ___/
# | |/ / / / / /_/ /_/ / /_/ / / __/ / / /_/ / / /_/ (__ )
diff --git a/apps/cpu/FlowAroundCylinder/cylinder.cpp b/apps/cpu/FlowAroundCylinder/cylinder.cpp
index 3e5be5e080a5702e9608bf037d492b8cbc809dfc..0599060816ff806e087680d892be5232b1558930 100644
--- a/apps/cpu/FlowAroundCylinder/cylinder.cpp
+++ b/apps/cpu/FlowAroundCylinder/cylinder.cpp
@@ -1,6 +1,7 @@
#include <iostream>
#include <string>
+#include "K17CompressibleNavierStokes.h"
#include "VirtualFluids.h"
using namespace std;
@@ -245,7 +246,7 @@ void run(string configname)
UBLOG(logINFO, "Available memory per process = "<<availMem<<" bytes");
}
- SPtr<LBMKernel> kernel(new CompressibleCumulantLBMKernel());
+ SPtr<LBMKernel> kernel(new K17CompressibleNavierStokes());
SPtr<BCSet> bcProc(new BCSet());
kernel->setBCSet(bcProc);
diff --git a/pythonbindings/src/gpu/submodules/turbulence_models.cpp b/pythonbindings/src/gpu/submodules/turbulence_models.cpp
index cfbb9e56127fee0cd90a482dde258d8b96389989..79b8ac00d66e2c32abf2e2b708404f7c265de013 100644
--- a/pythonbindings/src/gpu/submodules/turbulence_models.cpp
+++ b/pythonbindings/src/gpu/submodules/turbulence_models.cpp
@@ -37,6 +37,7 @@
namespace turbulence_model
{
namespace py = pybind11;
+ using namespace vf::lbm;
void makeModule(py::module_ &parentModule)
{
diff --git a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
index 11119ecbd59f1a4d145fb130b3d0fe2b2a671fac..ae213baab78c7c18912c2b33af684aad47e3b1dd 100644
--- a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
+++ b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
@@ -56,70 +56,72 @@ void D3Q27EsoTwist3DSplittedVector::swap() { std::swap(this->localDistributions,
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::getDistribution(real *const f, size_t x1, size_t x2, size_t x3)
{
- using namespace vf::lbm::dir;
+ const size_t x1p = x1 + 1;
+ const size_t x2p = x2 + 1;
+ const size_t x3p = x3 + 1;
- f[DIR_P00] = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
- f[DIR_0P0] = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
- f[DIR_00P] = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
- f[DIR_PP0] = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
- f[DIR_MP0] = (*this->localDistributions)(D3Q27System::ET_NW, x1 + 1, x2, x3);
- f[DIR_P0P] = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
- f[DIR_M0P] = (*this->localDistributions)(D3Q27System::ET_TW, x1 + 1, x2, x3);
- f[DIR_0PP] = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3);
- f[DIR_0MP] = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2 + 1, x3);
- f[DIR_PPP] = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3);
- f[DIR_MPP] = (*this->localDistributions)(D3Q27System::ET_TNW, x1 + 1, x2, x3);
- f[DIR_PMP] = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2 + 1, x3);
- f[DIR_MMP] = (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3);
+ f[vf::lbm::dir::DIR_P00] = (*this->localDistributions)(D3Q27System::ET_P00, x1, x2, x3);
+ f[vf::lbm::dir::DIR_0P0] = (*this->localDistributions)(D3Q27System::ET_0P0, x1, x2, x3);
+ f[vf::lbm::dir::DIR_00P] = (*this->localDistributions)(D3Q27System::ET_00P, x1, x2, x3);
+ f[vf::lbm::dir::DIR_PP0] = (*this->localDistributions)(D3Q27System::ET_PP0, x1, x2, x3);
+ f[vf::lbm::dir::DIR_MP0] = (*this->localDistributions)(D3Q27System::ET_MP0, x1p, x2, x3);
+ f[vf::lbm::dir::DIR_P0P] = (*this->localDistributions)(D3Q27System::ET_P0P, x1, x2, x3);
+ f[vf::lbm::dir::DIR_M0P] = (*this->localDistributions)(D3Q27System::ET_M0P, x1p, x2, x3);
+ f[vf::lbm::dir::DIR_0PP] = (*this->localDistributions)(D3Q27System::ET_0PP, x1, x2, x3);
+ f[vf::lbm::dir::DIR_0MP] = (*this->localDistributions)(D3Q27System::ET_0MP, x1, x2p, x3);
+ f[vf::lbm::dir::DIR_PPP] = (*this->localDistributions)(D3Q27System::ET_PPP, x1, x2, x3);
+ f[vf::lbm::dir::DIR_MPP] = (*this->localDistributions)(D3Q27System::ET_MPP, x1p, x2, x3);
+ f[vf::lbm::dir::DIR_PMP] = (*this->localDistributions)(D3Q27System::ET_PMP, x1, x2p, x3);
+ f[vf::lbm::dir::DIR_MMP] = (*this->localDistributions)(D3Q27System::ET_MMP, x1p, x2p, x3);
- f[DIR_M00] = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1 + 1, x2, x3);
- f[DIR_0M0] = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2 + 1, x3);
- f[DIR_00M] = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3 + 1);
- f[DIR_MM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1 + 1, x2 + 1, x3);
- f[DIR_PM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2 + 1, x3);
- f[DIR_M0M] = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1 + 1, x2, x3 + 1);
- f[DIR_P0M] = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3 + 1);
- f[DIR_0MM] = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2 + 1, x3 + 1);
- f[DIR_0PM] = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3 + 1);
- f[DIR_MMM] = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1 + 1, x2 + 1, x3 + 1);
- f[DIR_PMM] = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2 + 1, x3 + 1);
- f[DIR_MPM] = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1);
- f[DIR_PPM] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
+ f[vf::lbm::dir::DIR_M00] = (*this->nonLocalDistributions)(D3Q27System::ET_M00, x1p, x2, x3);
+ f[vf::lbm::dir::DIR_0M0] = (*this->nonLocalDistributions)(D3Q27System::ET_0M0, x1, x2p, x3);
+ f[vf::lbm::dir::DIR_00M] = (*this->nonLocalDistributions)(D3Q27System::ET_00M, x1, x2, x3p);
+ f[vf::lbm::dir::DIR_MM0] = (*this->nonLocalDistributions)(D3Q27System::ET_MM0, x1p, x2p, x3);
+ f[vf::lbm::dir::DIR_PM0] = (*this->nonLocalDistributions)(D3Q27System::ET_PM0, x1, x2p, x3);
+ f[vf::lbm::dir::DIR_M0M] = (*this->nonLocalDistributions)(D3Q27System::ET_M0M, x1p, x2, x3p);
+ f[vf::lbm::dir::DIR_P0M] = (*this->nonLocalDistributions)(D3Q27System::ET_P0M, x1, x2, x3p);
+ f[vf::lbm::dir::DIR_0MM] = (*this->nonLocalDistributions)(D3Q27System::ET_0MM, x1, x2p, x3p);
+ f[vf::lbm::dir::DIR_0PM] = (*this->nonLocalDistributions)(D3Q27System::ET_0PM, x1, x2, x3p);
+ f[vf::lbm::dir::DIR_MMM] = (*this->nonLocalDistributions)(D3Q27System::ET_MMM, x1p, x2p, x3p);
+ f[vf::lbm::dir::DIR_PMM] = (*this->nonLocalDistributions)(D3Q27System::ET_PMM, x1, x2p, x3p);
+ f[vf::lbm::dir::DIR_MPM] = (*this->nonLocalDistributions)(D3Q27System::ET_MPM, x1p, x2, x3p);
+ f[vf::lbm::dir::DIR_PPM] = (*this->nonLocalDistributions)(D3Q27System::ET_PPM, x1, x2, x3p);
- f[DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
+ f[vf::lbm::dir::DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setDistribution(const real *const f, size_t x1, size_t x2, size_t x3)
{
using namespace vf::lbm::dir;
- (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = f[INV_P00];
- (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = f[INV_0P0];
- (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = f[INV_00P];
- (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = f[INV_PP0];
- (*this->localDistributions)(D3Q27System::ET_NW, x1 + 1, x2, x3) = f[INV_MP0];
- (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = f[INV_P0P];
- (*this->localDistributions)(D3Q27System::ET_TW, x1 + 1, x2, x3) = f[INV_M0P];
- (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = f[INV_0PP];
- (*this->localDistributions)(D3Q27System::ET_TS, x1, x2 + 1, x3) = f[INV_0MP];
- (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = f[INV_PPP];
- (*this->localDistributions)(D3Q27System::ET_TNW, x1 + 1, x2, x3) = f[INV_MPP];
- (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2 + 1, x3) = f[INV_PMP];
- (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3) = f[INV_MMP];
+ (*this->localDistributions)(D3Q27System::ET_P00, x1, x2, x3) = f[INV_P00];
+ (*this->localDistributions)(D3Q27System::ET_0P0, x1, x2, x3) = f[INV_0P0];
+ (*this->localDistributions)(D3Q27System::ET_00P, x1, x2, x3) = f[INV_00P];
+ (*this->localDistributions)(D3Q27System::ET_PP0, x1, x2, x3) = f[INV_PP0];
+ (*this->localDistributions)(D3Q27System::ET_MP0, x1 + 1, x2, x3) = f[INV_MP0];
+ (*this->localDistributions)(D3Q27System::ET_P0P, x1, x2, x3) = f[INV_P0P];
+ (*this->localDistributions)(D3Q27System::ET_M0P, x1 + 1, x2, x3) = f[INV_M0P];
+ (*this->localDistributions)(D3Q27System::ET_0PP, x1, x2, x3) = f[INV_0PP];
+ (*this->localDistributions)(D3Q27System::ET_0MP, x1, x2 + 1, x3) = f[INV_0MP];
+ (*this->localDistributions)(D3Q27System::ET_PPP, x1, x2, x3) = f[INV_PPP];
+ (*this->localDistributions)(D3Q27System::ET_MPP, x1 + 1, x2, x3) = f[INV_MPP];
+ (*this->localDistributions)(D3Q27System::ET_PMP, x1, x2 + 1, x3) = f[INV_PMP];
+ (*this->localDistributions)(D3Q27System::ET_MMP, x1 + 1, x2 + 1, x3) = f[INV_MMP];
- (*this->nonLocalDistributions)(D3Q27System::ET_W, x1 + 1, x2, x3) = f[INV_M00];
- (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2 + 1, x3) = f[INV_0M0];
- (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3 + 1) = f[INV_00M];
- (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1 + 1, x2 + 1, x3) = f[INV_MM0];
- (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2 + 1, x3) = f[INV_PM0];
- (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1 + 1, x2, x3 + 1) = f[INV_M0M];
- (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3 + 1) = f[INV_P0M];
- (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2 + 1, x3 + 1) = f[INV_0MM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3 + 1) = f[INV_0PM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1 + 1, x2 + 1, x3 + 1) = f[INV_MMM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2 + 1, x3 + 1) = f[INV_PMM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1) = f[INV_MPM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f[INV_PPM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_M00, x1 + 1, x2, x3) = f[INV_M00];
+ (*this->nonLocalDistributions)(D3Q27System::ET_0M0, x1, x2 + 1, x3) = f[INV_0M0];
+ (*this->nonLocalDistributions)(D3Q27System::ET_00M, x1, x2, x3 + 1) = f[INV_00M];
+ (*this->nonLocalDistributions)(D3Q27System::ET_MM0, x1 + 1, x2 + 1, x3) = f[INV_MM0];
+ (*this->nonLocalDistributions)(D3Q27System::ET_PM0, x1, x2 + 1, x3) = f[INV_PM0];
+ (*this->nonLocalDistributions)(D3Q27System::ET_M0M, x1 + 1, x2, x3 + 1) = f[INV_M0M];
+ (*this->nonLocalDistributions)(D3Q27System::ET_P0M, x1, x2, x3 + 1) = f[INV_P0M];
+ (*this->nonLocalDistributions)(D3Q27System::ET_0MM, x1, x2 + 1, x3 + 1) = f[INV_0MM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_0PM, x1, x2, x3 + 1) = f[INV_0PM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_MMM, x1 + 1, x2 + 1, x3 + 1) = f[INV_MMM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_PMM, x1, x2 + 1, x3 + 1) = f[INV_PMM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_MPM, x1 + 1, x2, x3 + 1) = f[INV_MPM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_PPM, x1, x2, x3 + 1) = f[INV_PPM];
(*this->zeroDistributions)(x1, x2, x3) = f[DIR_000];
}
@@ -128,33 +130,33 @@ void D3Q27EsoTwist3DSplittedVector::getDistributionInv(real *const f, size_t x1,
{
using namespace vf::lbm::dir;
- f[INV_P00] = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
- f[INV_0P0] = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
- f[INV_00P] = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
- f[INV_PP0] = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
- f[INV_MP0] = (*this->localDistributions)(D3Q27System::ET_NW, x1 + 1, x2, x3);
- f[INV_P0P] = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
- f[INV_M0P] = (*this->localDistributions)(D3Q27System::ET_TW, x1 + 1, x2, x3);
- f[INV_0PP] = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3);
- f[INV_0MP] = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2 + 1, x3);
- f[INV_PPP] = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3);
- f[INV_MPP] = (*this->localDistributions)(D3Q27System::ET_TNW, x1 + 1, x2, x3);
- f[INV_PMP] = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2 + 1, x3);
- f[INV_MMP] = (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3);
+ f[INV_P00] = (*this->localDistributions)(D3Q27System::ET_P00, x1, x2, x3);
+ f[INV_0P0] = (*this->localDistributions)(D3Q27System::ET_0P0, x1, x2, x3);
+ f[INV_00P] = (*this->localDistributions)(D3Q27System::ET_00P, x1, x2, x3);
+ f[INV_PP0] = (*this->localDistributions)(D3Q27System::ET_PP0, x1, x2, x3);
+ f[INV_MP0] = (*this->localDistributions)(D3Q27System::ET_MP0, x1 + 1, x2, x3);
+ f[INV_P0P] = (*this->localDistributions)(D3Q27System::ET_P0P, x1, x2, x3);
+ f[INV_M0P] = (*this->localDistributions)(D3Q27System::ET_M0P, x1 + 1, x2, x3);
+ f[INV_0PP] = (*this->localDistributions)(D3Q27System::ET_0PP, x1, x2, x3);
+ f[INV_0MP] = (*this->localDistributions)(D3Q27System::ET_0MP, x1, x2 + 1, x3);
+ f[INV_PPP] = (*this->localDistributions)(D3Q27System::ET_PPP, x1, x2, x3);
+ f[INV_MPP] = (*this->localDistributions)(D3Q27System::ET_MPP, x1 + 1, x2, x3);
+ f[INV_PMP] = (*this->localDistributions)(D3Q27System::ET_PMP, x1, x2 + 1, x3);
+ f[INV_MMP] = (*this->localDistributions)(D3Q27System::ET_MMP, x1 + 1, x2 + 1, x3);
- f[INV_M00] = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1 + 1, x2, x3);
- f[INV_0M0] = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2 + 1, x3);
- f[INV_00M] = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3 + 1);
- f[INV_MM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1 + 1, x2 + 1, x3);
- f[INV_PM0] = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2 + 1, x3);
- f[INV_M0M] = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1 + 1, x2, x3 + 1);
- f[INV_P0M] = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3 + 1);
- f[INV_0MM] = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2 + 1, x3 + 1);
- f[INV_0PM] = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3 + 1);
- f[INV_MMM] = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1 + 1, x2 + 1, x3 + 1);
- f[INV_PMM] = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2 + 1, x3 + 1);
- f[INV_MPM] = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1);
- f[INV_PPM] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
+ f[INV_M00] = (*this->nonLocalDistributions)(D3Q27System::ET_M00, x1 + 1, x2, x3);
+ f[INV_0M0] = (*this->nonLocalDistributions)(D3Q27System::ET_0M0, x1, x2 + 1, x3);
+ f[INV_00M] = (*this->nonLocalDistributions)(D3Q27System::ET_00M, x1, x2, x3 + 1);
+ f[INV_MM0] = (*this->nonLocalDistributions)(D3Q27System::ET_MM0, x1 + 1, x2 + 1, x3);
+ f[INV_PM0] = (*this->nonLocalDistributions)(D3Q27System::ET_PM0, x1, x2 + 1, x3);
+ f[INV_M0M] = (*this->nonLocalDistributions)(D3Q27System::ET_M0M, x1 + 1, x2, x3 + 1);
+ f[INV_P0M] = (*this->nonLocalDistributions)(D3Q27System::ET_P0M, x1, x2, x3 + 1);
+ f[INV_0MM] = (*this->nonLocalDistributions)(D3Q27System::ET_0MM, x1, x2 + 1, x3 + 1);
+ f[INV_0PM] = (*this->nonLocalDistributions)(D3Q27System::ET_0PM, x1, x2, x3 + 1);
+ f[INV_MMM] = (*this->nonLocalDistributions)(D3Q27System::ET_MMM, x1 + 1, x2 + 1, x3 + 1);
+ f[INV_PMM] = (*this->nonLocalDistributions)(D3Q27System::ET_PMM, x1, x2 + 1, x3 + 1);
+ f[INV_MPM] = (*this->nonLocalDistributions)(D3Q27System::ET_MPM, x1 + 1, x2, x3 + 1);
+ f[INV_PPM] = (*this->nonLocalDistributions)(D3Q27System::ET_PPM, x1, x2, x3 + 1);
f[DIR_000] = (*this->zeroDistributions)(x1, x2, x3);
}
@@ -163,33 +165,33 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionInv(const real *const f, size
{
using namespace vf::lbm::dir;
- (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = f[DIR_P00];
- (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = f[DIR_0P0];
- (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = f[DIR_00P];
- (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = f[DIR_PP0];
- (*this->localDistributions)(D3Q27System::ET_NW, x1 + 1, x2, x3) = f[DIR_MP0];
- (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = f[DIR_P0P];
- (*this->localDistributions)(D3Q27System::ET_TW, x1 + 1, x2, x3) = f[DIR_M0P];
- (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = f[DIR_0PP];
- (*this->localDistributions)(D3Q27System::ET_TS, x1, x2 + 1, x3) = f[DIR_0MP];
- (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = f[DIR_PPP];
- (*this->localDistributions)(D3Q27System::ET_TNW, x1 + 1, x2, x3) = f[DIR_MPP];
- (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2 + 1, x3) = f[DIR_PMP];
- (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3) = f[DIR_MMP];
+ (*this->localDistributions)(D3Q27System::ET_P00, x1, x2, x3) = f[DIR_P00];
+ (*this->localDistributions)(D3Q27System::ET_0P0, x1, x2, x3) = f[DIR_0P0];
+ (*this->localDistributions)(D3Q27System::ET_00P, x1, x2, x3) = f[DIR_00P];
+ (*this->localDistributions)(D3Q27System::ET_PP0, x1, x2, x3) = f[DIR_PP0];
+ (*this->localDistributions)(D3Q27System::ET_MP0, x1 + 1, x2, x3) = f[DIR_MP0];
+ (*this->localDistributions)(D3Q27System::ET_P0P, x1, x2, x3) = f[DIR_P0P];
+ (*this->localDistributions)(D3Q27System::ET_M0P, x1 + 1, x2, x3) = f[DIR_M0P];
+ (*this->localDistributions)(D3Q27System::ET_0PP, x1, x2, x3) = f[DIR_0PP];
+ (*this->localDistributions)(D3Q27System::ET_0MP, x1, x2 + 1, x3) = f[DIR_0MP];
+ (*this->localDistributions)(D3Q27System::ET_PPP, x1, x2, x3) = f[DIR_PPP];
+ (*this->localDistributions)(D3Q27System::ET_MPP, x1 + 1, x2, x3) = f[DIR_MPP];
+ (*this->localDistributions)(D3Q27System::ET_PMP, x1, x2 + 1, x3) = f[DIR_PMP];
+ (*this->localDistributions)(D3Q27System::ET_MMP, x1 + 1, x2 + 1, x3) = f[DIR_MMP];
- (*this->nonLocalDistributions)(D3Q27System::ET_W, x1 + 1, x2, x3) = f[DIR_M00];
- (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2 + 1, x3) = f[DIR_0M0];
- (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3 + 1) = f[DIR_00M];
- (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1 + 1, x2 + 1, x3) = f[DIR_MM0];
- (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2 + 1, x3) = f[DIR_PM0];
- (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1 + 1, x2, x3 + 1) = f[DIR_M0M];
- (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3 + 1) = f[DIR_P0M];
- (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2 + 1, x3 + 1) = f[DIR_0MM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3 + 1) = f[DIR_0PM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1 + 1, x2 + 1, x3 + 1) = f[DIR_MMM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2 + 1, x3 + 1) = f[DIR_PMM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1) = f[DIR_MPM];
- (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f[DIR_PPM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_M00, x1 + 1, x2, x3) = f[DIR_M00];
+ (*this->nonLocalDistributions)(D3Q27System::ET_0M0, x1, x2 + 1, x3) = f[DIR_0M0];
+ (*this->nonLocalDistributions)(D3Q27System::ET_00M, x1, x2, x3 + 1) = f[DIR_00M];
+ (*this->nonLocalDistributions)(D3Q27System::ET_MM0, x1 + 1, x2 + 1, x3) = f[DIR_MM0];
+ (*this->nonLocalDistributions)(D3Q27System::ET_PM0, x1, x2 + 1, x3) = f[DIR_PM0];
+ (*this->nonLocalDistributions)(D3Q27System::ET_M0M, x1 + 1, x2, x3 + 1) = f[DIR_M0M];
+ (*this->nonLocalDistributions)(D3Q27System::ET_P0M, x1, x2, x3 + 1) = f[DIR_P0M];
+ (*this->nonLocalDistributions)(D3Q27System::ET_0MM, x1, x2 + 1, x3 + 1) = f[DIR_0MM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_0PM, x1, x2, x3 + 1) = f[DIR_0PM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_MMM, x1 + 1, x2 + 1, x3 + 1) = f[DIR_MMM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_PMM, x1, x2 + 1, x3 + 1) = f[DIR_PMM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_MPM, x1 + 1, x2, x3 + 1) = f[DIR_MPM];
+ (*this->nonLocalDistributions)(D3Q27System::ET_PPM, x1, x2, x3 + 1) = f[DIR_PPM];
(*this->zeroDistributions)(x1, x2, x3) = f[DIR_000];
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernelUnified.cpp b/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernelUnified.cpp
deleted file mode 100644
index 2a791636d42e8f17761873cd5e95c8dcf933e0b6..0000000000000000000000000000000000000000
--- a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernelUnified.cpp
+++ /dev/null
@@ -1,338 +0,0 @@
-//=======================================================================================
-// ____ ____ __ ______ __________ __ __ __ __
-// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
-// \ \ | | | | | |_) | | | | | | | / \ | |
-// \ \ | | | | | _ / | | | | | | / /\ \ | |
-// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
-// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
-// \ \ | | ________________________________________________________________
-// \ \ | | | ______________________________________________________________|
-// \ \| | | | __ __ __ __ ______ _______
-// \ | | |_____ | | | | | | | | | _ \ / _____)
-// \ | | _____| | | | | | | | | | | \ \ \_______
-// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
-// \ _____| |__| |________| \_______/ |__| |______/ (_______/
-//
-// This file is part of VirtualFluids. VirtualFluids is free software: you can
-// redistribute it and/or modify it under the terms of the GNU General Public
-// License as published by the Free Software Foundation, either version 3 of
-// the License, or (at your option) any later version.
-//
-// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
-// for more details.
-//
-// You should have received a copy of the GNU General Public License along
-// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
-//
-//! \file CumulantK17LBMKernel.cpp
-//! \ingroup LBM
-//! \author Konstantin Kutscher, Martin Geier
-//=======================================================================================
-#include <lbm/KernelParameter.h>
-#include <lbm/CumulantChimera.h>
-#include <lbm/constants/D3Q27.h>
-
-#include "CumulantK17LBMKernelUnified.h"
-#include "D3Q27System.h"
-#include "D3Q27EsoTwist3DSplittedVector.h"
-#include <cmath>
-#include "DataSet3D.h"
-#include "LBMKernel.h"
-#include "Block3D.h"
-#include "BCArray3D.h"
-
-
-//#define PROOF_CORRECTNESS
-
-//using namespace UbMath;
-using namespace vf::basics::constant;
-
-//////////////////////////////////////////////////////////////////////////
-CumulantK17LBMKernelUnified::CumulantK17LBMKernelUnified()
-{
- this->compressible = true;
-}
-//////////////////////////////////////////////////////////////////////////
-void CumulantK17LBMKernelUnified::initDataSet()
-{
- SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0] + 2, nx[1] + 2, nx[2] + 2, -999.9));
- dataSet->setFdistributions(d);
-}
-//////////////////////////////////////////////////////////////////////////
-SPtr<LBMKernel> CumulantK17LBMKernelUnified::clone()
-{
- SPtr<LBMKernel> kernel(new CumulantK17LBMKernelUnified());
- kernel->setNX(nx);
- std::dynamic_pointer_cast<CumulantK17LBMKernelUnified>(kernel)->initDataSet();
- kernel->setCollisionFactor(this->collFactor);
- kernel->setBCSet(bcSet->clone(kernel));
- kernel->setWithForcing(withForcing);
- kernel->setForcingX1(muForcingX1);
- kernel->setForcingX2(muForcingX2);
- kernel->setForcingX3(muForcingX3);
- kernel->setIndex(ix1, ix2, ix3);
- kernel->setDeltaT(deltaT);
- kernel->setBlock(block.lock());
-
- return kernel;
-}
-//////////////////////////////////////////////////////////////////////////
-void CumulantK17LBMKernelUnified::calculate(int step)
-{
- //////////////////////////////////////////////////////////////////////////
- //! Cumulant K17 Kernel is based on
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
- //! and
- //! <a href="https://doi.org/10.1016/j.jcp.2017.07.004"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.07.004 ]</b></a>
- //!
- //! The cumulant kernel is executed in the following steps
- //!
- ////////////////////////////////////////////////////////////////////////////////
- //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
- //!
-
- using namespace std;
-
- //initializing of forcing stuff
- if (withForcing)
- {
- muForcingX1.DefineVar("x1", &muX1); muForcingX1.DefineVar("x2", &muX2); muForcingX1.DefineVar("x3", &muX3);
- muForcingX2.DefineVar("x1", &muX1); muForcingX2.DefineVar("x2", &muX2); muForcingX2.DefineVar("x3", &muX3);
- muForcingX3.DefineVar("x1", &muX1); muForcingX3.DefineVar("x2", &muX2); muForcingX3.DefineVar("x3", &muX3);
-
- muDeltaT = deltaT;
-
- muForcingX1.DefineVar("dt", &muDeltaT);
- muForcingX2.DefineVar("dt", &muDeltaT);
- muForcingX3.DefineVar("dt", &muDeltaT);
-
- muNu = (c1o1 / c3o1) * (c1o1 / collFactor - c1o1 / c2o1);
-
- muForcingX1.DefineVar("nu", &muNu);
- muForcingX2.DefineVar("nu", &muNu);
- muForcingX3.DefineVar("nu", &muNu);
- }
- /////////////////////////////////////
-
- localDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
- nonLocalDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
- restDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
-
- SPtr<BCArray3D> bcArray = this->getBCSet()->getBCArray();
-
- const int bcArrayMaxX1 = (int)bcArray->getNX1();
- const int bcArrayMaxX2 = (int)bcArray->getNX2();
- const int bcArrayMaxX3 = (int)bcArray->getNX3();
-
- int minX1 = ghostLayerWidth;
- int minX2 = ghostLayerWidth;
- int minX3 = ghostLayerWidth;
- int maxX1 = bcArrayMaxX1 - ghostLayerWidth;
- int maxX2 = bcArrayMaxX2 - ghostLayerWidth;
- int maxX3 = bcArrayMaxX3 - ghostLayerWidth;
-
- real omega = collFactor;
-
- for (int x3 = minX3; x3 < maxX3; x3++)
- {
- for (int x2 = minX2; x2 < maxX2; x2++)
- {
- for (int x1 = minX1; x1 < maxX1; x1++)
- {
- if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3))
- {
- int x1p = x1 + 1;
- int x2p = x2 + 1;
- int x3p = x3 + 1;
- //////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////
- //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on timestep is based on the esoteric twist algorithm
- //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
- //!
- ////////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////
-
- //E N T
- //c c c
- //////////
- //W S B
- //a a a
-
- //Rest is b
-
- //mfxyz
- //a - negative
- //b - null
- //c - positive
-
- // a b c
- //-1 0 1
-
- real mfcbb = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
- real mfbcb = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
- real mfbbc = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
- real mfccb = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
- real mfacb = (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3);
- real mfcbc = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
- real mfabc = (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3);
- real mfbcc = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3);
- real mfbac = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3);
- real mfccc = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3);
- real mfacc = (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3);
- real mfcac = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3);
- real mfaac = (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3);
-
- real mfabb = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3);
- real mfbab = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3);
- real mfbba = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p);
- real mfaab = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3);
- real mfcab = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3);
- real mfaba = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p);
- real mfcba = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p);
- real mfbaa = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p);
- real mfbca = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p);
- real mfaaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p);
- real mfcaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p);
- real mfaca = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p);
- real mfcca = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p);
-
- real mfbbb = (*this->restDistributions)(x1, x2, x3);
-
-
- real forces[3] = { c0o1, c0o1, c0o1 };
- if (withForcing)
- {
- muX1 = static_cast<real>(x1 - 1 + ix1 * maxX1);
- muX2 = static_cast<real>(x2 - 1 + ix2 * maxX2);
- muX3 = static_cast<real>(x3 - 1 + ix3 * maxX3);
-
- forcingX1 = muForcingX1.Eval();
- forcingX2 = muForcingX2.Eval();
- forcingX3 = muForcingX3.Eval();
-
- forces[0] = forcingX1 * deltaT;
- forces[1] = forcingX2 * deltaT;
- forces[2] = forcingX3 * deltaT;
- }
-
- vf::lbm::Distribution27 distribution;
-
- distribution.f[vf::lbm::dir::PZZ] = mfcbb;
- distribution.f[vf::lbm::dir::MZZ] = mfabb;
- distribution.f[vf::lbm::dir::ZPZ] = mfbcb;
- distribution.f[vf::lbm::dir::ZMZ] = mfbab;
- distribution.f[vf::lbm::dir::ZZP] = mfbbc;
- distribution.f[vf::lbm::dir::ZZM] = mfbba;
- distribution.f[vf::lbm::dir::PPZ] = mfccb;
- distribution.f[vf::lbm::dir::MMZ] = mfaab;
- distribution.f[vf::lbm::dir::PMZ] = mfcab;
- distribution.f[vf::lbm::dir::MPZ] = mfacb;
- distribution.f[vf::lbm::dir::PZP] = mfcbc;
- distribution.f[vf::lbm::dir::MZM] = mfaba;
- distribution.f[vf::lbm::dir::PZM] = mfcba;
- distribution.f[vf::lbm::dir::MZP] = mfabc;
- distribution.f[vf::lbm::dir::ZPP] = mfbcc;
- distribution.f[vf::lbm::dir::ZMM] = mfbaa;
- distribution.f[vf::lbm::dir::ZPM] = mfbca;
- distribution.f[vf::lbm::dir::ZMP] = mfbac;
- distribution.f[vf::lbm::dir::PPP] = mfccc;
- distribution.f[vf::lbm::dir::MPP] = mfacc;
- distribution.f[vf::lbm::dir::PMP] = mfcac;
- distribution.f[vf::lbm::dir::MMP] = mfaac;
- distribution.f[vf::lbm::dir::PPM] = mfcca;
- distribution.f[vf::lbm::dir::MPM] = mfaca;
- distribution.f[vf::lbm::dir::PMM] = mfcaa;
- distribution.f[vf::lbm::dir::MMM] = mfaaa;
- distribution.f[vf::lbm::dir::ZZZ] = mfbbb;
-
- vf::lbm::KernelParameter parameter {distribution, omega, forces};
- vf::lbm::cumulantChimera(parameter, vf::lbm::setRelaxationRatesK17);
-
- mfcbb = distribution.f[vf::lbm::dir::PZZ];
- mfabb = distribution.f[vf::lbm::dir::MZZ];
- mfbcb = distribution.f[vf::lbm::dir::ZPZ];
- mfbab = distribution.f[vf::lbm::dir::ZMZ];
- mfbbc = distribution.f[vf::lbm::dir::ZZP];
- mfbba = distribution.f[vf::lbm::dir::ZZM];
- mfccb = distribution.f[vf::lbm::dir::PPZ];
- mfaab = distribution.f[vf::lbm::dir::MMZ];
- mfcab = distribution.f[vf::lbm::dir::PMZ];
- mfacb = distribution.f[vf::lbm::dir::MPZ];
- mfcbc = distribution.f[vf::lbm::dir::PZP];
- mfaba = distribution.f[vf::lbm::dir::MZM];
- mfcba = distribution.f[vf::lbm::dir::PZM];
- mfabc = distribution.f[vf::lbm::dir::MZP];
- mfbcc = distribution.f[vf::lbm::dir::ZPP];
- mfbaa = distribution.f[vf::lbm::dir::ZMM];
- mfbca = distribution.f[vf::lbm::dir::ZPM];
- mfbac = distribution.f[vf::lbm::dir::ZMP];
- mfccc = distribution.f[vf::lbm::dir::PPP];
- mfacc = distribution.f[vf::lbm::dir::MPP];
- mfcac = distribution.f[vf::lbm::dir::PMP];
- mfaac = distribution.f[vf::lbm::dir::MMP];
- mfcca = distribution.f[vf::lbm::dir::PPM];
- mfaca = distribution.f[vf::lbm::dir::MPM];
- mfcaa = distribution.f[vf::lbm::dir::PMM];
- mfaaa = distribution.f[vf::lbm::dir::MMM];
- mfbbb = distribution.f[vf::lbm::dir::ZZZ];
-
- //////////////////////////////////////////////////////////////////////////
- //proof correctness
- //////////////////////////////////////////////////////////////////////////
-#ifdef PROOF_CORRECTNESS
- real drho_post = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca)
- + (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
- + (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
- real dif = distribution.getDensity_() - drho_post;
-#ifdef SINGLEPRECISION
- if (dif > 10.0E-7 || dif < -10.0E-7)
-#else
- if (dif > 10.0E-15 || dif < -10.0E-15)
-#endif
- {
- UB_THROW(UbException(UB_EXARGS, "rho=" + UbSystem::toString(distribution.getDensity_()) + ", rho_post=" + UbSystem::toString(drho_post)
- + " dif=" + UbSystem::toString(dif)
- + " rho is not correct for node " + UbSystem::toString(x1) + "," + UbSystem::toString(x2) + "," + UbSystem::toString(x3)
- + " in " + block.lock()->toString() + " step = " + UbSystem::toString(step)));
- }
-#endif
-
- (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = mfcbb;
- (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = mfbcb;
- (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = mfbbc;
- (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = mfccb;
- (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = mfacb;
- (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = mfcbc;
- (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = mfabc;
- (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = mfbcc;
- (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = mfbac;
- (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = mfccc;
- (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = mfacc;
- (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = mfcac;
- (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfaac;
-
- (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = mfabb;
- (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = mfbab;
- (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = mfbba;
- (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = mfaab;
- (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = mfcab;
- (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = mfaba;
- (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = mfcba;
- (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbaa;
- (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = mfbca;
- (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfaaa;
- (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfcaa;
- (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfaca;
- (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = mfcca;
- (*this->restDistributions)(x1, x2, x3) = mfbbb;
- //////////////////////////////////////////////////////////////////////////
-
- }
- }
- }
- }
-}
-//////////////////////////////////////////////////////////////////////////
-
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
index 2f5cbdc30ea8ae1368b59b6888a464700c1d8af5..c63e9f554367808d95c87fd630922333cf17790b 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
@@ -215,6 +215,33 @@ static const int ET_BNW = 11;
static const int ET_TSW = 12;
static const int ET_BNE = 12;
+static const int ET_P00 = 0;
+static const int ET_M00 = 0;
+static const int ET_0P0 = 1;
+static const int ET_0M0 = 1;
+static const int ET_00P = 2;
+static const int ET_00M = 2;
+static const int ET_PP0 = 3;
+static const int ET_MM0 = 3;
+static const int ET_PM0 = 4;
+static const int ET_MP0 = 4;
+static const int ET_P0P = 5;
+static const int ET_M0M = 5;
+static const int ET_P0M = 6;
+static const int ET_M0P = 6;
+static const int ET_0PP = 7;
+static const int ET_0MM = 7;
+static const int ET_0PM = 8;
+static const int ET_0MP = 8;
+static const int ET_PPP = 9;
+static const int ET_MMM = 9;
+static const int ET_MPP = 10;
+static const int ET_PMM = 10;
+static const int ET_PMP = 11;
+static const int ET_MPM = 11;
+static const int ET_MMP = 12;
+static const int ET_PPM = 12;
+
//////////////////////////////////////////////////////////////////////////
inline std::string getDirectionString(int direction)
{
diff --git a/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.cpp b/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..186574807cab3cea5dc3d46edb8111a392cf19ed
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.cpp
@@ -0,0 +1,165 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file CumulantK17LBMKernel.cpp
+//! \ingroup LBM
+//! \author Konstantin Kutscher, Martin Geier
+//=======================================================================================
+#include "K17CompressibleNavierStokes.h"
+
+#include <lbm/collision/CollisionParameter.h>
+#include <lbm/collision/K17CompressibleNavierStokes.h>
+#include <lbm/collision/TurbulentViscosity.h>
+
+#include "BCArray3D.h"
+#include "BCSet.h"
+#include "Block3D.h"
+#include "D3Q27EsoTwist3DSplittedVector.h"
+#include "D3Q27System.h"
+#include "DataSet3D.h"
+#include "LBMKernel.h"
+
+K17CompressibleNavierStokes::K17CompressibleNavierStokes()
+{
+ this->compressible = true;
+}
+
+void K17CompressibleNavierStokes::initDataSet()
+{
+ SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0] + 2, nx[1] + 2, nx[2] + 2, -999.9));
+ dataSet->setFdistributions(d);
+}
+
+SPtr<LBMKernel> K17CompressibleNavierStokes::clone()
+{
+ SPtr<LBMKernel> kernel(new K17CompressibleNavierStokes());
+ kernel->setNX(nx);
+ std::dynamic_pointer_cast<K17CompressibleNavierStokes>(kernel)->initDataSet();
+ kernel->setCollisionFactor(this->collFactor);
+ kernel->setBCSet(bcSet->clone(kernel));
+ kernel->setWithForcing(withForcing);
+ kernel->setForcingX1(muForcingX1);
+ kernel->setForcingX2(muForcingX2);
+ kernel->setForcingX3(muForcingX3);
+ kernel->setIndex(ix1, ix2, ix3);
+ kernel->setDeltaT(deltaT);
+ kernel->setBlock(block.lock());
+
+ return kernel;
+}
+
+void K17CompressibleNavierStokes::calculate(int step)
+{
+ (void)step; // parameter unused
+
+ mu::value_type muX1, muX2, muX3;
+ mu::value_type muDeltaT;
+ mu::value_type muNu;
+
+ if (withForcing) {
+ muForcingX1.DefineVar("x1", &muX1); muForcingX1.DefineVar("x2", &muX2); muForcingX1.DefineVar("x3", &muX3);
+ muForcingX2.DefineVar("x1", &muX1); muForcingX2.DefineVar("x2", &muX2); muForcingX2.DefineVar("x3", &muX3);
+ muForcingX3.DefineVar("x1", &muX1); muForcingX3.DefineVar("x2", &muX2); muForcingX3.DefineVar("x3", &muX3);
+
+ muDeltaT = deltaT;
+
+ muForcingX1.DefineVar("dt", &muDeltaT);
+ muForcingX2.DefineVar("dt", &muDeltaT);
+ muForcingX3.DefineVar("dt", &muDeltaT);
+
+ muNu = (c1o1 / c3o1) * (c1o1 / collFactor - c1o1 / c2o1);
+
+ muForcingX1.DefineVar("nu", &muNu);
+ muForcingX2.DefineVar("nu", &muNu);
+ muForcingX3.DefineVar("nu", &muNu);
+ }
+
+ auto localDistributions = std::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
+ auto nonLocalDistributions = std::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
+ auto restDistributions = std::dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
+
+ SPtr<BCArray3D> bcArray = this->getBCSet()->getBCArray();
+
+ const int bcArrayMaxX1 = (int)bcArray->getNX1();
+ const int bcArrayMaxX2 = (int)bcArray->getNX2();
+ const int bcArrayMaxX3 = (int)bcArray->getNX3();
+
+ const int minX1 = ghostLayerWidth;
+ const int minX2 = ghostLayerWidth;
+ const int minX3 = ghostLayerWidth;
+ const int maxX1 = bcArrayMaxX1 - ghostLayerWidth;
+ const int maxX2 = bcArrayMaxX2 - ghostLayerWidth;
+ const int maxX3 = bcArrayMaxX3 - ghostLayerWidth;
+
+ const real omega = collFactor;
+
+ for (int x3 = minX3; x3 < maxX3; x3++) {
+ for (int x2 = minX2; x2 < maxX2; x2++) {
+ for (int x1 = minX1; x1 < maxX1; x1++) {
+
+ if (bcArray->isSolid(x1, x2, x3) || bcArray->isUndefined(x1, x2, x3)) {
+ continue;
+ }
+
+ vf::lbm::CollisionParameter parameter;
+ dataSet->getFdistributions()->getDistribution(parameter.distribution, x1, x2, x3);
+
+ real forces[3] = { c0o1, c0o1, c0o1 };
+ if (withForcing) // TODO: add level factor?
+ {
+ muX1 = static_cast<real>(x1 - 1 + ix1 * maxX1);
+ muX2 = static_cast<real>(x2 - 1 + ix2 * maxX2);
+ muX3 = static_cast<real>(x3 - 1 + ix3 * maxX3);
+
+ const real forcingX1 = muForcingX1.Eval();
+ const real forcingX2 = muForcingX2.Eval();
+ const real forcingX3 = muForcingX3.Eval();
+
+ forces[0] = forcingX1 * deltaT;
+ forces[1] = forcingX2 * deltaT;
+ forces[2] = forcingX3 * deltaT;
+ }
+ parameter.forceX = forces[0];
+ parameter.forceY = forces[1];
+ parameter.forceZ = forces[2];
+
+ parameter.omega = omega;
+
+ real quadricLimiter[3] = { 0.01, 0.01, 0.01 }; // TODO: Where do we configure the quadricLimiter?
+ parameter.quadricLimiter = quadricLimiter;
+
+ const bool writeMacroscopicVariables = false;
+ vf::lbm::MacroscopicValues mv; // not used
+ vf::lbm::TurbulentViscosity tv; // not used
+ vf::lbm::runK17CompressibleNavierStokes<vf::lbm::TurbulenceModel::None, writeMacroscopicVariables>(parameter,
+ mv, tv);
+ dataSet->getFdistributions()->setDistribution(parameter.distribution, x1, x2, x3);
+ }
+ }
+ }
+}
diff --git a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernelUnified.h b/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.h
similarity index 73%
rename from src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernelUnified.h
rename to src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.h
index 9c6876f4d4de855f931ea90e8764cb1e9a39202a..6c6268206b09398ab07150f3a03fd2cdeaf40e30 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernelUnified.h
+++ b/src/cpu/VirtualFluidsCore/LBM/K17CompressibleNavierStokes.h
@@ -31,15 +31,14 @@
//! \author Konstantin Kutscher, Martin Geier
//=======================================================================================
-#ifndef CumulantK17LBMKernelUnified_h__
-#define CumulantK17LBMKernelUnified_h__
+#ifndef CPU_K17COMPRESSIBLENAVIERSTOKES_H
+#define CPU_K17COMPRESSIBLENAVIERSTOKES_H
+
+#include <memory>
+
+#include <basics/DataTypes.h>
#include "LBMKernel.h"
-#include "BCSet.h"
-#include "D3Q27System.h"
-#include "basics/utilities/UbTiming.h"
-#include "basics/container/CbArray4D.h"
-#include "basics/container/CbArray3D.h"
//! \brief Compressible cumulant LBM kernel.
//! \details LBM implementation that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model
@@ -48,30 +47,19 @@
//! <a href="http://dx.doi.org/10.1016/j.jcp.2017.05.040"><b>[ Geier et al., (2017), 10.1016/j.jcp.2017.05.040]</b></a>,
//! <a href="http://dx.doi.org/10.1016/j.jcp.2017.07.004"><b>[ Geier et al., (2017), 10.1016/j.jcp.2017.07.004]</b></a>
//!
-class CumulantK17LBMKernelUnified : public LBMKernel
+class K17CompressibleNavierStokes : public LBMKernel
{
public:
- CumulantK17LBMKernelUnified();
- ~CumulantK17LBMKernelUnified() = default;
+ K17CompressibleNavierStokes();
void calculate(int step) override;
- SPtr<LBMKernel> clone() override;
- real getCalculationTime() override { return .0; }
-
-protected:
- virtual void initDataSet();
- real f[D3Q27System::ENDF + 1];
+ std::shared_ptr<LBMKernel> clone() override;
+ real getCalculationTime() override
+ {
+ return .0;
+ }
- CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr localDistributions;
- CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributions;
- CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr restDistributions;
-
- mu::value_type muX1, muX2, muX3;
- mu::value_type muDeltaT;
- mu::value_type muNu;
- real forcingX1;
- real forcingX2;
- real forcingX3;
+private:
+ void initDataSet();
};
-
-#endif // CumulantK17LBMKernel_h__
\ No newline at end of file
+#endif
diff --git a/src/gpu/VirtualFluids_GPU/CMakeLists.txt b/src/gpu/VirtualFluids_GPU/CMakeLists.txt
index 4cd0d1df5db24385c2b528294bb665d405e73937..ef2b2bed068f415799c7aad19c4c7a5ad6567f49 100644
--- a/src/gpu/VirtualFluids_GPU/CMakeLists.txt
+++ b/src/gpu/VirtualFluids_GPU/CMakeLists.txt
@@ -25,7 +25,7 @@ if(BUILD_VF_UNIT_TESTS)
set_source_files_properties(DataStructureInitializer/GridReaderGenerator/GridGeneratorTest.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(DataStructureInitializer/GridReaderGenerator/IndexRearrangementForStreamsTest.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(Kernel/Kernels/BasicKernels/FluidFlow/Compressible/CumulantK17/CumulantK17Test.cpp PROPERTIES LANGUAGE CUDA)
- set_source_files_properties(Kernel/Utilities/DistributionHelperTests.cpp PROPERTIES LANGUAGE CUDA)
+ set_source_files_properties(LBM/GPUHelperFunctions/KernelUtilitiesTests.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(Parameter/ParameterTest.cpp PROPERTIES LANGUAGE CUDA)
set_source_files_properties(PreCollisionInteractor/ActuatorFarmInlinesTest.cpp PROPERTIES LANGUAGE CUDA)
endif()
diff --git a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
index 8907e846757c8923c3aed46f9c90d6c67f465eee..1ae2c62d05eacad55e1101e19f8195b5dd937a42 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/CalcMac27.cu
@@ -30,12 +30,14 @@
//! \ingroup GPU
//! \author Martin Schoenherr, Soeren Peters
//======================================================================================
-#include "LBM/LB.h"
+#include "LBM/LB.h"
+
+#include "LBM/GPUHelperFunctions/KernelUtilities.h"
+
#include "lbm/constants/D3Q27.h"
#include "basics/constants/NumericConstants.h"
#include "lbm/MacroscopicQuantities.h"
-#include "Kernel/Utilities/DistributionHelper.cuh"
using namespace vf::basics::constant;
using namespace vf::lbm::dir;
@@ -79,13 +81,17 @@ __global__ void LBCalcMac27(
vyD[k] = c0o1;
vzD[k] = c0o1;
- DistributionWrapper distr_wrapper(distributions, numberOfLBnodes, isEvenTimestep, k, neighborX, neighborY, neighborZ);
- const auto& distribution = distr_wrapper.distribution;
+ Distributions27 dist;
+ vf::gpu::getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
+ vf::gpu::ListIndices listIndices(k, neighborX, neighborY, neighborZ);
+
+ real distribution[27];
+ vf::gpu::read(distribution, dist, listIndices);
- rhoD[k] = vf::lbm::getDensity(distribution.f);
- vxD[k] = vf::lbm::getIncompressibleVelocityX1(distribution.f);
- vyD[k] = vf::lbm::getIncompressibleVelocityX2(distribution.f);
- vzD[k] = vf::lbm::getIncompressibleVelocityX3(distribution.f);
+ rhoD[k] = vf::lbm::getDensity(distribution);
+ vxD[k] = vf::lbm::getIncompressibleVelocityX1(distribution);
+ vyD[k] = vf::lbm::getIncompressibleVelocityX2(distribution);
+ vzD[k] = vf::lbm::getIncompressibleVelocityX3(distribution);
}
@@ -297,14 +303,18 @@ __global__ void LBCalcMacCompSP27(
if (!isValidFluidNode(geoD[nodeIndex]))
return;
- DistributionWrapper distr_wrapper(distributions, numberOfLBnodes, isEvenTimestep, nodeIndex, neighborX, neighborY, neighborZ);
- const auto &distribution = distr_wrapper.distribution;
+ Distributions27 dist;
+ vf::gpu::getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
+ vf::gpu::ListIndices listIndices(nodeIndex, neighborX, neighborY, neighborZ);
+
+ real distribution[27];
+ vf::gpu::read(distribution, dist, listIndices);
- rhoD[nodeIndex] = vf::lbm::getDensity(distribution.f);
- vxD[nodeIndex] = vf::lbm::getCompressibleVelocityX1(distribution.f, rhoD[nodeIndex]);
- vyD[nodeIndex] = vf::lbm::getCompressibleVelocityX2(distribution.f, rhoD[nodeIndex]);
- vzD[nodeIndex] = vf::lbm::getCompressibleVelocityX3(distribution.f, rhoD[nodeIndex]);
- pressD[nodeIndex] = vf::lbm::getPressure(distribution.f, rhoD[nodeIndex], vxD[nodeIndex], vyD[nodeIndex], vzD[nodeIndex]);
+ rhoD[nodeIndex] = vf::lbm::getDensity(distribution);
+ vxD[nodeIndex] = vf::lbm::getCompressibleVelocityX1(distribution, rhoD[nodeIndex]);
+ vyD[nodeIndex] = vf::lbm::getCompressibleVelocityX2(distribution, rhoD[nodeIndex]);
+ vzD[nodeIndex] = vf::lbm::getCompressibleVelocityX3(distribution, rhoD[nodeIndex]);
+ pressD[nodeIndex] = vf::lbm::getPressure(distribution, rhoD[nodeIndex], vxD[nodeIndex], vyD[nodeIndex], vzD[nodeIndex]);
}
diff --git a/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu b/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu
index 2c482c3b0fc368c52cca1e74246c75210131b326..11cbb45c565e1346a9664db8a6c6164d8031aa3e 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/Cumulant27chim.cu
@@ -31,16 +31,16 @@
//! \author Martin Schoenherr
//=======================================================================================
/* Device code */
-#include "LBM/LB.h"
+#include "LBM/LB.h"
#include "lbm/constants/D3Q27.h"
#include <basics/constants/NumericConstants.h>
using namespace vf::basics::constant;
using namespace vf::lbm::dir;
-#include "math.h"
+#include <cmath>
-#include "lbm/Chimera.h"
+#include <lbm/ChimeraTransformation.h>
////////////////////////////////////////////////////////////////////////////////
diff --git a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
index bad30b31130347dd2543714318effb61dbd21972..2492cb2ac438c1276981b4ce019556b76b51173c 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/GridScaling/scaleFC_compressible.cu
@@ -35,7 +35,6 @@
#include "LBM/GPUHelperFunctions/KernelUtilities.h"
#include "LBM/GPUHelperFunctions/ScalingUtilities.h"
-#include <lbm/KernelParameter.h>
#include <lbm/refinement/InterpolationFC.h>
#include <lbm/interpolation/InterpolationCoefficients.h>
diff --git a/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu b/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu
index ecd7427665427376aaee290e918fd5c723576f73..3289d6543aa2e13c4beb1310985f40ad73445a4d 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/SlipBCs27.cu
@@ -32,7 +32,6 @@
//======================================================================================
#include "LBM/LB.h"
#include "lbm/constants/D3Q27.h"
-#include "Kernel/Utilities/DistributionHelper.cuh"
#include "basics/constants/NumericConstants.h"
#include "LBM/GPUHelperFunctions/KernelUtilities.h"
diff --git a/src/gpu/VirtualFluids_GPU/GPU/StressBCs27.cu b/src/gpu/VirtualFluids_GPU/GPU/StressBCs27.cu
index 1cc5017816aed29d52e74823a8c910bfed35ad42..4b9db6d8b05a0ed958b7416ae8f3d07d761fb71c 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/StressBCs27.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/StressBCs27.cu
@@ -42,7 +42,6 @@
#include "LBM/LB.h"
#include "lbm/constants/D3Q27.h"
-#include "Kernel/Utilities/DistributionHelper.cuh"
#include <basics/constants/NumericConstants.h>
#include "LBM/GPUHelperFunctions/KernelUtilities.h"
diff --git a/src/gpu/VirtualFluids_GPU/GPU/TurbulenceIntensity.cu b/src/gpu/VirtualFluids_GPU/GPU/TurbulenceIntensity.cu
index 82e5f98fda0086458f2bc937dbe33e7b66feb2c5..46911fa39f61e455aa46af2645d37f1ef3786f3f 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/TurbulenceIntensity.cu
+++ b/src/gpu/VirtualFluids_GPU/GPU/TurbulenceIntensity.cu
@@ -12,7 +12,6 @@
#include "basics/constants/NumericConstants.h"
#include "lbm/MacroscopicQuantities.h"
-#include "../Kernel/Utilities/DistributionHelper.cuh"
#include "LBM/GPUHelperFunctions/KernelUtilities.h"
@@ -44,34 +43,37 @@ __global__ void CalcTurbulenceIntensity(
//!
const unsigned nodeIndex = getNodeIndex();
- if (nodeIndex >= numberOfLBnodes)
- return;
+ if (nodeIndex >= numberOfLBnodes)
+ return;
- if (!isValidFluidNode(typeOfGridNode[nodeIndex]))
- return;
+ if (!isValidFluidNode(typeOfGridNode[nodeIndex]))
+ return;
- DistributionWrapper distr_wrapper(distributions, numberOfLBnodes, isEvenTimestep, nodeIndex, neighborX, neighborY, neighborZ);
- const auto &distribution = distr_wrapper.distribution;
+ Distributions27 dist;
+ vf::gpu::getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
+ vf::gpu::ListIndices listIndices(nodeIndex, neighborX, neighborY, neighborZ);
- // analogue to LBCalcMacCompSP27
- real rho = vf::lbm::getDensity(distribution.f);
- real vx = vf::lbm::getCompressibleVelocityX1(distribution.f, rho);
- real vy = vf::lbm::getCompressibleVelocityX2(distribution.f, rho);
- real vz = vf::lbm::getCompressibleVelocityX3(distribution.f, rho);
+ real distribution[27];
+ vf::gpu::read(distribution, dist, listIndices);
+ // analogue to LBCalcMacCompSP27
+ real rho = vf::lbm::getDensity(distribution);
+ real vx = vf::lbm::getCompressibleVelocityX1(distribution, rho);
+ real vy = vf::lbm::getCompressibleVelocityX2(distribution, rho);
+ real vz = vf::lbm::getCompressibleVelocityX3(distribution, rho);
- // compute subtotals:
- // fluctuations
- vxx[nodeIndex] = vxx[nodeIndex] + vx * vx;
- vyy[nodeIndex] = vyy[nodeIndex] + vy * vy;
- vzz[nodeIndex] = vzz[nodeIndex] + vz * vz;
- vxy[nodeIndex] = vxy[nodeIndex] + vx * vy;
- vxz[nodeIndex] = vxz[nodeIndex] + vx * vz;
- vyz[nodeIndex] = vyz[nodeIndex] + vy * vz;
+ // compute subtotals:
+ // fluctuations
+ vxx[nodeIndex] = vxx[nodeIndex] + vx * vx;
+ vyy[nodeIndex] = vyy[nodeIndex] + vy * vy;
+ vzz[nodeIndex] = vzz[nodeIndex] + vz * vz;
+ vxy[nodeIndex] = vxy[nodeIndex] + vx * vy;
+ vxz[nodeIndex] = vxz[nodeIndex] + vx * vz;
+ vyz[nodeIndex] = vyz[nodeIndex] + vy * vz;
- // velocity (for mean velocity)
- vx_mean[nodeIndex] = vx_mean[nodeIndex] + vx;
- vy_mean[nodeIndex] = vy_mean[nodeIndex] + vy;
- vz_mean[nodeIndex] = vz_mean[nodeIndex] + vz;
+ // velocity (for mean velocity)
+ vx_mean[nodeIndex] = vx_mean[nodeIndex] + vx;
+ vy_mean[nodeIndex] = vy_mean[nodeIndex] + vy;
+ vz_mean[nodeIndex] = vz_mean[nodeIndex] + vz;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplusUnified.cu b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplusUnified.cu
deleted file mode 100644
index bafd4477dc611b77f17c896dc85832b71acdccc0..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplusUnified.cu
+++ /dev/null
@@ -1,57 +0,0 @@
-#include "B15CompressibleNavierStokesBGKplusUnified.h"
-
-#include <stdexcept>
-
-#include "Parameter/Parameter.h"
-#include "../RunLBMKernel.cuh"
-
-#include <lbm/BGK.h>
-#include <lbm/KernelParameter.h>
-
-
-namespace vf
-{
-namespace gpu
-{
-
-
-B15CompressibleNavierStokesBGKplusUnified::B15CompressibleNavierStokesBGKplusUnified(std::shared_ptr<Parameter> para, int level)
- : KernelImp(para, level)
-{
-#ifndef BUILD_CUDA_LTO
- throw std::invalid_argument("To use the BKGUnified kernel, pass -DBUILD_CUDA_LTO=ON to cmake. Requires: CUDA 11.2 & cc 5.0");
-#endif
-
- myPreProcessorTypes.push_back(InitCompSP27);
-
-
-
- this->cudaGrid = cuda::CudaGrid(para->getParD(level)->numberofthreads, para->getParD(level)->numberOfNodes);
-}
-
-
-void B15CompressibleNavierStokesBGKplusUnified::run()
-{
- GPUKernelParameter kernelParameter{
- para->getParD(level)->omega,
- para->getParD(level)->typeOfGridNode,
- para->getParD(level)->neighborX,
- para->getParD(level)->neighborY,
- para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- (int)para->getParD(level)->numberOfNodes,
- nullptr, /* forces not used in bgk kernel */
- para->getParD(level)->isEvenTimestep };
-
- auto lambda = [] __device__(lbm::KernelParameter parameter) {
- return lbm::bgk(parameter);
- };
-
- runKernel<<<cudaGrid.grid, cudaGrid.threads>>>(lambda, kernelParameter);
-
- getLastCudaError("LB_Kernel_BGKUnified execution failed");
-}
-
-
-}
-}
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplusUnified.h b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplusUnified.h
deleted file mode 100644
index ffc8ce41d13369b78d19ff0d55bab988696cdb5d..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplusUnified.h
+++ /dev/null
@@ -1,22 +0,0 @@
-#ifndef B15CompressibleNavierStokesBGKplusUnified_H
-#define B15CompressibleNavierStokesBGKplusUnified_H
-
-#include "Kernel/KernelImp.h"
-
-namespace vf
-{
-namespace gpu
-{
-
-class B15CompressibleNavierStokesBGKplusUnified : public KernelImp
-{
-public:
- B15CompressibleNavierStokesBGKplusUnified(std::shared_ptr<Parameter> para, int level);
-
- void run();
-};
-
-}
-}
-
-#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesUnified.cu b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesUnified.cu
deleted file mode 100644
index 6e1fbfbe60808fda3980fff462f237096b6a4568..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesUnified.cu
+++ /dev/null
@@ -1,55 +0,0 @@
-#include "K15CompressibleNavierStokesUnified.h"
-
-#include <stdexcept>
-
-#include "../RunLBMKernel.cuh"
-
-#include "Parameter/Parameter.h"
-
-#include <lbm/CumulantChimera.h>
-
-namespace vf
-{
-namespace gpu
-{
-
-K15CompressibleNavierStokesUnified::K15CompressibleNavierStokesUnified(std::shared_ptr<Parameter> para, int level)
- : KernelImp(para, level)
-{
-#ifndef BUILD_CUDA_LTO
- throw std::invalid_argument(
- "To use the CumulantK15Unified kernel, pass -DBUILD_CUDA_LTO=ON to cmake. Requires: CUDA 11.2 & cc 5.0");
-#endif
-
- myPreProcessorTypes.push_back(InitCompSP27);
-
-
-
- this->cudaGrid = cuda::CudaGrid(para->getParD(level)->numberofthreads, para->getParD(level)->numberOfNodes);
-}
-
-void K15CompressibleNavierStokesUnified::run()
-{
- GPUKernelParameter kernelParameter{
- para->getParD(level)->omega,
- para->getParD(level)->typeOfGridNode,
- para->getParD(level)->neighborX,
- para->getParD(level)->neighborY,
- para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- (int)para->getParD(level)->numberOfNodes,
- para->getParD(level)->forcing,
- para->getParD(level)->isEvenTimestep };
-
- auto lambda = [] __device__(lbm::KernelParameter parameter) {
- return lbm::cumulantChimera(parameter, lbm::setRelaxationRatesK15);
- };
-
- vf::gpu::runKernel<<<cudaGrid.grid, cudaGrid.threads>>>(lambda, kernelParameter);
-
- getLastCudaError("LB_Kernel_CumulantK15Comp execution failed");
-}
-
-
-}
-}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesUnified.h b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesUnified.h
deleted file mode 100644
index e16211520a4e310d4ee1a73c00437a4d18441354..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesUnified.h
+++ /dev/null
@@ -1,21 +0,0 @@
-#ifndef K15CompressibleNavierStokesUnified_H
-#define K15CompressibleNavierStokesUnified_H
-
-#include "Kernel/KernelImp.h"
-
-namespace vf
-{
-namespace gpu
-{
-class K15CompressibleNavierStokesUnified : public KernelImp
-{
-public:
- K15CompressibleNavierStokesUnified(std::shared_ptr<Parameter> para, int level);
-
- void run();
-};
-
-}
-}
-
-#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu
index fccd9bfd246c9bd927bd9dcce6da8993e3d57374..02df4e149018668f830b161f86639f102ca7b12d 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.cu
@@ -1,118 +1,123 @@
#include "K17CompressibleNavierStokes.h"
+
+#include <optional>
+#include <stdexcept>
+
+#include <cuda.h>
+
#include <logger/Logger.h>
-#include "Parameter/Parameter.h"
+
+#include <lbm/collision/TurbulentViscosity.h>
+#include <lbm/collision/K17CompressibleNavierStokes.h>
+
+#include "Kernel/KernelImp.h"
#include "Parameter/CudaStreamManager.h"
-#include "K17CompressibleNavierStokes_Device.cuh"
+#include "Parameter/Parameter.h"
+#include "LBM/GPUHelperFunctions/KernelUtilities.h"
+#include "LBM/GPUHelperFunctions/RunCollision.cuh"
-#include <cuda.h>
+template <vf::lbm::TurbulenceModel turbulenceModel>
+std::shared_ptr<K17CompressibleNavierStokes<turbulenceModel>>
+K17CompressibleNavierStokes<turbulenceModel>::getNewInstance(std::shared_ptr<Parameter> para, int level)
+{
+ return std::shared_ptr<K17CompressibleNavierStokes<turbulenceModel>>(new K17CompressibleNavierStokes(para, level));
+}
-template<TurbulenceModel turbulenceModel>
-std::shared_ptr< K17CompressibleNavierStokes<turbulenceModel> > K17CompressibleNavierStokes<turbulenceModel>::getNewInstance(std::shared_ptr<Parameter> para, int level)
+vf::gpu::GPUCollisionParameter getCollisionParameter(const std::shared_ptr<Parameter>& para, int level, const uint* indices,
+ uint sizeIndices)
{
- return std::shared_ptr<K17CompressibleNavierStokes<turbulenceModel> >(new K17CompressibleNavierStokes<turbulenceModel>(para,level));
+ vf::gpu::GPUCollisionParameter collisionParameter { para->getParD(level)->omega,
+ para->getParD(level)->neighborX,
+ para->getParD(level)->neighborY,
+ para->getParD(level)->neighborZ,
+ para->getParD(level)->distributions.f[0],
+ para->getParD(level)->rho,
+ para->getParD(level)->velocityX,
+ para->getParD(level)->velocityY,
+ para->getParD(level)->velocityZ,
+ para->getParD(level)->turbViscosity,
+ para->getSGSConstant(),
+ (int)para->getParD(level)->numberOfNodes,
+ vf::gpu::getForceFactor(level),
+ para->getForcesDev(),
+ para->getParD(level)->forceX_SP,
+ para->getParD(level)->forceY_SP,
+ para->getParD(level)->forceZ_SP,
+ para->getQuadricLimitersDev(),
+ para->getParD(level)->isEvenTimestep,
+ indices,
+ sizeIndices };
+
+ return collisionParameter;
}
-template<TurbulenceModel turbulenceModel>
+template <vf::lbm::TurbulenceModel turbulenceModel>
void K17CompressibleNavierStokes<turbulenceModel>::run()
{
- K17CompressibleNavierStokes_Device < turbulenceModel, false, false > <<< cudaGrid.grid, cudaGrid.threads >>>( para->getParD(level)->omega,
- para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- para->getParD(level)->rho,
- para->getParD(level)->velocityX, para->getParD(level)->velocityY, para->getParD(level)->velocityZ,
- para->getParD(level)->turbViscosity,
- para->getSGSConstant(),
- para->getParD(level)->numberOfNodes,
- level,
- para->getForcesDev(),
- para->getParD(level)->forceX_SP, para->getParD(level)->forceY_SP, para->getParD(level)->forceZ_SP,
- para->getQuadricLimitersDev(),
- para->getParD(level)->isEvenTimestep,
- para->getParD(level)->taggedFluidNodeIndices[CollisionTemplate::Default],
- para->getParD(level)->numberOfTaggedFluidNodes[CollisionTemplate::Default]);
+ vf::gpu::GPUCollisionParameter kernelParameter =
+ getCollisionParameter(para, level, para->getParD(level)->taggedFluidNodeIndices[CollisionTemplate::Default],
+ para->getParD(level)->numberOfTaggedFluidNodes[CollisionTemplate::Default]);
- getLastCudaError("K17CompressibleNavierStokes_Device execution failed");
+ auto collision = [] __device__(vf::lbm::CollisionParameter & parameter, vf::lbm::MacroscopicValues & macroscopicValues,
+ vf::lbm::TurbulentViscosity & turbulentViscosity) {
+ return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, false>(parameter, macroscopicValues, turbulentViscosity);
+ };
+
+ vf::gpu::runCollision<decltype(collision), turbulenceModel, false, false><<<cudaGrid.grid, cudaGrid.threads>>>(collision, kernelParameter);
+
+ getLastCudaError("K17CompressibleNavierStokesUnified execution failed");
}
-template<TurbulenceModel turbulenceModel>
-void K17CompressibleNavierStokes<turbulenceModel>::runOnIndices( const unsigned int *indices, unsigned int size_indices, CollisionTemplate collisionTemplate, CudaStreamIndex streamIndex )
+template <vf::lbm::TurbulenceModel turbulenceModel>
+void K17CompressibleNavierStokes<turbulenceModel>::runOnIndices(const unsigned int* indices, unsigned int size_indices,
+ CollisionTemplate collisionTemplate,
+ CudaStreamIndex streamIndex)
{
cudaStream_t stream = para->getStreamManager()->getStream(streamIndex);
- switch (collisionTemplate)
- {
- case CollisionTemplate::Default:
- K17CompressibleNavierStokes_Device < turbulenceModel, false, false > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>(para->getParD(level)->omega,
- para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- para->getParD(level)->rho,
- para->getParD(level)->velocityX, para->getParD(level)->velocityY, para->getParD(level)->velocityZ,
- para->getParD(level)->turbViscosity,
- para->getSGSConstant(),
- para->getParD(level)->numberOfNodes,
- level,
- para->getForcesDev(),
- para->getParD(level)->forceX_SP, para->getParD(level)->forceY_SP, para->getParD(level)->forceZ_SP,
- para->getQuadricLimitersDev(),
- para->getParD(level)->isEvenTimestep,
- indices,
- size_indices);
- break;
+ switch (collisionTemplate) {
+ case CollisionTemplate::Default: {
+ vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
+
+ auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
+ return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, false>(parameter, macroscopicValues, turbulentViscosity);
+ };
+
+ vf::gpu::runCollision<decltype(collision), turbulenceModel, false, false><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
- case CollisionTemplate::WriteMacroVars:
- K17CompressibleNavierStokes_Device < turbulenceModel, true, false > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>( para->getParD(level)->omega,
- para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- para->getParD(level)->rho,
- para->getParD(level)->velocityX, para->getParD(level)->velocityY, para->getParD(level)->velocityZ,
- para->getParD(level)->turbViscosity,
- para->getSGSConstant(),
- para->getParD(level)->numberOfNodes,
- level,
- para->getForcesDev(),
- para->getParD(level)->forceX_SP, para->getParD(level)->forceY_SP, para->getParD(level)->forceZ_SP,
- para->getQuadricLimitersDev(),
- para->getParD(level)->isEvenTimestep,
- indices,
- size_indices);
break;
+ }
+ case CollisionTemplate::WriteMacroVars: {
+ vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
+
+ auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
+ return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, true>(parameter, macroscopicValues, turbulentViscosity);
+ };
+ vf::gpu::runCollision<decltype(collision), turbulenceModel, true, false><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
+ break;
+ }
case CollisionTemplate::SubDomainBorder:
- case CollisionTemplate::AllFeatures:
- K17CompressibleNavierStokes_Device < turbulenceModel, true, true > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>( para->getParD(level)->omega,
- para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- para->getParD(level)->rho,
- para->getParD(level)->velocityX, para->getParD(level)->velocityY, para->getParD(level)->velocityZ,
- para->getParD(level)->turbViscosity,
- para->getSGSConstant(),
- para->getParD(level)->numberOfNodes,
- level,
- para->getForcesDev(),
- para->getParD(level)->forceX_SP, para->getParD(level)->forceY_SP, para->getParD(level)->forceZ_SP,
- para->getQuadricLimitersDev(),
- para->getParD(level)->isEvenTimestep,
- indices,
- size_indices);
+ case CollisionTemplate::AllFeatures: {
+ vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
+
+ auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
+ return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, true>(parameter, macroscopicValues, turbulentViscosity);
+ };
+ vf::gpu::runCollision<decltype(collision), turbulenceModel, true, true><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
+
break;
+ }
+ case CollisionTemplate::ApplyBodyForce: {
+ vf::gpu::GPUCollisionParameter kernelParameter = getCollisionParameter(para, level, indices, size_indices);
+
+ auto collision = [] __device__(vf::lbm::CollisionParameter& parameter, vf::lbm::MacroscopicValues& macroscopicValues, vf::lbm::TurbulentViscosity& turbulentViscosity) {
+ return vf::lbm::runK17CompressibleNavierStokes<turbulenceModel, false>(parameter, macroscopicValues, turbulentViscosity);
+ };
+ vf::gpu::runCollision<decltype(collision), turbulenceModel, false, true><<<cudaGrid.grid, cudaGrid.threads, 0, stream>>>(collision, kernelParameter);
- case CollisionTemplate::ApplyBodyForce:
- K17CompressibleNavierStokes_Device < turbulenceModel, false, true > <<< cudaGrid.grid, cudaGrid.threads, 0, stream >>>( para->getParD(level)->omega,
- para->getParD(level)->neighborX, para->getParD(level)->neighborY, para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- para->getParD(level)->rho,
- para->getParD(level)->velocityX, para->getParD(level)->velocityY, para->getParD(level)->velocityZ,
- para->getParD(level)->turbViscosity,
- para->getSGSConstant(),
- para->getParD(level)->numberOfNodes,
- level,
- para->getForcesDev(),
- para->getParD(level)->forceX_SP, para->getParD(level)->forceY_SP, para->getParD(level)->forceZ_SP,
- para->getQuadricLimitersDev(),
- para->getParD(level)->isEvenTimestep,
- indices,
- size_indices);
break;
+ }
default:
throw std::runtime_error("Invalid CollisionTemplate in CumulantK17::runOnIndices()");
break;
@@ -121,23 +126,18 @@ void K17CompressibleNavierStokes<turbulenceModel>::runOnIndices( const unsigned
getLastCudaError("K17CompressibleNavierStokes_Device execution failed");
}
-template<TurbulenceModel turbulenceModel>
-K17CompressibleNavierStokes<turbulenceModel>::K17CompressibleNavierStokes(std::shared_ptr<Parameter> para, int level)
+template <vf::lbm::TurbulenceModel turbulenceModel>
+K17CompressibleNavierStokes<turbulenceModel>::K17CompressibleNavierStokes(std::shared_ptr<Parameter> para, int level) : KernelImp(para, level)
{
- this->para = para;
- this->level = level;
-
myPreProcessorTypes.push_back(InitCompSP27);
-
-
this->cudaGrid = vf::cuda::CudaGrid(para->getParD(level)->numberofthreads, para->getParD(level)->numberOfNodes);
this->kernelUsesFluidNodeIndices = true;
VF_LOG_INFO("Using turbulence model: {}", turbulenceModel);
}
-template class K17CompressibleNavierStokes<TurbulenceModel::AMD>;
-template class K17CompressibleNavierStokes<TurbulenceModel::Smagorinsky>;
-template class K17CompressibleNavierStokes<TurbulenceModel::QR>;
-template class K17CompressibleNavierStokes<TurbulenceModel::None>;
+template class K17CompressibleNavierStokes<vf::lbm::TurbulenceModel::AMD>;
+template class K17CompressibleNavierStokes<vf::lbm::TurbulenceModel::Smagorinsky>;
+template class K17CompressibleNavierStokes<vf::lbm::TurbulenceModel::QR>;
+template class K17CompressibleNavierStokes<vf::lbm::TurbulenceModel::None>;
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.h b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.h
index 543b2de3e2ef9ad7c4bf4881aae2afb6110ec552..376d6eab87760dcbd7d2c607c913c1a7bb8b911c 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.h
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.h
@@ -1,20 +1,24 @@
-#ifndef K17CompressibleNavierStokes_H
-#define K17CompressibleNavierStokes_H
+#ifndef GPU_K17CompressibleNavierStokes_H
+#define GPU_K17CompressibleNavierStokes_H
#include "Kernel/KernelImp.h"
-#include "Parameter/Parameter.h"
-template<TurbulenceModel turbulenceModel>
+#include <lbm/collision/TurbulentViscosity.h>
+
+class Parameter;
+
+template <vf::lbm::TurbulenceModel turbulenceModel>
class K17CompressibleNavierStokes : public KernelImp
{
public:
- static std::shared_ptr< K17CompressibleNavierStokes<turbulenceModel> > getNewInstance(std::shared_ptr< Parameter> para, int level);
- void run() override;
- void runOnIndices(const unsigned int *indices, unsigned int size_indices, CollisionTemplate collisionTemplate, CudaStreamIndex streamIndex) override;
+ static std::shared_ptr<K17CompressibleNavierStokes<turbulenceModel>> getNewInstance(std::shared_ptr<Parameter> para,
+ int level);
+ void run() override;
+ void runOnIndices(const unsigned int* indices, unsigned int size_indices, CollisionTemplate collisionTemplate,
+ CudaStreamIndex streamIndex) override;
private:
- K17CompressibleNavierStokes();
K17CompressibleNavierStokes(std::shared_ptr<Parameter> para, int level);
};
-#endif
+#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesUnified.cu b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesUnified.cu
deleted file mode 100644
index 45a998e215b999a9f920358132865fbfccc458aa..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesUnified.cu
+++ /dev/null
@@ -1,56 +0,0 @@
-#include "K17CompressibleNavierStokesUnified.h"
-
-#include <stdexcept>
-
-#include "Parameter/Parameter.h"
-#include "../RunLBMKernel.cuh"
-
-#include <lbm/CumulantChimera.h>
-
-namespace vf
-{
-namespace gpu
-{
-
-
-K17CompressibleNavierStokesUnified::K17CompressibleNavierStokesUnified(std::shared_ptr<Parameter> para, int level)
- : KernelImp(para, level)
-{
-#ifndef BUILD_CUDA_LTO
- throw std::invalid_argument("To use the CumulantK17Unified kernel, pass -DBUILD_CUDA_LTO=ON to cmake. Requires: CUDA 11.2 & cc 5.0");
-#endif
-
- myPreProcessorTypes.push_back(InitCompSP27);
-
-
-
- this->cudaGrid = cuda::CudaGrid(para->getParD(level)->numberofthreads, para->getParD(level)->numberOfNodes);
-}
-
-
-
-void K17CompressibleNavierStokesUnified::run()
-{
- GPUKernelParameter kernelParameter{
- para->getParD(level)->omega,
- para->getParD(level)->typeOfGridNode,
- para->getParD(level)->neighborX,
- para->getParD(level)->neighborY,
- para->getParD(level)->neighborZ,
- para->getParD(level)->distributions.f[0],
- (int)para->getParD(level)->numberOfNodes,
- para->getParD(level)->forcing,
- para->getParD(level)->isEvenTimestep };
-
- auto lambda = [] __device__(lbm::KernelParameter parameter) {
- return lbm::cumulantChimera(parameter, lbm::setRelaxationRatesK17);
- };
-
- runKernel<<<cudaGrid.grid, cudaGrid.threads>>>(lambda, kernelParameter);
-
- getLastCudaError("K17CompressibleNavierStokesUnified execution failed");
-}
-
-
-}
-}
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesUnified.h b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesUnified.h
deleted file mode 100644
index 2b416a6d80d6bbb961152e6d3c559922e280b997..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesUnified.h
+++ /dev/null
@@ -1,23 +0,0 @@
-#ifndef K17CompressibleNavierStokesUnified_H
-#define K17CompressibleNavierStokesUnified_H
-
-#include "Kernel/KernelImp.h"
-
-namespace vf
-{
-namespace gpu
-{
-
-
-class K17CompressibleNavierStokesUnified : public KernelImp
-{
-public:
- K17CompressibleNavierStokesUnified(std::shared_ptr<Parameter> para, int level);
-
- void run();
-};
-
-}
-}
-
-#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cuh
deleted file mode 100644
index 49f144422bf5803582edebf9826f34b2af3f929a..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cuh
+++ /dev/null
@@ -1,29 +0,0 @@
-#ifndef K17CompressibleNavierStokes_Device_H
-#define K17CompressibleNavierStokes_Device_H
-
-#include <DataTypes.h>
-#include <curand.h>
-
-template< TurbulenceModel turbulenceModel, bool writeMacroscopicVariables, bool applyBodyForce > __global__ void K17CompressibleNavierStokes_Device(
- real omega_in,
- uint* neighborX,
- uint* neighborY,
- uint* neighborZ,
- real* distributions,
- real* rho,
- real* vx,
- real* vy,
- real* vz,
- real* turbulentViscosity,
- real SGSconstant,
- unsigned long long numberOfLBnodes,
- int level,
- real* forces,
- real* bodyForceX,
- real* bodyForceY,
- real* bodyForceZ,
- real* quadricLimiters,
- bool isEvenTimestep,
- const uint *fluidNodeIndices,
- uint numberOfFluidNodes);
-#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/RunLBMKernel.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/RunLBMKernel.cuh
deleted file mode 100644
index 3be594e3e39a57cd71741cd060e9dddda15d6035..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/RunLBMKernel.cuh
+++ /dev/null
@@ -1,64 +0,0 @@
-#ifndef GPU_CUMULANT_KERNEL_H
-#define GPU_CUMULANT_KERNEL_H
-
-
-#include <DataTypes.h>
-#include <cuda_runtime.h>
-
-#include "lbm/KernelParameter.h"
-
-#include "Kernel/Utilities/DistributionHelper.cuh"
-
-namespace vf
-{
-namespace gpu
-{
-
-
-struct GPUKernelParameter
-{
- real omega;
- unsigned int* typeOfGridNode;
- unsigned int* neighborX;
- unsigned int* neighborY;
- unsigned int* neighborZ;
- real* distributions;
- int numberOfLBnodes;
- real* forces;
- bool isEvenTimestep;
-};
-
-template<typename KernelFunctor>
-__global__ void runKernel(KernelFunctor kernel, GPUKernelParameter kernelParameter)
-{
- ////////////////////////////////////////////////////////////////////////////////
- //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
- //!
- const unsigned nodeIndex = getNodeIndex();
-
- if(nodeIndex >= kernelParameter.numberOfLBnodes)
- return;
-
- if (!isValidFluidNode(kernelParameter.typeOfGridNode[nodeIndex]))
- return;
-
- DistributionWrapper distributionWrapper {
- kernelParameter.distributions,
- (unsigned int)kernelParameter.numberOfLBnodes,
- kernelParameter.isEvenTimestep,
- nodeIndex,
- kernelParameter.neighborX,
- kernelParameter.neighborY,
- kernelParameter.neighborZ
- };
-
- lbm::KernelParameter parameter {distributionWrapper.distribution, kernelParameter.omega, kernelParameter.forces};
- kernel(parameter);
-
- distributionWrapper.write();
-}
-
-}
-}
-
-#endif
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/CheckParameterStrategy/CheckParameterStrategy.h b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/CheckParameterStrategy/CheckParameterStrategy.h
index 2fcd9a081df520af8d5517747b021a3f0353df6f..63e69c57cfd7b0d9138672a3a27de3c918b48297 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/CheckParameterStrategy/CheckParameterStrategy.h
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/CheckParameterStrategy/CheckParameterStrategy.h
@@ -10,6 +10,5 @@ class CheckParameterStrategy
public:
virtual ~CheckParameterStrategy() = default;
virtual bool checkParameter(std::shared_ptr<Parameter> para) = 0;
-
};
-#endif
\ No newline at end of file
+#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cu b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cu
deleted file mode 100644
index 3cc4b33f0b752943c97b1b770c62189a6c17d9c8..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cu
+++ /dev/null
@@ -1,82 +0,0 @@
-#include "DistributionHelper.cuh"
-
-#include <cuda_runtime.h>
-
-#include "basics/constants/NumericConstants.h"
-#include "lbm/constants/D3Q27.h"
-using namespace vf::lbm::dir;
-
-namespace vf::gpu
-{
-
-__device__ DistributionWrapper::DistributionWrapper(real *distributions, unsigned int size_Mat, bool isEvenTimestep,
- uint k, uint *neighborX, uint *neighborY, uint *neighborZ)
- : distribution_references(getDistributionReferences27(distributions, size_Mat, isEvenTimestep)), k(k), kw(neighborX[k]), ks(neighborY[k]),
- kb(neighborZ[k]), ksw(neighborY[kw]), kbw(neighborZ[kw]), kbs(neighborZ[ks]), kbsw(neighborZ[ksw])
-{
- read();
-}
-
-__device__ void DistributionWrapper::read()
-{
- distribution.f[vf::lbm::dir::DIR_P00] = (distribution_references.f[DIR_P00])[k];
- distribution.f[vf::lbm::dir::DIR_M00] = (distribution_references.f[DIR_M00])[kw];
- distribution.f[vf::lbm::dir::DIR_0P0] = (distribution_references.f[DIR_0P0])[k];
- distribution.f[vf::lbm::dir::DIR_0M0] = (distribution_references.f[DIR_0M0])[ks];
- distribution.f[vf::lbm::dir::DIR_00P] = (distribution_references.f[DIR_00P])[k];
- distribution.f[vf::lbm::dir::DIR_00M] = (distribution_references.f[DIR_00M])[kb];
- distribution.f[vf::lbm::dir::DIR_PP0] = (distribution_references.f[DIR_PP0])[k];
- distribution.f[vf::lbm::dir::DIR_MM0] = (distribution_references.f[DIR_MM0])[ksw];
- distribution.f[vf::lbm::dir::DIR_PM0] = (distribution_references.f[DIR_PM0])[ks];
- distribution.f[vf::lbm::dir::DIR_MP0] = (distribution_references.f[DIR_MP0])[kw];
- distribution.f[vf::lbm::dir::DIR_P0P] = (distribution_references.f[DIR_P0P])[k];
- distribution.f[vf::lbm::dir::DIR_M0M] = (distribution_references.f[DIR_M0M])[kbw];
- distribution.f[vf::lbm::dir::DIR_P0M] = (distribution_references.f[DIR_P0M])[kb];
- distribution.f[vf::lbm::dir::DIR_M0P] = (distribution_references.f[DIR_M0P])[kw];
- distribution.f[vf::lbm::dir::DIR_0PP] = (distribution_references.f[DIR_0PP])[k];
- distribution.f[vf::lbm::dir::DIR_0MM] = (distribution_references.f[DIR_0MM])[kbs];
- distribution.f[vf::lbm::dir::DIR_0PM] = (distribution_references.f[DIR_0PM])[kb];
- distribution.f[vf::lbm::dir::DIR_0MP] = (distribution_references.f[DIR_0MP])[ks];
- distribution.f[vf::lbm::dir::DIR_PPP] = (distribution_references.f[DIR_PPP])[k];
- distribution.f[vf::lbm::dir::DIR_MPP] = (distribution_references.f[DIR_MPP])[kw];
- distribution.f[vf::lbm::dir::DIR_PMP] = (distribution_references.f[DIR_PMP])[ks];
- distribution.f[vf::lbm::dir::DIR_MMP] = (distribution_references.f[DIR_MMP])[ksw];
- distribution.f[vf::lbm::dir::DIR_PPM] = (distribution_references.f[DIR_PPM])[kb];
- distribution.f[vf::lbm::dir::DIR_MPM] = (distribution_references.f[DIR_MPM])[kbw];
- distribution.f[vf::lbm::dir::DIR_PMM] = (distribution_references.f[DIR_PMM])[kbs];
- distribution.f[vf::lbm::dir::DIR_MMM] = (distribution_references.f[DIR_MMM])[kbsw];
- distribution.f[vf::lbm::dir::DIR_000] = (distribution_references.f[DIR_000])[k];
-}
-
-__device__ void DistributionWrapper::write()
-{
- (distribution_references.f[DIR_P00])[k] = distribution.f[vf::lbm::dir::PZZ];
- (distribution_references.f[DIR_M00])[kw] = distribution.f[vf::lbm::dir::MZZ];
- (distribution_references.f[DIR_0P0])[k] = distribution.f[vf::lbm::dir::ZPZ];
- (distribution_references.f[DIR_0M0])[ks] = distribution.f[vf::lbm::dir::ZMZ];
- (distribution_references.f[DIR_00P])[k] = distribution.f[vf::lbm::dir::ZZP];
- (distribution_references.f[DIR_00M])[kb] = distribution.f[vf::lbm::dir::ZZM];
- (distribution_references.f[DIR_PP0])[k] = distribution.f[vf::lbm::dir::PPZ];
- (distribution_references.f[DIR_MM0])[ksw] = distribution.f[vf::lbm::dir::MMZ];
- (distribution_references.f[DIR_PM0])[ks] = distribution.f[vf::lbm::dir::PMZ];
- (distribution_references.f[DIR_MP0])[kw] = distribution.f[vf::lbm::dir::MPZ];
- (distribution_references.f[DIR_P0P])[k] = distribution.f[vf::lbm::dir::PZP];
- (distribution_references.f[DIR_M0M])[kbw] = distribution.f[vf::lbm::dir::MZM];
- (distribution_references.f[DIR_P0M])[kb] = distribution.f[vf::lbm::dir::PZM];
- (distribution_references.f[DIR_M0P])[kw] = distribution.f[vf::lbm::dir::MZP];
- (distribution_references.f[DIR_0PP])[k] = distribution.f[vf::lbm::dir::ZPP];
- (distribution_references.f[DIR_0MM])[kbs] = distribution.f[vf::lbm::dir::ZMM];
- (distribution_references.f[DIR_0PM])[kb] = distribution.f[vf::lbm::dir::ZPM];
- (distribution_references.f[DIR_0MP])[ks] = distribution.f[vf::lbm::dir::ZMP];
- (distribution_references.f[DIR_PPP])[k] = distribution.f[vf::lbm::dir::PPP];
- (distribution_references.f[DIR_MPP])[kw] = distribution.f[vf::lbm::dir::MPP];
- (distribution_references.f[DIR_PMP])[ks] = distribution.f[vf::lbm::dir::PMP];
- (distribution_references.f[DIR_MMP])[ksw] = distribution.f[vf::lbm::dir::MMP];
- (distribution_references.f[DIR_PPM])[kb] = distribution.f[vf::lbm::dir::PPM];
- (distribution_references.f[DIR_MPM])[kbw] = distribution.f[vf::lbm::dir::MPM];
- (distribution_references.f[DIR_PMM])[kbs] = distribution.f[vf::lbm::dir::PMM];
- (distribution_references.f[DIR_MMM])[kbsw] = distribution.f[vf::lbm::dir::MMM];
- (distribution_references.f[DIR_000])[k] = distribution.f[vf::lbm::dir::ZZZ];
-}
-
-}
\ No newline at end of file
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh
deleted file mode 100644
index 599f3f46668c07da49725770177d77239f8ef9df..0000000000000000000000000000000000000000
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelper.cuh
+++ /dev/null
@@ -1,99 +0,0 @@
-//=======================================================================================
-// ____ ____ __ ______ __________ __ __ __ __
-// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
-// \ \ | | | | | |_) | | | | | | | / \ | |
-// \ \ | | | | | _ / | | | | | | / /\ \ | |
-// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
-// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
-// \ \ | | ________________________________________________________________
-// \ \ | | | ______________________________________________________________|
-// \ \| | | | __ __ __ __ ______ _______
-// \ | | |_____ | | | | | | | | | _ \ / _____)
-// \ | | _____| | | | | | | | | | | \ \ \_______
-// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
-// \ _____| |__| |________| \_______/ |__| |______/ (_______/
-//
-// This file is part of VirtualFluids. VirtualFluids is free software: you can
-// redistribute it and/or modify it under the terms of the GNU General Public
-// License as published by the Free Software Foundation, either version 3 of
-// the License, or (at your option) any later version.
-//
-// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
-// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
-// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
-// for more details.
-//
-// You should have received a copy of the GNU General Public License along
-// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
-//
-//! \file Cumulant27chim.cu
-//! \ingroup GPU
-//! \author Martin Schoenherr, Soeren Peters
-//=======================================================================================
-#ifndef DISTRIBUTUION_HELPER_CUH
-#define DISTRIBUTUION_HELPER_CUH
-
-#include "LBM/LB.h"
-
-#include "lbm/KernelParameter.h"
-#include "lbm/constants/D3Q27.h"
-#include "LBM/GPUHelperFunctions/KernelUtilities.h"
-
-using namespace vf::lbm::dir;
-
-namespace vf::gpu
-{
-
-/**
-* Getting references to the 27 directions.
-* @params distributions 1D real* array containing all data (number of elements = 27 * matrix_size)
-* @params matrix_size number of discretizations nodes
-* @params isEvenTimestep: stored data dependent on timestep is based on the esoteric twist algorithm
-* @return a data struct containing the addresses to the 27 directions within the 1D distribution array
-*/
-__inline__ __device__ __host__ DistributionReferences27 getDistributionReferences27(real* distributions, const unsigned long long numberOfLBnodes, const bool isEvenTimestep){
- DistributionReferences27 distribution_references;
- getPointersToDistributions(distribution_references, distributions, numberOfLBnodes, isEvenTimestep);
- return distribution_references;
-}
-
-
-/**
-* Holds the references to all directions and the concrete distributions for a single node.
-* After instantiation the distributions are read to the member "distribution" from "distribution_references".
-* After computation the data can be written back to "distribution_references".
-*/
-struct DistributionWrapper
-{
- __device__ DistributionWrapper(
- real* distributions,
- unsigned int size_Mat,
- bool isEvenTimestep,
- uint k,
- uint* neighborX,
- uint* neighborY,
- uint* neighborZ);
-
- __device__ void read();
-
- __device__ void write();
-
- // origin distributions to read from and write to after computation
- DistributionReferences27 distribution_references;
-
- // distribution pass to kernel computation
- vf::lbm::Distribution27 distribution;
-
- const uint k;
- const uint kw;
- const uint ks;
- const uint kb;
- const uint ksw;
- const uint kbw;
- const uint kbs;
- const uint kbsw;
-};
-
-}
-
-#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp
index 337456db79d372acef014d78e3588f7aa59fa74e..1ae8519670a93dcf6c3e3b1ca7a160e98167d27d 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp
+++ b/src/gpu/VirtualFluids_GPU/Kernel/Utilities/KernelFactory/KernelFactoryImp.cpp
@@ -8,7 +8,6 @@
//LBM kernel (compressible)
#include "Kernel/Compressible/FluidFlow/B92/B92CompressibleNavierStokes.h"
-#include "Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplusUnified.h"
#include "Kernel/Compressible/FluidFlow/B15/B15CompressibleNavierStokesBGKplus.h"
#include "Kernel/Compressible/FluidFlow/M02/M02CompressibleNavierStokes.h"
#include "Kernel/Compressible/FluidFlow/C06/C06CompressibleNavierStokes.h"
@@ -16,11 +15,9 @@
#include "Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokes.h"
#include "Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesBulkViscosity.h"
#include "Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesSponge.h"
-#include "Kernel/Compressible/FluidFlow/K15/K15CompressibleNavierStokesUnified.h"
#include "Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes.h"
#include "Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesChimeraLegacy.h"
#include "Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesBulkViscosity.h"
-#include "Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesUnified.h"
#include "Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokesSecondDerivatesFrom5thCumulants.h"
#include "Kernel/Compressible/FluidFlow/K18/K18CompressibleNavierStokes.h"
#include "Kernel/Compressible/FluidFlow/K20/K20CompressibleNavierStokes.h"
@@ -61,6 +58,8 @@
#include "Kernel/Kernels/PorousMediaKernels/FluidFlow/Compressible/PMFluidFlowCompStrategy.h"
#include "Kernel/Kernels/WaleKernels/FluidFlow/Compressible/WaleFluidFlowCompStrategy.h"
+#include <lbm/collision/TurbulentViscosity.h>
+
using namespace vf;
std::vector<std::shared_ptr<Kernel>> KernelFactoryImp::makeKernels(std::shared_ptr<Parameter> para)
@@ -103,9 +102,6 @@ std::shared_ptr<Kernel> KernelFactoryImp::makeKernel(std::shared_ptr<Parameter>
if (kernel == collisionKernel::compressible::BGK) {
newKernel = B92CompressibleNavierStokes::getNewInstance(para, level); // compressible
checkStrategy = FluidFlowCompStrategy::getInstance(); // ||
- } else if (kernel == collisionKernel::compressible::BGKUnified) { // \/
- newKernel = std::make_shared<vf::gpu::B15CompressibleNavierStokesBGKplusUnified>(para, level);
- checkStrategy = FluidFlowCompStrategy::getInstance();
} else if (kernel == collisionKernel::compressible::BGKPlus) {
newKernel = B15CompressibleNavierStokesBGKplus::getNewInstance(para, level);
checkStrategy = FluidFlowCompStrategy::getInstance();
@@ -118,12 +114,6 @@ std::shared_ptr<Kernel> KernelFactoryImp::makeKernel(std::shared_ptr<Parameter>
} else if (kernel == collisionKernel::compressible::CumulantClassic) {
newKernel = K08CompressibleNavierStokes::getNewInstance(para, level);
checkStrategy = FluidFlowCompStrategy::getInstance();
- } else if (kernel == collisionKernel::compressible::CumulantK15Unified) {
- newKernel = std::make_shared<vf::gpu::K15CompressibleNavierStokesUnified>(para, level);
- checkStrategy = FluidFlowCompStrategy::getInstance();
- } else if (kernel == collisionKernel::compressible::K17CompressibleNavierStokesUnified) {
- newKernel = std::make_shared<vf::gpu::K17CompressibleNavierStokesUnified>(para, level);
- checkStrategy = FluidFlowCompStrategy::getInstance();
} else if (kernel == collisionKernel::compressible::K17CompressibleNavierStokesBulkViscosity) {
newKernel = K17CompressibleNavierStokesBulkViscosity::getNewInstance(para, level);
checkStrategy = FluidFlowCompStrategy::getInstance();
@@ -133,17 +123,17 @@ std::shared_ptr<Kernel> KernelFactoryImp::makeKernel(std::shared_ptr<Parameter>
} else if (kernel == collisionKernel::compressible::K17CompressibleNavierStokes){
switch(para->getTurbulenceModel())
{
- case TurbulenceModel::AMD:
- newKernel = K17CompressibleNavierStokes<TurbulenceModel::AMD>::getNewInstance(para, level);
+ case lbm::TurbulenceModel::AMD:
+ newKernel = K17CompressibleNavierStokes<lbm::TurbulenceModel::AMD>::getNewInstance(para, level);
break;
- case TurbulenceModel::Smagorinsky:
- newKernel = K17CompressibleNavierStokes<TurbulenceModel::Smagorinsky>::getNewInstance(para, level);
+ case lbm::TurbulenceModel::Smagorinsky:
+ newKernel = K17CompressibleNavierStokes<lbm::TurbulenceModel::Smagorinsky>::getNewInstance(para, level);
break;
- case TurbulenceModel::QR:
- newKernel = K17CompressibleNavierStokes<TurbulenceModel::QR>::getNewInstance(para, level);
+ case lbm::TurbulenceModel::QR:
+ newKernel = K17CompressibleNavierStokes<lbm::TurbulenceModel::QR>::getNewInstance(para, level);
break;
- case TurbulenceModel::None:
- newKernel = K17CompressibleNavierStokes<TurbulenceModel::None>::getNewInstance(para, level);
+ case lbm::TurbulenceModel::None:
+ newKernel = K17CompressibleNavierStokes<lbm::TurbulenceModel::None>::getNewInstance(para, level);
break;
default:
throw std::runtime_error("Unknown turbulence model!");
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h
index ebd1b8ab367b31c62fc1608b39215622819d35d9..da13a0db463acb23f032779b6b40caed3310d184 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilities.h
@@ -30,14 +30,19 @@
//! \ingroup LBM/GPUHelperFunctions
//! \author Martin Schoenherr, Anna Wellmann, Soeren Peters
//=======================================================================================
-#ifndef KERNEL_UTILITIES_H
-#define KERNEL_UTILITIES_H
+#ifndef GPU_DISTRIBUTION_HELPER_H
+#define GPU_DISTRIBUTION_HELPER_H
#include "LBM/LB.h"
+#include <array>
#include <cassert>
+#include <optional>
+
+#include <cuda_runtime.h>
#include <lbm/constants/D3Q27.h>
+
#include <basics/constants/NumericConstants.h>
using namespace vf::basics::constant;
@@ -46,10 +51,19 @@ using namespace vf::lbm::dir;
namespace vf::gpu
{
+inline real getForceFactor(int level)
+{
+ real factor = c1o1;
+ for (size_t i = 1; i <= level; i++) {
+ factor *= c2o1;
+ }
+ factor = 1 / factor;
+ return factor;
+}
+
__inline__ __device__ __host__ void getPointersToDistributions(Distributions27 &dist, real *distributionArray, const unsigned long long numberOfLBnodes, const bool isEvenTimestep)
{
- if (isEvenTimestep)
- {
+ if (isEvenTimestep) {
dist.f[DIR_000] = &distributionArray[DIR_000 * numberOfLBnodes];
dist.f[DIR_P00] = &distributionArray[DIR_P00 * numberOfLBnodes];
dist.f[DIR_M00] = &distributionArray[DIR_M00 * numberOfLBnodes];
@@ -77,37 +91,48 @@ __inline__ __device__ __host__ void getPointersToDistributions(Distributions27 &
dist.f[DIR_MMM] = &distributionArray[DIR_MMM * numberOfLBnodes];
dist.f[DIR_PMM] = &distributionArray[DIR_PMM * numberOfLBnodes];
dist.f[DIR_MPM] = &distributionArray[DIR_MPM * numberOfLBnodes];
+ } else {
+ dist.f[DIR_M00] = &distributionArray[DIR_P00 * numberOfLBnodes];
+ dist.f[DIR_P00] = &distributionArray[DIR_M00 * numberOfLBnodes];
+ dist.f[DIR_0M0] = &distributionArray[DIR_0P0 * numberOfLBnodes];
+ dist.f[DIR_0P0] = &distributionArray[DIR_0M0 * numberOfLBnodes];
+ dist.f[DIR_00M] = &distributionArray[DIR_00P * numberOfLBnodes];
+ dist.f[DIR_00P] = &distributionArray[DIR_00M * numberOfLBnodes];
+ dist.f[DIR_MM0] = &distributionArray[DIR_PP0 * numberOfLBnodes];
+ dist.f[DIR_PP0] = &distributionArray[DIR_MM0 * numberOfLBnodes];
+ dist.f[DIR_MP0] = &distributionArray[DIR_PM0 * numberOfLBnodes];
+ dist.f[DIR_PM0] = &distributionArray[DIR_MP0 * numberOfLBnodes];
+ dist.f[DIR_M0M] = &distributionArray[DIR_P0P * numberOfLBnodes];
+ dist.f[DIR_P0P] = &distributionArray[DIR_M0M * numberOfLBnodes];
+ dist.f[DIR_M0P] = &distributionArray[DIR_P0M * numberOfLBnodes];
+ dist.f[DIR_P0M] = &distributionArray[DIR_M0P * numberOfLBnodes];
+ dist.f[DIR_0MM] = &distributionArray[DIR_0PP * numberOfLBnodes];
+ dist.f[DIR_0PP] = &distributionArray[DIR_0MM * numberOfLBnodes];
+ dist.f[DIR_0MP] = &distributionArray[DIR_0PM * numberOfLBnodes];
+ dist.f[DIR_0PM] = &distributionArray[DIR_0MP * numberOfLBnodes];
+ dist.f[DIR_000] = &distributionArray[DIR_000 * numberOfLBnodes];
+ dist.f[DIR_PPP] = &distributionArray[DIR_MMM * numberOfLBnodes];
+ dist.f[DIR_MMP] = &distributionArray[DIR_PPM * numberOfLBnodes];
+ dist.f[DIR_PMP] = &distributionArray[DIR_MPM * numberOfLBnodes];
+ dist.f[DIR_MPP] = &distributionArray[DIR_PMM * numberOfLBnodes];
+ dist.f[DIR_PPM] = &distributionArray[DIR_MMP * numberOfLBnodes];
+ dist.f[DIR_MMM] = &distributionArray[DIR_PPP * numberOfLBnodes];
+ dist.f[DIR_PMM] = &distributionArray[DIR_MPP * numberOfLBnodes];
+ dist.f[DIR_MPM] = &distributionArray[DIR_PMP * numberOfLBnodes];
}
- else
- {
- dist.f[DIR_M00] = &distributionArray[DIR_P00 * numberOfLBnodes];
- dist.f[DIR_P00] = &distributionArray[DIR_M00 * numberOfLBnodes];
- dist.f[DIR_0M0] = &distributionArray[DIR_0P0 * numberOfLBnodes];
- dist.f[DIR_0P0] = &distributionArray[DIR_0M0 * numberOfLBnodes];
- dist.f[DIR_00M] = &distributionArray[DIR_00P * numberOfLBnodes];
- dist.f[DIR_00P] = &distributionArray[DIR_00M * numberOfLBnodes];
- dist.f[DIR_MM0] = &distributionArray[DIR_PP0 * numberOfLBnodes];
- dist.f[DIR_PP0] = &distributionArray[DIR_MM0 * numberOfLBnodes];
- dist.f[DIR_MP0] = &distributionArray[DIR_PM0 * numberOfLBnodes];
- dist.f[DIR_PM0] = &distributionArray[DIR_MP0 * numberOfLBnodes];
- dist.f[DIR_M0M] = &distributionArray[DIR_P0P * numberOfLBnodes];
- dist.f[DIR_P0P] = &distributionArray[DIR_M0M * numberOfLBnodes];
- dist.f[DIR_M0P] = &distributionArray[DIR_P0M * numberOfLBnodes];
- dist.f[DIR_P0M] = &distributionArray[DIR_M0P * numberOfLBnodes];
- dist.f[DIR_0MM] = &distributionArray[DIR_0PP * numberOfLBnodes];
- dist.f[DIR_0PP] = &distributionArray[DIR_0MM * numberOfLBnodes];
- dist.f[DIR_0MP] = &distributionArray[DIR_0PM * numberOfLBnodes];
- dist.f[DIR_0PM] = &distributionArray[DIR_0MP * numberOfLBnodes];
- dist.f[DIR_000] = &distributionArray[DIR_000 * numberOfLBnodes];
- dist.f[DIR_PPP] = &distributionArray[DIR_MMM * numberOfLBnodes];
- dist.f[DIR_MMP] = &distributionArray[DIR_PPM * numberOfLBnodes];
- dist.f[DIR_PMP] = &distributionArray[DIR_MPM * numberOfLBnodes];
- dist.f[DIR_MPP] = &distributionArray[DIR_PMM * numberOfLBnodes];
- dist.f[DIR_PPM] = &distributionArray[DIR_MMP * numberOfLBnodes];
- dist.f[DIR_MMM] = &distributionArray[DIR_PPP * numberOfLBnodes];
- dist.f[DIR_PMM] = &distributionArray[DIR_MPP * numberOfLBnodes];
- dist.f[DIR_MPM] = &distributionArray[DIR_PMP * numberOfLBnodes];
- }
+}
+
+/**
+* Getting references to the 27 directions.
+* @params distributions 1D real* array containing all data (number of elements = 27 * matrix_size)
+* @params matrix_size number of discretizations nodes
+* @params isEvenTimestep: stored data dependent on timestep is based on the esoteric twist algorithm
+* @return a data struct containing the addresses to the 27 directions within the 1D distribution array
+*/
+__inline__ __device__ __host__ DistributionReferences27 getDistributionReferences27(real* distributions, const unsigned long long numberOfLBnodes, const bool isEvenTimestep){
+ DistributionReferences27 distribution_references;
+ getPointersToDistributions(distribution_references, distributions, numberOfLBnodes, isEvenTimestep);
+ return distribution_references;
}
__inline__ __device__ void getPointersToSubgridDistances(SubgridDistances27& subgridD, real* subgridDistances, const unsigned int numberOfSubgridIndices)
@@ -268,8 +293,33 @@ __device__ __inline__ void read(real *f, const Distributions27 &dist, const List
__device__ __inline__ void readInverse(real *f, const Distributions27 &dist, const ListIndices &indices)
{
- //TODO: https://git.rz.tu-bs.de/irmb/VirtualFluids_dev/-/issues/101
- assert((false) || !fprintf(stderr, "Not implemented yet.\n"));
+ f[DIR_000] = (dist.f[DIR_000])[indices.k_000];
+ f[DIR_P00] = (dist.f[DIR_P00])[indices.k_000];
+ f[DIR_M00] = (dist.f[DIR_M00])[indices.k_M00];
+ f[DIR_0P0] = (dist.f[DIR_0P0])[indices.k_000];
+ f[DIR_0M0] = (dist.f[DIR_0M0])[indices.k_0M0];
+ f[DIR_00P] = (dist.f[DIR_00P])[indices.k_000];
+ f[DIR_00M] = (dist.f[DIR_00M])[indices.k_00M];
+ f[DIR_PP0] = (dist.f[DIR_PP0])[indices.k_000];
+ f[DIR_MM0] = (dist.f[DIR_MM0])[indices.k_MM0];
+ f[DIR_PM0] = (dist.f[DIR_PM0])[indices.k_0M0];
+ f[DIR_MP0] = (dist.f[DIR_MP0])[indices.k_M00];
+ f[DIR_P0P] = (dist.f[DIR_P0P])[indices.k_000];
+ f[DIR_M0M] = (dist.f[DIR_M0M])[indices.k_M0M];
+ f[DIR_P0M] = (dist.f[DIR_P0M])[indices.k_00M];
+ f[DIR_M0P] = (dist.f[DIR_M0P])[indices.k_M00];
+ f[DIR_0PP] = (dist.f[DIR_0PP])[indices.k_000];
+ f[DIR_0MM] = (dist.f[DIR_0MM])[indices.k_0MM];
+ f[DIR_0PM] = (dist.f[DIR_0PM])[indices.k_00M];
+ f[DIR_0MP] = (dist.f[DIR_0MP])[indices.k_0M0];
+ f[DIR_PPP] = (dist.f[DIR_PPP])[indices.k_000];
+ f[DIR_MPP] = (dist.f[DIR_MPP])[indices.k_M00];
+ f[DIR_PMP] = (dist.f[DIR_PMP])[indices.k_0M0];
+ f[DIR_MMP] = (dist.f[DIR_MMP])[indices.k_MM0];
+ f[DIR_PPM] = (dist.f[DIR_PPM])[indices.k_00M];
+ f[DIR_MPM] = (dist.f[DIR_MPM])[indices.k_M0M];
+ f[DIR_PMM] = (dist.f[DIR_PMM])[indices.k_0MM];
+ f[DIR_MMM] = (dist.f[DIR_MMM])[indices.k_MMM];
}
@@ -278,7 +328,7 @@ __device__ __inline__ void readInverse(real *f, const Distributions27 &dist, con
//! stored arrays dependent on timestep is based on the esoteric twist algorithm
//! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017),
//! DOI:10.3390/computation5020019 ]</b></a>
-__inline__ __device__ void write(Distributions27 &destination, const ListIndices &indices, const real* f)
+__inline__ __device__ void write(Distributions27& destination, const ListIndices& indices, const real* f)
{
(destination.f[DIR_000])[indices.k_000] = f[DIR_000];
(destination.f[DIR_P00])[indices.k_000] = f[DIR_P00];
@@ -309,15 +359,37 @@ __inline__ __device__ void write(Distributions27 &destination, const ListIndices
(destination.f[DIR_MMM])[indices.k_MMM] = f[DIR_MMM];
}
-__inline__ __device__ void writeInverse(Distributions27 &destination, const ListIndices &indices, const real* f)
+__inline__ __device__ void writeInverse(Distributions27& destination, const ListIndices& indices, const real* f)
{
- //TODO: https://git.rz.tu-bs.de/irmb/VirtualFluids_dev/-/issues/101
- assert((false) || !fprintf(stderr, "Not implemented yet.\n"));
+ (destination.f[DIR_000])[indices.k_000] = f[DIR_000];
+ (destination.f[DIR_P00])[indices.k_000] = f[DIR_M00];
+ (destination.f[DIR_M00])[indices.k_M00] = f[DIR_P00];
+ (destination.f[DIR_0P0])[indices.k_000] = f[DIR_0M0];
+ (destination.f[DIR_0M0])[indices.k_0M0] = f[DIR_0P0];
+ (destination.f[DIR_00P])[indices.k_000] = f[DIR_00M];
+ (destination.f[DIR_00M])[indices.k_00M] = f[DIR_00P];
+ (destination.f[DIR_PP0])[indices.k_000] = f[DIR_MM0];
+ (destination.f[DIR_MM0])[indices.k_MM0] = f[DIR_PP0];
+ (destination.f[DIR_PM0])[indices.k_0M0] = f[DIR_MP0];
+ (destination.f[DIR_MP0])[indices.k_M00] = f[DIR_PM0];
+ (destination.f[DIR_P0P])[indices.k_000] = f[DIR_M0M];
+ (destination.f[DIR_M0M])[indices.k_M0M] = f[DIR_P0P];
+ (destination.f[DIR_P0M])[indices.k_00M] = f[DIR_M0P];
+ (destination.f[DIR_M0P])[indices.k_M00] = f[DIR_P0M];
+ (destination.f[DIR_0PP])[indices.k_000] = f[DIR_0MM];
+ (destination.f[DIR_0MM])[indices.k_0MM] = f[DIR_0PP];
+ (destination.f[DIR_0PM])[indices.k_00M] = f[DIR_0MP];
+ (destination.f[DIR_0MP])[indices.k_0M0] = f[DIR_0PM];
+ (destination.f[DIR_PPP])[indices.k_000] = f[DIR_MMM];
+ (destination.f[DIR_MPP])[indices.k_M00] = f[DIR_PMM];
+ (destination.f[DIR_PMP])[indices.k_0M0] = f[DIR_MPM];
+ (destination.f[DIR_MMP])[indices.k_MM0] = f[DIR_PPM];
+ (destination.f[DIR_PPM])[indices.k_00M] = f[DIR_MMP];
+ (destination.f[DIR_MPM])[indices.k_M0M] = f[DIR_PMP];
+ (destination.f[DIR_PMM])[indices.k_0MM] = f[DIR_MPP];
+ (destination.f[DIR_MMM])[indices.k_MMM] = f[DIR_PPP];
}
-
-
-
__inline__ __device__ real calculateOmega(const real omega_old, real turbulenceViscosity)
{
return omega_old / (c1o1 + c3o1 * omega_old * turbulenceViscosity);
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelperTests.cpp b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilitiesTests.cpp
similarity index 98%
rename from src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelperTests.cpp
rename to src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilitiesTests.cpp
index 9d6529404284f786f85a7445e494487fe75d8ad7..87f7979bedd8995b075c2adf164a68c1a5b95022 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Utilities/DistributionHelperTests.cpp
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/KernelUtilitiesTests.cpp
@@ -1,9 +1,11 @@
#include <gmock/gmock.h>
-#include "basics/tests/testUtilities.h"
-#include "DistributionHelper.cuh"
+#include <basics/tests/testUtilities.h>
+
+#include "KernelUtilities.h"
+
+#include <lbm/constants/D3Q27.h>
-#include "lbm/constants/D3Q27.h"
using namespace vf::lbm::dir;
TEST(DistributionHelperTests, getPointerToDistribution_WhenEvenTimeStep_ShouldBeEqualToInput)
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/RunCollision.cuh b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/RunCollision.cuh
new file mode 100644
index 0000000000000000000000000000000000000000..ba1034ac59b52522c8afd881d9f91f90e701d5ae
--- /dev/null
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/RunCollision.cuh
@@ -0,0 +1,114 @@
+#ifndef GPU_RUNCOLLISIONKERNEL_CUH
+#define GPU_RUNCOLLISIONKERNEL_CUH
+
+#include <cuda_runtime.h>
+
+#include <basics/DataTypes.h>
+
+#include <lbm/collision/CollisionParameter.h>
+#include <lbm/collision/TurbulentViscosity.h>
+
+#include "LBM/GPUHelperFunctions/KernelUtilities.h"
+
+namespace vf::gpu
+{
+
+struct GPUCollisionParameter
+{
+ real omega;
+ unsigned int* neighborX;
+ unsigned int* neighborY;
+ unsigned int* neighborZ;
+ real* distributions;
+ real* rho;
+ real* vx;
+ real* vy;
+ real* vz;
+ real* turbulentViscosity;
+ real SGSconstant;
+ int numberOfLBnodes;
+ real forceFactor;
+ real* forces;
+ real* bodyForceX;
+ real* bodyForceY;
+ real* bodyForceZ;
+ real* quadricLimiters;
+ bool isEvenTimestep;
+ const uint* fluidNodeIndices;
+ uint numberOfFluidNodes;
+};
+
+template <typename CollisionFunctor, vf::lbm::TurbulenceModel turbulenceModel, bool writeMacroscopicVariables, bool applyBodyForce>
+__global__ void runCollision(CollisionFunctor collision, GPUCollisionParameter collisionParameter)
+{
+ const unsigned nodeIndex = getNodeIndex();
+
+ if (nodeIndex >= collisionParameter.numberOfFluidNodes)
+ return;
+
+ const unsigned k_000 = collisionParameter.fluidNodeIndices[nodeIndex];
+
+ vf::lbm::CollisionParameter para;
+ para.omega = collisionParameter.omega;
+ para.quadricLimiter = collisionParameter.quadricLimiters;
+
+ if (applyBodyForce) {
+ para.forceX = (collisionParameter.forces[0] + collisionParameter.bodyForceX[k_000]) * c1o2 * collisionParameter.forceFactor;
+ para.forceY = (collisionParameter.forces[1] + collisionParameter.bodyForceY[k_000]) * c1o2 * collisionParameter.forceFactor;
+ para.forceZ = (collisionParameter.forces[2] + collisionParameter.bodyForceZ[k_000]) * c1o2 * collisionParameter.forceFactor;
+
+ // Reset body force. To be used when not using round-off correction.
+ collisionParameter.bodyForceX[k_000] = 0.0f;
+ collisionParameter.bodyForceX[k_000] = 0.0f;
+ collisionParameter.bodyForceX[k_000] = 0.0f;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //!> Round-off correction
+ //!
+ //!> Similar to Kahan summation algorithm (https://en.wikipedia.org/wiki/Kahan_summation_algorithm)
+ //!> Essentially computes the round-off error of the applied force and adds it in the next time step as a compensation.
+ //!> Seems to be necesseary at very high Re boundary layers, where the forcing and velocity can
+ //!> differ by several orders of magnitude.
+ //!> \note 16/05/2022: Testing, still ongoing!
+ //!
+ // bodyForceX[k_000] = (acc_x-(vvx-vx))*factor*c2o1;
+ // bodyForceY[k_000] = (acc_y-(vvy-vy))*factor*c2o1;
+ // bodyForceZ[k_000] = (acc_z-(vvz-vz))*factor*c2o1;
+ } else {
+ para.forceX = collisionParameter.forces[0] * c1o2 * collisionParameter.forceFactor;
+ para.forceY = collisionParameter.forces[1] * c1o2 * collisionParameter.forceFactor;
+ para.forceZ = collisionParameter.forces[2] * c1o2 * collisionParameter.forceFactor;
+ }
+
+ vf::lbm::TurbulentViscosity turbulentViscosity;
+ if (turbulenceModel != vf::lbm::TurbulenceModel::None) {
+ turbulentViscosity.value = collisionParameter.turbulentViscosity[k_000];
+ turbulentViscosity.SGSconstant = collisionParameter.SGSconstant;
+ }
+
+ Distributions27 dist;
+ getPointersToDistributions(dist, collisionParameter.distributions, collisionParameter.numberOfLBnodes,
+ collisionParameter.isEvenTimestep);
+
+ ListIndices listIndices(k_000, collisionParameter.neighborX, collisionParameter.neighborY, collisionParameter.neighborZ);
+
+ read(para.distribution, dist, listIndices);
+
+ vf::lbm::MacroscopicValues macroscopicValues;
+ collision(para, macroscopicValues, turbulentViscosity);
+
+ if (writeMacroscopicVariables) {
+ collisionParameter.vx[k_000] = macroscopicValues.vx;
+ collisionParameter.vy[k_000] = macroscopicValues.vy;
+ collisionParameter.vz[k_000] = macroscopicValues.vz;
+ collisionParameter.rho[k_000] = macroscopicValues.rho;
+ }
+ if (turbulenceModel != vf::lbm::TurbulenceModel::None)
+ collisionParameter.turbulentViscosity[k_000] = turbulentViscosity.value;
+
+ writeInverse(dist, listIndices, para.distribution);
+}
+
+} // namespace vf::gpu
+
+#endif
diff --git a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
index c825ae108dac5aed46b6741606026b0ba9e52a08..2d71eaf428803147056e1ea5ee03758f2b1d11d8 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/GPUHelperFunctions/ScalingUtilities.h
@@ -38,7 +38,6 @@
#include <basics/constants/NumericConstants.h>
#include <lbm/constants/D3Q27.h>
-#include <lbm/KernelParameter.h>
#include <basics/DataTypes.h>
#include <lbm/interpolation/InterpolationCoefficients.h>
diff --git a/src/gpu/VirtualFluids_GPU/LBM/LB.h b/src/gpu/VirtualFluids_GPU/LBM/LB.h
index 5da2a2a82c59d9a70e26849d4ee413b7c6e8f1d3..979e08d87e170207451e7ba80304cca8d45d1dd4 100644
--- a/src/gpu/VirtualFluids_GPU/LBM/LB.h
+++ b/src/gpu/VirtualFluids_GPU/LBM/LB.h
@@ -75,19 +75,6 @@
#include <string>
#include <vector>
-//! \brief An enumeration for selecting a turbulence model
-enum class TurbulenceModel {
- //! - Smagorinsky
- Smagorinsky,
- //! - AMD (Anisotropic Minimum Dissipation) model, see e.g. Rozema et al., Phys. Fluids 27, 085107 (2015), https://doi.org/10.1063/1.4928700
- AMD,
- //! - QR model by Verstappen
- QR,
- //! - TODO: move the WALE model here from the old kernels
- //WALE
- //! - No turbulence model
- None
-};
//! \brief An enumeration for selecting a template of the collision kernel (CumulantK17)
enum class CollisionTemplate {
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
index ede6344270cf1f119ec30b699ae2c72b9778dc81..d11088f041251d1048a1548f4d8a04f4c00c1feb 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.cpp
@@ -1008,7 +1008,7 @@ void Parameter::setUseWale(bool useWale)
if (useWale)
setUseTurbulentViscosity(true);
}
-void Parameter::setTurbulenceModel(TurbulenceModel turbulenceModel)
+void Parameter::setTurbulenceModel(vf::lbm::TurbulenceModel turbulenceModel)
{
this->turbulenceModel = turbulenceModel;
}
@@ -2383,7 +2383,7 @@ bool Parameter::getUseWale()
{
return this->isWale;
}
-TurbulenceModel Parameter::getTurbulenceModel()
+vf::lbm::TurbulenceModel Parameter::getTurbulenceModel()
{
return this->turbulenceModel;
}
diff --git a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
index 9140f43fa5487d6db9959fd299b836272f1ac4b2..78d00850aee5f9ce0d6800d8c4301efb3348052c 100644
--- a/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
+++ b/src/gpu/VirtualFluids_GPU/Parameter/Parameter.h
@@ -45,6 +45,8 @@
#include "TurbulenceModels/TurbulenceModelFactory.h"
#include "VirtualFluids_GPU/Kernel/Utilities/KernelTypes.h"
+#include <lbm/collision/TurbulentViscosity.h>
+
#include "VirtualFluids_GPU_export.h"
struct curandStateXORWOW;
@@ -602,7 +604,7 @@ public:
void setConcFile(bool concFile);
void setUseMeasurePoints(bool useMeasurePoints);
void setUseWale(bool useWale);
- void setTurbulenceModel(TurbulenceModel turbulenceModel);
+ void setTurbulenceModel(vf::lbm::TurbulenceModel turbulenceModel);
void setUseTurbulentViscosity(bool useTurbulentViscosity);
void setSGSConstant(real SGSConstant);
void setHasWallModelMonitor(bool hasWallModelMonitor);
@@ -896,7 +898,7 @@ public:
bool getConcFile();
bool getUseMeasurePoints();
bool getUseWale();
- TurbulenceModel getTurbulenceModel();
+ vf::lbm::TurbulenceModel getTurbulenceModel();
bool getUseTurbulentViscosity();
real getSGSConstant();
bool getHasWallModelMonitor();
@@ -1082,7 +1084,7 @@ private:
std::string concentration;
std::string geomNormalX, geomNormalY, geomNormalZ, inflowNormalX, inflowNormalY, inflowNormalZ, outflowNormalX, outflowNormalY, outflowNormalZ;
- TurbulenceModel turbulenceModel{ TurbulenceModel::None };
+ vf::lbm::TurbulenceModel turbulenceModel{ vf::lbm::TurbulenceModel::None };
// Kernel
diff --git a/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.cpp b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.cpp
index 980a199b03710ac91310d26424c20c7371f1097d..b17c05890ea440a3c6bad39e4770bbf4ddd73c96 100644
--- a/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.cpp
+++ b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.cpp
@@ -38,14 +38,14 @@
#include <variant>
-void TurbulenceModelFactory::setTurbulenceModel(TurbulenceModel _turbulenceModel)
+void TurbulenceModelFactory::setTurbulenceModel(vf::lbm::TurbulenceModel _turbulenceModel)
{
this->turbulenceModel = _turbulenceModel;
para->setTurbulenceModel(_turbulenceModel);
- if(this->turbulenceModel != TurbulenceModel::None) para->setUseTurbulentViscosity(true);
+ if(this->turbulenceModel != vf::lbm::TurbulenceModel::None) para->setUseTurbulentViscosity(true);
switch (this->turbulenceModel) {
- case TurbulenceModel::AMD:
+ case vf::lbm::TurbulenceModel::AMD:
this->turbulenceModelKernel = calcTurbulentViscosityAMD;
break;
default:
@@ -64,10 +64,10 @@ void TurbulenceModelFactory::readConfigFile(const vf::basics::ConfigurationFile
{
std::string config = configData.getValue<std::string>("TurbulenceModel");
- if (config == "Smagorinsky") this->setTurbulenceModel( TurbulenceModel::Smagorinsky );
- else if (config == "AMD") this->setTurbulenceModel( TurbulenceModel::AMD );
- else if (config == "QR" ) this->setTurbulenceModel( TurbulenceModel::QR );
- else if (config == "None") this->setTurbulenceModel( TurbulenceModel::None );
+ if (config == "Smagorinsky") this->setTurbulenceModel( vf::lbm::TurbulenceModel::Smagorinsky );
+ else if (config == "AMD") this->setTurbulenceModel( vf::lbm::TurbulenceModel::AMD );
+ else if (config == "QR" ) this->setTurbulenceModel( vf::lbm::TurbulenceModel::QR );
+ else if (config == "None") this->setTurbulenceModel( vf::lbm::TurbulenceModel::None );
else std::runtime_error("TurbulenceModelFactory: Invalid turbulence model!");
VF_LOG_INFO("Turbulence model: {}", config);
diff --git a/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h
index e71c8ed5f7be016a4a800b83cfb3252ee6b8246e..63ab3570734a4f9967fb0a8ad60aa811a7f50515 100644
--- a/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h
+++ b/src/gpu/VirtualFluids_GPU/TurbulenceModels/TurbulenceModelFactory.h
@@ -37,12 +37,13 @@
#include <map>
#include <string>
#include <variant>
+#include <memory>
#include "LBM/LB.h"
-#include "Parameter/Parameter.h"
-
#include <basics/config/ConfigurationFile.h>
+#include <lbm/collision/TurbulentViscosity.h>
+
class Parameter;
using TurbulenceModelKernel = std::function<void(Parameter *, int )>;
@@ -51,9 +52,9 @@ class TurbulenceModelFactory
{
public:
- TurbulenceModelFactory(SPtr<Parameter> parameter): para(parameter) {}
+ TurbulenceModelFactory(std::shared_ptr<Parameter> parameter): para(parameter) {}
- void setTurbulenceModel(TurbulenceModel _turbulenceModel);
+ void setTurbulenceModel(vf::lbm::TurbulenceModel _turbulenceModel);
void setModelConstant(real modelConstant);
@@ -62,9 +63,9 @@ public:
void runTurbulenceModelKernel(const int level) const;
private:
- TurbulenceModel turbulenceModel = TurbulenceModel::None;
+ vf::lbm::TurbulenceModel turbulenceModel = vf::lbm::TurbulenceModel::None;
TurbulenceModelKernel turbulenceModelKernel = nullptr;
- SPtr<Parameter> para;
+ std::shared_ptr<Parameter> para;
};
diff --git a/src/lbm/BGK.cpp b/src/lbm/BGK.cpp
deleted file mode 100644
index b59572293e3ffd8f15e8edc864b2993e48b00296..0000000000000000000000000000000000000000
--- a/src/lbm/BGK.cpp
+++ /dev/null
@@ -1,135 +0,0 @@
-#include "BGK.h"
-
-
-#include <basics/DataTypes.h>
-#include <basics/constants/NumericConstants.h>
-
-#include "constants/D3Q27.h"
-#include "MacroscopicQuantities.h"
-
-namespace vf::lbm
-{
-
-using namespace vf::basics::constant;
-
-
-
-__host__ __device__ void bgk(KernelParameter parameter)
-{
- auto& distribution = parameter.distribution;
- const auto omega = parameter.omega;
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Read distributions: style of reading and writing the distributions from/to
- //! stored arrays dependent on timestep is based on the esoteric twist algorithm
- //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
- //!
- real mfcbb = distribution.f[dir::PZZ];
- real mfabb = distribution.f[dir::MZZ];
- real mfbcb = distribution.f[dir::ZPZ];
- real mfbab = distribution.f[dir::ZMZ];
- real mfbbc = distribution.f[dir::ZZP];
- real mfbba = distribution.f[dir::ZZM];
- real mfccb = distribution.f[dir::PPZ];
- real mfaab = distribution.f[dir::MMZ];
- real mfcab = distribution.f[dir::PMZ];
- real mfacb = distribution.f[dir::MPZ];
- real mfcbc = distribution.f[dir::PZP];
- real mfaba = distribution.f[dir::MZM];
- real mfcba = distribution.f[dir::PZM];
- real mfabc = distribution.f[dir::MZP];
- real mfbcc = distribution.f[dir::ZPP];
- real mfbaa = distribution.f[dir::ZMM];
- real mfbca = distribution.f[dir::ZPM];
- real mfbac = distribution.f[dir::ZMP];
- real mfccc = distribution.f[dir::PPP];
- real mfacc = distribution.f[dir::MPP];
- real mfcac = distribution.f[dir::PMP];
- real mfaac = distribution.f[dir::MMP];
- real mfcca = distribution.f[dir::PPM];
- real mfaca = distribution.f[dir::MPM];
- real mfcaa = distribution.f[dir::PMM];
- real mfaaa = distribution.f[dir::MMM];
- real mfbbb = distribution.f[dir::ZZZ];
-
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Acquire macroscopic quantities
- const real drho = getDensity(distribution.f);
- const real rho = c1o1 + drho;
- const real OOrho = c1o1 / (c1o1 + drho);
-
- const real vvx = getIncompressibleVelocityX1(distribution.f) * OOrho;
- const real vvy = getIncompressibleVelocityX2(distribution.f) * OOrho;
- const real vvz = getIncompressibleVelocityX3(distribution.f) * OOrho;
-
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - BGK computation
- const real cusq = c3o2*(vvx*vvx + vvy*vvy + vvz*vvz);
-
- mfbbb = mfbbb *(c1o1 + (-omega)) - (-omega)* c8o27* (drho - rho * cusq);
- mfcbb = mfcbb *(c1o1 + (-omega)) - (-omega)* c2o27* (drho + rho * (c3o1*(vvx)+c9o2*(vvx)*(vvx)-cusq));
- mfabb = mfabb *(c1o1 + (-omega)) - (-omega)* c2o27* (drho + rho * (c3o1*(-vvx) + c9o2*(-vvx)*(-vvx) - cusq));
- mfbcb = mfbcb *(c1o1 + (-omega)) - (-omega)* c2o27* (drho + rho * (c3o1*(vvy)+c9o2*(vvy)*(vvy)-cusq));
- mfbab = mfbab *(c1o1 + (-omega)) - (-omega)* c2o27* (drho + rho * (c3o1*(-vvy) + c9o2*(-vvy)*(-vvy) - cusq));
- mfbbc = mfbbc *(c1o1 + (-omega)) - (-omega)* c2o27* (drho + rho * (c3o1*(vvz)+c9o2*(vvz)*(vvz)-cusq));
- mfbba = mfbba *(c1o1 + (-omega)) - (-omega)* c2o27* (drho + rho * (c3o1*(-vvz) + c9o2*(-vvz)*(-vvz) - cusq));
- mfccb = mfccb *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(vvx + vvy) + c9o2*(vvx + vvy)*(vvx + vvy) - cusq));
- mfaab = mfaab *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(-vvx - vvy) + c9o2*(-vvx - vvy)*(-vvx - vvy) - cusq));
- mfcab = mfcab *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(vvx - vvy) + c9o2*(vvx - vvy)*(vvx - vvy) - cusq));
- mfacb = mfacb *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(-vvx + vvy) + c9o2*(-vvx + vvy)*(-vvx + vvy) - cusq));
- mfcbc = mfcbc *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(vvx + vvz) + c9o2*(vvx + vvz)*(vvx + vvz) - cusq));
- mfaba = mfaba *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(-vvx - vvz) + c9o2*(-vvx - vvz)*(-vvx - vvz) - cusq));
- mfcba = mfcba *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(vvx - vvz) + c9o2*(vvx - vvz)*(vvx - vvz) - cusq));
- mfabc = mfabc *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(-vvx + vvz) + c9o2*(-vvx + vvz)*(-vvx + vvz) - cusq));
- mfbcc = mfbcc *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(vvy + vvz) + c9o2*(vvy + vvz)*(vvy + vvz) - cusq));
- mfbaa = mfbaa *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(-vvy - vvz) + c9o2*(-vvy - vvz)*(-vvy - vvz) - cusq));
- mfbca = mfbca *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(vvy - vvz) + c9o2*(vvy - vvz)*(vvy - vvz) - cusq));
- mfbac = mfbac *(c1o1 + (-omega)) - (-omega)* c1o54* (drho + rho * (c3o1*(-vvy + vvz) + c9o2*(-vvy + vvz)*(-vvy + vvz) - cusq));
- mfccc = mfccc *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(vvx + vvy + vvz) + c9o2*(vvx + vvy + vvz)*(vvx + vvy + vvz) - cusq));
- mfaaa = mfaaa *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(-vvx - vvy - vvz) + c9o2*(-vvx - vvy - vvz)*(-vvx - vvy - vvz) - cusq));
- mfcca = mfcca *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(vvx + vvy - vvz) + c9o2*(vvx + vvy - vvz)*(vvx + vvy - vvz) - cusq));
- mfaac = mfaac *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(-vvx - vvy + vvz) + c9o2*(-vvx - vvy + vvz)*(-vvx - vvy + vvz) - cusq));
- mfcac = mfcac *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(vvx - vvy + vvz) + c9o2*(vvx - vvy + vvz)*(vvx - vvy + vvz) - cusq));
- mfaca = mfaca *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(-vvx + vvy - vvz) + c9o2*(-vvx + vvy - vvz)*(-vvx + vvy - vvz) - cusq));
- mfcaa = mfcaa *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(vvx - vvy - vvz) + c9o2*(vvx - vvy - vvz)*(vvx - vvy - vvz) - cusq));
- mfacc = mfacc *(c1o1 + (-omega)) - (-omega)* c1o216*(drho + rho * (c3o1*(-vvx + vvy + vvz) + c9o2*(-vvx + vvy + vvz)*(-vvx + vvy + vvz) - cusq));
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Write distributions: style of reading and writing the distributions from/to
- //! stored arrays dependent on timestep is based on the esoteric twist algorithm
- //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
- //!
- distribution.f[dir::MZZ] = mfcbb;
- distribution.f[dir::PZZ] = mfabb;
- distribution.f[dir::ZMZ] = mfbcb;
- distribution.f[dir::ZPZ] = mfbab;
- distribution.f[dir::ZZM] = mfbbc;
- distribution.f[dir::ZZP] = mfbba;
- distribution.f[dir::MMZ] = mfccb;
- distribution.f[dir::PPZ] = mfaab;
- distribution.f[dir::MPZ] = mfcab;
- distribution.f[dir::PMZ] = mfacb;
- distribution.f[dir::MZM] = mfcbc;
- distribution.f[dir::PZP] = mfaba;
- distribution.f[dir::MZP] = mfcba;
- distribution.f[dir::PZM] = mfabc;
- distribution.f[dir::ZMM] = mfbcc;
- distribution.f[dir::ZPP] = mfbaa;
- distribution.f[dir::ZMP] = mfbca;
- distribution.f[dir::ZPM] = mfbac;
- distribution.f[dir::MMM] = mfccc;
- distribution.f[dir::PMM] = mfacc;
- distribution.f[dir::MPM] = mfcac;
- distribution.f[dir::PPM] = mfaac;
- distribution.f[dir::MMP] = mfcca;
- distribution.f[dir::PMP] = mfaca;
- distribution.f[dir::MPP] = mfcaa;
- distribution.f[dir::PPP] = mfaaa;
- distribution.f[dir::ZZZ] = mfbbb;
-}
-
-
-}
-
diff --git a/src/lbm/BGK.h b/src/lbm/BGK.h
deleted file mode 100644
index 6cde85013dd92472022bbf7b93bc73e7940049a1..0000000000000000000000000000000000000000
--- a/src/lbm/BGK.h
+++ /dev/null
@@ -1,24 +0,0 @@
-#ifndef LBM_BGK_H
-#define LBM_BGK_H
-
-#ifndef __host__
-#define __host__
-#endif
-#ifndef __device__
-#define __device__
-#endif
-
-#include <basics/DataTypes.h>
-
-#include "KernelParameter.h"
-
-namespace vf
-{
-namespace lbm
-{
-
-__host__ __device__ void bgk(KernelParameter parameter);
-
-}
-}
-#endif
diff --git a/src/lbm/CMakeLists.txt b/src/lbm/CMakeLists.txt
index 52ab3d78710c8551475307463334c9d1d0baf36f..73a98aed61ab05bfaee065e60480bac6f41d262f 100644
--- a/src/lbm/CMakeLists.txt
+++ b/src/lbm/CMakeLists.txt
@@ -3,10 +3,6 @@ vf_add_library(PUBLIC_LINK basics)
if(BUILD_VF_GPU)
set_target_properties(lbm PROPERTIES CUDA_SEPARABLE_COMPILATION ON POSITION_INDEPENDENT_CODE ON)
-
- set_source_files_properties(KernelParameter.cpp PROPERTIES LANGUAGE CUDA)
- set_source_files_properties(CumulantChimera.cpp PROPERTIES LANGUAGE CUDA)
- set_source_files_properties(BGK.cpp PROPERTIES LANGUAGE CUDA)
endif()
vf_add_tests()
\ No newline at end of file
diff --git a/src/lbm/ChimeraTests.cpp b/src/lbm/ChimeraTests.cpp
index 01abbe82764c276f53a4cdd479d333536199c3a9..4e0370c60a5d2ead900b562390acec0646cf8400 100644
--- a/src/lbm/ChimeraTests.cpp
+++ b/src/lbm/ChimeraTests.cpp
@@ -1,6 +1,6 @@
#include <gmock/gmock.h>
-#include "Chimera.h"
+#include "ChimeraTransformation.h"
#ifdef VF_DOUBLE_ACCURACY
#define REAL_EQ(a) testing::DoubleEq(a)
diff --git a/src/lbm/Chimera.h b/src/lbm/ChimeraTransformation.h
similarity index 100%
rename from src/lbm/Chimera.h
rename to src/lbm/ChimeraTransformation.h
diff --git a/src/lbm/CumulantChimera.cpp b/src/lbm/CumulantChimera.cpp
deleted file mode 100644
index f1257c4e579b36f868f3476282e79d44291750b1..0000000000000000000000000000000000000000
--- a/src/lbm/CumulantChimera.cpp
+++ /dev/null
@@ -1,448 +0,0 @@
-
-#include <cmath>
-
-#include <basics/DataTypes.h>
-#include <basics/constants/NumericConstants.h>
-
-#include "constants/D3Q27.h"
-#include "CumulantChimera.h"
-#include "Chimera.h"
-#include "MacroscopicQuantities.h"
-
-namespace vf::lbm
-{
-
-using namespace vf::basics::constant;
-
-
-////////////////////////////////////////////////////////////////////////////////////
-//! - Setting relaxation rates for non-hydrodynamic cumulants (default values). Variable names and equations according to
-//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
-//! => [NAME IN PAPER]=[NAME IN CODE]=[DEFAULT VALUE].
-//! - Trace of second order cumulants \f$ C_{200}+C_{020}+C_{002} \f$ used to adjust bulk viscosity:\f$\omega_2=OxxPyyPzz=1.0 \f$.
-//! - Third order cumulants \f$ C_{120}+C_{102}, C_{210}+C_{012}, C_{201}+C_{021} \f$: \f$ \omega_3=OxyyPxzz \f$ set according to Eq. (111) with simplifications assuming \f$ \omega_2=1.0\f$.
-//! - Third order cumulants \f$ C_{120}-C_{102}, C_{210}-C_{012}, C_{201}-C_{021} \f$: \f$ \omega_4 = OxyyMxzz \f$ set according to Eq. (112) with simplifications assuming \f$ \omega_2 = 1.0\f$.
-//! - Third order cumulants \f$ C_{111} \f$: \f$ \omega_5 = Oxyz \f$ set according to Eq. (113) with simplifications assuming \f$ \omega_2 = 1.0\f$ (modify for different bulk viscosity).
-//! - Fourth order cumulants \f$ C_{220}, C_{202}, C_{022}, C_{211}, C_{121}, C_{112} \f$: for simplification all set to the same default value \f$ \omega_6=\omega_7=\omega_8=O4=1.0 \f$.
-//! - Fifth order cumulants \f$ C_{221}, C_{212}, C_{122}\f$: \f$\omega_9=O5=1.0\f$.
-//! - Sixth order cumulant \f$ C_{222}\f$: \f$\omega_{10}=O6=1.0\f$.
-//////////////////////////////////////////////////////////////////////////
-__host__ __device__ void setRelaxationRatesK17(real omega, real &OxxPyyPzz, real &OxyyPxzz, real &OxyyMxzz, real &Oxyz,
- real &O4, real &O5, real &O6)
-{
- OxxPyyPzz = c1o1;
-
- OxyyPxzz = c8o1 * (-c2o1 + omega) * (c1o1 + c2o1 * omega) / (-c8o1 - c14o1 * omega + c7o1 * omega * omega);
- OxyyMxzz = c8o1 * (-c2o1 + omega) * (-c7o1 + c4o1 * omega) / (c56o1 - c50o1 * omega + c9o1 * omega * omega);
- Oxyz = c24o1 * (-c2o1 + omega) * (-c2o1 - c7o1 * omega + c3o1 * omega * omega) /
- (c48o1 + c152o1 * omega - c130o1 * omega * omega + c29o1 * omega * omega * omega);
-
- O4 = c1o1;
-
- O5 = c1o1;
-
- O6 = c1o1;
-}
-
-
-__host__ __device__ void setRelaxationRatesK15(real omega, real &OxxPyyPzz, real &OxyyPxzz, real &OxyyMxzz, real &Oxyz,
- real &O4, real &O5, real &O6)
-{
- OxxPyyPzz = c1o1;
-
- OxyyPxzz = c1o1;
- OxyyMxzz = c1o1;
- Oxyz = c1o1;
-
- O4 = c1o1;
-
- O5 = c1o1;
-
- O6 = c1o1;
-}
-
-//////////////////////////////////////////////////////////////////////////
-//! Cumulant K17 Kernel is based on \ref
-//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
-//! and \ref
-//! <a href="https://doi.org/10.1016/j.jcp.2017.07.004"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.07.004 ]</b></a>
-//////////////////////////////////////////////////////////////////////////
-__host__ __device__ void cumulantChimera(KernelParameter parameter, RelaxationRatesFunctor setRelaxationRates)
-{
- auto& distribution = parameter.distribution;
- const auto omega = parameter.omega;
- const auto* forces = parameter.forces;
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Read distributions: style of reading and writing the distributions from/to
- //! stored arrays dependent on timestep is based on the esoteric twist algorithm
- //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
- //!
- real mfcbb = distribution.f[dir::PZZ];
- real mfabb = distribution.f[dir::MZZ];
- real mfbcb = distribution.f[dir::ZPZ];
- real mfbab = distribution.f[dir::ZMZ];
- real mfbbc = distribution.f[dir::ZZP];
- real mfbba = distribution.f[dir::ZZM];
- real mfccb = distribution.f[dir::PPZ];
- real mfaab = distribution.f[dir::MMZ];
- real mfcab = distribution.f[dir::PMZ];
- real mfacb = distribution.f[dir::MPZ];
- real mfcbc = distribution.f[dir::PZP];
- real mfaba = distribution.f[dir::MZM];
- real mfcba = distribution.f[dir::PZM];
- real mfabc = distribution.f[dir::MZP];
- real mfbcc = distribution.f[dir::ZPP];
- real mfbaa = distribution.f[dir::ZMM];
- real mfbca = distribution.f[dir::ZPM];
- real mfbac = distribution.f[dir::ZMP];
- real mfccc = distribution.f[dir::PPP];
- real mfacc = distribution.f[dir::MPP];
- real mfcac = distribution.f[dir::PMP];
- real mfaac = distribution.f[dir::MMP];
- real mfcca = distribution.f[dir::PPM];
- real mfaca = distribution.f[dir::MPM];
- real mfcaa = distribution.f[dir::PMM];
- real mfaaa = distribution.f[dir::MMM];
- real mfbbb = distribution.f[dir::ZZZ];
-
-
- const real drho = getDensity(distribution.f);
- const real OOrho = c1o1 / (c1o1 + drho);
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref
- //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
- //!
- const real vvx = getIncompressibleVelocityX1(distribution.f) * OOrho + forces[0] * c1o2;
- const real vvy = getIncompressibleVelocityX2(distribution.f) * OOrho + forces[0] * c1o2;
- const real vvz = getIncompressibleVelocityX3(distribution.f) * OOrho + forces[0] * c1o2;
-
- ////////////////////////////////////////////////////////////////////////////////////
- // calculate the square of velocities for this lattice node
- const real vx2 = vvx*vvx;
- const real vy2 = vvy*vvy;
- const real vz2 = vvz*vvz;
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Chimera transform from well conditioned distributions to central moments as defined in Appendix J in \ref
- //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
- //! see also Eq. (6)-(14) in \ref
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //!
- ////////////////////////////////////////////////////////////////////////////////////
- // Z - Dir
- vf::lbm::forwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36o1, c1o36);
- vf::lbm::forwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9o1, c1o9);
- vf::lbm::forwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36o1, c1o36);
- vf::lbm::forwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9o1, c1o9);
- vf::lbm::forwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
- vf::lbm::forwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9o1, c1o9);
- vf::lbm::forwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36o1, c1o36);
- vf::lbm::forwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9o1, c1o9);
- vf::lbm::forwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36o1, c1o36);
- ////////////////////////////////////////////////////////////////////////////////////
- // Y - Dir
- vf::lbm::forwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6o1, c1o6);
- vf::lbm::forwardChimera( mfaab, mfabb, mfacb, vvy, vy2);
- vf::lbm::forwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18o1, c1o18);
- vf::lbm::forwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
- vf::lbm::forwardChimera( mfbab, mfbbb, mfbcb, vvy, vy2);
- vf::lbm::forwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
- vf::lbm::forwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6o1, c1o6);
- vf::lbm::forwardChimera( mfcab, mfcbb, mfccb, vvy, vy2);
- vf::lbm::forwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18o1, c1o18);
- ////////////////////////////////////////////////////////////////////////////////////
- // X - Dir
- vf::lbm::forwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1o1, c1o1);
- vf::lbm::forwardChimera( mfaba, mfbba, mfcba, vvx, vx2);
- vf::lbm::forwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3o1, c1o3);
- vf::lbm::forwardChimera( mfaab, mfbab, mfcab, vvx, vx2);
- vf::lbm::forwardChimera( mfabb, mfbbb, mfcbb, vvx, vx2);
- vf::lbm::forwardChimera( mfacb, mfbcb, mfccb, vvx, vx2);
- vf::lbm::forwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3o1, c1o3);
- vf::lbm::forwardChimera( mfabc, mfbbc, mfcbc, vvx, vx2);
- vf::lbm::forwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c3o1, c1o9);
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Setting relaxation rates for non-hydrodynamic cumulants (default values). Variable names and equations
- real OxxPyyPzz;
- real OxyyPxzz;
- real OxyyMxzz;
- real Oxyz;
- real O4;
- real O5;
- real O6;
-
- setRelaxationRates(omega, OxxPyyPzz, OxyyPxzz, OxyyMxzz, Oxyz, O4, O5, O6);
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - A and B: parameters for fourth order convergence of the diffusion term according to Eq. (115) and (116)
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //! with simplifications assuming \f$ \omega_2 = 1.0 \f$ (modify for different bulk viscosity).
- //!
- const real A = (c4o1 + c2o1*omega - c3o1*omega*omega) / (c2o1 - c7o1*omega + c5o1*omega*omega);
- const real B = (c4o1 + c28o1*omega - c14o1*omega*omega) / (c6o1 - c21o1*omega + c15o1*omega*omega);
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Compute cumulants from central moments according to Eq. (20)-(23) in
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //!
- ////////////////////////////////////////////////////////////
- //4.
- real CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + c2o1 * mfbba * mfbab) * OOrho;
- real CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + c2o1 * mfbba * mfabb) * OOrho;
- real CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + c2o1 * mfbab * mfabb) * OOrho;
- real CUMcca = mfcca - (((mfcaa * mfaca + c2o1 * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) * OOrho - c1o9*(drho * OOrho));
- real CUMcac = mfcac - (((mfcaa * mfaac + c2o1 * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) * OOrho - c1o9*(drho * OOrho));
- real CUMacc = mfacc - (((mfaac * mfaca + c2o1 * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) * OOrho - c1o9*(drho * OOrho));
- ////////////////////////////////////////////////////////////
- //5.
- real CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + c4o1 * mfabb * mfbbb + c2o1 * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) * OOrho;
- real CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + c4o1 * mfbab * mfbbb + c2o1 * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) * OOrho;
- real CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + c4o1 * mfbba * mfbbb + c2o1 * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) * OOrho;
- ////////////////////////////////////////////////////////////
- //6.
- real CUMccc = mfccc + ((-c4o1 * mfbbb * mfbbb
- - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
- - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
- - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
- + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
- + c2o1 * (mfcaa * mfaca * mfaac)
- + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
- - c1o3 * (mfacc + mfcac + mfcca) * OOrho
- - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
- + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
- + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
- + c1o27*((drho * drho - drho) * OOrho * OOrho));
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Compute linear combinations of second and third order cumulants
- //!
- ////////////////////////////////////////////////////////////
- //2.
- real mxxPyyPzz = mfcaa + mfaca + mfaac;
- real mxxMyy = mfcaa - mfaca;
- real mxxMzz = mfcaa - mfaac;
- ////////////////////////////////////////////////////////////
- //3.
- real mxxyPyzz = mfcba + mfabc;
- real mxxyMyzz = mfcba - mfabc;
- real mxxzPyyz = mfcab + mfacb;
- real mxxzMyyz = mfcab - mfacb;
- real mxyyPxzz = mfbca + mfbac;
- real mxyyMxzz = mfbca - mfbac;
- ////////////////////////////////////////////////////////////////////////////////////
- //incl. correction
- ////////////////////////////////////////////////////////////
- //! - Compute velocity gradients from second order cumulants according to Eq. (27)-(32)
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //! Further explanations of the correction in viscosity in Appendix H of
- //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
- //! Note that the division by rho is omitted here as we need rho times the gradients later.
- //!
- const real Dxy = -c3o1*omega*mfbba;
- const real Dxz = -c3o1*omega*mfbab;
- const real Dyz = -c3o1*omega*mfabb;
- const real dxux = c1o2 * (-omega) *(mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
- const real dyuy = dxux + omega * c3o2 * mxxMyy;
- const real dzuz = dxux + omega * c3o2 * mxxMzz;
- ////////////////////////////////////////////////////////////
- //! - Relaxation of second order cumulants with correction terms according to Eq. (33)-(35) in
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //!
- mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz) - c3o1 * (c1o1 - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
- mxxMyy += omega * (-mxxMyy) - c3o1 * (c1o1 + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
- mxxMzz += omega * (-mxxMzz) - c3o1 * (c1o1 + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
- ////////////////////////////////////////////////////////////////////////////////////
- ////no correction
- //mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz);
- //mxxMyy += -(-omega) * (-mxxMyy);
- //mxxMzz += -(-omega) * (-mxxMzz);
- //////////////////////////////////////////////////////////////////////////
- mfabb += omega * (-mfabb);
- mfbab += omega * (-mfbab);
- mfbba += omega * (-mfbba);
- ////////////////////////////////////////////////////////////////////////////////////
- //relax
- //////////////////////////////////////////////////////////////////////////
- // incl. limiter
- //! Set relaxation limiters for third order cumulants to default value \f$ \lambda=0.001 \f$ according to section 6 in \ref
- //! - Relaxation of third order cumulants including limiter according to Eq. (116)-(123)
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //!
-
- const real qudricLimitP = c1o100;
- const real qudricLimitM = c1o100;
- const real qudricLimitD = c1o100;
-
- real wadjust = Oxyz + (c1o1 - Oxyz)*abs_internal(mfbbb) / (abs_internal(mfbbb) + qudricLimitD);
- mfbbb += wadjust * (-mfbbb);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs_internal(mxxyPyzz) / (abs_internal(mxxyPyzz) + qudricLimitP);
- mxxyPyzz += wadjust * (-mxxyPyzz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs_internal(mxxyMyzz) / (abs_internal(mxxyMyzz) + qudricLimitM);
- mxxyMyzz += wadjust * (-mxxyMyzz);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs_internal(mxxzPyyz) / (abs_internal(mxxzPyyz) + qudricLimitP);
- mxxzPyyz += wadjust * (-mxxzPyyz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs_internal(mxxzMyyz) / (abs_internal(mxxzMyyz) + qudricLimitM);
- mxxzMyyz += wadjust * (-mxxzMyyz);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz)*abs_internal(mxyyPxzz) / (abs_internal(mxyyPxzz) + qudricLimitP);
- mxyyPxzz += wadjust * (-mxyyPxzz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz)*abs_internal(mxyyMxzz) / (abs_internal(mxyyMxzz) + qudricLimitM);
- mxyyMxzz += wadjust * (-mxyyMxzz);
- //////////////////////////////////////////////////////////////////////////
- // no limiter
- //mfbbb += OxyyMxzz * (-mfbbb);
- //mxxyPyzz += OxyyPxzz * (-mxxyPyzz);
- //mxxyMyzz += OxyyMxzz * (-mxxyMyzz);
- //mxxzPyyz += OxyyPxzz * (-mxxzPyyz);
- //mxxzMyyz += OxyyMxzz * (-mxxzMyyz);
- //mxyyPxzz += OxyyPxzz * (-mxyyPxzz);
- //mxyyMxzz += OxyyMxzz * (-mxyyMxzz);
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Compute inverse linear combinations of second and third order cumulants
- //!
- mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
- mfaca = c1o3 * (-c2o1* mxxMyy + mxxMzz + mxxPyyPzz);
- mfaac = c1o3 * (mxxMyy - c2o1* mxxMzz + mxxPyyPzz);
- mfcba = ( mxxyMyzz + mxxyPyzz) * c1o2;
- mfabc = (-mxxyMyzz + mxxyPyzz) * c1o2;
- mfcab = ( mxxzMyyz + mxxzPyyz) * c1o2;
- mfacb = (-mxxzMyyz + mxxzPyyz) * c1o2;
- mfbca = ( mxyyMxzz + mxyyPxzz) * c1o2;
- mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2;
- //////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////
- //4.
- // no limiter
- //! - Relax fourth order cumulants to modified equilibrium for fourth order convergence of diffusion according to Eq. (43)-(48)
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //!
- CUMacc = -O4*(c1o1 / omega - c1o2) * (dyuy + dzuz) * c2o3 * A + (c1o1 - O4) * (CUMacc);
- CUMcac = -O4*(c1o1 / omega - c1o2) * (dxux + dzuz) * c2o3 * A + (c1o1 - O4) * (CUMcac);
- CUMcca = -O4*(c1o1 / omega - c1o2) * (dyuy + dxux) * c2o3 * A + (c1o1 - O4) * (CUMcca);
- CUMbbc = -O4*(c1o1 / omega - c1o2) * Dxy * c1o3 * B + (c1o1 - O4) * (CUMbbc);
- CUMbcb = -O4*(c1o1 / omega - c1o2) * Dxz * c1o3 * B + (c1o1 - O4) * (CUMbcb);
- CUMcbb = -O4*(c1o1 / omega - c1o2) * Dyz * c1o3 * B + (c1o1 - O4) * (CUMcbb);
- //////////////////////////////////////////////////////////////////////////
- //5.
- CUMbcc += O5 * (-CUMbcc);
- CUMcbc += O5 * (-CUMcbc);
- CUMccb += O5 * (-CUMccb);
- //////////////////////////////////////////////////////////////////////////
- //6.
- CUMccc += O6 * (-CUMccc);
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Compute central moments from post collision cumulants according to Eq. (53)-(56) in
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //!
- //////////////////////////////////////////////////////////////////////////
- //4.
- mfcbb = CUMcbb + c1o3*((c3o1*mfcaa + c1o1) * mfabb + c6o1 * mfbba * mfbab) * OOrho;
- mfbcb = CUMbcb + c1o3*((c3o1*mfaca + c1o1) * mfbab + c6o1 * mfbba * mfabb) * OOrho;
- mfbbc = CUMbbc + c1o3*((c3o1*mfaac + c1o1) * mfbba + c6o1 * mfbab * mfabb) * OOrho;
- mfcca = CUMcca + (((mfcaa * mfaca + c2o1 * mfbba * mfbba)*c9o1 + c3o1 * (mfcaa + mfaca)) * OOrho - (drho * OOrho))*c1o9;
- mfcac = CUMcac + (((mfcaa * mfaac + c2o1 * mfbab * mfbab)*c9o1 + c3o1 * (mfcaa + mfaac)) * OOrho - (drho * OOrho))*c1o9;
- mfacc = CUMacc + (((mfaac * mfaca + c2o1 * mfabb * mfabb)*c9o1 + c3o1 * (mfaac + mfaca)) * OOrho - (drho * OOrho))*c1o9;
- //////////////////////////////////////////////////////////////////////////
- //5.
- mfbcc = CUMbcc + c1o3 *(c3o1*(mfaac * mfbca + mfaca * mfbac + c4o1 * mfabb * mfbbb + c2o1 * (mfbab * mfacb + mfbba * mfabc)) + (mfbca + mfbac)) * OOrho;
- mfcbc = CUMcbc + c1o3 *(c3o1*(mfaac * mfcba + mfcaa * mfabc + c4o1 * mfbab * mfbbb + c2o1 * (mfabb * mfcab + mfbba * mfbac)) + (mfcba + mfabc)) * OOrho;
- mfccb = CUMccb + c1o3 *(c3o1*(mfcaa * mfacb + mfaca * mfcab + c4o1 * mfbba * mfbbb + c2o1 * (mfbab * mfbca + mfabb * mfcba)) + (mfacb + mfcab)) * OOrho;
- //////////////////////////////////////////////////////////////////////////
- //6.
- mfccc = CUMccc - ((-c4o1 * mfbbb * mfbbb
- - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
- - c4o1 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
- - c2o1 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
- + (c4o1 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
- + c2o1 * (mfcaa * mfaca * mfaac)
- + c16o1 * mfbba * mfbab * mfabb) * OOrho * OOrho
- - c1o3 * (mfacc + mfcac + mfcca) * OOrho
- - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
- + (c2o1 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
- + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 *(mfaac + mfaca + mfcaa)) * OOrho * OOrho * c2o3
- + c1o27*((drho * drho - drho) * OOrho * OOrho));
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Add acceleration (body force) to first order cumulants according to Eq. (85)-(87) in
- //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
- //!
- mfbaa = -mfbaa;
- mfaba = -mfaba;
- mfaab = -mfaab;
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Chimera transform from central moments to well conditioned distributions as defined in Appendix J in
- //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015), DOI:10.1016/j.camwa 2015.05.001 ]</b></a>
- //! see also Eq. (88)-(96) in
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05 040 ]</b></a>
- //!
- ////////////////////////////////////////////////////////////////////////////////////
- // X - Dir
- vf::lbm::backwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1o1, c1o1);
- vf::lbm::backwardChimera( mfaba, mfbba, mfcba, vvx, vx2);
- vf::lbm::backwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3o1, c1o3);
- vf::lbm::backwardChimera( mfaab, mfbab, mfcab, vvx, vx2);
- vf::lbm::backwardChimera( mfabb, mfbbb, mfcbb, vvx, vx2);
- vf::lbm::backwardChimera( mfacb, mfbcb, mfccb, vvx, vx2);
- vf::lbm::backwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3o1, c1o3);
- vf::lbm::backwardChimera( mfabc, mfbbc, mfcbc, vvx, vx2);
- vf::lbm::backwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9o1, c1o9);
- ////////////////////////////////////////////////////////////////////////////////////
- // Y - Dir
- vf::lbm::backwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6o1, c1o6);
- vf::lbm::backwardChimera( mfaab, mfabb, mfacb, vvy, vy2);
- vf::lbm::backwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18o1, c1o18);
- vf::lbm::backwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
- vf::lbm::backwardChimera( mfbab, mfbbb, mfbcb, vvy, vy2);
- vf::lbm::backwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
- vf::lbm::backwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6o1, c1o6);
- vf::lbm::backwardChimera( mfcab, mfcbb, mfccb, vvy, vy2);
- vf::lbm::backwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18o1, c1o18);
- ////////////////////////////////////////////////////////////////////////////////////
- // Z - Dir
- vf::lbm::backwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36o1, c1o36);
- vf::lbm::backwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9o1, c1o9);
- vf::lbm::backwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36o1, c1o36);
- vf::lbm::backwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9o1, c1o9);
- vf::lbm::backwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
- vf::lbm::backwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9o1, c1o9);
- vf::lbm::backwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36o1, c1o36);
- vf::lbm::backwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9o1, c1o9);
- vf::lbm::backwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36o1, c1o36);
-
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Write distributions: style of reading and writing the distributions from/to
- //! stored arrays dependent on timestep is based on the esoteric twist algorithm
- //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017), DOI:10.3390/computation5020019 ]</b></a>
- //!
- distribution.f[dir::MZZ] = mfcbb;
- distribution.f[dir::PZZ] = mfabb;
- distribution.f[dir::ZMZ] = mfbcb;
- distribution.f[dir::ZPZ] = mfbab;
- distribution.f[dir::ZZM] = mfbbc;
- distribution.f[dir::ZZP] = mfbba;
- distribution.f[dir::MMZ] = mfccb;
- distribution.f[dir::PPZ] = mfaab;
- distribution.f[dir::MPZ] = mfcab;
- distribution.f[dir::PMZ] = mfacb;
- distribution.f[dir::MZM] = mfcbc;
- distribution.f[dir::PZP] = mfaba;
- distribution.f[dir::MZP] = mfcba;
- distribution.f[dir::PZM] = mfabc;
- distribution.f[dir::ZMM] = mfbcc;
- distribution.f[dir::ZPP] = mfbaa;
- distribution.f[dir::ZMP] = mfbca;
- distribution.f[dir::ZPM] = mfbac;
- distribution.f[dir::MMM] = mfccc;
- distribution.f[dir::PMM] = mfacc;
- distribution.f[dir::MPM] = mfcac;
- distribution.f[dir::PPM] = mfaac;
- distribution.f[dir::MMP] = mfcca;
- distribution.f[dir::PMP] = mfaca;
- distribution.f[dir::MPP] = mfcaa;
- distribution.f[dir::PPP] = mfaaa;
- distribution.f[dir::ZZZ] = mfbbb;
-}
-
-
-}
-
diff --git a/src/lbm/CumulantChimera.h b/src/lbm/CumulantChimera.h
deleted file mode 100644
index c30a0c07912953cff45e3734c0b60fa7a03acd53..0000000000000000000000000000000000000000
--- a/src/lbm/CumulantChimera.h
+++ /dev/null
@@ -1,34 +0,0 @@
-#ifndef LBM_CUMULANT_CHIMERA_H
-#define LBM_CUMULANT_CHIMERA_H
-
-#ifndef __host__
-#define __host__
-#endif
-#ifndef __device__
-#define __device__
-#endif
-
-#include <basics/DataTypes.h>
-
-#include "KernelParameter.h"
-
-namespace vf
-{
-namespace lbm
-{
-
-__host__ __device__ void setRelaxationRatesK17(real omega, real &OxxPyyPzz, real &OxyyPxzz, real &OxyyMxzz, real &Oxyz,
- real &O4, real &O5, real &O6);
-
-__host__ __device__ void setRelaxationRatesK15(real omega, real &OxxPyyPzz, real &OxyyPxzz, real &OxyyMxzz, real &Oxyz,
- real &O4, real &O5, real &O6);
-
-using RelaxationRatesFunctor = void(*)(real omega, real &OxxPyyPzz, real &OxyyPxzz, real &OxyyMxzz, real &Oxyz,
- real &O4, real &O5, real &O6);
-
-
-__host__ __device__ void cumulantChimera(KernelParameter parameter, RelaxationRatesFunctor setRelaxationRates);
-
-}
-}
-#endif
diff --git a/src/lbm/KernelParameter.cpp b/src/lbm/KernelParameter.cpp
deleted file mode 100644
index 7bf5a369d0e5d4e673d79dcb30bc22fc2c330e68..0000000000000000000000000000000000000000
--- a/src/lbm/KernelParameter.cpp
+++ /dev/null
@@ -1,27 +0,0 @@
-#include "KernelParameter.h"
-
-#include <cmath>
-
-#include "MacroscopicQuantities.h"
-
-
-namespace vf::lbm
-{
-
-
-inline __host__ __device__ real Distribution27::getDensity_() const
-{
- return getDensity(f);
-}
-
-__host__ __device__ real abs_internal(real value)
-{
-#ifdef __CUDA_ARCH__
- return ::abs(value);
-#else
- return std::abs(value);
-#endif
-}
-
-
-}
diff --git a/src/lbm/KernelParameter.h b/src/lbm/KernelParameter.h
deleted file mode 100644
index 9c07524226a40aaa9e2c65e7ab028b07aec62ddc..0000000000000000000000000000000000000000
--- a/src/lbm/KernelParameter.h
+++ /dev/null
@@ -1,38 +0,0 @@
-#ifndef LBM_KERNEL_PARAMETER_H
-#define LBM_KERNEL_PARAMETER_H
-
-#ifndef __host__
-#define __host__
-#endif
-#ifndef __device__
-#define __device__
-#endif
-
-#include <basics/DataTypes.h>
-
-
-namespace vf::lbm
-{
-
-struct Distribution27
-{
- real f[27];
-
- __host__ __device__ real getDensity_() const;
-};
-
-
-__host__ __device__ real abs_internal(real value);
-
-
-struct KernelParameter
-{
- Distribution27& distribution;
- real omega;
- real* forces;
-};
-
-
-}
-
-#endif
diff --git a/src/lbm/collision/CollisionParameter.h b/src/lbm/collision/CollisionParameter.h
new file mode 100644
index 0000000000000000000000000000000000000000..1f358f32f5bcf6116ade856442e2a6c1fe1627a1
--- /dev/null
+++ b/src/lbm/collision/CollisionParameter.h
@@ -0,0 +1,35 @@
+#ifndef LBM_KERNEL_PARAMETER_H
+#define LBM_KERNEL_PARAMETER_H
+
+#include <basics/DataTypes.h>
+
+namespace vf::lbm
+{
+
+struct CollisionParameter
+{
+ real distribution[27];
+ real omega;
+ real* quadricLimiter;
+ real forceX;
+ real forceY;
+ real forceZ;
+};
+
+struct TurbulentViscosity
+{
+ real SGSconstant;
+ real value { 0. };
+};
+
+struct MacroscopicValues
+{
+ real vx { 0. };
+ real vy { 0. };
+ real vz { 0. };
+ real rho { 0. };
+};
+
+} // namespace vf::lbm
+
+#endif
diff --git a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cu b/src/lbm/collision/K17CompressibleNavierStokes.h
similarity index 60%
rename from src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cu
rename to src/lbm/collision/K17CompressibleNavierStokes.h
index 732e96ce44e64d2f86bf899bd87bb33de4725874..af0d55e419212265cd4760eb6d82198f80e61499 100644
--- a/src/gpu/VirtualFluids_GPU/Kernel/Compressible/FluidFlow/K17/K17CompressibleNavierStokes_Device.cu
+++ b/src/lbm/collision/K17CompressibleNavierStokes.h
@@ -41,113 +41,76 @@
//! are provided by the utilized PostCollisionInteractiors depending on they specific requirements (e.g. writeMacroscopicVariables for probes).
//=======================================================================================
-#include "LBM/LB.h"
+#include <basics/constants/NumericConstants.h>
+
#include "lbm/constants/D3Q27.h"
-#include "basics/constants/NumericConstants.h"
-#include "LBM/GPUHelperFunctions/KernelUtilities.h"
-#include "LBM/GPUHelperFunctions/ChimeraTransformation.h"
-#include "GPU/TurbulentViscosityInlines.cuh"
+#include "ChimeraTransformation.h"
+
+#include "TurbulentViscosity.h"
+
+#include "MacroscopicQuantities.h"
+
+#include "CollisionParameter.h"
+
+#ifndef __host__
+#define __host__
+#endif
+#ifndef __device__
+#define __device__
+#endif
+
+#ifdef __CUDACC__
+#define KERNEL_ABS abs
+#else
+#include <cmath>
+#define KERNEL_ABS std::abs
+#endif
using namespace vf::basics::constant;
using namespace vf::lbm::dir;
-using namespace vf::gpu;
-////////////////////////////////////////////////////////////////////////////////
-template<TurbulenceModel turbulenceModel, bool writeMacroscopicVariables, bool applyBodyForce>
-__global__ void K17CompressibleNavierStokes_Device(
- real omega_in,
- uint* neighborX,
- uint* neighborY,
- uint* neighborZ,
- real* distributions,
- real* rho,
- real* vx,
- real* vy,
- real* vz,
- real* turbulentViscosity,
- real SGSconstant,
- unsigned long long numberOfLBnodes,
- int level,
- real* forces,
- real* bodyForceX,
- real* bodyForceY,
- real* bodyForceZ,
- real* quadricLimiters,
- bool isEvenTimestep,
- const uint *fluidNodeIndices,
- uint numberOfFluidNodes)
+namespace vf::lbm
{
- //////////////////////////////////////////////////////////////////////////
- //! Cumulant K17 Kernel is based on \ref
- //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040
- //! ]</b></a> and \ref <a href="https://doi.org/10.1016/j.jcp.2017.07.004"><b>[ M. Geier et al. (2017),
- //! DOI:10.1016/j.jcp.2017.07.004 ]</b></a>
- //!
- //! The cumulant kernel is executed in the following steps
- //!
- ////////////////////////////////////////////////////////////////////////////////
- //! - Get node index coordinates from threadIdx, blockIdx, blockDim and gridDim.
- //!
- const unsigned nodeIndex = getNodeIndex();
- //////////////////////////////////////////////////////////////////////////
- // run for all indices in size_Mat and fluid nodes
- if (nodeIndex >= numberOfFluidNodes)
- return;
- ////////////////////////////////////////////////////////////////////////////////
- //! - Get the node index from the array containing all indices of fluid nodes
- //!
- const unsigned k_000 = fluidNodeIndices[nodeIndex];
-
- //////////////////////////////////////////////////////////////////////////
- //! - Read distributions: style of reading and writing the distributions from/to stored arrays dependent on
- //! timestep is based on the esoteric twist algorithm \ref <a
- //! href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017),
- //! DOI:10.3390/computation5020019 ]</b></a>
- //!
- Distributions27 dist;
- getPointersToDistributions(dist, distributions, numberOfLBnodes, isEvenTimestep);
-
- ////////////////////////////////////////////////////////////////////////////////
- //! - Set neighbor indices (necessary for indirect addressing)
- uint k_M00 = neighborX[k_000];
- uint k_0M0 = neighborY[k_000];
- uint k_00M = neighborZ[k_000];
- uint k_MM0 = neighborY[k_M00];
- uint k_M0M = neighborZ[k_M00];
- uint k_0MM = neighborZ[k_0M0];
- uint k_MMM = neighborZ[k_MM0];
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Set local distributions
- //!
- real f_000 = (dist.f[DIR_000])[k_000];
- real f_P00 = (dist.f[DIR_P00])[k_000];
- real f_M00 = (dist.f[DIR_M00])[k_M00];
- real f_0P0 = (dist.f[DIR_0P0])[k_000];
- real f_0M0 = (dist.f[DIR_0M0])[k_0M0];
- real f_00P = (dist.f[DIR_00P])[k_000];
- real f_00M = (dist.f[DIR_00M])[k_00M];
- real f_PP0 = (dist.f[DIR_PP0])[k_000];
- real f_MM0 = (dist.f[DIR_MM0])[k_MM0];
- real f_PM0 = (dist.f[DIR_PM0])[k_0M0];
- real f_MP0 = (dist.f[DIR_MP0])[k_M00];
- real f_P0P = (dist.f[DIR_P0P])[k_000];
- real f_M0M = (dist.f[DIR_M0M])[k_M0M];
- real f_P0M = (dist.f[DIR_P0M])[k_00M];
- real f_M0P = (dist.f[DIR_M0P])[k_M00];
- real f_0PP = (dist.f[DIR_0PP])[k_000];
- real f_0MM = (dist.f[DIR_0MM])[k_0MM];
- real f_0PM = (dist.f[DIR_0PM])[k_00M];
- real f_0MP = (dist.f[DIR_0MP])[k_0M0];
- real f_PPP = (dist.f[DIR_PPP])[k_000];
- real f_MPP = (dist.f[DIR_MPP])[k_M00];
- real f_PMP = (dist.f[DIR_PMP])[k_0M0];
- real f_MMP = (dist.f[DIR_MMP])[k_MM0];
- real f_PPM = (dist.f[DIR_PPM])[k_00M];
- real f_MPM = (dist.f[DIR_MPM])[k_M0M];
- real f_PMM = (dist.f[DIR_PMM])[k_0MM];
- real f_MMM = (dist.f[DIR_MMM])[k_MMM];
+//////////////////////////////////////////////////////////////////////////
+//! Cumulant K17 Kernel is based on \ref
+//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040
+//! ]</b></a> and \ref <a href="https://doi.org/10.1016/j.jcp.2017.07.004"><b>[ M. Geier et al. (2017),
+//! DOI:10.1016/j.jcp.2017.07.004 ]</b></a>
+////////////////////////////////////////////////////////////////////////////////
+template <TurbulenceModel turbulenceModel, bool writeMacroscopicVariables>
+__host__ __device__ void runK17CompressibleNavierStokes(CollisionParameter& parameter, MacroscopicValues& macroscopicValues, TurbulentViscosity& turbulentViscosity)
+{
+ auto& distribution = parameter.distribution;
+
+ real f_000 = distribution[DIR_000];
+ real f_P00 = distribution[DIR_P00];
+ real f_M00 = distribution[DIR_M00];
+ real f_0P0 = distribution[DIR_0P0];
+ real f_0M0 = distribution[DIR_0M0];
+ real f_00P = distribution[DIR_00P];
+ real f_00M = distribution[DIR_00M];
+ real f_PP0 = distribution[DIR_PP0];
+ real f_MM0 = distribution[DIR_MM0];
+ real f_PM0 = distribution[DIR_PM0];
+ real f_MP0 = distribution[DIR_MP0];
+ real f_P0P = distribution[DIR_P0P];
+ real f_M0M = distribution[DIR_M0M];
+ real f_P0M = distribution[DIR_P0M];
+ real f_M0P = distribution[DIR_M0P];
+ real f_0PP = distribution[DIR_0PP];
+ real f_0MM = distribution[DIR_0MM];
+ real f_0PM = distribution[DIR_0PM];
+ real f_0MP = distribution[DIR_0MP];
+ real f_PPP = distribution[DIR_PPP];
+ real f_MPP = distribution[DIR_MPP];
+ real f_PMP = distribution[DIR_PMP];
+ real f_MMP = distribution[DIR_MMP];
+ real f_PPM = distribution[DIR_PPM];
+ real f_MPM = distribution[DIR_MPM];
+ real f_PMM = distribution[DIR_PMM];
+ real f_MMM = distribution[DIR_MMM];
////////////////////////////////////////////////////////////////////////////////////
//! - Define aliases to use the same variable for the moments (m's):
@@ -185,78 +148,17 @@ __global__ void K17CompressibleNavierStokes_Device(
//! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
//! DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
//!
- real drho = ((((f_PPP + f_MMM) + (f_MPM + f_PMP)) + ((f_MPP + f_PMM) + (f_MMP + f_PPM))) +
- (((f_0MP + f_0PM) + (f_0MM + f_0PP)) + ((f_M0P + f_P0M) + (f_M0M + f_P0P)) +
- ((f_MP0 + f_PM0) + (f_MM0 + f_PP0))) +
- ((f_M00 + f_P00) + (f_0M0 + f_0P0) + (f_00M + f_00P))) +
- f_000;
-
- real oneOverRho = c1o1 / (c1o1 + drho);
-
- real vvx = ((((f_PPP - f_MMM) + (f_PMP - f_MPM)) + ((f_PMM - f_MPP) + (f_PPM - f_MMP))) +
- (((f_P0M - f_M0P) + (f_P0P - f_M0M)) + ((f_PM0 - f_MP0) + (f_PP0 - f_MM0))) + (f_P00 - f_M00)) *
- oneOverRho;
- real vvy = ((((f_PPP - f_MMM) + (f_MPM - f_PMP)) + ((f_MPP - f_PMM) + (f_PPM - f_MMP))) +
- (((f_0PM - f_0MP) + (f_0PP - f_0MM)) + ((f_MP0 - f_PM0) + (f_PP0 - f_MM0))) + (f_0P0 - f_0M0)) *
- oneOverRho;
- real vvz = ((((f_PPP - f_MMM) + (f_PMP - f_MPM)) + ((f_MPP - f_PMM) + (f_MMP - f_PPM))) +
- (((f_0MP - f_0PM) + (f_0PP - f_0MM)) + ((f_M0P - f_P0M) + (f_P0P - f_M0M))) + (f_00P - f_00M)) *
- oneOverRho;
+ real drho, oneOverRho, vvx, vvy, vvz;
+ getCompressibleMacroscopicValues(distribution, drho, oneOverRho, vvx, vvy, vvz);
////////////////////////////////////////////////////////////////////////////////////
//! - Add half of the acceleration (body force) to the velocity as in Eq. (42) \ref
//! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
//! DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
//!
- real factor = c1o1;
- for (size_t i = 1; i <= level; i++) {
- factor *= c2o1;
- }
-
- real fx = forces[0];
- real fy = forces[1];
- real fz = forces[2];
-
- if( applyBodyForce ){
- fx += bodyForceX[k_000];
- fy += bodyForceY[k_000];
- fz += bodyForceZ[k_000];
-
- // real vx = vvx;
- // real vy = vvy;
- // real vz = vvz;
- real acc_x = fx * c1o2 / factor;
- real acc_y = fy * c1o2 / factor;
- real acc_z = fz * c1o2 / factor;
-
- vvx += acc_x;
- vvy += acc_y;
- vvz += acc_z;
-
- // Reset body force. To be used when not using round-off correction.
- bodyForceX[k_000] = 0.0f;
- bodyForceY[k_000] = 0.0f;
- bodyForceZ[k_000] = 0.0f;
-
- ////////////////////////////////////////////////////////////////////////////////////
- //!> Round-off correction
- //!
- //!> Similar to Kahan summation algorithm (https://en.wikipedia.org/wiki/Kahan_summation_algorithm)
- //!> Essentially computes the round-off error of the applied force and adds it in the next time step as a compensation.
- //!> Seems to be necesseary at very high Re boundary layers, where the forcing and velocity can
- //!> differ by several orders of magnitude.
- //!> \note 16/05/2022: Testing, still ongoing!
- //!
- // bodyForceX[k_000] = (acc_x-(vvx-vx))*factor*c2o1;
- // bodyForceY[k_000] = (acc_y-(vvy-vy))*factor*c2o1;
- // bodyForceZ[k_000] = (acc_z-(vvz-vz))*factor*c2o1;
- }
- else{
- vvx += fx * c1o2 / factor;
- vvy += fy * c1o2 / factor;
- vvz += fz * c1o2 / factor;
- }
-
+ vvx += parameter.forceX;
+ vvy += parameter.forceY;
+ vvz += parameter.forceZ;
////////////////////////////////////////////////////////////////////////////////////
// calculate the square of velocities for this lattice node
@@ -268,9 +170,9 @@ __global__ void K17CompressibleNavierStokes_Device(
//! section 6 in \ref <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
//! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
//!
- real quadricLimitP = quadricLimiters[0];
- real quadricLimitM = quadricLimiters[1];
- real quadricLimitD = quadricLimiters[2];
+ real quadricLimitP = parameter.quadricLimiter[0];
+ real quadricLimitM = parameter.quadricLimiter[1];
+ real quadricLimitD = parameter.quadricLimiter[2];
////////////////////////////////////////////////////////////////////////////////////
//! - Chimera transform from well conditioned distributions to central moments as defined in Appendix J in \ref
//! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
@@ -335,8 +237,10 @@ __global__ void K17CompressibleNavierStokes_Device(
////////////////////////////////////////////////////////////////////////////////////
//! - Calculate modified omega with turbulent viscosity
//!
- real omega = omega_in;
- if(turbulenceModel != TurbulenceModel::None){ omega /= (c1o1 + c3o1*omega_in*turbulentViscosity[k_000]); }
+ real omega = parameter.omega;
+ if (turbulenceModel != TurbulenceModel::None) {
+ omega /= (c1o1 + c3o1 * parameter.omega * turbulentViscosity.value);
+ }
////////////////////////////////////////////////////////////
// 2.
real OxxPyyPzz = c1o1;
@@ -452,10 +356,10 @@ __global__ void K17CompressibleNavierStokes_Device(
case TurbulenceModel::AMD: //AMD is computed in separate kernel
break;
case TurbulenceModel::Smagorinsky:
- turbulentViscosity[k_000] = calcTurbulentViscositySmagorinsky(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
+ turbulentViscosity.value = calcTurbulentViscositySmagorinsky(turbulentViscosity.SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
break;
case TurbulenceModel::QR:
- turbulentViscosity[k_000] = calcTurbulentViscosityQR(SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
+ turbulentViscosity.value = calcTurbulentViscosityQR(turbulentViscosity.SGSconstant, dxux, dyuy, dzuz, Dxy, Dxz , Dyz);
break;
default:
break;
@@ -487,19 +391,19 @@ __global__ void K17CompressibleNavierStokes_Device(
//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
//! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
//!
- real wadjust = Oxyz + (c1o1 - Oxyz) * abs(m_111) / (abs(m_111) + quadricLimitD);
+ real wadjust = Oxyz + (c1o1 - Oxyz) * KERNEL_ABS(m_111) / (KERNEL_ABS(m_111) + quadricLimitD);
m_111 += wadjust * (-m_111);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxxyPyzz) / (abs(mxxyPyzz) + quadricLimitP);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * KERNEL_ABS(mxxyPyzz) / (KERNEL_ABS(mxxyPyzz) + quadricLimitP);
mxxyPyzz += wadjust * (-mxxyPyzz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxxyMyzz) / (abs(mxxyMyzz) + quadricLimitM);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * KERNEL_ABS(mxxyMyzz) / (KERNEL_ABS(mxxyMyzz) + quadricLimitM);
mxxyMyzz += wadjust * (-mxxyMyzz);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxxzPyyz) / (abs(mxxzPyyz) + quadricLimitP);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * KERNEL_ABS(mxxzPyyz) / (KERNEL_ABS(mxxzPyyz) + quadricLimitP);
mxxzPyyz += wadjust * (-mxxzPyyz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxxzMyyz) / (abs(mxxzMyyz) + quadricLimitM);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * KERNEL_ABS(mxxzMyyz) / (KERNEL_ABS(mxxzMyyz) + quadricLimitM);
mxxzMyyz += wadjust * (-mxxzMyyz);
- wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * abs(mxyyPxzz) / (abs(mxyyPxzz) + quadricLimitP);
+ wadjust = OxyyPxzz + (c1o1 - OxyyPxzz) * KERNEL_ABS(mxyyPxzz) / (KERNEL_ABS(mxyyPxzz) + quadricLimitP);
mxyyPxzz += wadjust * (-mxyyPxzz);
- wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * abs(mxyyMxzz) / (abs(mxyyMxzz) + quadricLimitM);
+ wadjust = OxyyMxzz + (c1o1 - OxyyMxzz) * KERNEL_ABS(mxyyMxzz) / (KERNEL_ABS(mxyyMxzz) + quadricLimitM);
mxyyMxzz += wadjust * (-mxyyMxzz);
//////////////////////////////////////////////////////////////////////////
// no limiter
@@ -606,13 +510,12 @@ __global__ void K17CompressibleNavierStokes_Device(
m_010 = -m_010;
m_001 = -m_001;
- //Write to array here to distribute read/write
- if(writeMacroscopicVariables || turbulenceModel==TurbulenceModel::AMD)
- {
- rho[k_000] = drho;
- vx[k_000] = vvx;
- vy[k_000] = vvy;
- vz[k_000] = vvz;
+ // Write to array here to distribute read/write
+ if (writeMacroscopicVariables || turbulenceModel == TurbulenceModel::AMD) {
+ macroscopicValues.rho = drho;
+ macroscopicValues.vx = vvx;
+ macroscopicValues.vy = vvy;
+ macroscopicValues.vz = vvz;
}
////////////////////////////////////////////////////////////////////////////////////
@@ -658,69 +561,33 @@ __global__ void K17CompressibleNavierStokes_Device(
backwardInverseChimeraWithK(m_210, m_211, m_212, vvz, vz2, c9o1, c1o9);
backwardInverseChimeraWithK(m_220, m_221, m_222, vvz, vz2, c36o1, c1o36);
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Write distributions: style of reading and writing the distributions from/to
- //! stored arrays dependent on timestep is based on the esoteric twist algorithm
- //! <a href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017),
- //! DOI:10.3390/computation5020019 ]</b></a>
- //!
- (dist.f[DIR_P00])[k_000] = f_M00;
- (dist.f[DIR_M00])[k_M00] = f_P00;
- (dist.f[DIR_0P0])[k_000] = f_0M0;
- (dist.f[DIR_0M0])[k_0M0] = f_0P0;
- (dist.f[DIR_00P])[k_000] = f_00M;
- (dist.f[DIR_00M])[k_00M] = f_00P;
- (dist.f[DIR_PP0])[k_000] = f_MM0;
- (dist.f[DIR_MM0])[k_MM0] = f_PP0;
- (dist.f[DIR_PM0])[k_0M0] = f_MP0;
- (dist.f[DIR_MP0])[k_M00] = f_PM0;
- (dist.f[DIR_P0P])[k_000] = f_M0M;
- (dist.f[DIR_M0M])[k_M0M] = f_P0P;
- (dist.f[DIR_P0M])[k_00M] = f_M0P;
- (dist.f[DIR_M0P])[k_M00] = f_P0M;
- (dist.f[DIR_0PP])[k_000] = f_0MM;
- (dist.f[DIR_0MM])[k_0MM] = f_0PP;
- (dist.f[DIR_0PM])[k_00M] = f_0MP;
- (dist.f[DIR_0MP])[k_0M0] = f_0PM;
- (dist.f[DIR_000])[k_000] = f_000;
- (dist.f[DIR_PPP])[k_000] = f_MMM;
- (dist.f[DIR_PMP])[k_0M0] = f_MPM;
- (dist.f[DIR_PPM])[k_00M] = f_MMP;
- (dist.f[DIR_PMM])[k_0MM] = f_MPP;
- (dist.f[DIR_MPP])[k_M00] = f_PMM;
- (dist.f[DIR_MMP])[k_MM0] = f_PPM;
- (dist.f[DIR_MPM])[k_M0M] = f_PMP;
- (dist.f[DIR_MMM])[k_MMM] = f_PPP;
+ distribution[DIR_P00] = f_P00;
+ distribution[DIR_M00] = f_M00;
+ distribution[DIR_0P0] = f_0P0;
+ distribution[DIR_0M0] = f_0M0;
+ distribution[DIR_00P] = f_00P;
+ distribution[DIR_00M] = f_00M;
+ distribution[DIR_PP0] = f_PP0;
+ distribution[DIR_MM0] = f_MM0;
+ distribution[DIR_PM0] = f_PM0;
+ distribution[DIR_MP0] = f_MP0;
+ distribution[DIR_P0P] = f_P0P;
+ distribution[DIR_M0M] = f_M0M;
+ distribution[DIR_P0M] = f_P0M;
+ distribution[DIR_M0P] = f_M0P;
+ distribution[DIR_0PP] = f_0PP;
+ distribution[DIR_0MM] = f_0MM;
+ distribution[DIR_0PM] = f_0PM;
+ distribution[DIR_0MP] = f_0MP;
+ distribution[DIR_000] = f_000;
+ distribution[DIR_PPP] = f_PPP;
+ distribution[DIR_PMP] = f_PMP;
+ distribution[DIR_PPM] = f_PPM;
+ distribution[DIR_PMM] = f_PMM;
+ distribution[DIR_MPP] = f_MPP;
+ distribution[DIR_MMP] = f_MMP;
+ distribution[DIR_MPM] = f_MPM;
+ distribution[DIR_MMM] = f_MMM;
}
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::AMD, true, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::Smagorinsky, true, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::QR, true, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::None, true, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::AMD, true, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::Smagorinsky, true, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::QR, true, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::None, true, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::AMD, false, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::Smagorinsky, false, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::QR, false, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::None, false, true > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::AMD, false, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::Smagorinsky, false, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::QR, false, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
-
-template __global__ void K17CompressibleNavierStokes_Device < TurbulenceModel::None, false, false > ( real omega_in, uint* neighborX, uint* neighborY, uint* neighborZ, real* distributions, real* rho, real* vx, real* vy, real* vz, real* turbulentViscosity, real SGSconstant, unsigned long long numberOfLBnodes, int level, real* forces, real* bodyForceX, real* bodyForceY, real* bodyForceZ, real* quadricLimiters, bool isEvenTimestep, const uint *fluidNodeIndices, uint numberOfFluidNodes);
+} // namespace vf::lbm
diff --git a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityInlines.cuh b/src/lbm/collision/TurbulentViscosity.h
similarity index 58%
rename from src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityInlines.cuh
rename to src/lbm/collision/TurbulentViscosity.h
index ebf67339b65782a5c10c1b756c3fe5e06c3977d1..0122dda7710f0832511afb83103c97ebcda9a384 100644
--- a/src/gpu/VirtualFluids_GPU/GPU/TurbulentViscosityInlines.cuh
+++ b/src/lbm/collision/TurbulentViscosity.h
@@ -34,28 +34,59 @@
#ifndef TURBULENT_VISCOSITY_INLINES_CUH_
#define TURBULENT_VISCOSITY_INLINES_CUH_
-#include <cuda.h>
-#include <cuda_runtime.h>
+#ifndef __host__
+#define __host__
+#endif
+#ifndef __device__
+#define __device__
+#endif
-#include "LBM/LB.h"
+#include <basics/DataTypes.h>
#include <basics/constants/NumericConstants.h>
using namespace vf::basics::constant;
-__inline__ __device__ real calcTurbulentViscositySmagorinsky(real Cs, real dxux, real dyuy, real dzuz, real Dxy, real Dxz , real Dyz)
+namespace vf::lbm
{
- return Cs*Cs * sqrt( c2o1 * ( dxux*dxux + dyuy*dyuy + dzuz*dzuz ) + Dxy*Dxy + Dxz*Dxz + Dyz*Dyz );
+
+//! \brief An enumeration for selecting a turbulence model
+enum class TurbulenceModel {
+ //! - Smagorinsky
+ Smagorinsky,
+ //! - AMD (Anisotropic Minimum Dissipation) model, see e.g. Rozema et al., Phys. Fluids 27, 085107 (2015),
+ //! https://doi.org/10.1063/1.4928700
+ AMD,
+ //! - QR model by Verstappen
+ QR,
+ //! - TODO: move the WALE model here from the old kernels
+ // WALE
+ //! - No turbulence model
+ None
+};
+
+inline __host__ __device__ real calcTurbulentViscositySmagorinsky(real Cs, real dxux, real dyuy, real dzuz, real Dxy, real Dxz,
+ real Dyz)
+{
+ return Cs * Cs * sqrt(c2o1 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz);
}
-__inline__ __device__ real calcTurbulentViscosityQR(real C, real dxux, real dyuy, real dzuz, real Dxy, real Dxz , real Dyz)
+template <typename T>
+__host__ __device__ int max( T a, T b )
{
- // ! Verstappen's QR model
- //! Second invariant of the strain-rate tensor
- real Q = c1o2*( dxux*dxux + dyuy*dyuy + dzuz*dzuz ) + c1o4*( Dxy*Dxy + Dxz*Dxz + Dyz*Dyz);
- //! Third invariant of the strain-rate tensor (determinant)
- // real R = - dxux*dyuy*dzuz - c1o4*( Dxy*Dxz*Dyz + dxux*Dyz*Dyz + dyuy*Dxz*Dxz + dzuz*Dxy*Dxy );
- real R = - dxux*dyuy*dzuz + c1o4*( -Dxy*Dxz*Dyz + dxux*Dyz*Dyz + dyuy*Dxz*Dxz + dzuz*Dxy*Dxy );
- return C * max(R, c0o1) / Q;
+ return ( a > b ) ? a : b;
}
-#endif //TURBULENT_VISCOSITY_H_e
\ No newline at end of file
+inline __host__ __device__ real calcTurbulentViscosityQR(real C, real dxux, real dyuy, real dzuz, real Dxy, real Dxz, real Dyz)
+{
+ // ! Verstappen's QR model
+ //! Second invariant of the strain-rate tensor
+ real Q = c1o2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + c1o4 * (Dxy * Dxy + Dxz * Dxz + Dyz * Dyz);
+ //! Third invariant of the strain-rate tensor (determinant)
+ // real R = - dxux*dyuy*dzuz - c1o4*( Dxy*Dxz*Dyz + dxux*Dyz*Dyz + dyuy*Dxz*Dxz + dzuz*Dxy*Dxy );
+ real R = -dxux * dyuy * dzuz + c1o4 * (-Dxy * Dxz * Dyz + dxux * Dyz * Dyz + dyuy * Dxz * Dxz + dzuz * Dxy * Dxy);
+ return C * max(R, c0o1) / Q;
+}
+
+} // namespace vf::lbm
+
+#endif //TURBULENT_VISCOSITY_H
diff --git a/src/lbm/constants/D3Q27.h b/src/lbm/constants/D3Q27.h
index b6c05eae921ae66b43999ff01977f7a674ce505f..17b663169997bc384ab879a6a6fde6421ec0aee5 100644
--- a/src/lbm/constants/D3Q27.h
+++ b/src/lbm/constants/D3Q27.h
@@ -2,7 +2,8 @@
#define LBM_D3Q27_H
#include <map>
-#include "basics/DataTypes.h"
+
+#include <basics/DataTypes.h>
namespace vf::lbm::dir
{
@@ -93,7 +94,8 @@ static constexpr size_t SGD_MPM = 23;
static constexpr size_t SGD_PMM = 24;
static constexpr size_t SGD_MMM = 25;
-struct countersForPointerChasing{
+struct countersForPointerChasing
+{
uint counterInverse;
uint counterX;
uint counterY;
diff --git a/src/lbm/dummy.cpp b/src/lbm/dummy.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..e69de29bb2d1d6434b8b29ae775ad8c2e48c5391
diff --git a/src/lbm/interpolation/InterpolationCoefficients.h b/src/lbm/interpolation/InterpolationCoefficients.h
index acc1127b7d9c78e74d76ef9187442efdb47dfa40..2244568426c0e0e837f9260ba941accf3fd707fc 100644
--- a/src/lbm/interpolation/InterpolationCoefficients.h
+++ b/src/lbm/interpolation/InterpolationCoefficients.h
@@ -41,7 +41,6 @@
#include "lbm/constants/D3Q27.h"
-#include "lbm/KernelParameter.h"
#include "lbm/MacroscopicQuantities.h"
using namespace vf::basics::constant;
diff --git a/src/lbm/refinement/InterpolationCF.h b/src/lbm/refinement/InterpolationCF.h
index bf59a5d41208997220551e65c8e24feaf49ad38e..fb08d3c12502717c2d8fec1ac81064405bd89cd6 100644
--- a/src/lbm/refinement/InterpolationCF.h
+++ b/src/lbm/refinement/InterpolationCF.h
@@ -41,8 +41,7 @@
#include "lbm/constants/D3Q27.h"
-#include "lbm/KernelParameter.h"
-#include "lbm/Chimera.h"
+#include "lbm/ChimeraTransformation.h"
#include "lbm/interpolation/InterpolationCoefficients.h"
diff --git a/src/lbm/refinement/InterpolationFC.h b/src/lbm/refinement/InterpolationFC.h
index e88c96bac6d7509e139f1e4eb1c4c4eb6017febd..63ed388dbe076136702491d345bc2987703795c9 100644
--- a/src/lbm/refinement/InterpolationFC.h
+++ b/src/lbm/refinement/InterpolationFC.h
@@ -40,8 +40,7 @@
#include <basics/constants/NumericConstants.h>
#include "lbm/constants/D3Q27.h"
-#include "lbm/KernelParameter.h"
-#include "lbm/Chimera.h"
+#include "lbm/ChimeraTransformation.h"
#include "lbm/interpolation/InterpolationCoefficients.h"