diff --git a/.gitlab-ci.yml b/.gitlab-ci.yml
index 1d61280b2bab3e74b5dafc26dec8c62d4b0a5564..62f5f38df39256f74d205f13ff760ea698b57da5 100644
--- a/.gitlab-ci.yml
+++ b/.gitlab-ci.yml
@@ -1,4 +1,4 @@
-image: irmb/virtualfluids-python-deps
+image: irmb/virtualfluids-python-deps:latest
stages:
- build
@@ -97,6 +97,37 @@ gcc_9_rebuild:
paths:
- $CI_PROJECT_DIR/cache
+
+###############################################################################
+gcc_9_cpu_warning_like_errors:
+ stage: build
+
+ image: irmb/virtualfluids-deps-ubuntu20.04
+
+ tags:
+ - gpu
+ - linux
+
+ before_script:
+ - export CCACHE_BASEDIR=$CI_PROJECT_DIR
+ - export CCACHE_DIR=$CI_PROJECT_DIR/cache
+ - ccache -s
+
+ script:
+ - mkdir -p $CI_PROJECT_DIR/build
+ - cd $CI_PROJECT_DIR/build
+ - rm -r -f ./*
+ - cmake ..
+ --preset=cpu_make_ccache
+ -DBUILD_WARNINGS_AS_ERRORS=ON
+ - make -j4
+ - ccache -s
+
+ cache:
+ key: "$CI_JOB_NAME-$CI_COMMIT_REF_SLUG"
+ paths:
+ - $CI_PROJECT_DIR/cache
+
###############################################################################
msvc_16:
stage: build
diff --git a/3rdParty/MuParser/include/muParserDLL.h b/3rdParty/MuParser/include/muParserDLL.h
index b71111c285fbf86fb85092fd405d617364e2f148..14c65b48a463c6bfb52fb92cca2daf62988df937 100644
--- a/3rdParty/MuParser/include/muParserDLL.h
+++ b/3rdParty/MuParser/include/muParserDLL.h
@@ -32,6 +32,7 @@
#ifdef __clang__
#pragma clang system_header
#endif
+
#include "muParserFixes.h"
#ifdef __cplusplus
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
index 7cef8205ce13d89802c9bbd7eebb5d6eb3759b3f..734aa9342df011846ce345ba7c37a181383e8828 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
@@ -47,18 +47,30 @@ class BoundaryConditions;
class BCAlgorithm
{
public:
- static const char VelocityBCAlgorithm = 0;
- static const char EqDensityBCAlgorithm = 1;
- static const char NonEqDensityBCAlgorithm = 2;
- static const char NoSlipBCAlgorithm = 3;
- static const char SlipBCAlgorithm = 4;
- static const char HighViscosityNoSlipBCAlgorithm = 5;
- static const char ThinWallNoSlipBCAlgorithm = 6;
- static const char VelocityWithDensityBCAlgorithm = 7;
+ static const char VelocityBCAlgorithm = 0;
+ static const char EqDensityBCAlgorithm = 1;
+ static const char NonEqDensityBCAlgorithm = 2;
+ static const char NoSlipBCAlgorithm = 3;
+ static const char SlipBCAlgorithm = 4;
+ static const char HighViscosityNoSlipBCAlgorithm = 5;
+ static const char ThinWallNoSlipBCAlgorithm = 6;
+ static const char VelocityWithDensityBCAlgorithm = 7;
static const char NonReflectingOutflowBCAlgorithm = 8;
+ static const char VelocityAndThixotropyBCAlgorithm = 9;
+ static const char DensityAndThixotropyBCAlgorithm = 10;
+ static const char NoSlipAndThixotropyBCAlgorithm = 11;
+ static const char NonReflectingOutflowAndThixotropyBCAlgorithm = 12;
+ static const char VelocityWithDensityAndThixotropyBCAlgorithm = 13;
+ static const char BinghamModelNoSlipBCAlgorithm = 14;
+ static const char HerschelBulkleyModelNoSlipBCAlgorithm = 15;
+ static const char SimpleVelocityBCAlgorithm = 16;
+ static const char SimpleSlipBCAlgorithm = 17;
+ static const char PowellEyringModelNoSlipBCAlgorithm = 18;
+ static const char BinghamModelVelocityBCAlgorithm = 19;
+
public:
- BCAlgorithm() = default;
+ BCAlgorithm() = default;
virtual ~BCAlgorithm() = default;
virtual void addDistributions(SPtr<DistributionArray3D> distributions) = 0;
@@ -72,11 +84,13 @@ public:
SPtr<BCArray3D> getBcArray();
void setBcArray(SPtr<BCArray3D> bcarray);
virtual void applyBC() = 0;
+ bool getThixotropy(){ return thixotropy; };
protected:
- bool compressible{ false };
+ bool compressible { false };
char type;
bool preCollision;
+ bool thixotropy { false };
SPtr<BoundaryConditions> bcPtr;
SPtr<DistributionArray3D> distributions;
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.cpp
index 87606eecf03943259dfec89a805336d2a3190bfa..88f4a52b2ff0445838af5aade25d9c78ce6809a8 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.cpp
@@ -63,11 +63,11 @@ void BCArray3D::resize(std::size_t nx1, std::size_t nx2, std::size_t nx3, int va
//////////////////////////////////////////////////////////////////////////
bool BCArray3D::validIndices(std::size_t x1, std::size_t x2, std::size_t x3) const
{
- if (x1 < 0 || x1 >= this->bcindexmatrix.getNX1())
+ if (x1 >= this->bcindexmatrix.getNX1())
return false;
- if (x2 < 0 || x2 >= this->bcindexmatrix.getNX2())
+ if (x2 >= this->bcindexmatrix.getNX2())
return false;
- if (x3 < 0 || x3 >= this->bcindexmatrix.getNX3())
+ if (x3 >= this->bcindexmatrix.getNX3())
return false;
return true;
}
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.h
index 835e5b1c95454a9fbe8186d6942c3936a5e0e2cc..b9d08f7117d9dc41c008c9d92a5780aceedad21c 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCArray3D.h
@@ -35,7 +35,7 @@
#define BCArray_H
#include "BoundaryConditions.h"
-#include "basics/container/CbArray3D.h"
+#include "CbArray3D.h"
#include <typeinfo>
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BoundaryConditions.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/BoundaryConditions.h
index 1ab3b4e284de49fd1b80bad13d65fee98d221e57..84ba7a6041d38546da03323f35a23f7084df9809 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/BoundaryConditions.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BoundaryConditions.h
@@ -48,7 +48,6 @@ class BoundaryConditions
{
public:
BoundaryConditions()
-
{
UB_STATIC_ASSERT(sizeof(long long) >= 8);
UB_STATIC_ASSERT((sizeof(long long) * 8) >= (D3Q27System::FENDDIR + 1) * BoundaryConditions::optionDigits);
@@ -150,13 +149,13 @@ public:
{
return (short)(((slipBoundaryFlags >> (optionDigits * direction)) & maxOptionVal) - 1);
}
- void setNormalVector(const LBMReal &nx1, const LBMReal &nx2, const LBMReal &nx3)
+ void setNormalVector(const float &nx1, const float &nx2, const float &nx3)
{
this->nx1 = nx1;
this->nx2 = nx2;
this->nx3 = nx3;
}
- UbTupleDouble3 getNormalVector() { return makeUbTuple(nx1, nx2, nx3); }
+ UbTupleFloat3 getNormalVector() { return makeUbTuple(nx1, nx2, nx3); }
/*============== Velocity Boundary ========================*/
void setVelocityBoundaryFlag(const int &direction, const short &secOpt = 0)
@@ -181,72 +180,73 @@ public:
void setBoundaryVelocity(const Vector3D &vx)
{
- setBoundaryVelocityX1((LBMReal)vx[0]);
- setBoundaryVelocityX2((LBMReal)vx[1]);
- setBoundaryVelocityX3((LBMReal)vx[2]);
+ setBoundaryVelocityX1((float)vx[0]);
+ setBoundaryVelocityX2((float)vx[1]);
+ setBoundaryVelocityX3((float)vx[2]);
}
- void setBoundaryVelocityX1(const LBMReal &vx1) { this->bcVelocityX1 = vx1; }
- void setBoundaryVelocityX2(const LBMReal &vx2) { this->bcVelocityX2 = vx2; }
- void setBoundaryVelocityX3(const LBMReal &vx3) { this->bcVelocityX3 = vx3; }
- LBMReal getBoundaryVelocityX1() { return this->bcVelocityX1; }
- LBMReal getBoundaryVelocityX2() { return this->bcVelocityX2; }
- LBMReal getBoundaryVelocityX3() { return this->bcVelocityX3; }
- LBMReal getBoundaryVelocity(const int &direction)
+ void setBoundaryVelocityX1(const float &vx1) { this->bcVelocityX1 = vx1; }
+ void setBoundaryVelocityX2(const float &vx2) { this->bcVelocityX2 = vx2; }
+ void setBoundaryVelocityX3(const float &vx3) { this->bcVelocityX3 = vx3; }
+ float getBoundaryVelocityX1() { return this->bcVelocityX1; }
+ float getBoundaryVelocityX2() { return this->bcVelocityX2; }
+ float getBoundaryVelocityX3() { return this->bcVelocityX3; }
+ float getBoundaryVelocity(const int &direction)
{
switch (direction) {
case D3Q27System::E:
- return (LBMReal)(UbMath::c4o9 *
- (+bcVelocityX1)); //(2/cs^2)(=6)*rho_0(=1 for incompressible)*wi*u*ei with cs=1/sqrt(3)
+ return (float)(UbMath::c4o9 *
+ (+bcVelocityX1)); //(2/cs^2)(=6)*rho_0(=1 bei inkompr)*wi*u*ei mit cs=1/sqrt(3)
case D3Q27System::W:
- return (LBMReal)(UbMath::c4o9 * (-bcVelocityX1));
+ return (float)(UbMath::c4o9 *
+ (-bcVelocityX1)); // z.B. aus paper manfred MRT LB models in three dimensions (2002)
case D3Q27System::N:
- return (LBMReal)(UbMath::c4o9 * (+bcVelocityX2));
+ return (float)(UbMath::c4o9 * (+bcVelocityX2));
case D3Q27System::S:
- return (LBMReal)(UbMath::c4o9 * (-bcVelocityX2));
+ return (float)(UbMath::c4o9 * (-bcVelocityX2));
case D3Q27System::T:
- return (LBMReal)(UbMath::c4o9 * (+bcVelocityX3));
+ return (float)(UbMath::c4o9 * (+bcVelocityX3));
case D3Q27System::B:
- return (LBMReal)(UbMath::c4o9 * (-bcVelocityX3));
+ return (float)(UbMath::c4o9 * (-bcVelocityX3));
case D3Q27System::NE:
- return (LBMReal)(UbMath::c1o9 * (+bcVelocityX1 + bcVelocityX2));
+ return (float)(UbMath::c1o9 * (+bcVelocityX1 + bcVelocityX2));
case D3Q27System::SW:
- return (LBMReal)(UbMath::c1o9 * (-bcVelocityX1 - bcVelocityX2));
+ return (float)(UbMath::c1o9 * (-bcVelocityX1 - bcVelocityX2));
case D3Q27System::SE:
- return (LBMReal)(UbMath::c1o9 * (+bcVelocityX1 - bcVelocityX2));
+ return (float)(UbMath::c1o9 * (+bcVelocityX1 - bcVelocityX2));
case D3Q27System::NW:
- return (LBMReal)(UbMath::c1o9 * (-bcVelocityX1 + bcVelocityX2));
+ return (float)(UbMath::c1o9 * (-bcVelocityX1 + bcVelocityX2));
case D3Q27System::TE:
- return (LBMReal)(UbMath::c1o9 * (+bcVelocityX1 + bcVelocityX3));
+ return (float)(UbMath::c1o9 * (+bcVelocityX1 + bcVelocityX3));
case D3Q27System::BW:
- return (LBMReal)(UbMath::c1o9 * (-bcVelocityX1 - bcVelocityX3));
+ return (float)(UbMath::c1o9 * (-bcVelocityX1 - bcVelocityX3));
case D3Q27System::BE:
- return (LBMReal)(UbMath::c1o9 * (+bcVelocityX1 - bcVelocityX3));
+ return (float)(UbMath::c1o9 * (+bcVelocityX1 - bcVelocityX3));
case D3Q27System::TW:
- return (LBMReal)(UbMath::c1o9 * (-bcVelocityX1 + bcVelocityX3));
+ return (float)(UbMath::c1o9 * (-bcVelocityX1 + bcVelocityX3));
case D3Q27System::TN:
- return (LBMReal)(UbMath::c1o9 * (+bcVelocityX2 + bcVelocityX3));
+ return (float)(UbMath::c1o9 * (+bcVelocityX2 + bcVelocityX3));
case D3Q27System::BS:
- return (LBMReal)(UbMath::c1o9 * (-bcVelocityX2 - bcVelocityX3));
+ return (float)(UbMath::c1o9 * (-bcVelocityX2 - bcVelocityX3));
case D3Q27System::BN:
- return (LBMReal)(UbMath::c1o9 * (+bcVelocityX2 - bcVelocityX3));
+ return (float)(UbMath::c1o9 * (+bcVelocityX2 - bcVelocityX3));
case D3Q27System::TS:
- return (LBMReal)(UbMath::c1o9 * (-bcVelocityX2 + bcVelocityX3));
+ return (float)(UbMath::c1o9 * (-bcVelocityX2 + bcVelocityX3));
case D3Q27System::TNE:
- return (LBMReal)(UbMath::c1o36 * (+bcVelocityX1 + bcVelocityX2 + bcVelocityX3));
+ return (float)(UbMath::c1o36 * (+bcVelocityX1 + bcVelocityX2 + bcVelocityX3));
case D3Q27System::BSW:
- return (LBMReal)(UbMath::c1o36 * (-bcVelocityX1 - bcVelocityX2 - bcVelocityX3));
+ return (float)(UbMath::c1o36 * (-bcVelocityX1 - bcVelocityX2 - bcVelocityX3));
case D3Q27System::BNE:
- return (LBMReal)(UbMath::c1o36 * (+bcVelocityX1 + bcVelocityX2 - bcVelocityX3));
+ return (float)(UbMath::c1o36 * (+bcVelocityX1 + bcVelocityX2 - bcVelocityX3));
case D3Q27System::TSW:
- return (LBMReal)(UbMath::c1o36 * (-bcVelocityX1 - bcVelocityX2 + bcVelocityX3));
+ return (float)(UbMath::c1o36 * (-bcVelocityX1 - bcVelocityX2 + bcVelocityX3));
case D3Q27System::TSE:
- return (LBMReal)(UbMath::c1o36 * (+bcVelocityX1 - bcVelocityX2 + bcVelocityX3));
+ return (float)(UbMath::c1o36 * (+bcVelocityX1 - bcVelocityX2 + bcVelocityX3));
case D3Q27System::BNW:
- return (LBMReal)(UbMath::c1o36 * (-bcVelocityX1 + bcVelocityX2 - bcVelocityX3));
+ return (float)(UbMath::c1o36 * (-bcVelocityX1 + bcVelocityX2 - bcVelocityX3));
case D3Q27System::BSE:
- return (LBMReal)(UbMath::c1o36 * (+bcVelocityX1 - bcVelocityX2 - bcVelocityX3));
+ return (float)(UbMath::c1o36 * (+bcVelocityX1 - bcVelocityX2 - bcVelocityX3));
case D3Q27System::TNW:
- return (LBMReal)(UbMath::c1o36 * (-bcVelocityX1 + bcVelocityX2 + bcVelocityX3));
+ return (float)(UbMath::c1o36 * (-bcVelocityX1 + bcVelocityX2 + bcVelocityX3));
default:
throw UbException(UB_EXARGS, "unknown error");
}
@@ -273,44 +273,44 @@ public:
return (short)(((densityBoundaryFlags >> (optionDigits * direction)) & maxOptionVal) - 1);
}
- void setBoundaryDensity(LBMReal density) { this->bcDensity = density; }
- LBMReal getBoundaryDensity() { return this->bcDensity; }
+ void setBoundaryDensity(float density) { this->bcDensity = density; }
+ float getBoundaryDensity() { return this->bcDensity; }
- // Lodi extension
- void setDensityLodiDensity(const LBMReal &bcLodiDensity) { this->bcLodiDensity = bcLodiDensity; }
- void setDensityLodiVelocityX1(const LBMReal &bcLodiVelocityX1) { this->bcLodiVelocityX1 = bcLodiVelocityX1; }
- void setDensityLodiVelocityX2(const LBMReal &bcLodiVelocityX2) { this->bcLodiVelocityX2 = bcLodiVelocityX2; }
- void setDensityLodiVelocityX3(const LBMReal &bcLodiVelocityX3) { this->bcLodiVelocityX3 = bcLodiVelocityX3; }
- void setDensityLodiLength(const LBMReal &bcLodiLentgh) { this->bcLodiLentgh = bcLodiLentgh; }
- LBMReal getDensityLodiDensity() const { return this->bcLodiDensity; }
- LBMReal getDensityLodiVelocityX1() const { return this->bcLodiVelocityX1; }
- LBMReal getDensityLodiVelocityX2() const { return this->bcLodiVelocityX2; }
- LBMReal getDensityLodiVelocityX3() const { return this->bcLodiVelocityX3; }
- LBMReal getDensityLodiLength() const { return this->bcLodiLentgh; }
+ ////Lodi extension
+ void setDensityLodiDensity(const float &bcLodiDensity) { this->bcLodiDensity = bcLodiDensity; }
+ void setDensityLodiVelocityX1(const float &bcLodiVelocityX1) { this->bcLodiVelocityX1 = bcLodiVelocityX1; }
+ void setDensityLodiVelocityX2(const float &bcLodiVelocityX2) { this->bcLodiVelocityX2 = bcLodiVelocityX2; }
+ void setDensityLodiVelocityX3(const float &bcLodiVelocityX3) { this->bcLodiVelocityX3 = bcLodiVelocityX3; }
+ void setDensityLodiLength(const float &bcLodiLentgh) { this->bcLodiLentgh = bcLodiLentgh; }
+ float getDensityLodiDensity() const { return this->bcLodiDensity; }
+ float getDensityLodiVelocityX1() const { return this->bcLodiVelocityX1; }
+ float getDensityLodiVelocityX2() const { return this->bcLodiVelocityX2; }
+ float getDensityLodiVelocityX3() const { return this->bcLodiVelocityX3; }
+ float getDensityLodiLength() const { return this->bcLodiLentgh; }
- LBMReal &densityLodiDensity() { return this->bcLodiDensity; }
- LBMReal &densityLodiVelocityX1() { return this->bcLodiVelocityX1; }
- LBMReal &densityLodiVelocityX2() { return this->bcLodiVelocityX2; }
- LBMReal &densityLodiVelocityX3() { return this->bcLodiVelocityX3; }
- LBMReal &densityLodiLentgh() { return this->bcLodiLentgh; }
+ float &densityLodiDensity() { return this->bcLodiDensity; }
+ float &densityLodiVelocityX1() { return this->bcLodiVelocityX1; }
+ float &densityLodiVelocityX2() { return this->bcLodiVelocityX2; }
+ float &densityLodiVelocityX3() { return this->bcLodiVelocityX3; }
+ float &densityLodiLentgh() { return this->bcLodiLentgh; }
- const LBMReal &densityLodiDensity() const { return this->bcLodiDensity; }
- const LBMReal &densityLodiVelocityX1() const { return this->bcLodiVelocityX1; }
- const LBMReal &densityLodiVelocityX2() const { return this->bcLodiVelocityX2; }
- const LBMReal &densityLodiVelocityX3() const { return this->bcLodiVelocityX3; }
- const LBMReal &densityLodiLentgh() const { return this->bcLodiLentgh; }
+ const float &densityLodiDensity() const { return this->bcLodiDensity; }
+ const float &densityLodiVelocityX1() const { return this->bcLodiVelocityX1; }
+ const float &densityLodiVelocityX2() const { return this->bcLodiVelocityX2; }
+ const float &densityLodiVelocityX3() const { return this->bcLodiVelocityX3; }
+ const float &densityLodiLentgh() const { return this->bcLodiLentgh; }
/*======================= Qs =============================*/
- void setQ(const LBMReal &val, const int &direction) { q[direction] = val; }
- LBMReal getQ(const int &direction) { return q[direction]; }
+ void setQ(const float &val, const int &direction) { q[direction] = val; }
+ float getQ(const int &direction) { return q[direction]; }
virtual std::vector<std::string> getBCNames()
{
std::vector<std::string> tmp;
- tmp.emplace_back("NoSlipBC");
- tmp.emplace_back("SlipBC");
- tmp.emplace_back("VelocityBC");
- tmp.emplace_back("DensityBC");
+ tmp.push_back("NoSlipBC");
+ tmp.push_back("SlipBC");
+ tmp.push_back("VelocityBC");
+ tmp.push_back("DensityBC");
return tmp;
}
virtual std::vector<long long> getBCFlags()
@@ -336,7 +336,7 @@ public:
static const long long maxOptionVal; // = ( 1<<optionDigits ) - 1; //2^3-1 -> 7
protected:
- LBMReal q[D3Q27System::FENDDIR + 1];
+ float q[D3Q27System::FENDDIR + 1];
long long noslipBoundaryFlags{ 0 };
long long slipBoundaryFlags{ 0 };
@@ -344,20 +344,26 @@ protected:
long long densityBoundaryFlags{ 0 };
long long wallModelBoundaryFlags{ 0 };
- LBMReal bcVelocityX1{ 0.0f };
- LBMReal bcVelocityX2{ 0.0f };
- LBMReal bcVelocityX3{ 0.0f };
- LBMReal bcDensity{ 0.0f };
+ float bcVelocityX1{ 0.0f };
+ float bcVelocityX2{ 0.0f };
+ float bcVelocityX3{ 0.0f };
+ float bcDensity{ 0.0f };
+ // float bcThixotropy{ 0.0f };
+
+ float bcLodiDensity{ 0.0f };
+ float bcLodiVelocityX1{ 0.0f };
+ float bcLodiVelocityX2{ 0.0f };
+ float bcLodiVelocityX3{ 0.0f };
+ float bcLodiLentgh{ 0.0f };
- LBMReal bcLodiDensity{ 0.0f };
- LBMReal bcLodiVelocityX1{ 0.0f };
- LBMReal bcLodiVelocityX2{ 0.0f };
- LBMReal bcLodiVelocityX3{ 0.0f };
- LBMReal bcLodiLentgh{ 0.0f };
+ float nx1{ 0.0f }, nx2{ 0.0f }, nx3{ 0.0f };
- LBMReal nx1{ 0.0f }, nx2{ 0.0f }, nx3{ 0.0f };
+ char algorithmType { -1 };
- char algorithmType{ -1 };
+private:
+ friend class MPIIORestartCoProcessor;
+ friend class MPIIOMigrationCoProcessor;
+ friend class MPIIOMigrationBECoProcessor;
};
#endif
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
index a7ba3a84505f1695a0c59175d2d2d8adf3d4dbdb..e85f2806df40e11c7f30cca1c86bcb5dc639ee73 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteBoundaryConditionsCoProcessor.cpp
@@ -151,8 +151,6 @@ void WriteBoundaryConditionsCoProcessor::addDataGeo(SPtr<Block3D> block)
SPtr<ILBMKernel> kernel = block->getKernel();
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
- // knotennummerierung faengt immer bei 0 an!
- unsigned int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
int minX1 = 0;
int minX2 = 0;
@@ -171,9 +169,9 @@ void WriteBoundaryConditionsCoProcessor::addDataGeo(SPtr<Block3D> block)
maxX2 -= 1;
maxX3 -= 1;
- for (size_t ix3 = minX3; ix3 <= maxX3; ix3++) {
- for (size_t ix2 = minX2; ix2 <= maxX2; ix2++) {
- for (size_t ix1 = minX1; ix1 <= maxX1; ix1++) {
+ for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
+ for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
+ for (int ix1 = minX1; ix1 <= maxX1; ix1++) {
if (!bcArray->isUndefined(ix1, ix2, ix3)) {
int index = 0;
nodeNumbers(ix1, ix2, ix3) = nr++;
@@ -239,6 +237,9 @@ void WriteBoundaryConditionsCoProcessor::addDataGeo(SPtr<Block3D> block)
maxX2 -= 1;
maxX3 -= 1;
+ // knotennummerierung faengt immer bei 0 an!
+ int SWB = 0, SEB = 0, NEB = 0, NWB = 0, SWT = 0, SET = 0, NET = 0, NWT = 0;
+
// cell vector erstellen
for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
for (int ix2 = minX2; ix2 <= maxX2; ix2++) {
@@ -248,7 +249,9 @@ void WriteBoundaryConditionsCoProcessor::addDataGeo(SPtr<Block3D> block)
(SWT = nodeNumbers(ix1, ix2, ix3 + 1)) >= 0 && (SET = nodeNumbers(ix1 + 1, ix2, ix3 + 1)) >= 0 &&
(NET = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0 &&
(NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0) {
- cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
+ cells.push_back(makeUbTuple((unsigned int)SWB, (unsigned int)SEB, (unsigned int)NEB,
+ (unsigned int)NWB, (unsigned int)SWT, (unsigned int)SET,
+ (unsigned int)NET, (unsigned int)NWT));
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
index 27cf92056c88f1481f256ced4080003ed511f241..eef7d3bf6a0d0ad80aa1bdd7d83dfb469b044584 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
@@ -53,20 +53,24 @@ WriteMacroscopicQuantitiesCoProcessor::WriteMacroscopicQuantitiesCoProcessor(SPt
WbWriter *const writer,
SPtr<LBMUnitConverter> conv,
SPtr<Communicator> comm)
- : CoProcessor(grid, s), path(path), writer(writer), conv(conv), comm(comm)
+ : CoProcessor(grid, s), path(path), writer(writer), conv(conv), comm(comm)
{
- gridRank = comm->getProcessID();
+ gridRank = comm->getProcessID();
minInitLevel = this->grid->getCoarsestInitializedLevel();
maxInitLevel = this->grid->getFinestInitializedLevel();
blockVector.resize(maxInitLevel + 1);
- for (int level = minInitLevel; level <= maxInitLevel; level++) {
+ for (int level = minInitLevel; level <= maxInitLevel; level++)
+ {
grid->getBlocks(level, gridRank, true, blockVector[level]);
}
}
+
//////////////////////////////////////////////////////////////////////////
-void WriteMacroscopicQuantitiesCoProcessor::init() {}
+void WriteMacroscopicQuantitiesCoProcessor::init()
+{}
+
//////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::process(double step)
{
@@ -75,14 +79,18 @@ void WriteMacroscopicQuantitiesCoProcessor::process(double step)
UBLOG(logDEBUG3, "WriteMacroscopicQuantitiesCoProcessor::update:" << step);
}
+
//////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
{
int istep = static_cast<int>(step);
- for (int level = minInitLevel; level <= maxInitLevel; level++) {
- for (SPtr<Block3D> block : blockVector[level]) {
- if (block) {
+ for (int level = minInitLevel; level <= maxInitLevel; level++)
+ {
+ for (SPtr<Block3D> block : blockVector[level])
+ {
+ if (block)
+ {
addDataMQ(block);
}
}
@@ -93,27 +101,29 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
subfolder = "mq" + UbSystem::toString(istep);
pfilePath = path + "/mq/" + subfolder;
cfilePath = path + "/mq/mq_collection";
- partPath = pfilePath + "/mq" + UbSystem::toString(gridRank) + "_" + UbSystem::toString(istep);
+ partPath = pfilePath + "/mq" + UbSystem::toString(gridRank) + "_" + UbSystem::toString(istep);
std::string partName = writer->writeOctsWithNodeData(partPath, nodes, cells, datanames, data);
- size_t found = partName.find_last_of("/");
- std::string piece = partName.substr(found + 1);
- piece = subfolder + "/" + piece;
+ size_t found = partName.find_last_of("/");
+ std::string piece = partName.substr(found + 1);
+ piece = subfolder + "/" + piece;
std::vector<std::string> cellDataNames;
std::vector<std::string> pieces;
pieces.push_back(piece);
if (comm->getProcessID() == comm->getRoot()) {
std::string pname =
- WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
+ WbWriterVtkXmlASCII::getInstance()->writeParallelFile(pfilePath, pieces, datanames, cellDataNames);
found = pname.find_last_of("/");
piece = pname.substr(found + 1);
std::vector<std::string> filenames;
filenames.push_back(piece);
- if (step == CoProcessor::scheduler->getMinBegin()) {
+ if (step == CoProcessor::scheduler->getMinBegin())
+ {
WbWriterVtkXmlASCII::getInstance()->writeCollection(cfilePath, filenames, istep, false);
- } else {
+ } else
+ {
WbWriterVtkXmlASCII::getInstance()->addFilesToCollection(cfilePath, filenames, istep, false);
}
UBLOG(logINFO, "WriteMacroscopicQuantitiesCoProcessor step: " << istep);
@@ -121,6 +131,7 @@ void WriteMacroscopicQuantitiesCoProcessor::collectData(double step)
clearData();
}
+
//////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::clearData()
{
@@ -129,16 +140,23 @@ void WriteMacroscopicQuantitiesCoProcessor::clearData()
datanames.clear();
data.clear();
}
+
//////////////////////////////////////////////////////////////////////////
void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
{
- // This data is written:
+ double level = (double)block->getLevel();
+
+ // Diese Daten werden geschrieben:
datanames.resize(0);
- datanames.emplace_back("DRho");
- datanames.emplace_back("Press");
- datanames.emplace_back("Vx");
- datanames.emplace_back("Vy");
- datanames.emplace_back("Vz");
+ datanames.push_back("Rho");
+ datanames.push_back("Vx");
+ datanames.push_back("Vy");
+ datanames.push_back("Vz");
+ // datanames.push_back("Press");
+ datanames.push_back("Level");
+ // datanames.push_back("BlockID");
+ // datanames.push_back("gamma");
+ // datanames.push_back("collFactor");
data.resize(datanames.size());
@@ -146,10 +164,10 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
LBMReal f[D3Q27System::ENDF + 1];
- LBMReal vx1, vx2, vx3, drho, press;
+ LBMReal vx1, vx2, vx3, rho;
- // node numbering always starts at 0!
- unsigned int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
+ // knotennummerierung faengt immer bei 0 an!
+ int SWB, SEB, NEB, NWB, SWT, SET, NET, NWT;
if (block->getKernel()->getCompressible()) {
calcMacros = &D3Q27System::calcCompMacroscopicValues;
@@ -165,12 +183,21 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
int maxX2 = (int)(distributions->getNX2());
int maxX3 = (int)(distributions->getNX3());
- // assign numbers and create node vector + collect data
+ // int minX1 = 1;
+ // int minX2 = 1;
+ // int minX3 = 1;
+
+ // int maxX1 = (int)(distributions->getNX1());
+ // int maxX2 = (int)(distributions->getNX2());
+ // int maxX3 = (int)(distributions->getNX3());
+
+ // nummern vergeben und node vector erstellen + daten sammeln
CbArray3D<int> nodeNumbers((int)maxX1, (int)maxX2, (int)maxX3, -1);
maxX1 -= 2;
maxX2 -= 2;
maxX3 -= 2;
+ // D3Q27BoundaryConditionPtr bcPtr;
int nr = (int)nodes.size();
for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
@@ -180,44 +207,60 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
int index = 0;
nodeNumbers(ix1, ix2, ix3) = nr++;
Vector3D worldCoordinates = grid->getNodeCoordinates(block, ix1, ix2, ix3);
- nodes.emplace_back(float(worldCoordinates[0]), float(worldCoordinates[1]),
- float(worldCoordinates[2]));
+ nodes.push_back(UbTupleFloat3(float(worldCoordinates[0]), float(worldCoordinates[1]),
+ float(worldCoordinates[2])));
distributions->getDistribution(f, ix1, ix2, ix3);
- calcMacros(f, drho, vx1, vx2, vx3);
- press = D3Q27System::calcPress(f, drho, vx1, vx2, vx3);
-
- if (UbMath::isNaN(drho) || UbMath::isInfinity(drho))
- UB_THROW(UbException(
- UB_EXARGS, "drho is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
- block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
- UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+ calcMacros(f, rho, vx1, vx2, vx3);
+ double press = D3Q27System::getPressure(f); // D3Q27System::calcPress(f,rho,vx1,vx2,vx3);
+
+ if (UbMath::isNaN(rho) || UbMath::isInfinity(rho))
+ // UB_THROW( UbException(UB_EXARGS,"rho is not a number (nan or -1.#IND) or infinity number
+ // -1.#INF in block="+block->toString()+",
+ // node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ rho = 999.0;
if (UbMath::isNaN(press) || UbMath::isInfinity(press))
- UB_THROW(UbException(
- UB_EXARGS, "press is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
- block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
- UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+ // UB_THROW( UbException(UB_EXARGS,"press is not a number (nan or -1.#IND) or infinity number
+ // -1.#INF in block="+block->toString()+",
+ // node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ press = 999.0;
if (UbMath::isNaN(vx1) || UbMath::isInfinity(vx1))
- UB_THROW(UbException(
- UB_EXARGS, "vx1 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
- block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
- UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+ // UB_THROW( UbException(UB_EXARGS,"vx1 is not a number (nan or -1.#IND) or infinity number
+ // -1.#INF in block="+block->toString()+",
+ // node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ vx1 = 999.0;
if (UbMath::isNaN(vx2) || UbMath::isInfinity(vx2))
- UB_THROW(UbException(
- UB_EXARGS, "vx2 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
- block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
- UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
+ // UB_THROW( UbException(UB_EXARGS,"vx2 is not a number (nan or -1.#IND) or infinity number
+ // -1.#INF in block="+block->toString()+",
+ // node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ vx2 = 999.0;
if (UbMath::isNaN(vx3) || UbMath::isInfinity(vx3))
- UB_THROW(UbException(
- UB_EXARGS, "vx3 is not a number (nan or -1.#IND) or infinity number -1.#INF in block=" +
- block->toString() + ", node=" + UbSystem::toString(ix1) + "," +
- UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
-
- data[index++].push_back(drho * conv->getFactorDensityLbToW());
- data[index++].push_back(press * conv->getFactorPressureLbToW());
- data[index++].push_back(vx1 * conv->getFactorVelocityLbToW());
- data[index++].push_back(vx2 * conv->getFactorVelocityLbToW());
- data[index++].push_back(vx3 * conv->getFactorVelocityLbToW());
+ // UB_THROW( UbException(UB_EXARGS,"vx3 is not a number (nan or -1.#IND) or infinity number
+ // -1.#INF in block="+block->toString()+",
+ // node="+UbSystem::toString(ix1)+","+UbSystem::toString(ix2)+","+UbSystem::toString(ix3)));
+ vx3 = 999.0;
+
+ data[index++].push_back(rho);
+ data[index++].push_back(vx1);
+ data[index++].push_back(vx2);
+ data[index++].push_back(vx3);
+
+ // shearRate = D3Q27System::getShearRate(f, collFactor);
+
+ // LBMReal collFactorF = BinghamModelLBMKernel::getBinghamCollFactor(collFactor, yieldStress,
+ // shearRate, rho);
+
+ // data[index++].push_back(shearRate);
+ // data[index++].push_back(collFactorF);
+
+ // data[index++].push_back((rho+1.0) * conv->getFactorDensityLbToW() );
+ // data[index++].push_back(vx1 * conv->getFactorVelocityLbToW());
+ // data[index++].push_back(vx2 * conv->getFactorVelocityLbToW());
+ // data[index++].push_back(vx3 * conv->getFactorVelocityLbToW());
+ // data[index++].push_back((press * conv->getFactorPressureLbToW()) / ((rho+1.0) *
+ // conv->getFactorDensityLbToW()));
+ data[index++].push_back(level);
+ // data[index++].push_back(blockID);
}
}
}
@@ -234,7 +277,9 @@ void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
(SWT = nodeNumbers(ix1, ix2, ix3 + 1)) >= 0 && (SET = nodeNumbers(ix1 + 1, ix2, ix3 + 1)) >= 0 &&
(NET = nodeNumbers(ix1 + 1, ix2 + 1, ix3 + 1)) >= 0 &&
(NWT = nodeNumbers(ix1, ix2 + 1, ix3 + 1)) >= 0) {
- cells.push_back(makeUbTuple(SWB, SEB, NEB, NWB, SWT, SET, NET, NWT));
+ cells.push_back(makeUbTuple((unsigned int)SWB, (unsigned int)SEB, (unsigned int)NEB,
+ (unsigned int)NWB, (unsigned int)SWT, (unsigned int)SET,
+ (unsigned int)NET, (unsigned int)NWT));
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
index d7e2286d4599aabf4e7b6e4c9d6a824b38e6d873..6f8a6e74664cf82a550b9000071d4f6beb9ebac2 100644
--- a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
+++ b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.cpp
@@ -47,7 +47,7 @@ D3Q27EsoTwist3DSplittedVector::D3Q27EsoTwist3DSplittedVector(size_t nx1, size_t
this->nonLocalDistributions =
std::make_shared<CbArray4D<LBMReal, IndexerX4X3X2X1>>(13, nx1 + 1, nx2 + 1, nx3 + 1, value);
- this->restDistributions = std::make_shared<CbArray3D<LBMReal, IndexerX3X2X1>>(nx1, nx2, nx3, value);
+ this->zeroDistributions = std::make_shared<CbArray3D<LBMReal, IndexerX3X2X1>>(nx1, nx2, nx3, value);
}
//////////////////////////////////////////////////////////////////////////
D3Q27EsoTwist3DSplittedVector::~D3Q27EsoTwist3DSplittedVector() = default;
@@ -84,7 +84,7 @@ void D3Q27EsoTwist3DSplittedVector::getDistribution(LBMReal *const f, size_t x1,
f[D3Q27System::BNW] = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1);
f[D3Q27System::BNE] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
- f[D3Q27System::REST] = (*this->restDistributions)(x1, x2, x3);
+ f[D3Q27System::ZERO] = (*this->zeroDistributions)(x1, x2, x3);
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setDistribution(const LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -117,7 +117,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistribution(const LBMReal *const f, size
(*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1) = f[D3Q27System::INV_BNW];
(*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f[D3Q27System::INV_BNE];
- (*this->restDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+ (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::ZERO];
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::getDistributionInv(LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -150,7 +150,7 @@ void D3Q27EsoTwist3DSplittedVector::getDistributionInv(LBMReal *const f, size_t
f[D3Q27System::INV_BNW] = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1);
f[D3Q27System::INV_BNE] = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
- f[D3Q27System::REST] = (*this->restDistributions)(x1, x2, x3);
+ f[D3Q27System::ZERO] = (*this->zeroDistributions)(x1, x2, x3);
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setDistributionInv(const LBMReal *const f, size_t x1, size_t x2, size_t x3)
@@ -183,7 +183,7 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionInv(const LBMReal *const f, s
(*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1 + 1, x2, x3 + 1) = f[D3Q27System::BNW];
(*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f[D3Q27System::BNE];
- (*this->restDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+ (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::ZERO];
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(const LBMReal *const f, size_t x1, size_t x2, size_t x3,
@@ -241,8 +241,8 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(const LBMReal *c
(*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3) = f[D3Q27System::BNE];
if ((direction & EsoTwistD3Q27System::etTSW) == EsoTwistD3Q27System::etTSW)
(*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f[D3Q27System::TSW];
- if ((direction & EsoTwistD3Q27System::REST) == EsoTwistD3Q27System::REST)
- (*this->restDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+ if ((direction & EsoTwistD3Q27System::ZERO) == EsoTwistD3Q27System::ZERO)
+ (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::ZERO];
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(LBMReal f, size_t x1, size_t x2, size_t x3,
@@ -327,8 +327,8 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionForDirection(LBMReal f, size_
case D3Q27System::TSW:
(*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f;
break;
- case D3Q27System::REST:
- (*this->restDistributions)(x1, x2, x3) = f;
+ case D3Q27System::ZERO:
+ (*this->zeroDistributions)(x1, x2, x3) = f;
break;
default:
UB_THROW(UbException(UB_EXARGS, "Direction didn't find"));
@@ -390,8 +390,8 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionInvForDirection(const LBMReal
(*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1) = f[D3Q27System::BNE];
if ((direction & EsoTwistD3Q27System::etTSW) == EsoTwistD3Q27System::etTSW)
(*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3) = f[D3Q27System::TSW];
- if ((direction & EsoTwistD3Q27System::REST) == EsoTwistD3Q27System::REST)
- (*this->restDistributions)(x1, x2, x3) = f[D3Q27System::REST];
+ if ((direction & EsoTwistD3Q27System::ZERO) == EsoTwistD3Q27System::ZERO)
+ (*this->zeroDistributions)(x1, x2, x3) = f[D3Q27System::ZERO];
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setDistributionInvForDirection(LBMReal f, size_t x1, size_t x2, size_t x3,
@@ -476,8 +476,8 @@ void D3Q27EsoTwist3DSplittedVector::setDistributionInvForDirection(LBMReal f, si
case D3Q27System::TSW:
(*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3) = f;
break;
- case D3Q27System::REST:
- (*this->restDistributions)(x1, x2, x3) = f;
+ case D3Q27System::ZERO:
+ (*this->zeroDistributions)(x1, x2, x3) = f;
break;
default:
UB_THROW(UbException(UB_EXARGS, "Direction didn't find"));
@@ -539,8 +539,8 @@ LBMReal D3Q27EsoTwist3DSplittedVector::getDistributionForDirection(size_t x1, si
return (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3);
case D3Q27System::BNE:
return (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
- case D3Q27System::REST:
- return (*this->restDistributions)(x1, x2, x3);
+ case D3Q27System::ZERO:
+ return (*this->zeroDistributions)(x1, x2, x3);
default:
UB_THROW(UbException(UB_EXARGS, "Direction didn't find"));
}
@@ -601,8 +601,8 @@ LBMReal D3Q27EsoTwist3DSplittedVector::getDistributionInvForDirection(size_t x1,
return (*this->localDistributions)(D3Q27System::ET_TSW, x1 + 1, x2 + 1, x3);
case D3Q27System::TSW:
return (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3 + 1);
- case D3Q27System::REST:
- return (*this->restDistributions)(x1, x2, x3);
+ case D3Q27System::ZERO:
+ return (*this->zeroDistributions)(x1, x2, x3);
default:
UB_THROW(UbException(UB_EXARGS, "Direction didn't find"));
}
@@ -626,7 +626,7 @@ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr D3Q27EsoTwist3DSplittedVector:
//////////////////////////////////////////////////////////////////////////
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr D3Q27EsoTwist3DSplittedVector::getZeroDistributions()
{
- return this->restDistributions;
+ return this->zeroDistributions;
}
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setNX1(size_t newNX1) { NX1 = newNX1; }
@@ -647,7 +647,7 @@ void D3Q27EsoTwist3DSplittedVector::setNonLocalDistributions(CbArray4D<LBMReal,
//////////////////////////////////////////////////////////////////////////
void D3Q27EsoTwist3DSplittedVector::setZeroDistributions(CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr array)
{
- restDistributions = array;
+ zeroDistributions = array;
}
//////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.h b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.h
index df44457b04a36918643400d60ec5f514e32982ea..1c0d7d05f1392c8c116863e9e0b41000c90ed15e 100644
--- a/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.h
+++ b/src/cpu/VirtualFluidsCore/Data/D3Q27EsoTwist3DSplittedVector.h
@@ -100,7 +100,7 @@ public:
protected:
CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr localDistributions;
CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributions;
- CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr restDistributions;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr zeroDistributions;
size_t NX1, NX2, NX3;
friend class MPIIORestartCoProcessor;
diff --git a/src/cpu/VirtualFluidsCore/Data/DataSet3D.h b/src/cpu/VirtualFluidsCore/Data/DataSet3D.h
index c0171588bcd4a74680326c44db062dba63d4c41c..4930beeab491caf541c436a2e60b326b7cd54c64 100644
--- a/src/cpu/VirtualFluidsCore/Data/DataSet3D.h
+++ b/src/cpu/VirtualFluidsCore/Data/DataSet3D.h
@@ -51,6 +51,9 @@ public:
SPtr<DistributionArray3D> getFdistributions() const;
void setFdistributions(SPtr<DistributionArray3D> distributions);
+ SPtr<DistributionArray3D> getHdistributions() const;
+ void setHdistributions(SPtr<DistributionArray3D> distributions);
+
SPtr<AverageValuesArray3D> getAverageDensity() const;
void setAverageDensity(SPtr<AverageValuesArray3D> values);
@@ -71,10 +74,12 @@ public:
SPtr<RelaxationFactorArray3D> getRelaxationFactor() const;
void setRelaxationFactor(SPtr<RelaxationFactorArray3D> values);
-
protected:
private:
SPtr<DistributionArray3D> fdistributions;
+
+ SPtr<DistributionArray3D> hdistributions;
+
SPtr<AverageValuesArray3D> averageValues;
SPtr<AverageValuesArray3D> averageDensity;
@@ -85,40 +90,96 @@ private:
SPtr<ShearStressValuesArray3D> shearStressValues;
SPtr<RelaxationFactorArray3D> relaxationFactor;
+
};
-inline SPtr<DistributionArray3D> DataSet3D::getFdistributions() const { return fdistributions; }
+inline SPtr<DistributionArray3D> DataSet3D::getFdistributions() const
+{
+ return fdistributions;
+}
-inline void DataSet3D::setFdistributions(SPtr<DistributionArray3D> distributions) { fdistributions = distributions; }
+inline void DataSet3D::setFdistributions(SPtr<DistributionArray3D> distributions)
+{
+ fdistributions = distributions;
+}
-inline SPtr<AverageValuesArray3D> DataSet3D::getAverageValues() const { return averageValues; }
+inline SPtr<DistributionArray3D> DataSet3D::getHdistributions() const
+{
+ return hdistributions;
+}
-inline void DataSet3D::setAverageValues(SPtr<AverageValuesArray3D> values) { averageValues = values; }
+inline void DataSet3D::setHdistributions(SPtr<DistributionArray3D> distributions)
+{
+ hdistributions = distributions;
+}
-inline SPtr<AverageValuesArray3D> DataSet3D::getAverageDensity() const { return averageDensity; }
+inline SPtr<AverageValuesArray3D> DataSet3D::getAverageValues() const
+{
+ return averageValues;
+}
-inline void DataSet3D::setAverageDensity(SPtr<AverageValuesArray3D> values) { averageDensity = values; }
+inline void DataSet3D::setAverageValues(SPtr<AverageValuesArray3D> values)
+{
+ averageValues = values;
+}
-inline SPtr<AverageValuesArray3D> DataSet3D::getAverageVelocity() const { return averageVelocity; }
+inline SPtr<AverageValuesArray3D> DataSet3D::getAverageDensity() const
+{
+ return averageDensity;
+}
-inline void DataSet3D::setAverageVelocity(SPtr<AverageValuesArray3D> values) { averageVelocity = values; }
+inline void DataSet3D::setAverageDensity(SPtr<AverageValuesArray3D> values)
+{
+ averageDensity = values;
+}
-inline SPtr<AverageValuesArray3D> DataSet3D::getAverageFluctuations() const { return averageFluktuations; }
+inline SPtr<AverageValuesArray3D> DataSet3D::getAverageVelocity() const
+{
+ return averageVelocity;
+}
+
+inline void DataSet3D::setAverageVelocity(SPtr<AverageValuesArray3D> values)
+{
+ averageVelocity = values;
+}
-inline void DataSet3D::setAverageFluctuations(SPtr<AverageValuesArray3D> values) { averageFluktuations = values; }
+inline SPtr<AverageValuesArray3D> DataSet3D::getAverageFluctuations() const
+{
+ return averageFluktuations;
+}
+
+inline void DataSet3D::setAverageFluctuations(SPtr<AverageValuesArray3D> values)
+{
+ averageFluktuations = values;
+}
-inline SPtr<AverageValuesArray3D> DataSet3D::getAverageTriplecorrelations() const { return averageTriplecorrelations; }
+inline SPtr<AverageValuesArray3D> DataSet3D::getAverageTriplecorrelations() const
+{
+ return averageTriplecorrelations;
+}
inline void DataSet3D::setAverageTriplecorrelations(SPtr<AverageValuesArray3D> values)
{
averageTriplecorrelations = values;
}
-inline SPtr<ShearStressValuesArray3D> DataSet3D::getShearStressValues() const { return shearStressValues; }
+inline SPtr<ShearStressValuesArray3D> DataSet3D::getShearStressValues() const
+{
+ return shearStressValues;
+}
-inline void DataSet3D::setShearStressValues(SPtr<ShearStressValuesArray3D> values) { shearStressValues = values; }
+inline void DataSet3D::setShearStressValues(SPtr<ShearStressValuesArray3D> values)
+{
+ shearStressValues = values;
+}
-inline SPtr<RelaxationFactorArray3D> DataSet3D::getRelaxationFactor() const { return relaxationFactor; }
+inline SPtr<RelaxationFactorArray3D> DataSet3D::getRelaxationFactor() const
+{
+ return relaxationFactor;
+}
-inline void DataSet3D::setRelaxationFactor(SPtr<RelaxationFactorArray3D> values) { relaxationFactor = values; }
+inline void DataSet3D::setRelaxationFactor(SPtr<RelaxationFactorArray3D> values)
+{
+ relaxationFactor = values;
+}
#endif
diff --git a/src/cpu/VirtualFluidsCore/Data/DistributionArray3D.h b/src/cpu/VirtualFluidsCore/Data/DistributionArray3D.h
index 242c6b7a1c156216a88d4c06c4945040383e39af..8fe4dccea1b53da0513a093e8a741cd0071caf48 100644
--- a/src/cpu/VirtualFluidsCore/Data/DistributionArray3D.h
+++ b/src/cpu/VirtualFluidsCore/Data/DistributionArray3D.h
@@ -41,9 +41,9 @@ class DistributionArray3D
{
public:
DistributionArray3D() = default;
- ;
+
virtual ~DistributionArray3D() = default;
- ;
+
//! get number of nodes for x1 direction
virtual size_t getNX1() const = 0;
//! get number of nodes for x2 direction
diff --git a/src/cpu/VirtualFluidsCore/Data/EsoTwist3D.h b/src/cpu/VirtualFluidsCore/Data/EsoTwist3D.h
index 689589222634edfc52dcfbf358ffc8d32dd1186b..319a9200cc204b0f9b869b2e52353e717a89d783 100644
--- a/src/cpu/VirtualFluidsCore/Data/EsoTwist3D.h
+++ b/src/cpu/VirtualFluidsCore/Data/EsoTwist3D.h
@@ -43,13 +43,23 @@
// Geier, M., & Schönherr, M. (2017). Esoteric twist: an efficient in-place streaming algorithmus for the lattice
// Boltzmann method on massively parallel hardware. Computation, 5(2), 19.
+class EsoTwistD3Q27UnrollArray
+{
+};
+class EsoTwistPlusD3Q27UnrollArray
+{
+};
+class EsoTwistPlusD3Q19UnrollArray
+{
+};
+
class EsoTwist3D : public DistributionArray3D
{
public:
EsoTwist3D() = default;
- ;
+
~EsoTwist3D() override = default;
- ;
+
//////////////////////////////////////////////////////////////////////////
void swap() override = 0;
//////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp
index c456be678449744475a0ac6932850dceb0ee6f1c..1a13aa008ab49a48f1d16c7a2a71ea39dfb191ab 100644
--- a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp
+++ b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.cpp
@@ -35,7 +35,7 @@
// index 0 1 2 3 4 5 6 7 8 9 10 11 12
// 13 14 15 16 17 18 19 20 21 22 23 24 25 26 f: E, W, N, S, T, B, NE, SW, SE, NW, TE, BW, BE, TW, TN,
-// BS, BN, TS, TNE TNW TSE TSW BNE BNW BSE BSW REST
+// BS, BN, TS, TNE TNW TSE TSW BNE BNW BSE BSW ZERO
const int EsoTwistD3Q27System::ETX1[EsoTwistD3Q27System::ENDF + 1] = { 0, 1, 0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1,
0, 0, 0, 0, 0, 1, 0, 1, 0, 1, 0, 1, 0 };
const int EsoTwistD3Q27System::ETX2[EsoTwistD3Q27System::ENDF + 1] = { 0, 0, 0, 1, 0, 0, 0, 1, 0, -1, 0, 0, 0, 0,
@@ -49,7 +49,7 @@ const int EsoTwistD3Q27System::etINVDIR[EsoTwistD3Q27System::ENDF + 1] = {
D3Q27System::INV_TE, D3Q27System::INV_BW, D3Q27System::INV_BE, D3Q27System::INV_TW, D3Q27System::INV_TN,
D3Q27System::INV_BS, D3Q27System::INV_BN, D3Q27System::INV_TS, D3Q27System::INV_TNE, D3Q27System::INV_TNW,
D3Q27System::INV_TSE, D3Q27System::INV_TSW, D3Q27System::INV_BNE, D3Q27System::INV_BNW, D3Q27System::INV_BSE,
- D3Q27System::INV_BSW, D3Q27System::REST
+ D3Q27System::INV_BSW, D3Q27System::ZERO
};
const unsigned long int EsoTwistD3Q27System::etDIR[EsoTwistD3Q27System::ENDF + 1] = {
diff --git a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h
index 21752cc48a84b02bc24cb7efe9e3c5912f476dfd..a9214673ec4b4a66a52fa53f9b625ead0180768b 100644
--- a/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/Data/EsoTwistD3Q27System.h
@@ -47,7 +47,7 @@ struct EsoTwistD3Q27System {
const static int STARTDIR = D3Q27System::STARTDIR;
const static int ENDDIR = D3Q27System::ENDDIR;
- static const int REST = D3Q27System::REST; /*f0 */
+ static const int ZERO = D3Q27System::ZERO; /*f0 */
static const int E = D3Q27System::E; /*f1 */
static const int W = D3Q27System::W; /*f2 */
static const int N = D3Q27System::N; /*f3 */
diff --git a/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp b/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp
index 724d855e3dbfecbe3388ad4071a2a1b1666c1010..f1b2e5ad8c62babaea92cea50a646434f9757cd9 100644
--- a/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp
+++ b/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp
@@ -127,7 +127,7 @@ void BasicCalculator::calculate()
if (refinement) {
if (straightStartLevel < maxInitLevel)
exchangeBlockData(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
time[4] = timer.stop();
UBLOG(logINFO, "refinement exchangeBlockData time = " << time[4]);
@@ -155,14 +155,16 @@ void BasicCalculator::calculate()
} catch (std::exception &e) {
UBLOG(logERROR, e.what());
UBLOG(logERROR, " step = " << calcStep);
- // throw;
- exit(EXIT_FAILURE);
+ // throw e;
+ // exit(EXIT_FAILURE);
} catch (std::string &s) {
UBLOG(logERROR, s);
- exit(EXIT_FAILURE);
+ // exit(EXIT_FAILURE);
+ // throw s;
} catch (...) {
UBLOG(logERROR, "unknown exception");
- exit(EXIT_FAILURE);
+ // exit(EXIT_FAILURE);
+ // throw;
}
}
//////////////////////////////////////////////////////////////////////////
@@ -173,28 +175,27 @@ void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calc
#endif
{
SPtr<Block3D> blockTemp;
- try {
- // startLevel bis maxInitLevel
- for (int level = startLevel; level <= maxInitLevel; level++) {
- // timer.resetAndStart();
- // call LBM kernel
- int size = (int)blocks[level].size();
+ // startLevel bis maxInitLevel
+ for (int level = startLevel; level <= maxInitLevel; level++) {
+ // timer.resetAndStart();
+ // call LBM kernel
+ int size = (int)blocks[level].size();
#ifdef _OPENMP
#pragma omp for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
+ for (int i = 0; i < size; i++) {
+ try {
blockTemp = blocks[level][i];
blockTemp->getKernel()->calculate(calcStep);
+ } catch (std::exception &e) {
+ UBLOG(logERROR, e.what());
+ UBLOG(logERROR, blockTemp->toString() << " step = " << calcStep);
+ std::exit(EXIT_FAILURE);
}
- // timer.stop();
- // UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " <<
- // timer.getTotalTime());
}
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- // UBLOG(logERROR, blockTemp->toString()<<" step = "<<calcStep);
- // throw;
- exit(EXIT_FAILURE);
+ // timer.stop();
+ // UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " <<
+ // timer.getTotalTime());
}
}
}
@@ -239,8 +240,13 @@ void BasicCalculator::connectorsPrepareLocal(std::vector<SPtr<Block3DConnector>>
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
for (int i = 0; i < size; i++) {
- connectors[i]->prepareForReceive();
- connectors[i]->prepareForSend();
+ try {
+ connectors[i]->prepareForReceive();
+ connectors[i]->prepareForSend();
+ } catch (std::exception &e) {
+ UBLOG(logERROR, e.what());
+ std::exit(EXIT_FAILURE);
+ }
}
}
//////////////////////////////////////////////////////////////////////////
@@ -251,8 +257,13 @@ void BasicCalculator::connectorsSendLocal(std::vector<SPtr<Block3DConnector>> &c
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
for (int i = 0; i < size; i++) {
- connectors[i]->fillSendVectors();
- connectors[i]->sendVectors();
+ try {
+ connectors[i]->fillSendVectors();
+ connectors[i]->sendVectors();
+ } catch (std::exception &e) {
+ UBLOG(logERROR, e.what());
+ std::exit(EXIT_FAILURE);
+ }
}
}
//////////////////////////////////////////////////////////////////////////
@@ -321,36 +332,43 @@ void BasicCalculator::applyPreCollisionBC(int startLevel, int maxInitLevel)
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
for (int i = 0; i < size; i++) {
- blocks[level][i]->getKernel()->getBCProcessor()->applyPreCollisionBC();
+ try {
+ blocks[level][i]->getKernel()->getBCProcessor()->applyPreCollisionBC();
+ } catch (std::exception &e) {
+ UBLOG(logERROR, e.what());
+ exit(EXIT_FAILURE);
+ } catch (std::string &s) {
+ UBLOG(logERROR, s);
+ exit(EXIT_FAILURE);
+ } catch (...) {
+ UBLOG(logERROR, "unknown exception");
+ exit(EXIT_FAILURE);
+ }
}
}
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::applyPostCollisionBC(int startLevel, int maxInitLevel)
{
- try {
- // from startLevel to maxInitLevel
- for (int level = startLevel; level <= maxInitLevel; level++) {
- int size = (int)blocks[level].size();
+ // from startLevel to maxInitLevel
+ for (int level = startLevel; level <= maxInitLevel; level++) {
+ int size = (int)blocks[level].size();
#ifdef _OPENMP
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
+ for (int i = 0; i < size; i++) {
+ try {
blocks[level][i]->getKernel()->getBCProcessor()->applyPostCollisionBC();
+ } catch (std::exception &e) {
+ UBLOG(logERROR, e.what());
+ exit(EXIT_FAILURE);
+ } catch (std::string &s) {
+ UBLOG(logERROR, s);
+ exit(EXIT_FAILURE);
+ } catch (...) {
+ UBLOG(logERROR, "unknown exception");
+ exit(EXIT_FAILURE);
}
}
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- // UBLOG(logERROR, " step = "<<calcStep);
- // throw;
- exit(EXIT_FAILURE);
- } catch (std::string &s) {
- UBLOG(logERROR, s);
- // throw;
- exit(EXIT_FAILURE);
- } catch (...) {
- UBLOG(logERROR, "unknown exception");
- // throw;
- exit(EXIT_FAILURE);
}
}
diff --git a/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.h b/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.h
index fb6e68ce0a56ed7e407ab23605a1facd83eace52..3ef1f4c712e552ea5d5b5e82306e2bd94d74d7ab 100644
--- a/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.h
+++ b/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.h
@@ -39,7 +39,8 @@
class Block3DConnector;
//! \class BasicCalculator
-//! \brief Class implements basic functionality with OpenMP parallelization for main calculation LBM loop
+//! \brief Class implements basic functionality with MPI + OpenMP parallelization for main calculation loop
+//! \author Konstantin Kutscher
class BasicCalculator : public Calculator
{
diff --git a/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp b/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp
index 5e921c4f20390da51724f8e53caaab23e0270549..79753c144f5cfff831f1d0415e9434c50b11bcea 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp
+++ b/src/cpu/VirtualFluidsCore/Grid/Block3D.cpp
@@ -318,6 +318,13 @@ void Block3D::deleteInterpolationFlag()
interpolationFlagCF = 0;
}
//////////////////////////////////////////////////////////////////////////
+double Block3D::getWorkLoad()
+{
+ double l = kernel->getCalculationTime();
+ l *= static_cast<double>(1 << level);
+ return l;
+}
+//////////////////////////////////////////////////////////////////////////
std::string Block3D::toString()
{
std::stringstream ss;
diff --git a/src/cpu/VirtualFluidsCore/Grid/Block3D.h b/src/cpu/VirtualFluidsCore/Grid/Block3D.h
index 7a4a2aad75825d6e2aabd252f9c480ca3a96c5a5..b2279b069e6ee322023d30419f8eed5c587f95e8 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Block3D.h
+++ b/src/cpu/VirtualFluidsCore/Grid/Block3D.h
@@ -134,6 +134,8 @@ public:
bool hasInterpolationFlagFC(int dir);
bool hasInterpolationFlagFC();
+ double getWorkLoad();
+
std::string toString();
static int getMaxGlobalID() { return counter; }
diff --git a/src/cpu/VirtualFluidsCore/Grid/Calculator.cpp b/src/cpu/VirtualFluidsCore/Grid/Calculator.cpp
index 083b30b9a4080183b2543dfe14f33d3551a093ec..fbeb2de979bb31dfb87441b5cfcfdf3393f0043c 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Calculator.cpp
+++ b/src/cpu/VirtualFluidsCore/Grid/Calculator.cpp
@@ -124,11 +124,11 @@ void Calculator::initRemoteConnectors()
// grid->getBlocks(level, gridRank, true, blockVector);
grid->getBlocks(l, blockVector);
for (SPtr<Block3D> block : blockVector) {
- int l = block->getLevel();
- block->pushBackRemoteSameLevelConnectors(remoteConns[l]);
+ int block_level = block->getLevel();
+ block->pushBackRemoteSameLevelConnectors(remoteConns[block_level]);
- block->pushBackRemoteInterpolationConnectorsCF(remoteInterConnsCF[l]);
- block->pushBackRemoteInterpolationConnectorsFC(remoteInterConnsFC[l]);
+ block->pushBackRemoteInterpolationConnectorsCF(remoteInterConnsCF[block_level]);
+ block->pushBackRemoteInterpolationConnectorsFC(remoteInterConnsFC[block_level]);
}
}
diff --git a/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp b/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp
index c78fda22d5f1b2fe95fcac5a94435f72cea1c6a5..92be5ed5a06e1909a34144cdd0d1b31000309281 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp
+++ b/src/cpu/VirtualFluidsCore/Grid/Grid3D.cpp
@@ -153,6 +153,17 @@ bool Grid3D::deleteBlock(int ix1, int ix2, int ix3, int level)
return false;
}
}
+void Grid3D::deleteBlocks()
+{
+ std::vector<std::vector<SPtr<Block3D>>> blocksVector(25);
+ int minInitLevel = Grid3DSystem::MINLEVEL;
+ int maxInitLevel = Grid3DSystem::MAXLEVEL;
+ for (int level = minInitLevel; level < maxInitLevel; level++) {
+ getBlocks(level, blocksVector[level]);
+ for (SPtr<Block3D> block : blocksVector[level]) // blocks of the current level
+ deleteBlock(block);
+ }
+}
//////////////////////////////////////////////////////////////////////////
void Grid3D::replaceBlock(SPtr<Block3D> block)
{
@@ -1328,7 +1339,7 @@ void Grid3D::getNeighborBlocksForDirectionWithDirZero(int dir, int ix1, int ix2,
case Grid3DSystem::BSW:
this->getNeighborsBottomSouthWest(ix1, ix2, ix3, level, levelDepth, blocks);
break;
- case Grid3DSystem::REST:
+ case Grid3DSystem::ZERO:
this->getNeighborsZero(ix1, ix2, ix3, level, levelDepth, blocks);
break;
default:
diff --git a/src/cpu/VirtualFluidsCore/Grid/Grid3D.h b/src/cpu/VirtualFluidsCore/Grid/Grid3D.h
index 69c5106847bba69b651ee2c9ed84c0616798b1c3..84c821e84b8c98f17e39814a211de3262e75f804 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Grid3D.h
+++ b/src/cpu/VirtualFluidsCore/Grid/Grid3D.h
@@ -74,6 +74,7 @@ public:
void addBlock(SPtr<Block3D> block);
bool deleteBlock(SPtr<Block3D> block);
bool deleteBlock(int ix1, int ix2, int ix3, int level);
+ void deleteBlocks();
void deleteBlocks(const std::vector<int> &ids);
void replaceBlock(SPtr<Block3D> block);
SPtr<Block3D> getBlock(int ix1, int ix2, int ix3, int level) const;
diff --git a/src/cpu/VirtualFluidsCore/Grid/Grid3DSystem.h b/src/cpu/VirtualFluidsCore/Grid/Grid3DSystem.h
index 008e38b88aec6be8695411d51c263084816fd2e1..ee61b8f7327e76a9393d4d3caa13c3a796470c08 100644
--- a/src/cpu/VirtualFluidsCore/Grid/Grid3DSystem.h
+++ b/src/cpu/VirtualFluidsCore/Grid/Grid3DSystem.h
@@ -71,7 +71,7 @@ static const int BNE = 22;
static const int BNW = 23;
static const int BSE = 24;
static const int BSW = 25;
-static const int REST /*f0 */ = 26;
+static const int ZERO /*f0 */ = 26;
static const int ENDDIR = 25;
@@ -104,6 +104,7 @@ static const int INV_BSW = TNE;
extern const int INVDIR[ENDDIR + 1];
+static const int MINLEVEL = 0;
static const int MAXLEVEL = 25;
extern const int EX1[ENDDIR + 1];
diff --git a/src/cpu/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp b/src/cpu/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
index 55ce8d1dd332c51d1acd8e1285ea0e81ec9ab00f..bf1895b930f1c61d36d537319b53fe4b0abcd960 100644
--- a/src/cpu/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
+++ b/src/cpu/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
@@ -219,7 +219,7 @@ void D3Q27Interactor::initInteractor(const double &timeStep)
else
this->unsetTimeDependent();
- Interactor3D::initInteractor(timeStep);
+ updateBlocks();
}
//////////////////////////////////////////////////////////////////////////
void D3Q27Interactor::updateInteractor(const double ×tep)
diff --git a/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.cpp b/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.cpp
index 0127c9c880f03d574f657ebc43e53ccaa4b67c7e..84526c62598b1d718b1f179228ae2a3f51839856 100644
--- a/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.cpp
+++ b/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.cpp
@@ -290,11 +290,10 @@ void Interactor3D::setInactive() { active = false; }
//////////////////////////////////////////////////////////////////////////
bool Interactor3D::isActive() { return active; }
//////////////////////////////////////////////////////////////////////////
-void Interactor3D::initInteractor(const double & /*timeStep*/)
+void Interactor3D::updateBlocks()
{
- // UBLOG(logINFO, "transBlocks.size = "<<transBlocks.size());
-
- for (SPtr<Block3D> block : bcBlocks) {
+ for (SPtr<Block3D> block : bcBlocks)
+ {
this->setDifferencesToGbObject3D(block);
}
}
diff --git a/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.h b/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.h
index 9bf3a03ba179a1da7fa932a209d78e39b6622bcf..74627b76addaf6badaea678d1c4a20b274234b3a 100644
--- a/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.h
+++ b/src/cpu/VirtualFluidsCore/Interactors/Interactor3D.h
@@ -57,7 +57,7 @@ public:
Interactor3D(SPtr<GbObject3D> geoObject3D, SPtr<Grid3D> grid, int type, Interactor3D::Accuracy a);
virtual ~Interactor3D();
- virtual void initInteractor(const double ×tep = 0);
+ virtual void initInteractor(const double ×tep = 0) = 0;
virtual void updateInteractor(const double ×tep = 0) = 0;
void setSolidBlock(SPtr<Block3D> block);
@@ -76,7 +76,7 @@ public:
SPtr<Grid3D> getGrid3D() const { return grid.lock(); }
void setGrid3D(SPtr<Grid3D> grid) { this->grid = grid; }
virtual SPtr<GbObject3D> getGbObject3D() const { return geoObject3D; }
- virtual bool setDifferencesToGbObject3D(const SPtr<Block3D> /*block*//*, const double& x1, const double& x2, const double& x3, const double& blockLengthX1, const double& blockLengthX2, const double& blockLengthX3, const double& timestep=0*/)
+ virtual bool setDifferencesToGbObject3D(const SPtr<Block3D>)
{
// UBLOG(logINFO, "Interactor3D::setDifferencesToGbObject3D()");
return false;
@@ -123,6 +123,8 @@ protected:
bool isBlockCuttingGeoObject(double minX1, double minX2, double minX3, double maxX1, double maxX2, double maxX3,
double delta);
+ void updateBlocks();
+
SPtr<GbObject3D> geoObject3D;
WPtr<Grid3D> grid;
int type;
diff --git a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.cpp
index 9cecabe4653b78df759db47fe3f48a43a090dad5..daed493b9cc1afddbd92acabcd551da0f463ea26 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.cpp
@@ -30,24 +30,24 @@
//! \ingroup LBM
//! \author Konstantin Kutscher, Martin Geier
//=======================================================================================
-
#include "CumulantK17LBMKernel.h"
-#include "BCArray3D.h"
-#include "Block3D.h"
-#include "D3Q27EsoTwist3DSplittedVector.h"
#include "D3Q27System.h"
+#include "D3Q27EsoTwist3DSplittedVector.h"
+#include <cmath>
#include "DataSet3D.h"
#include "LBMKernel.h"
-#include <cmath>
+#include "Block3D.h"
+#include "BCArray3D.h"
#define PROOF_CORRECTNESS
using namespace UbMath;
//////////////////////////////////////////////////////////////////////////
-CumulantK17LBMKernel::CumulantK17LBMKernel() { this->compressible = true; }
-//////////////////////////////////////////////////////////////////////////
-CumulantK17LBMKernel::~CumulantK17LBMKernel(void) = default;
+CumulantK17LBMKernel::CumulantK17LBMKernel()
+{
+ this->compressible = true;
+}
//////////////////////////////////////////////////////////////////////////
void CumulantK17LBMKernel::initDataSet()
{
@@ -75,561 +75,561 @@ SPtr<LBMKernel> CumulantK17LBMKernel::clone()
//////////////////////////////////////////////////////////////////////////
void CumulantK17LBMKernel::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 thredIdx, blockIdx, blockDim and gridDim.
- //!
- using namespace D3Q27System;
- 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 = (1.0 / 3.0) * (1.0 / collFactor - 1.0 / 2.0);
-
- 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->getBCProcessor()->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;
-
- LBMReal 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))
+ //////////////////////////////////////////////////////////////////////////
+ //! 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 thredIdx, 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 = (1.0 / 3.0) * (1.0 / collFactor - 1.0 / 2.0);
+
+ 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->getBCProcessor()->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;
+
+ LBMReal omega = collFactor;
+
+ for (int x3 = minX3; x3 < maxX3; x3++)
+ {
+ for (int x2 = minX2; x2 < maxX2; x2++)
+ {
+ for (int x1 = minX1; x1 < maxX1; x1++)
{
- 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
-
- LBMReal mfcbb = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
- LBMReal mfbcb = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
- LBMReal mfbbc = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
- LBMReal mfccb = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
- LBMReal mfacb = (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3);
- LBMReal mfcbc = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
- LBMReal mfabc = (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3);
- LBMReal mfbcc = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3);
- LBMReal mfbac = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3);
- LBMReal mfccc = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3);
- LBMReal mfacc = (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3);
- LBMReal mfcac = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3);
- LBMReal mfaac = (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3);
-
- LBMReal mfabb = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3);
- LBMReal mfbab = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3);
- LBMReal mfbba = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p);
- LBMReal mfaab = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3);
- LBMReal mfcab = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3);
- LBMReal mfaba = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p);
- LBMReal mfcba = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p);
- LBMReal mfbaa = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p);
- LBMReal mfbca = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p);
- LBMReal mfaaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p);
- LBMReal mfcaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p);
- LBMReal mfaca = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p);
- LBMReal mfcca = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p);
-
- LBMReal mfbbb = (*this->restDistributions)(x1, x2, x3);
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq. (J1)-(J3)
- //! <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>
- //!
- LBMReal drho = ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
- (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
- ((mfabb + mfcbb) + (mfbab + mfbcb)) + (mfbba + mfbbc)) + mfbbb;
-
- LBMReal rho = c1 + drho;
- LBMReal OOrho = c1 / rho;
- ////////////////////////////////////////////////////////////////////////////////////
- LBMReal vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
- (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
- (mfcbb - mfabb)) / rho;
- LBMReal vvy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
- (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
- (mfbcb - mfbab)) / rho;
- LBMReal vvz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
- (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
- (mfbbc - mfbba)) / rho;
- ////////////////////////////////////////////////////////////////////////////////////
- //forcing
- ///////////////////////////////////////////////////////////////////////////////////////////
- if (withForcing)
- {
- muX1 = static_cast<double>(x1 - 1 + ix1 * maxX1);
- muX2 = static_cast<double>(x2 - 1 + ix2 * maxX2);
- muX3 = static_cast<double>(x3 - 1 + ix3 * maxX3);
-
- forcingX1 = muForcingX1.Eval();
- forcingX2 = muForcingX2.Eval();
- forcingX3 = muForcingX3.Eval();
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Add half of the acceleration (body force) to the velocity as in Eq. (42)
- //! <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>
- //!
- vvx += forcingX1 * deltaT * c1o2; // X
- vvy += forcingX2 * deltaT * c1o2; // Y
- vvz += forcingX3 * deltaT * c1o2; // Z
- }
- ////////////////////////////////////////////////////////////////////////////////////
- // calculate the square of velocities for this lattice node
- LBMReal vx2 = vvx * vvx;
- LBMReal vy2 = vvy * vvy;
- LBMReal vz2 = vvz * vvz;
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Set relaxation limiters for third order cumulants to default value \f$ \lambda=0.001 \f$ according to section 6 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>
- //!
- LBMReal wadjust;
- LBMReal qudricLimitP = c1o100;
- LBMReal qudricLimitM = c1o100;
- LBMReal qudricLimitD = c1o100;
- ////////////////////////////////////////////////////////////////////////////////////
- //! - Chimera transform from well conditioned distributions to central moments 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. (6)-(14) 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>
- //!
- ////////////////////////////////////////////////////////////////////////////////////
- // Z - Dir
- forwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36, c1o36);
- forwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9, c1o9);
- forwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36, c1o36);
- forwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9, c1o9);
- forwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
- forwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9, c1o9);
- forwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36, c1o36);
- forwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9, c1o9);
- forwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36, c1o36);
-
- ////////////////////////////////////////////////////////////////////////////////////
- // Y - Dir
- forwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6, c1o6);
- forwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
- forwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18, c1o18);
- forwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
- forwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
- forwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
- forwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6, c1o6);
- forwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
- forwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18, c1o18);
-
- ////////////////////////////////////////////////////////////////////////////////////
- // X - Dir
- forwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1, c1);
- forwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
- forwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3, c1o3);
- forwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
- forwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
- forwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
- forwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3, c1o3);
- forwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
- forwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9, c1o9);
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - 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} \f$, \f$ C_{210}+C_{012} \f$, \f$ 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} \f$, \f$ C_{210}-C_{012} \f$, \f$ 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} \f$, \f$ C_{202} \f$, \f$ C_{022} \f$, \f$ C_{211} \f$, \f$ C_{121} \f$, \f$ 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}\f$, \f$C_{212}\f$, \f$C_{122}\f$: \f$\omega_9=O5=1.0\f$.
- //! - Sixth order cumulant \f$ C_{222}\f$: \f$\omega_{10}=O6=1.0\f$.
- //!
- ////////////////////////////////////////////////////////////
- //2.
- LBMReal OxxPyyPzz = c1;
- ////////////////////////////////////////////////////////////
- //3.
- LBMReal OxyyPxzz = c8 * (-c2 + omega) * ( c1 + c2*omega) / (-c8 - c14*omega + c7*omega*omega);
- LBMReal OxyyMxzz = c8 * (-c2 + omega) * (-c7 + c4*omega) / (c56 - c50*omega + c9*omega*omega);
- LBMReal Oxyz = c24 * (-c2 + omega) * (-c2 - c7*omega + c3*omega*omega) / (c48 + c152*omega - c130*omega*omega + c29*omega*omega*omega);
- ////////////////////////////////////////////////////////////
- //4.
- LBMReal O4 = c1;
- ////////////////////////////////////////////////////////////
- //5.
- LBMReal O5 = c1;
- ////////////////////////////////////////////////////////////
- //6.
- LBMReal O6 = c1;
-
- ////////////////////////////////////////////////////////////////////////////////////
- //! - A and B: parameters for fourth order convergence of the diffusion term according to Eq. (114) and (115)
- //! <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).
- //!
- LBMReal A = (c4 + c2*omega - c3*omega*omega) / (c2 - c7*omega + c5*omega*omega);
- LBMReal B = (c4 + c28*omega - c14*omega*omega) / (c6 - c21*omega + c15*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.
- LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + c2 * mfbba * mfbab) * OOrho;
- LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + c2 * mfbba * mfabb) * OOrho;
- LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + c2 * mfbab * mfabb) * OOrho;
-
- LBMReal CUMcca = mfcca - (((mfcaa * mfaca + c2 * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) * OOrho - c1o9 * (drho * OOrho));
- LBMReal CUMcac = mfcac - (((mfcaa * mfaac + c2 * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) * OOrho - c1o9 * (drho * OOrho));
- LBMReal CUMacc = mfacc - (((mfaac * mfaca + c2 * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) * OOrho - c1o9 * (drho * OOrho));
- ////////////////////////////////////////////////////////////
- //5.
- LBMReal CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + c4 * mfabb * mfbbb + c2 * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) * OOrho;
- LBMReal CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + c4 * mfbab * mfbbb + c2 * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) * OOrho;
- LBMReal CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + c4 * mfbba * mfbbb + c2 * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) * OOrho;
- ////////////////////////////////////////////////////////////
- //6.
- LBMReal CUMccc = mfccc + ((-c4 * mfbbb * mfbbb
- - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
- - c4 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
- - c2 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
- + (c4 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
- + c2 * (mfcaa * mfaca * mfaac)
- + c16 * mfbba * mfbab * mfabb) * OOrho * OOrho
- - c1o3 * (mfacc + mfcac + mfcca) * OOrho
- - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
- + (c2 * (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.
- LBMReal mxxPyyPzz = mfcaa + mfaca + mfaac;
- LBMReal mxxMyy = mfcaa - mfaca;
- LBMReal mxxMzz = mfcaa - mfaac;
- ////////////////////////////////////////////////////////////
- //3.
- LBMReal mxxyPyzz = mfcba + mfabc;
- LBMReal mxxyMyzz = mfcba - mfabc;
-
- LBMReal mxxzPyyz = mfcab + mfacb;
- LBMReal mxxzMyyz = mfcab - mfacb;
-
- LBMReal mxyyPxzz = mfbca + mfbac;
- LBMReal 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.
- //!
- LBMReal Dxy = -c3 * omega * mfbba;
- LBMReal Dxz = -c3 * omega * mfbab;
- LBMReal Dyz = -c3 * omega * mfabb;
- LBMReal dxux = c1o2 * (-omega) * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
- LBMReal dyuy = dxux + omega * c3o2 * mxxMyy;
- LBMReal 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) - c3 * (c1 - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
- mxxMyy += omega * (-mxxMyy) - c3 * (c1 + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
- mxxMzz += omega * (-mxxMzz) - c3 * (c1 + 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
- //! - 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>
- //!
- wadjust = Oxyz + (c1 - Oxyz) * abs(mfbbb) / (abs(mfbbb) + qudricLimitD);
- mfbbb += wadjust * (-mfbbb);
- wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxxyPyzz) / (abs(mxxyPyzz) + qudricLimitP);
- mxxyPyzz += wadjust * (-mxxyPyzz);
- wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxxyMyzz) / (abs(mxxyMyzz) + qudricLimitM);
- mxxyMyzz += wadjust * (-mxxyMyzz);
- wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxxzPyyz) / (abs(mxxzPyyz) + qudricLimitP);
- mxxzPyyz += wadjust * (-mxxzPyyz);
- wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxxzMyyz) / (abs(mxxzMyyz) + qudricLimitM);
- mxxzMyyz += wadjust * (-mxxzMyyz);
- wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxyyPxzz) / (abs(mxyyPxzz) + qudricLimitP);
- mxyyPxzz += wadjust * (-mxyyPxzz);
- wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxyyMxzz) / (abs(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 * (-c2 * mxxMyy + mxxMzz + mxxPyyPzz);
- mfaac = c1o3 * (mxxMyy - c2 * 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 * (c1 / omega - c1o2) * (dyuy + dzuz) * c2o3 * A + (c1 - O4) * (CUMacc);
- CUMcac = -O4 * (c1 / omega - c1o2) * (dxux + dzuz) * c2o3 * A + (c1 - O4) * (CUMcac);
- CUMcca = -O4 * (c1 / omega - c1o2) * (dyuy + dxux) * c2o3 * A + (c1 - O4) * (CUMcca);
- CUMbbc = -O4 * (c1 / omega - c1o2) * Dxy * c1o3 * B + (c1 - O4) * (CUMbbc);
- CUMbcb = -O4 * (c1 / omega - c1o2) * Dxz * c1o3 * B + (c1 - O4) * (CUMbcb);
- CUMcbb = -O4 * (c1 / omega - c1o2) * Dyz * c1o3 * B + (c1 - 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 * ((c3 * mfcaa + c1) * mfabb + c6 * mfbba * mfbab) * OOrho;
- mfbcb = CUMbcb + c1o3 * ((c3 * mfaca + c1) * mfbab + c6 * mfbba * mfabb) * OOrho;
- mfbbc = CUMbbc + c1o3 * ((c3 * mfaac + c1) * mfbba + c6 * mfbab * mfabb) * OOrho;
-
- mfcca = CUMcca + (((mfcaa * mfaca + c2 * mfbba * mfbba) * c9 + c3 * (mfcaa + mfaca)) * OOrho - (drho * OOrho)) * c1o9;
- mfcac = CUMcac + (((mfcaa * mfaac + c2 * mfbab * mfbab) * c9 + c3 * (mfcaa + mfaac)) * OOrho - (drho * OOrho)) * c1o9;
- mfacc = CUMacc + (((mfaac * mfaca + c2 * mfabb * mfabb) * c9 + c3 * (mfaac + mfaca)) * OOrho - (drho * OOrho)) * c1o9;
-
- //////////////////////////////////////////////////////////////////////////
- //5.
- mfbcc = CUMbcc + c1o3 * (c3 * (mfaac * mfbca + mfaca * mfbac + c4 * mfabb * mfbbb + c2 * (mfbab * mfacb + mfbba * mfabc)) + (mfbca + mfbac)) * OOrho;
- mfcbc = CUMcbc + c1o3 * (c3 * (mfaac * mfcba + mfcaa * mfabc + c4 * mfbab * mfbbb + c2 * (mfabb * mfcab + mfbba * mfbac)) + (mfcba + mfabc)) * OOrho;
- mfccb = CUMccb + c1o3 * (c3 * (mfcaa * mfacb + mfaca * mfcab + c4 * mfbba * mfbbb + c2 * (mfbab * mfbca + mfabb * mfcba)) + (mfacb + mfcab)) * OOrho;
-
- //////////////////////////////////////////////////////////////////////////
- //6.
- mfccc = CUMccc - ((-c4 * mfbbb * mfbbb
- - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
- - c4 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
- - c2 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
- + (c4 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
- + c2 * (mfcaa * mfaca * mfaac)
- + c16 * mfbba * mfbab * mfabb) * OOrho * OOrho
- - c1o3 * (mfacc + mfcac + mfcca) * OOrho
- - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
- + (c2 * (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
- backwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1, c1);
- backwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
- backwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3, c1o3);
- backwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
- backwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
- backwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
- backwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3, c1o3);
- backwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
- backwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9, c1o9);
-
- ////////////////////////////////////////////////////////////////////////////////////
- // Y - Dir
- backwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6, c1o6);
- backwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
- backwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18, c1o18);
- backwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
- backwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
- backwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
- backwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6, c1o6);
- backwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
- backwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18, c1o18);
-
- ////////////////////////////////////////////////////////////////////////////////////
- // Z - Dir
- backwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36, c1o36);
- backwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9, c1o9);
- backwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36, c1o36);
- backwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9, c1o9);
- backwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
- backwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9, c1o9);
- backwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36, c1o36);
- backwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9, c1o9);
- backwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36, c1o36);
- ////////////////////////////////////////////////////////////////////////////////////
-
- //////////////////////////////////////////////////////////////////////////
- //proof correctness
- //////////////////////////////////////////////////////////////////////////
+ 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
+
+ LBMReal mfcbb = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
+ LBMReal mfbcb = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
+ LBMReal mfbbc = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
+ LBMReal mfccb = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
+ LBMReal mfacb = (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3);
+ LBMReal mfcbc = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
+ LBMReal mfabc = (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3);
+ LBMReal mfbcc = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3);
+ LBMReal mfbac = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3);
+ LBMReal mfccc = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3);
+ LBMReal mfacc = (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3);
+ LBMReal mfcac = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3);
+ LBMReal mfaac = (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3);
+
+ LBMReal mfabb = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3);
+ LBMReal mfbab = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3);
+ LBMReal mfbba = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p);
+ LBMReal mfaab = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3);
+ LBMReal mfcab = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3);
+ LBMReal mfaba = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p);
+ LBMReal mfcba = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p);
+ LBMReal mfbaa = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p);
+ LBMReal mfbca = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p);
+ LBMReal mfaaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p);
+ LBMReal mfcaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p);
+ LBMReal mfaca = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p);
+ LBMReal mfcca = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p);
+
+ LBMReal mfbbb = (*this->restDistributions)(x1, x2, x3);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq. (J1)-(J3)
+ //! <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>
+ //!
+ LBMReal drho = ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
+ (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
+ ((mfabb + mfcbb) + (mfbab + mfbcb)) + (mfbba + mfbbc)) + mfbbb;
+
+ LBMReal rho = c1 + drho;
+ LBMReal OOrho = c1 / rho;
+ ////////////////////////////////////////////////////////////////////////////////////
+ LBMReal vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ (mfcbb - mfabb)) / rho;
+ LBMReal vvy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ (mfbcb - mfbab)) / rho;
+ LBMReal vvz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ (mfbbc - mfbba)) / rho;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //forcing
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ if (withForcing)
+ {
+ muX1 = static_cast<double>(x1 - 1 + ix1 * maxX1);
+ muX2 = static_cast<double>(x2 - 1 + ix2 * maxX2);
+ muX3 = static_cast<double>(x3 - 1 + ix3 * maxX3);
+
+ forcingX1 = muForcingX1.Eval();
+ forcingX2 = muForcingX2.Eval();
+ forcingX3 = muForcingX3.Eval();
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Add half of the acceleration (body force) to the velocity as in Eq. (42)
+ //! <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>
+ //!
+ vvx += forcingX1 * deltaT * c1o2; // X
+ vvy += forcingX2 * deltaT * c1o2; // Y
+ vvz += forcingX3 * deltaT * c1o2; // Z
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+ // calculate the square of velocities for this lattice node
+ LBMReal vx2 = vvx * vvx;
+ LBMReal vy2 = vvy * vvy;
+ LBMReal vz2 = vvz * vvz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Set relaxation limiters for third order cumulants to default value \f$ \lambda=0.001 \f$ according to section 6 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>
+ //!
+ LBMReal wadjust;
+ LBMReal qudricLimitP = c1o100;
+ LBMReal qudricLimitM = c1o100;
+ LBMReal qudricLimitD = c1o100;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from well conditioned distributions to central moments 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. (6)-(14) 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>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ forwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36, c1o36);
+ forwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36, c1o36);
+ forwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
+ forwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36, c1o36);
+ forwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36, c1o36);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ forwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6, c1o6);
+ forwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
+ forwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18, c1o18);
+ forwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
+ forwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
+ forwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
+ forwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6, c1o6);
+ forwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
+ forwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18, c1o18);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ forwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1, c1);
+ forwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
+ forwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3, c1o3);
+ forwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
+ forwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
+ forwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
+ forwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3, c1o3);
+ forwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
+ forwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9, c1o9);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - 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} \f$, \f$ C_{210}+C_{012} \f$, \f$ 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} \f$, \f$ C_{210}-C_{012} \f$, \f$ 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} \f$, \f$ C_{202} \f$, \f$ C_{022} \f$, \f$ C_{211} \f$, \f$ C_{121} \f$, \f$ 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}\f$, \f$C_{212}\f$, \f$C_{122}\f$: \f$\omega_9=O5=1.0\f$.
+ //! - Sixth order cumulant \f$ C_{222}\f$: \f$\omega_{10}=O6=1.0\f$.
+ //!
+ ////////////////////////////////////////////////////////////
+ //2.
+ LBMReal OxxPyyPzz = c1;
+ ////////////////////////////////////////////////////////////
+ //3.
+ LBMReal OxyyPxzz = c8 * (-c2 + omega) * ( c1 + c2*omega) / (-c8 - c14*omega + c7*omega*omega);
+ LBMReal OxyyMxzz = c8 * (-c2 + omega) * (-c7 + c4*omega) / (c56 - c50*omega + c9*omega*omega);
+ LBMReal Oxyz = c24 * (-c2 + omega) * (-c2 - c7*omega + c3*omega*omega) / (c48 + c152*omega - c130*omega*omega + c29*omega*omega*omega);
+ ////////////////////////////////////////////////////////////
+ //4.
+ LBMReal O4 = c1;
+ ////////////////////////////////////////////////////////////
+ //5.
+ LBMReal O5 = c1;
+ ////////////////////////////////////////////////////////////
+ //6.
+ LBMReal O6 = c1;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - A and B: parameters for fourth order convergence of the diffusion term according to Eq. (114) and (115)
+ //! <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).
+ //!
+ LBMReal A = (c4 + c2*omega - c3*omega*omega) / (c2 - c7*omega + c5*omega*omega);
+ LBMReal B = (c4 + c28*omega - c14*omega*omega) / (c6 - c21*omega + c15*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.
+ LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + c2 * mfbba * mfbab) * OOrho;
+ LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + c2 * mfbba * mfabb) * OOrho;
+ LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + c2 * mfbab * mfabb) * OOrho;
+
+ LBMReal CUMcca = mfcca - (((mfcaa * mfaca + c2 * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) * OOrho - c1o9 * (drho * OOrho));
+ LBMReal CUMcac = mfcac - (((mfcaa * mfaac + c2 * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) * OOrho - c1o9 * (drho * OOrho));
+ LBMReal CUMacc = mfacc - (((mfaac * mfaca + c2 * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) * OOrho - c1o9 * (drho * OOrho));
+ ////////////////////////////////////////////////////////////
+ //5.
+ LBMReal CUMbcc = mfbcc - ((mfaac * mfbca + mfaca * mfbac + c4 * mfabb * mfbbb + c2 * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac)) * OOrho;
+ LBMReal CUMcbc = mfcbc - ((mfaac * mfcba + mfcaa * mfabc + c4 * mfbab * mfbbb + c2 * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc)) * OOrho;
+ LBMReal CUMccb = mfccb - ((mfcaa * mfacb + mfaca * mfcab + c4 * mfbba * mfbbb + c2 * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab)) * OOrho;
+ ////////////////////////////////////////////////////////////
+ //6.
+ LBMReal CUMccc = mfccc + ((-c4 * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - c4 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - c2 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
+ + (c4 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + c2 * (mfcaa * mfaca * mfaac)
+ + c16 * mfbba * mfbab * mfabb) * OOrho * OOrho
+ - c1o3 * (mfacc + mfcac + mfcca) * OOrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
+ + (c2 * (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.
+ LBMReal mxxPyyPzz = mfcaa + mfaca + mfaac;
+ LBMReal mxxMyy = mfcaa - mfaca;
+ LBMReal mxxMzz = mfcaa - mfaac;
+ ////////////////////////////////////////////////////////////
+ //3.
+ LBMReal mxxyPyzz = mfcba + mfabc;
+ LBMReal mxxyMyzz = mfcba - mfabc;
+
+ LBMReal mxxzPyyz = mfcab + mfacb;
+ LBMReal mxxzMyyz = mfcab - mfacb;
+
+ LBMReal mxyyPxzz = mfbca + mfbac;
+ LBMReal 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.
+ //!
+ LBMReal Dxy = -c3 * omega * mfbba;
+ LBMReal Dxz = -c3 * omega * mfbab;
+ LBMReal Dyz = -c3 * omega * mfabb;
+ LBMReal dxux = c1o2 * (-omega) * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
+ LBMReal dyuy = dxux + omega * c3o2 * mxxMyy;
+ LBMReal 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) - c3 * (c1 - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+ mxxMyy += omega * (-mxxMyy) - c3 * (c1 + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+ mxxMzz += omega * (-mxxMzz) - c3 * (c1 + 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
+ //! - 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>
+ //!
+ wadjust = Oxyz + (c1 - Oxyz) * abs(mfbbb) / (abs(mfbbb) + qudricLimitD);
+ mfbbb += wadjust * (-mfbbb);
+ wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxxyPyzz) / (abs(mxxyPyzz) + qudricLimitP);
+ mxxyPyzz += wadjust * (-mxxyPyzz);
+ wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxxyMyzz) / (abs(mxxyMyzz) + qudricLimitM);
+ mxxyMyzz += wadjust * (-mxxyMyzz);
+ wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxxzPyyz) / (abs(mxxzPyyz) + qudricLimitP);
+ mxxzPyyz += wadjust * (-mxxzPyyz);
+ wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxxzMyyz) / (abs(mxxzMyyz) + qudricLimitM);
+ mxxzMyyz += wadjust * (-mxxzMyyz);
+ wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxyyPxzz) / (abs(mxyyPxzz) + qudricLimitP);
+ mxyyPxzz += wadjust * (-mxyyPxzz);
+ wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxyyMxzz) / (abs(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 * (-c2 * mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaac = c1o3 * (mxxMyy - c2 * 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 * (c1 / omega - c1o2) * (dyuy + dzuz) * c2o3 * A + (c1 - O4) * (CUMacc);
+ CUMcac = -O4 * (c1 / omega - c1o2) * (dxux + dzuz) * c2o3 * A + (c1 - O4) * (CUMcac);
+ CUMcca = -O4 * (c1 / omega - c1o2) * (dyuy + dxux) * c2o3 * A + (c1 - O4) * (CUMcca);
+ CUMbbc = -O4 * (c1 / omega - c1o2) * Dxy * c1o3 * B + (c1 - O4) * (CUMbbc);
+ CUMbcb = -O4 * (c1 / omega - c1o2) * Dxz * c1o3 * B + (c1 - O4) * (CUMbcb);
+ CUMcbb = -O4 * (c1 / omega - c1o2) * Dyz * c1o3 * B + (c1 - 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 * ((c3 * mfcaa + c1) * mfabb + c6 * mfbba * mfbab) * OOrho;
+ mfbcb = CUMbcb + c1o3 * ((c3 * mfaca + c1) * mfbab + c6 * mfbba * mfabb) * OOrho;
+ mfbbc = CUMbbc + c1o3 * ((c3 * mfaac + c1) * mfbba + c6 * mfbab * mfabb) * OOrho;
+
+ mfcca = CUMcca + (((mfcaa * mfaca + c2 * mfbba * mfbba) * c9 + c3 * (mfcaa + mfaca)) * OOrho - (drho * OOrho)) * c1o9;
+ mfcac = CUMcac + (((mfcaa * mfaac + c2 * mfbab * mfbab) * c9 + c3 * (mfcaa + mfaac)) * OOrho - (drho * OOrho)) * c1o9;
+ mfacc = CUMacc + (((mfaac * mfaca + c2 * mfabb * mfabb) * c9 + c3 * (mfaac + mfaca)) * OOrho - (drho * OOrho)) * c1o9;
+
+ //////////////////////////////////////////////////////////////////////////
+ //5.
+ mfbcc = CUMbcc + c1o3 * (c3 * (mfaac * mfbca + mfaca * mfbac + c4 * mfabb * mfbbb + c2 * (mfbab * mfacb + mfbba * mfabc)) + (mfbca + mfbac)) * OOrho;
+ mfcbc = CUMcbc + c1o3 * (c3 * (mfaac * mfcba + mfcaa * mfabc + c4 * mfbab * mfbbb + c2 * (mfabb * mfcab + mfbba * mfbac)) + (mfcba + mfabc)) * OOrho;
+ mfccb = CUMccb + c1o3 * (c3 * (mfcaa * mfacb + mfaca * mfcab + c4 * mfbba * mfbbb + c2 * (mfbab * mfbca + mfabb * mfcba)) + (mfacb + mfcab)) * OOrho;
+
+ //////////////////////////////////////////////////////////////////////////
+ //6.
+ mfccc = CUMccc - ((-c4 * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - c4 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - c2 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) * OOrho
+ + (c4 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + c2 * (mfcaa * mfaca * mfaac)
+ + c16 * mfbba * mfbab * mfabb) * OOrho * OOrho
+ - c1o3 * (mfacc + mfcac + mfcca) * OOrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * OOrho
+ + (c2 * (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
+ backwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1, c1);
+ backwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
+ backwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3, c1o3);
+ backwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
+ backwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
+ backwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
+ backwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3, c1o3);
+ backwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
+ backwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9, c1o9);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ backwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6, c1o6);
+ backwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
+ backwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18, c1o18);
+ backwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
+ backwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
+ backwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
+ backwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6, c1o6);
+ backwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
+ backwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18, c1o18);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ backwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36, c1o36);
+ backwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36, c1o36);
+ backwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
+ backwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36, c1o36);
+ backwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36, c1o36);
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ //////////////////////////////////////////////////////////////////////////
+ //proof correctness
+ //////////////////////////////////////////////////////////////////////////
#ifdef PROOF_CORRECTNESS
- LBMReal 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;
- LBMReal dif = drho - drho_post;
+ LBMReal 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;
+ LBMReal dif = drho - drho_post;
#ifdef SINGLEPRECISION
- if (dif > 10.0E-7 || dif < -10.0E-7)
+ if (dif > 10.0E-7 || dif < -10.0E-7)
#else
- if (dif > 10.0E-15 || dif < -10.0E-15)
+ if (dif > 10.0E-15 || dif < -10.0E-15)
#endif
- {
- UB_THROW(UbException(UB_EXARGS, "rho=" + UbSystem::toString(drho) + ", 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)));
- }
+ {
+ UB_THROW(UbException(UB_EXARGS, "rho=" + UbSystem::toString(drho) + ", 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
- ////////////////////////////////////////////////////////////////////////////////////
- //! - 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>
- //!
- (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = mfabb;
- (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = mfbab;
- (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = mfbba;
- (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = mfaab;
- (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab;
- (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = mfaba;
- (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba;
- (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa;
- (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca;
- (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa;
- (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa;
- (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca;
- (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca;
-
- (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb;
- (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb;
- (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc;
- (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb;
- (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb;
- (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac;
- (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac;
- (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
-
- (*this->restDistributions)(x1, x2, x3) = mfbbb;
- //////////////////////////////////////////////////////////////////////////
-
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - 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>
+ //!
+ (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = mfabb;
+ (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = mfbab;
+ (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = mfbba;
+ (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = mfaab;
+ (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab;
+ (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = mfaba;
+ (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba;
+ (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa;
+ (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca;
+ (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa;
+ (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa;
+ (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca;
+ (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca;
+
+ (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb;
+ (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb;
+ (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc;
+ (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb;
+ (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
+
+ (*this->restDistributions)(x1, x2, x3) = mfbbb;
+ //////////////////////////////////////////////////////////////////////////
+
+ }
}
- }
- }
- }
+ }
+ }
}
//////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.h
index ca652cb37953f530d8a49a6917026652e4382a87..10cfd49264bb829eac1fc6b9bedeee3b6eace265 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/CumulantK17LBMKernel.h
@@ -34,12 +34,12 @@
#ifndef CumulantK17LBMKernel_h__
#define CumulantK17LBMKernel_h__
+#include "LBMKernel.h"
#include "BCProcessor.h"
#include "D3Q27System.h"
-#include "LBMKernel.h"
-#include "basics/container/CbArray3D.h"
-#include "basics/container/CbArray4D.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
@@ -52,17 +52,16 @@ class CumulantK17LBMKernel : public LBMKernel
{
public:
CumulantK17LBMKernel();
- ~CumulantK17LBMKernel() override;
+ ~CumulantK17LBMKernel() = default;
void calculate(int step) override;
SPtr<LBMKernel> clone() override;
+ double getCalculationTime() override { return .0; }
protected:
- inline void forwardInverseChimeraWithK(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv, LBMReal v2,
- LBMReal Kinverse, LBMReal K);
- inline void backwardInverseChimeraWithK(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv, LBMReal v2,
- LBMReal Kinverse, LBMReal K);
- inline void forwardChimera(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv, LBMReal v2);
- inline void backwardChimera(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv, LBMReal v2);
+ inline void forwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K);
+ inline void backwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K);
+ inline void forwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
+ inline void backwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
virtual void initDataSet();
LBMReal f[D3Q27System::ENDF + 1];
@@ -80,19 +79,18 @@ protected:
};
////////////////////////////////////////////////////////////////////////////////
-//! \brief forward chimera transformation \ref forwardInverseChimeraWithK
+//! \brief forward chimera transformation \ref forwardInverseChimeraWithK
//! Transformation from distributions to central moments according to Eq. (6)-(14) 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> Modified for lower round-off errors.
+//! <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>
+//! Modified for lower round-off errors.
////////////////////////////////////////////////////////////////////////////////
-inline void CumulantK17LBMKernel::forwardInverseChimeraWithK(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv,
- LBMReal v2, LBMReal Kinverse, LBMReal K)
+inline void CumulantK17LBMKernel::forwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K)
{
using namespace UbMath;
LBMReal m2 = mfa + mfc;
LBMReal m1 = mfc - mfa;
LBMReal m0 = m2 + mfb;
- mfa = m0;
+ mfa = m0;
m0 *= Kinverse;
m0 += c1;
mfb = (m1 * Kinverse - m0 * vv) * K;
@@ -101,50 +99,49 @@ inline void CumulantK17LBMKernel::forwardInverseChimeraWithK(LBMReal &mfa, LBMRe
////////////////////////////////////////////////////////////////////////////////
//! \brief backward chimera transformation \ref backwardInverseChimeraWithK
//! Transformation from central moments to distributions according to Eq. (57)-(65) 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> Modified for lower round-off errors.
+//! <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>
+//! ] Modified for lower round-off errors.
////////////////////////////////////////////////////////////////////////////////
-inline void CumulantK17LBMKernel::backwardInverseChimeraWithK(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv,
- LBMReal v2, LBMReal Kinverse, LBMReal K)
+inline void CumulantK17LBMKernel::backwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K)
{
using namespace UbMath;
LBMReal m0 = (((mfc - mfb) * c1o2 + mfb * vv) * Kinverse + (mfa * Kinverse + c1) * (v2 - vv) * c1o2) * K;
LBMReal m1 = (((mfa - mfc) - c2 * mfb * vv) * Kinverse + (mfa * Kinverse + c1) * (-v2)) * K;
- mfc = (((mfc + mfb) * c1o2 + mfb * vv) * Kinverse + (mfa * Kinverse + c1) * (v2 + vv) * c1o2) * K;
- mfa = m0;
- mfb = m1;
+ mfc = (((mfc + mfb) * c1o2 + mfb * vv) * Kinverse + (mfa * Kinverse + c1) * (v2 + vv) * c1o2) * K;
+ mfa = m0;
+ mfb = m1;
}
////////////////////////////////////////////////////////////////////////////////
-//! \brief forward chimera transformation \ref forwardChimera
+//! \brief forward chimera transformation \ref forwardChimera
//! Transformation from distributions to central moments according to Eq. (6)-(14) 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> for \f$ K_{abc}=0 \f$. This is to avoid unnessary floating point operations. Modified for lower round-off
-//! errors.
+//! <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>
+//! for \f$ K_{abc}=0 \f$. This is to avoid unnessary floating point operations.
+//! Modified for lower round-off errors.
////////////////////////////////////////////////////////////////////////////////
-inline void CumulantK17LBMKernel::forwardChimera(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv, LBMReal v2)
+inline void CumulantK17LBMKernel::forwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2)
{
using namespace UbMath;
LBMReal m1 = (mfa + mfc) + mfb;
LBMReal m2 = mfc - mfa;
- mfc = (mfc + mfa) + (v2 * m1 - c2 * vv * m2);
- mfb = m2 - vv * m1;
- mfa = m1;
+ mfc = (mfc + mfa) + (v2 * m1 - c2 * vv * m2);
+ mfb = m2 - vv * m1;
+ mfa = m1;
}
////////////////////////////////////////////////////////////////////////////////
-//! \brief backward chimera transformation \ref backwardChimera
+//! \brief backward chimera transformation \ref backwardChimera
//! Transformation from central moments to distributions according to Eq. (57)-(65) 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> for \f$ K_{abc}=0 \f$. This is to avoid unnessary floating point operations. Modified for lower round-off
-//! errors.
+//! <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>
+//! for \f$ K_{abc}=0 \f$. This is to avoid unnessary floating point operations.
+//! Modified for lower round-off errors.
////////////////////////////////////////////////////////////////////////////////
-inline void CumulantK17LBMKernel::backwardChimera(LBMReal &mfa, LBMReal &mfb, LBMReal &mfc, LBMReal vv, LBMReal v2)
+inline void CumulantK17LBMKernel::backwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2)
{
using namespace UbMath;
LBMReal ma = (mfc + mfa * (v2 - vv)) * c1o2 + mfb * (vv - c1o2);
LBMReal mb = ((mfa - mfc) - mfa * v2) - c2 * mfb * vv;
- mfc = (mfc + mfa * (v2 + vv)) * c1o2 + mfb * (vv + c1o2);
- mfb = mb;
- mfa = ma;
+ mfc = (mfc + mfa * (v2 + vv)) * c1o2 + mfb * (vv + c1o2);
+ mfb = mb;
+ mfa = ma;
}
#endif // CumulantK17LBMKernel_h__
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
index e6362563e31ad994ccc7b42c26c5e81a56a6e101..e4bea8735887c3ee9237e4fc368554e4e0002b55 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
@@ -1,9 +1,11 @@
#include "D3Q27System.h"
+
namespace D3Q27System
{
using namespace UbMath;
+
// index 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18//falsch
-// f: REST, E, W, N, S, T, B, NE, SW, SE, NW, TE, BW, BE, TW, TN, BS, BN, TS, TNE TNW TSE TSW BNE BNW
+// f: ZERO, E, W, N, S, T, B, NE, SW, SE, NW, TE, BW, BE, TW, TN, BS, BN, TS, TNE TNW TSE TSW BNE BNW
// BSE BSW const int EX1[] = { 0, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1, 1, -1, 1, -1, 0, 0, 0, 0, 1, -1, 1, -1,
// 1, -1, 1, -1 }; const int EX2[] = { 0, 0, 0, 1, -1, 0, 0, 1, -1, -1, 1, 0, 0, 0, 0, 1, -1, 1, -1, 1, 1,
// -1, -1, 1, 1, -1, -1 }; const int EX3[] = { 0, 0, 0, 0, 0, 1, -1, 0, 0, 0, 0, 1, -1, -1, 1, 1, -1, -1,
@@ -152,7 +154,7 @@ const double cNorm[3][ENDDIR] = { { double(DX1[0]),
// const int BNW = 23;
// const int BSE = 24;
// const int BSW = 25;
-// const int REST /*f0 */ = 26;
+// const int ZERO /*f0 */ = 26;
// const int INV_E = W;
// const int INV_W = E;
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
index b3da79948c60e986b1c96f6f4266012dd7bffd23..b5d88d6c3791d716cd0dca567d7aaa803e863536 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
@@ -35,8 +35,8 @@
#define D3Q27SYSTEM_H
#include <cmath>
-#include <iostream>
#include <string>
+#include <iostream>
#include "LBMSystem.h"
#include "UbException.h"
@@ -89,7 +89,7 @@ static const int BNE = 22;
static const int BNW = 23;
static const int BSE = 24;
static const int BSW = 25;
-static const int REST = 26;
+static const int ZERO = 26;
static const int INV_E = W;
static const int INV_W = E;
@@ -155,7 +155,7 @@ static LBMReal getDensity(const LBMReal *const &f /*[27]*/)
return ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
(((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
- ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+ ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[ZERO];
}
/*=====================================================================*/
// ATTENTION: does not apply to all models -> use certificate instead of static! to do
@@ -184,7 +184,7 @@ static void calcDensity(const LBMReal *const &f /*[27]*/, LBMReal &rho)
rho = ((f[TNE] + f[BSW]) + (f[TSE] + f[BNW])) + ((f[BSE] + f[TNW]) + (f[TSW] + f[BNE])) +
(((f[NE] + f[SW]) + (f[SE] + f[NW])) + ((f[TE] + f[BW]) + (f[BE] + f[TW])) +
((f[BN] + f[TS]) + (f[TN] + f[BS]))) +
- ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[REST];
+ ((f[E] + f[W]) + (f[N] + f[S]) + (f[T] + f[B])) + f[ZERO];
}
/*=====================================================================*/
static void calcIncompVelocityX1(const LBMReal *const &f /*[27]*/, LBMReal &vx1)
@@ -279,7 +279,7 @@ static LBMReal getCompFeqForDirection(const int &direction, const LBMReal &drho,
////-----
LBMReal rho = drho + c1;
switch (direction) {
- case REST:
+ case ZERO:
return REAL_CAST(c8o27 * (drho + rho * (-cu_sq)));
case E:
return REAL_CAST(c2o27 * (drho + rho * (3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq)));
@@ -354,7 +354,7 @@ static void calcCompFeq(LBMReal *const &feq /*[27]*/, const LBMReal &drho, const
LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
LBMReal rho = drho + c1;
- feq[REST] = c8o27 * (drho + rho * (-cu_sq));
+ feq[ZERO] = c8o27 * (drho + rho * (-cu_sq));
feq[E] = c2o27 * (drho + rho * (3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq));
feq[W] = c2o27 * (drho + rho * (3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq));
feq[N] = c2o27 * (drho + rho * (3.0 * (vx2) + c9o2 * (vx2) * (vx2)-cu_sq));
@@ -395,7 +395,7 @@ static LBMReal getIncompFeqForDirection(const int &direction, const LBMReal &drh
LBMReal cu_sq = 1.5f * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
switch (direction) {
- case REST:
+ case ZERO:
return REAL_CAST(c8o27 * (drho - cu_sq));
case E:
return REAL_CAST(c2o27 * (drho + 3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq));
@@ -435,25 +435,25 @@ static LBMReal getIncompFeqForDirection(const int &direction, const LBMReal &drh
return REAL_CAST(c1o54 * (drho + 3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
case TNE:
return REAL_CAST(c1o216 *
- (drho + 3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
+ (drho + 3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
case BSW:
return REAL_CAST(
c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
case BNE:
return REAL_CAST(c1o216 *
- (drho + 3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
+ (drho + 3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
case TSW:
return REAL_CAST(
c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
case TSE:
return REAL_CAST(c1o216 *
- (drho + 3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
+ (drho + 3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
case BNW:
return REAL_CAST(
c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
case BSE:
return REAL_CAST(c1o216 *
- (drho + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
+ (drho + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
case TNW:
return REAL_CAST(
c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
@@ -469,7 +469,7 @@ static void calcIncompFeq(LBMReal *const &feq /*[27]*/, const LBMReal &drho, con
LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
- feq[REST] = c8o27 * (drho - cu_sq);
+ feq[ZERO] = c8o27 * (drho - cu_sq);
feq[E] = c2o27 * (drho + 3.0 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq);
feq[W] = c2o27 * (drho + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq);
feq[N] = c2o27 * (drho + 3.0 * (vx2) + c9o2 * (vx2) * (vx2)-cu_sq);
@@ -758,7 +758,7 @@ static inline LBMReal calcPress(const LBMReal *const f, LBMReal rho, LBMReal vx1
c2 * (f[NE] + f[SW] + f[SE] + f[NW] + f[TE] + f[BW] + f[BE] + f[TW] + f[TN] + f[BS] + f[BN] + f[TS]) +
c3 * (f[TNE] + f[TSW] + f[TSE] + f[TNW] + f[BNE] + f[BSW] + f[BSE] + f[BNW]) -
(vx1 * vx1 + vx2 * vx2 + vx3 * vx3)) *
- (c1 - c1o2 * OxxPyyPzz) +
+ (c1 - c1o2 * OxxPyyPzz) +
OxxPyyPzz * c1o2 * (rho)) *
c1o3;
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/ILBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/ILBMKernel.h
index 44d8d3273d6bfb555acb2c994b97dfbeb676f1c3..4dbe8eee09a37c0c220f47619b72bade2e6ec527 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ILBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/ILBMKernel.h
@@ -45,8 +45,9 @@ class ILBMKernel
public:
virtual ~ILBMKernel() = default;
- virtual void calculate(int step) = 0;
- virtual void swapDistributions() = 0;
+ virtual void calculate(int step) = 0;
+ virtual double getCalculationTime() = 0;
+ virtual void swapDistributions() = 0;
virtual bool getCompressible() const = 0;
virtual SPtr<BCProcessor> getBCProcessor() const = 0;
diff --git a/src/cpu/VirtualFluidsCore/LBM/LBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/LBMKernel.h
index bfaf9d31275d12ac2d4795c46787af864329980d..be29589b9b7ab239bece700126ae906795c83977 100644
--- a/src/cpu/VirtualFluidsCore/LBM/LBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/LBMKernel.h
@@ -56,6 +56,9 @@ public:
virtual SPtr<LBMKernel> clone() = 0;
+ void calculate(int step) override = 0;
+ double getCalculationTime() override = 0;
+
void setBCProcessor(SPtr<BCProcessor> bcp) override;
SPtr<BCProcessor> getBCProcessor() const override;
diff --git a/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h b/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h
index 376330314497224f2ac92e61de7ae4bf81589c93..40570cc3847f71a1942791afa7e95145daafb53b 100644
--- a/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h
+++ b/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h
@@ -97,6 +97,14 @@ public:
this->init(refLengthWorld, csWorld, rhoWorld, csWorld, refLengthLb, rhoLb, csLb);
}
+ LBMUnitConverter(int /*dummy*/, double uReal, double uLB, double nuReal, double nuLB)
+ {
+ factorVelocityLbToW = uReal / uLB;
+ factorViscosityLbToW = nuReal / nuLB;
+ factorDensityLbToW = factorViscosityLbToW * factorVelocityLbToW * factorVelocityLbToW;
+ factorPressureLbToW = factorDensityLbToW;
+ }
+
virtual ~LBMUnitConverter() = default;
double getRefRhoLb() { return refRhoLb; }
@@ -132,6 +140,10 @@ public:
double getFactorAccWToLb() { return 1.0 / this->getFactorAccLbToW(); }
double getFactorTimeLbToW(double deltaX) const { return factorTimeWithoutDx * deltaX; }
+ //////////////////////////////////////////////////////////////////////////
+ double getFactorVelocityLbToW2() { return factorVelocityLbToW; }
+ double getFactorDensityLbToW2() { return factorDensityLbToW; }
+ double getFactorPressureLbToW2() { return factorPressureLbToW; }
/*==========================================================*/
friend inline std::ostream &operator<<(std::ostream &os, LBMUnitConverter c)
@@ -199,7 +211,12 @@ protected:
double factorTimeLbToW{ 1.0 };
double factorMassLbToW{ 1.0 };
double refRhoLb{ 1.0 };
- double factorTimeWithoutDx;
+ double factorTimeWithoutDx{ 0.0 };
+
+ double factorVelocityLbToW{ 1.0 };
+ double factorViscosityLbToW{ 1.0 };
+ double factorDensityLbToW{ 1.0 };
+ double factorPressureLbToW{ 1.0 };
};
#endif // LBMUNITCONVERTER_H
diff --git a/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h b/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h
index b476fd7750fb5ef950d13b6157349607b6d4d7c9..670a597cb84bd4e98450dad2743a8100f04497ea 100644
--- a/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h
+++ b/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h
@@ -36,8 +36,8 @@
#if defined(_WIN32) || defined(_WIN64)
#define MEMORYUTIL_WINDOWS
-#include "psapi.h"
#include "windows.h"
+#include "psapi.h"
#pragma comment(lib, "psapi.lib")
#elif defined __APPLE__
#define MEMORYUTIL_APPLE
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
index 8cbf1801f94111f26f5874753271e24873420a1e..c424a6376a62159da2e4b9f73ccf01858fbde521 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
@@ -294,7 +294,7 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
f[BNW] = f_TSE + feq[BNW];
f[BSE] = f_TNW + feq[BSE];
f[BSW] = f_TNE + feq[BSW];
- f[REST] = f_ZERO + feq[REST];
+ f[ZERO] = f_ZERO + feq[ZERO];
// calcFeqsFct(f,rho,vx1,vx2,vx3);
// distributions->setDistribution(f, ix1, ix2, ix3);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp
index 9814ab036621d94833ad9b4baff61866b22eac4c..c0efdcc6135b1e06f766621dd4528f96fa32247d 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp
@@ -37,13 +37,11 @@
#include "Grid3D.h"
#include "Grid3DSystem.h"
-SetConnectorsBlockVisitor::SetConnectorsBlockVisitor(SPtr<Communicator> comm, bool fullConnector, int dirs, LBMReal nu)
- : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), comm(comm), fullConnector(fullConnector), dirs(dirs), nu(nu)
+SetConnectorsBlockVisitor::SetConnectorsBlockVisitor(SPtr<Communicator> comm, bool fullConnector, int dirs, LBMReal nue)
+ : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), comm(comm), fullConnector(fullConnector), dirs(dirs), nue(nue)
{
}
//////////////////////////////////////////////////////////////////////////
-SetConnectorsBlockVisitor::~SetConnectorsBlockVisitor(void) = default;
-//////////////////////////////////////////////////////////////////////////
void SetConnectorsBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
{
if (!block)
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
index 7d209e0524a13ab89a48dfc2179203eb17b95162..f6eb15206371af2ff6106a5c82c6c71eba26fb34 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
@@ -39,6 +39,7 @@
#include "Block3DVisitor.h"
#include "D3Q27System.h"
+
class Grid3D;
class Block3D;
class Communicator;
@@ -48,8 +49,8 @@ class InterpolationProcessor;
class SetConnectorsBlockVisitor : public Block3DVisitor
{
public:
- SetConnectorsBlockVisitor(SPtr<Communicator> comm, bool fullConnector, int dirs, LBMReal nu);
- ~SetConnectorsBlockVisitor() override;
+ SetConnectorsBlockVisitor(SPtr<Communicator> comm, bool fullConnector, int dirs, LBMReal nue);
+ ~SetConnectorsBlockVisitor() = default;
void visit(SPtr<Grid3D> grid, SPtr<Block3D> block) override;
//////////////////////////////////////////////////////////////////////////
protected:
@@ -58,7 +59,7 @@ protected:
bool fullConnector;
int dirs;
int gridRank;
- LBMReal nu;
+ LBMReal nue;
SPtr<InterpolationProcessor> iProcessor;
};
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
index bc4b3d25701789cfb62b107880680a3287486f8b..7dde0c34edf5c16dff22361a6bbc36394a9783ed 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
@@ -1,3 +1,4 @@
+#include "MemoryUtil.h"
//=======================================================================================
// ____ ____ __ ______ __________ __ __ __ __
// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
@@ -43,13 +44,24 @@
#include <utility>
//////////////////////////////////////////////////////////////////////////
-SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, SetKernelBlockVisitor::Action action)
+SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue,
+ SetKernelBlockVisitor::Action action)
: Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(std::move(kernel)), nue(nue), action(action), dataSetFlag(true)
{
}
+
+SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, int &numberOfProcesses,
+ SetKernelBlockVisitor::Action action)
+ : Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(std::move(kernel)), nue(nue), action(action), dataSetFlag(true),
+ numberOfProcesses(numberOfProcesses)
+{
+}
+
//////////////////////////////////////////////////////////////////////////
void SetKernelBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
{
+ throwExceptionIfNotEnoughMemory(grid);
+
if (kernel && (block->getRank() == grid->getRank())) {
LBMReal collFactor = LBMSystem::calcCollisionFactor(nue, block->getLevel());
kernel->setCollisionFactor(collFactor);
@@ -98,3 +110,27 @@ void SetKernelBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
}
void SetKernelBlockVisitor::setNoDataSetFlag(bool flag) { dataSetFlag = flag; }
+
+void SetKernelBlockVisitor::throwExceptionIfNotEnoughMemory(const SPtr<Grid3D> &grid)
+{
+ auto availableMemory = Utilities::getTotalPhysMem();
+ auto requiredMemory = getRequiredPhysicalMemory(grid);
+ if (requiredMemory > availableMemory)
+ throw UbException(UB_EXARGS, "SetKernelBlockVisitor: Not enough memory!!!");
+}
+
+double SetKernelBlockVisitor::getRequiredPhysicalMemory(const SPtr<Grid3D> &grid) const
+{
+ unsigned long long numberOfNodesPerBlockWithGhostLayer;
+ auto numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
+ auto blockNx = grid->getBlockNX();
+ int ghostLayer = 3;
+
+ numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (val<1>(blockNx) + ghostLayer) *
+ (val<2>(blockNx) + ghostLayer) * (val<3>(blockNx) + ghostLayer);
+
+ auto needMemAll =
+ double(numberOfNodesPerBlockWithGhostLayer * (27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
+
+ return needMemAll / double(numberOfProcesses);
+}
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
index c2863555c0490b24e23ae2d99933e62ce94ac467..7ce7c852e2a815bdf0a37f3ef2960e1b5b76e4b4 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
@@ -51,6 +51,10 @@ public:
SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue,
SetKernelBlockVisitor::Action action = SetKernelBlockVisitor::NewKernel);
+
+ SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, int &numberOfProcesses,
+ SetKernelBlockVisitor::Action action = SetKernelBlockVisitor::NewKernel);
+
~SetKernelBlockVisitor() override = default;
void visit(SPtr<Grid3D> grid, SPtr<Block3D> block) override;
@@ -62,6 +66,12 @@ private:
LBMReal nue;
Action action;
bool dataSetFlag;
+
+ int numberOfProcesses{ 1 };
+
+ double getRequiredPhysicalMemory(const SPtr<Grid3D> &grid) const;
+
+ void throwExceptionIfNotEnoughMemory(const SPtr<Grid3D> &grid);
};
#endif
diff --git a/src/gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.h b/src/gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.h
index 4cac8b6b4384265b7c611dae5bd4c279ea254345..b082e6a7402a606b72d08bc28e9b612fa6661974 100644
--- a/src/gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.h
+++ b/src/gpu/GridGenerator/grid/BoundaryConditions/BoundaryCondition.h
@@ -96,7 +96,7 @@ public:
//////////////////////////////////////////////////////////////////////////
-class VelocityBoundaryCondition : public gg::BoundaryCondition
+class VelocityBoundaryCondition : public gg ::BoundaryCondition
{
public:
static SPtr<VelocityBoundaryCondition> make(real vx, real vy, real vz)