From 7e39b5e92503a5d0bbdbd2aac04c7bb3663702a7 Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Thu, 3 Dec 2020 16:14:51 +0100
Subject: [PATCH] add new implementation of Bingham model / fix rheometer setup
---
apps/cpu/Applications.cmake | 1 +
apps/cpu/CouetteFlow/CMakeLists.txt | 10 +
apps/cpu/CouetteFlow/cf1.cfg | 32 +
apps/cpu/CouetteFlow/cflow.cpp | 307 ++++++
apps/cpu/HerschelBulkleyModel/hbf1.cfg | 6 +-
apps/cpu/HerschelBulkleyModel/hbflow.cpp | 10 +-
apps/cpu/rheometer/rheometer.cfg | 23 +-
apps/cpu/rheometer/rheometer.cpp | 87 +-
src/cpu/VirtualFluids.h | 5 +-
.../CalculateTorqueCoProcessor.cpp | 255 +++++
.../CoProcessors/CalculateTorqueCoProcessor.h | 54 ++
.../WriteThixotropyQuantitiesCoProcessor.cpp | 11 +-
src/cpu/VirtualFluidsCore/LBM/D3Q27System.h | 2 +
src/cpu/VirtualFluidsCore/LBM/Thixotropy.h | 36 +-
.../LBM/ThixotropyModelLBMKernel2.cpp | 899 ++++++++++++++++++
.../LBM/ThixotropyModelLBMKernel2.h | 51 +
16 files changed, 1736 insertions(+), 53 deletions(-)
create mode 100644 apps/cpu/CouetteFlow/CMakeLists.txt
create mode 100644 apps/cpu/CouetteFlow/cf1.cfg
create mode 100644 apps/cpu/CouetteFlow/cflow.cpp
create mode 100644 src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
create mode 100644 src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h
create mode 100644 src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
create mode 100644 src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.h
diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index d09703197..fb5b9cdf5 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -67,4 +67,5 @@ add_subdirectory(${APPS_ROOT_CPU}/LaminarTubeFlow)
add_subdirectory(${APPS_ROOT_CPU}/HerschelBulkleySphere)
add_subdirectory(${APPS_ROOT_CPU}/HerschelBulkleyModel)
add_subdirectory(${APPS_ROOT_CPU}/rheometer)
+add_subdirectory(${APPS_ROOT_CPU}/CouetteFlow)
diff --git a/apps/cpu/CouetteFlow/CMakeLists.txt b/apps/cpu/CouetteFlow/CMakeLists.txt
new file mode 100644
index 000000000..45293b4e2
--- /dev/null
+++ b/apps/cpu/CouetteFlow/CMakeLists.txt
@@ -0,0 +1,10 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
+
+########################################################
+## C++ PROJECT ###
+########################################################
+PROJECT(CouetteFlow)
+
+INCLUDE(${APPS_ROOT_CPU}/IncludsList.cmake)
+
+vf_add_library(BUILDTYPE binary DEPENDS VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} FILES cflow.cpp )
\ No newline at end of file
diff --git a/apps/cpu/CouetteFlow/cf1.cfg b/apps/cpu/CouetteFlow/cf1.cfg
new file mode 100644
index 000000000..0fa997611
--- /dev/null
+++ b/apps/cpu/CouetteFlow/cf1.cfg
@@ -0,0 +1,32 @@
+pathname = d:/temp/BinghamNewTest2
+
+numOfThreads = 4
+availMem = 8e9
+logToFile = false
+
+# blocknx = 2 32 32
+# boundingBox = 0.25 32 32
+# deltax = 0.125 #0.25 #0.5 #1
+
+# nuLB = 0.05
+# forcing = 1e-7 #2e-7 #4e-7 #8e-7
+
+# outTime = 10000
+# endTime = 1920000 #480000 #120000 #30000
+
+blocknx = 2 32 2
+boundingBox = 2 32 2
+deltax = 1
+
+#nuLB = 0.05
+forcing = 8e-7
+#n = 0.8
+#tau0 = 3e-6
+
+velocity = 1e-3
+n = 0.4
+Re = 1
+Bn = 0.01
+
+outTime = 1000
+endTime = 50000000
\ No newline at end of file
diff --git a/apps/cpu/CouetteFlow/cflow.cpp b/apps/cpu/CouetteFlow/cflow.cpp
new file mode 100644
index 000000000..f5a3abca7
--- /dev/null
+++ b/apps/cpu/CouetteFlow/cflow.cpp
@@ -0,0 +1,307 @@
+#include <iostream>
+#include <string>
+
+#include <VirtualFluids.h>
+
+using namespace std;
+
+
+void bflow(string configname)
+{
+ try
+ {
+ ConfigurationFile config;
+ config.load(configname);
+
+ string pathname = config.getValue<string>("pathname");
+ int numOfThreads = config.getValue<int>("numOfThreads");
+ vector<int> blocknx = config.getVector<int>("blocknx");
+ vector<double> boundingBox = config.getVector<double>("boundingBox");
+ //double nuLB = config.getValue<double>("nuLB");
+ double endTime = config.getValue<double>("endTime");
+ double outTime = config.getValue<double>("outTime");
+ double availMem = config.getValue<double>("availMem");
+ //int refineLevel = config.getValue<int>("refineLevel");
+ bool logToFile = config.getValue<bool>("logToFile");
+ //double restartStep = config.getValue<double>("restartStep");
+ double deltax = config.getValue<double>("deltax");
+ //double cpStep = config.getValue<double>("cpStep");
+ //double cpStepStart = config.getValue<double>("cpStepStart");
+ //bool newStart = config.getValue<bool>("newStart");
+ double forcing = config.getValue<double>("forcing");
+ //double n = config.getValue<double>("n");
+ //double k = config.getValue<double>("k");
+ //double tau0 = config.getValue<double>("tau0");
+ double velocity = config.getValue<double>("velocity");
+ double n = config.getValue<double>("n");
+ double Re = config.getValue<double>("Re");
+ double Bn = config.getValue<double>("Bn");
+
+ SPtr<Communicator> comm = MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ if (logToFile)
+ {
+#if defined(__unix__)
+ if (myid == 0)
+ {
+ const char* str = pathname.c_str();
+ mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ }
+#endif
+
+ if (myid == 0)
+ {
+ stringstream logFilename;
+ logFilename << pathname + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt";
+ UbLog::output_policy::setStream(logFilename.str());
+ }
+ }
+
+ LBMReal rhoLB = 0.0;
+
+ SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
+
+ //bounding box
+ //double g_minX1 = 0;
+ //double g_minX2 = 0;
+ //double g_minX3 = 0;
+
+ //double g_maxX1 = boundingBox[0];
+ //double g_maxX2 = boundingBox[1];
+ //double g_maxX3 = boundingBox[2];
+
+ double g_minX1 = 0.0;
+ double g_minX2 = -boundingBox[1]/2.0;
+ double g_minX3 = -boundingBox[2]/2.0;
+
+ double g_maxX1 = boundingBox[0];
+ double g_maxX2 = boundingBox[1]/2.0;
+ double g_maxX3 = boundingBox[2]/2.0;
+
+ double blockLength = 3.0 * deltax;
+
+ double h = (g_maxX2) / 2.0;
+ double dex = g_maxX1;
+
+ //LBMReal tau0 = 1.2e-7;
+ //LBMReal k = nuLB;
+ //LBMReal n = 0.4;
+
+
+ double d = boundingBox[1];
+ double U = velocity;
+ double Gamma = U / d;
+
+ double k = 0.05; // (U * d) / (Re * std::pow(Gamma, n - 1));
+ double tau0 = 1e-6;// Bn* k* std::pow(Gamma, n);
+
+ double beta = 14;
+ double c = 10; // 1.0 / 6.0;
+ double mu0 = 1e-4;
+
+ SPtr<Thixotropy> thix = Thixotropy::getInstance();
+ //Herschel-Bulkley
+ thix->setPowerIndex(n);
+ thix->setViscosityParameter(k);
+ thix->setYieldStress(tau0);
+ //Powell-Eyring
+ thix->setBeta(beta);
+ thix->setC(c);
+ thix->setMu0(mu0);
+
+ SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
+ //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new HerschelBulkleyModelNoSlipBCAlgorithm()));
+ //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new PowellEyringModelNoSlipBCAlgorithm()));
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelNoSlipBCAlgorithm()));
+
+ //BS visitor
+ BoundaryConditionsBlockVisitor bcVisitor;
+ bcVisitor.addBC(noSlipBCAdapter);
+
+ SPtr<BCProcessor> bcProc;
+ bcProc = SPtr<BCProcessor>(new BCProcessor());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new PowellEyringModelLBMKernel());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
+ SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
+ kernel->setForcingX1(forcing);
+ kernel->setWithForcing(true);
+ kernel->setBCProcessor(bcProc);
+
+
+ SPtr<Grid3D> grid(new Grid3D(comm));
+ grid->setPeriodicX1(true);
+ grid->setPeriodicX2(false);
+ grid->setPeriodicX3(true);
+ grid->setDeltaX(deltax);
+ grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
+
+ SPtr<GbObject3D> gridCube(new GbCuboid3D(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3));
+ if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+
+
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Parameters:");
+ UBLOG(logINFO, "forcing = " << forcing);
+ UBLOG(logINFO, "rho = " << rhoLB);
+ //UBLOG(logINFO, "nu = " << nuLB);
+ UBLOG(logINFO, "U = " << U);
+ UBLOG(logINFO, "Re = " << (U * d) / (k * std::pow(Gamma, n - 1)));
+ UBLOG(logINFO, "Bn = " << tau0 / (k * std::pow(Gamma, n)));
+ UBLOG(logINFO, "k = " << k);
+ UBLOG(logINFO, "n = " << n);
+ UBLOG(logINFO, "tau0 = " << tau0);
+ UBLOG(logINFO, "deltax = " << deltax);
+ //UBLOG(logINFO, "number of levels = " << refineLevel + 1);
+ UBLOG(logINFO, "numOfThreads = " << numOfThreads);
+ UBLOG(logINFO, "Preprozess - start");
+ }
+
+ //walls
+ //GbCuboid3DPtr addWallZmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_minX3));
+ //if (myid == 0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname + "/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance());
+
+ //GbCuboid3DPtr addWallZmax(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_maxX3, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ //if (myid == 0) GbSystem3D::writeGeoObject(addWallZmax.get(), pathname + "/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr addWallYmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_minX2, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallYmin.get(), pathname + "/geo/addWallYmin", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr addWallYmax(new GbCuboid3D(g_minX1 - blockLength, g_maxX2, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallYmax.get(), pathname + "/geo/addWallYmax", WbWriterVtkXmlASCII::getInstance());
+
+ //wall interactors
+ //SPtr<D3Q27Interactor> addWallZminInt(new D3Q27Interactor(addWallZmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ //SPtr<D3Q27Interactor> addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ SPtr<D3Q27Interactor> addWallYminInt(new D3Q27Interactor(addWallYmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> addWallYmaxInt(new D3Q27Interactor(addWallYmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ ////////////////////////////////////////////
+ //METIS
+ SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
+ ////////////////////////////////////////////
+ /////delete solid blocks
+ if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - start");
+ InteractorsHelper intHelper(grid, metisVisitor);
+ //intHelper.addInteractor(addWallZminInt);
+ //intHelper.addInteractor(addWallZmaxInt);
+ intHelper.addInteractor(addWallYminInt);
+ intHelper.addInteractor(addWallYmaxInt);
+ intHelper.selectBlocks();
+ if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
+ //////////////////////////////////////
+
+ SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppblocks->process(0);
+
+ unsigned long nob = grid->getNumberOfBlocks();
+ int gl = 3;
+ unsigned long nodb = (blocknx[0]) * (blocknx[1]) * (blocknx[2]);
+ unsigned long nod = nob * (blocknx[0]) * (blocknx[1]) * (blocknx[2]);
+ unsigned long nodg = nob * (blocknx[0] + gl) * (blocknx[1] + gl) * (blocknx[1] + gl);
+ double needMemAll = double(nodg * (27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
+ double needMem = needMemAll / double(comm->getNumberOfProcesses());
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Number of blocks = " << nob);
+ UBLOG(logINFO, "Number of nodes = " << nod);
+ int minInitLevel = grid->getCoarsestInitializedLevel();
+ int maxInitLevel = grid->getFinestInitializedLevel();
+ for (int level = minInitLevel; level <= maxInitLevel; level++)
+ {
+ int nobl = grid->getNumberOfBlocks(level);
+ UBLOG(logINFO, "Number of blocks for level " << level << " = " << nobl);
+ UBLOG(logINFO, "Number of nodes for level " << level << " = " << nobl * nodb);
+ }
+ UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
+ UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
+ UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
+ }
+
+ SetKernelBlockVisitor kernelVisitor(kernel, k, availMem, needMem);
+ grid->accept(kernelVisitor);
+
+ //BC
+ intHelper.setBC();
+
+ //initialization of distributions
+ InitDistributionsBlockVisitor initVisitor;
+ grid->accept(initVisitor);
+
+
+ if (myid == 0) UBLOG(logINFO, "Preprozess - end");
+
+
+ //set connectors
+ InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
+ SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, k, iProcessor);
+ grid->accept(setConnsVisitor);
+
+ grid->accept(bcVisitor);
+
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ WriteBoundaryConditionsCoProcessor ppgeo = WriteBoundaryConditionsCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm);
+ ppgeo.process(1);
+
+ SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
+ SPtr<CoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
+
+ //write data for visualization of macroscopic quantities
+ SPtr<UbScheduler> visSch(new UbScheduler(outTime));
+ SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, pathname,
+ WbWriterVtkXmlASCII::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+
+ SPtr<WriteThixotropyQuantitiesCoProcessor> writeThixotropicMQCoProcessor(new WriteThixotropyQuantitiesCoProcessor(grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+
+ SPtr<UbScheduler> stepGhostLayer(new UbScheduler(outTime));
+ SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
+ calculator->addCoProcessor(npr);
+ calculator->addCoProcessor(writeMQCoProcessor);
+ calculator->addCoProcessor(writeThixotropicMQCoProcessor);
+ //calculator->addCoProcessor(migCoProcessor);
+ //calculator->addCoProcessor(restartCoProcessor);
+
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
+ calculator->calculate();
+ if (myid == 0) UBLOG(logINFO, "Simulation-end");
+ }
+ catch (std::exception& e)
+ {
+ cerr << e.what() << endl << flush;
+ }
+ catch (std::string& s)
+ {
+ cerr << s << endl;
+ }
+ catch (...)
+ {
+ cerr << "unknown exception" << endl;
+ }
+
+}
+
+//////////////////////////////////////////////////////////////////////////
+int main(int argc, char* argv[])
+{
+ if (argv != NULL)
+ {
+ if (argv[1] != NULL)
+ {
+ //pflowForcing(string(argv[1]));
+ bflow(string(argv[1]));
+ }
+ else
+ {
+ cout << "Configuration file is missing!" << endl;
+ }
+ }
+
+ return 0;
+}
diff --git a/apps/cpu/HerschelBulkleyModel/hbf1.cfg b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
index 830fe2ae9..2cf4375f0 100644
--- a/apps/cpu/HerschelBulkleyModel/hbf1.cfg
+++ b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
@@ -1,4 +1,4 @@
-pathname = d:/temp/Herschel-BulkleyConvergenceStudy/HerschelBulkley32_n_0.4_Re_1_Bn_0.01_K15
+pathname = d:/temp/BinghamNewOmegaTest
numOfThreads = 4
availMem = 8e9
@@ -28,5 +28,5 @@ n = 0.4
Re = 1
Bn = 0.01
-outTime = 100000
-endTime = 50000000
\ No newline at end of file
+outTime = 1000
+endTime = 30000
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleyModel/hbflow.cpp b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
index 0e5897fa6..bd4879c96 100644
--- a/apps/cpu/HerschelBulkleyModel/hbflow.cpp
+++ b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
@@ -93,8 +93,8 @@ void bflow(string configname)
double U = velocity;
double Gamma = U / d;
- double k = (U * d) / (Re * std::pow(Gamma, n - 1));
- double tau0 = 0;// 1e-12;// Bn* k* std::pow(Gamma, n);
+ double k = 0.05; // (U * d) / (Re * std::pow(Gamma, n - 1));
+ double tau0 = 3e-6;// Bn* k* std::pow(Gamma, n);
double beta = 14;
double c = 10; // 1.0 / 6.0;
@@ -111,8 +111,9 @@ void bflow(string configname)
thix->setMu0(mu0);
SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
- noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new HerschelBulkleyModelNoSlipBCAlgorithm()));
+ //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new HerschelBulkleyModelNoSlipBCAlgorithm()));
//noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new PowellEyringModelNoSlipBCAlgorithm()));
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelNoSlipBCAlgorithm()));
//BS visitor
BoundaryConditionsBlockVisitor bcVisitor;
@@ -121,8 +122,9 @@ void bflow(string configname)
SPtr<BCProcessor> bcProc;
bcProc = SPtr<BCProcessor>(new BCProcessor());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new PowellEyringModelLBMKernel());
- SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
+ SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
kernel->setForcingX1(forcing);
kernel->setWithForcing(true);
kernel->setBCProcessor(bcProc);
diff --git a/apps/cpu/rheometer/rheometer.cfg b/apps/cpu/rheometer/rheometer.cfg
index 5f95613e5..530ece5ff 100644
--- a/apps/cpu/rheometer/rheometer.cfg
+++ b/apps/cpu/rheometer/rheometer.cfg
@@ -1,20 +1,29 @@
-outputPath = d:/temp/rheometer
+outputPath = d:/temp/rheometer/rheometerBingham128_8_test_New_Omega
numOfThreads = 4
availMem = 8e9
logToFile = false
-blocknx = 8 8 2
-boundingBox = 32 32 2
+blocknx = 8 8 1
+boundingBox = 128 128 1
deltax = 1
+#boundingBox = 0.02 0.02 0.00125
+#deltax = 0.000625
+
refineLevel = 0
radius = 5
sphereCenter = 25 25 25
-velocity = 1e-3
-n = 0.3
+#R0=0.02 #m
+#Ri=0.0075 #m
+
+uLB = 8e-4
+N = 20 #rpm
+
+
+n = 1
Re = 10
Bn = 0.01
@@ -25,5 +34,5 @@ restartStep = 100000
cpStart = 100000
cpStep = 100000
-outTime = 1000
-endTime = 100000
\ No newline at end of file
+outTime = 10000
+endTime = 1600000
\ No newline at end of file
diff --git a/apps/cpu/rheometer/rheometer.cpp b/apps/cpu/rheometer/rheometer.cpp
index 69de586fb..2cafcc7d3 100644
--- a/apps/cpu/rheometer/rheometer.cpp
+++ b/apps/cpu/rheometer/rheometer.cpp
@@ -29,10 +29,11 @@ void bflow(string configname)
double cpStep = config.getValue<double>("cpStep");
double cpStart = config.getValue<double>("cpStart");
bool newStart = config.getValue<bool>("newStart");
- double velocity = config.getValue<double>("velocity");
+ double uLB = config.getValue<double>("uLB");
double n = config.getValue<double>("n");
double Re = config.getValue<double>("Re");
double Bn = config.getValue<double>("Bn");
+ double N = config.getValue<double>("N");
vector<double> sphereCenter = config.getVector<double>("sphereCenter");
SPtr<Communicator> comm = MPICommunicator::getInstance();
@@ -58,7 +59,13 @@ void bflow(string configname)
LBMReal rhoLB = 0.0;
- SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
+ //SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
+ double uWorld = (N * PI) / 30.0; //0.0037699111843
+ double rhoWorld = 2350.0; //kg/m^3
+ double R0 = boundingBox[0] * 0.5;
+
+ SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter(deltax, uWorld*R0, rhoWorld, 1.0, uLB));
+ UBLOG(logINFO, conv->toString());
//bounding box
@@ -81,13 +88,15 @@ void bflow(string configname)
double blockLength = 3.0 * deltax;
double d = 2.0 * radius;
- double U = velocity;
+ double U = uLB;
double Gamma = U / d;
- double k = (U * d) / (Re);
+ double muWorld = 20; //Pa*s
+ double k = 0.0015; // muWorld / rhoWorld * conv->getFactorViscosityWToLb(); //(U * d) / (Re);
//double k = (U * d) / (Re * std::pow(Gamma, n - 1));
- double tau0 = Bn * k * std::pow(Gamma, n);
+ double yielStressWorld = 20; //Pa
+ double tau0 = 1e-6;// 3e-6;//yielStressWorld * conv->getFactorPressureWToLb(); //Bn * k * std::pow(Gamma, n);
//double k = 0.05; // (U * d) / (Re * std::pow(Gamma, n - 1));
//double tau0 = 3e-6; //Bn * k * std::pow(Gamma, n);
@@ -103,27 +112,31 @@ void bflow(string configname)
thix->setOmegaMin(omegaMin);
SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
- noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
+ //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
//noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new HerschelBulkleyModelNoSlipBCAlgorithm()));
- //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelNoSlipBCAlgorithm()));
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelNoSlipBCAlgorithm()));
//SPtr<BCAdapter> slipBCAdapter(new SlipBCAdapter());
//slipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new SimpleSlipBCAlgorithm()));
mu::Parser fctVx;
- fctVx.SetExpr("omega*(r-x2)");
- fctVx.DefineConst("omega", velocity);
+ //fctVx.SetExpr("omega*(r-x2)");
+ fctVx.SetExpr("-Omega*(x2-r)");
+ fctVx.DefineConst("Omega", uLB);
+ //fctVx.DefineConst("r", R0);
fctVx.DefineConst("r", g_maxX1*0.5);
mu::Parser fctVy;
- fctVy.SetExpr("omega*(x1-k)");
- fctVy.DefineConst("omega", velocity);
- fctVy.DefineConst("k", g_maxX1 * 0.5);
+ fctVy.SetExpr("Omega*(x1-r)");
+ fctVy.DefineConst("Omega", uLB);
+ //fctVy.DefineConst("r", R0);
+ fctVy.DefineConst("r", g_maxX2 * 0.5);
mu::Parser fctVz;
fctVz.SetExpr("0.0");
SPtr<BCAdapter> velocityBCAdapter(new VelocityBCAdapter(true, true, true, fctVx, fctVy, fctVz, 0, BCFunction::INFCONST));
+ //velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new SimpleVelocityBCAlgorithm()));
//velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
@@ -142,11 +155,11 @@ void bflow(string configname)
SPtr<BCProcessor> bcProc;
bcProc = SPtr<BCProcessor>(new BCProcessor());
- SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CumulantLBMKernel());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CumulantLBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulant4thOrderViscosityLBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
- //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
+ SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
kernel->setBCProcessor(bcProc);
//kernel->setForcingX1(forcing);
//kernel->setWithForcing(true);
@@ -175,6 +188,18 @@ void bflow(string configname)
//restartCoProcessor->setNu(k);
//////////////////////////////////////////////////////////////////////////
+ ////stator
+ SPtr<GbObject3D> stator(new GbCylinder3D(0.5 * g_maxX1, 0.5 * g_maxX2, g_minX3-2.0*deltax, 0.5 * g_maxX1, 0.5 * g_maxX2, g_maxX3+ 2.0 * deltax, 0.5 * g_maxX1));
+ GbSystem3D::writeGeoObject(stator.get(), outputPath + "/geo/stator", WbWriterVtkXmlBinary::getInstance());
+
+ SPtr<D3Q27Interactor> statorInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(stator, grid, velocityBCAdapter, Interactor3D::INVERSESOLID));
+
+ ////rotor (cylinder)
+ SPtr<GbObject3D> rotor(new GbCylinder3D(0.5 * g_maxX1, 0.5 * g_maxX2, g_minX3- 2.0 * deltax, 0.5 * g_maxX1, 0.5 * g_maxX2, g_maxX3+ 2.0 * deltax, 0.25 * g_maxX1));
+ GbSystem3D::writeGeoObject(rotor.get(), outputPath + "/geo/rotor", WbWriterVtkXmlBinary::getInstance());
+
+ SPtr<D3Q27Interactor> rotorInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(rotor, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
if (myid == 0)
{
UBLOG(logINFO, "Parameters:");
@@ -215,17 +240,6 @@ void bflow(string configname)
if (myid == 0) UBLOG(logINFO, "Refinement - end");
}
- ////stator
- SPtr<GbObject3D> stator(new GbCylinder3D(0.5*g_maxX1, 0.5 * g_maxX2, g_minX3-deltax, 0.5 * g_maxX1, 0.5 * g_maxX2, g_maxX3+deltax, 0.5 * g_maxX1));
- GbSystem3D::writeGeoObject(stator.get(), outputPath + "/geo/stator", WbWriterVtkXmlBinary::getInstance());
-
- SPtr<D3Q27Interactor> statorInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(stator, grid, velocityBCAdapter, Interactor3D::INVERSESOLID));
-
- ////rotor (cylinder)
- SPtr<GbObject3D> rotor(new GbCylinder3D(0.5 * g_maxX1, 0.5 * g_maxX2, g_minX3 - deltax, 0.5 * g_maxX1, 0.5 * g_maxX2, g_maxX3 + deltax, 0.25 * g_maxX1));
- GbSystem3D::writeGeoObject(rotor.get(), outputPath + "/geo/rotor", WbWriterVtkXmlBinary::getInstance());
-
- SPtr<D3Q27Interactor> rotorInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(rotor, grid, velocityBCAdapter, Interactor3D::SOLID));
//walls
GbCuboid3DPtr wallZmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_minX3));
@@ -326,9 +340,12 @@ void bflow(string configname)
{
restartCoProcessor->restart((int)restartStep);
grid->setTimeStep(restartStep);
- //SetBcBlocksBlockVisitor v(sphereInt);
- //grid->accept(v);
- //sphereInt->initInteractor();
+ SetBcBlocksBlockVisitor v1(rotorInt);
+ grid->accept(v1);
+ rotorInt->initInteractor();
+ SetBcBlocksBlockVisitor v2(statorInt);
+ grid->accept(v2);
+ statorInt->initInteractor();
}
omp_set_num_threads(numOfThreads);
@@ -354,19 +371,21 @@ void bflow(string configname)
SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
//writeMQCoProcessor->process(0);
- //double area = UbMath::PI*radius*radius;
- //SPtr<UbScheduler> forceSch(new UbScheduler(100));
- //SPtr<CalculateForcesCoProcessor> fp = make_shared<CalculateForcesCoProcessor>(grid, forceSch, outputPath + "/forces/forces.txt", comm, velocity, area);
- //fp->addInteractor(sphereInt);
+ SPtr<UbScheduler> forceSch(new UbScheduler(100));
+ SPtr<CalculateTorqueCoProcessor> fp = make_shared<CalculateTorqueCoProcessor>(grid, forceSch, outputPath + "/torque/TorqueRotor.txt", comm);
+ fp->addInteractor(rotorInt);
+ SPtr<CalculateTorqueCoProcessor> fp2 = make_shared<CalculateTorqueCoProcessor>(grid, forceSch, outputPath + "/torque/TorqueStator.txt", comm);
+ fp2->addInteractor(statorInt);
SPtr<WriteThixotropyQuantitiesCoProcessor> writeThixotropicMQCoProcessor(new WriteThixotropyQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
calculator->addCoProcessor(npr);
- //calculator->addCoProcessor(fp);
+ calculator->addCoProcessor(fp);
+ calculator->addCoProcessor(fp2);
calculator->addCoProcessor(writeMQCoProcessor);
- //calculator->addCoProcessor(writeThixotropicMQCoProcessor);
+ calculator->addCoProcessor(writeThixotropicMQCoProcessor);
calculator->addCoProcessor(restartCoProcessor);
if (myid == 0) UBLOG(logINFO, "Simulation-start");
diff --git a/src/cpu/VirtualFluids.h b/src/cpu/VirtualFluids.h
index 77a8a551e..e02a5cdf1 100644
--- a/src/cpu/VirtualFluids.h
+++ b/src/cpu/VirtualFluids.h
@@ -168,7 +168,9 @@
#include <CoProcessors/AdjustForcingCoProcessor.h>
#include <CoProcessors/CalculateForcesCoProcessor.h>
-#include <CoProcessors/WriteBlocksCoProcessor.h>
+#include <CoProcessors/CalculateTorqueCoProcessor.h>
+#include <CoProcessors/WriteMacroscopicQuantitiesCoProcessor.h>
+#include <CoProcessors/WriteMQFromSelectionCoProcessor.h>
#include <CoProcessors/WriteBoundaryConditionsCoProcessor.h>
#include <CoProcessors/WriteMQFromSelectionCoProcessor.h>
#include <CoProcessors/WriteMacroscopicQuantitiesCoProcessor.h>
@@ -221,6 +223,7 @@
#include <LBM/ThixotropyExpLBMKernel.h>
#include <LBM/CumulantLBMKernel.h>
#include <LBM/ThixotropyModelLBMKernel.h>
+#include <LBM/ThixotropyModelLBMKernel2.h>
#include <LBM/BinghamModelLBMKernel.h>
#include <LBM/HerschelBulkleyModelLBMKernel.h>
#include <LBM/ThixotropyInterpolationProcessor.h>
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
new file mode 100644
index 000000000..af7533ef5
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
@@ -0,0 +1,255 @@
+#include "CalculateTorqueCoProcessor.h"
+#include "BCProcessor.h"
+
+#include "Communicator.h"
+#include "D3Q27Interactor.h"
+#include "UbScheduler.h"
+#include "Grid3D.h"
+#include "BoundaryConditions.h"
+#include "DataSet3D.h"
+#include "Block3D.h"
+#include "LBMKernel.h"
+#include "BCArray3D.h"
+#include "EsoTwist3D.h"
+#include "DistributionArray3D.h"
+
+CalculateTorqueCoProcessor::CalculateTorqueCoProcessor( SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path, SPtr<Communicator> comm) : CoProcessor(grid, s), path(path), comm(comm), forceX1global(0), forceX2global(0), forceX3global(0)
+{
+ if (comm->getProcessID() == comm->getRoot())
+ {
+ std::ofstream ostr;
+ std::string fname = path;
+ ostr.open(fname.c_str(), std::ios_base::out | std::ios_base::app);
+ if(!ostr)
+ {
+ ostr.clear();
+ std::string path = UbSystem::getPathFromString(fname);
+ if(path.size()>0){ UbSystem::makeDirectory(path); ostr.open(fname.c_str(), std::ios_base::out | std::ios_base::app);}
+ if(!ostr) throw UbException(UB_EXARGS,"couldn't open file "+fname);
+ }
+ ostr.width(12);
+ ostr << "step" << "\t";
+ ostr.width(10);
+ ostr << "Tx" << "\t";
+ ostr.width(18);
+ ostr << "Ty" << "\t";
+ ostr.width(18);
+ ostr << "Tz" << std::endl;
+ ostr.close();
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+CalculateTorqueCoProcessor::~CalculateTorqueCoProcessor()
+{
+
+}
+//////////////////////////////////////////////////////////////////////////
+void CalculateTorqueCoProcessor::process( double step )
+{
+ if(scheduler->isDue(step) )
+ collectData(step);
+
+ UBLOG(logDEBUG3, "D3Q27ForcesCoProcessor::update:" << step);
+}
+//////////////////////////////////////////////////////////////////////////
+void CalculateTorqueCoProcessor::collectData( double step )
+{
+ calculateForces();
+
+ if (comm->getProcessID() == comm->getRoot())
+ {
+ int istep = static_cast<int>(step);
+ std::ofstream ostr;
+ std::string fname = path;
+ ostr.open(fname.c_str(), std::ios_base::out | std::ios_base::app);
+ if(!ostr)
+ {
+ ostr.clear();
+ std::string path = UbSystem::getPathFromString(fname);
+ if(path.size()>0){ UbSystem::makeDirectory(path); ostr.open(fname.c_str(), std::ios_base::out | std::ios_base::app);}
+ if(!ostr) throw UbException(UB_EXARGS,"couldn't open file "+fname);
+ }
+
+ ostr.width(12);
+ ostr.setf(std::ios::fixed);
+ ostr << istep << "\t";
+ write(&ostr, forceX1global, (char*)"\t");
+ write(&ostr, forceX2global, (char*)"\t");
+ write(&ostr, forceX3global, (char*)"\t");
+ ostr << std::endl;
+ ostr.close();
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+void CalculateTorqueCoProcessor::calculateForces()
+{
+ forceX1global = 0.0;
+ forceX2global = 0.0;
+ forceX3global = 0.0;
+ int counter = 0;
+
+ for(SPtr<D3Q27Interactor> interactor : interactors)
+ {
+ double x1Centre = interactor->getGbObject3D()->getX1Centroid();
+ double x2Centre = interactor->getGbObject3D()->getX2Centroid();
+ double x3Centre = interactor->getGbObject3D()->getX3Centroid();
+
+ for(BcNodeIndicesMap::value_type t : interactor->getBcNodeIndicesMap())
+ {
+ double torqueX1 = 0.0;
+ double torqueX2 = 0.0;
+ double torqueX3 = 0.0;
+
+ SPtr<Block3D> block = t.first;
+ std::set< std::vector<int> >& transNodeIndicesSet = t.second;
+
+ SPtr<ILBMKernel> kernel = block->getKernel();
+
+ if (kernel->getCompressible())
+ {
+ calcMacrosFct = &D3Q27System::calcCompMacroscopicValues;
+ compressibleFactor = 1.0;
+ }
+ else
+ {
+ calcMacrosFct = &D3Q27System::calcIncompMacroscopicValues;
+ compressibleFactor = 0.0;
+ }
+
+ SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
+ SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
+ distributions->swap();
+
+ int ghostLayerWidth = kernel->getGhostLayerWidth();
+ int minX1 = ghostLayerWidth;
+ int maxX1 = (int)bcArray->getNX1() - 1 - ghostLayerWidth;
+ int minX2 = ghostLayerWidth;
+ int maxX2 = (int)bcArray->getNX2() - 1 - ghostLayerWidth;
+ int minX3 = ghostLayerWidth;
+ int maxX3 = (int)bcArray->getNX3() - 1 - ghostLayerWidth;
+
+ for(std::vector<int> node : transNodeIndicesSet)
+ {
+ int x1 = node[0];
+ int x2 = node[1];
+ int x3 = node[2];
+
+ Vector3D worldCoordinates = grid->getNodeCoordinates(block, x1, x2, x3);
+
+ //without ghost nodes
+ if (x1 < minX1 || x1 > maxX1 || x2 < minX2 || x2 > maxX2 ||x3 < minX3 || x3 > maxX3 ) continue;
+
+ if(bcArray->isFluid(x1,x2,x3)) //es kann sein, dass der node von einem anderen interactor z.B. als solid gemarkt wurde!!!
+ {
+ SPtr<BoundaryConditions> bc = bcArray->getBC(x1,x2,x3);
+ UbTupleDouble3 forceVec = getForces(x1,x2,x3,distributions,bc);
+ torqueX1 += (worldCoordinates[1] - x2Centre) * val<3>(forceVec) - (worldCoordinates[2] - x3Centre) * val<2>(forceVec);
+ torqueX2 += (worldCoordinates[2] - x3Centre) * val<1>(forceVec) - (worldCoordinates[0] - x1Centre) * val<3>(forceVec);
+ torqueX3 += (worldCoordinates[0] - x1Centre) * val<2>(forceVec) - (worldCoordinates[1] - x2Centre) * val<1>(forceVec);
+ //counter++;
+ //UBLOG(logINFO, "x1="<<(worldCoordinates[1] - x2Centre)<<",x2=" << (worldCoordinates[2] - x3Centre)<< ",x3=" << (worldCoordinates[0] - x1Centre) <<" forceX3 = " << forceX3);
+ }
+ }
+ //if we have got discretization with more level
+ // deltaX is LBM deltaX and equal LBM deltaT
+ double deltaX = LBMSystem::getDeltaT(block->getLevel()); //grid->getDeltaT(block);
+ double deltaXquadrat = deltaX*deltaX;
+ torqueX1 *= deltaXquadrat;
+ torqueX2 *= deltaXquadrat;
+ torqueX3 *= deltaXquadrat;
+
+ distributions->swap();
+
+ forceX1global += torqueX1;
+ forceX2global += torqueX2;
+ forceX3global += torqueX3;
+
+ //UBLOG(logINFO, "forceX3global = " << forceX3global);
+ }
+ }
+ std::vector<double> values;
+ std::vector<double> rvalues;
+ values.push_back(forceX1global);
+ values.push_back(forceX2global);
+ values.push_back(forceX3global);
+
+ rvalues = comm->gather(values);
+ if (comm->getProcessID() == comm->getRoot())
+ {
+ forceX1global = 0.0;
+ forceX2global = 0.0;
+ forceX3global = 0.0;
+
+ for (int i = 0; i < (int)rvalues.size(); i+=3)
+ {
+ forceX1global += rvalues[i];
+ forceX2global += rvalues[i+1];
+ forceX3global += rvalues[i+2];
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+UbTupleDouble3 CalculateTorqueCoProcessor::getForces(int x1, int x2, int x3, SPtr<DistributionArray3D> distributions, SPtr<BoundaryConditions> bc)
+{
+ UbTupleDouble3 force(0.0,0.0,0.0);
+
+ LBMReal fs[D3Q27System::ENDF + 1];
+ distributions->getDistributionInv(fs, x1, x2, x3);
+ LBMReal rho = 0.0, vx1 = 0.0, vx2 = 0.0, vx3 = 0.0, drho = 0.0;
+ calcMacrosFct(fs, drho, vx1, vx2, vx3);
+ rho = 1.0 + drho * compressibleFactor;
+
+ if(bc)
+ {
+ //references to tuple "force"
+ double& forceX1 = val<1>(force);
+ double& forceX2 = val<2>(force);
+ double& forceX3 = val<3>(force);
+ double f, fnbr;
+
+ for(int fdir=D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++)
+ {
+ if(bc->hasNoSlipBoundaryFlag(fdir) || bc->hasVelocityBoundaryFlag(fdir))
+ {
+ const int invDir = D3Q27System::INVDIR[fdir];
+ f = dynamicPointerCast<EsoTwist3D>(distributions)->getDistributionInvForDirection(x1, x2, x3, invDir);
+ fnbr = dynamicPointerCast<EsoTwist3D>(distributions)->getDistributionInvForDirection(x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
+
+ Vector3D boundaryVelocity;
+ boundaryVelocity[0] = bc->getBoundaryVelocityX1();
+ boundaryVelocity[1] = bc->getBoundaryVelocityX2();
+ boundaryVelocity[2] = bc->getBoundaryVelocityX3();
+ double correction[3] = { 0.0, 0.0, 0.0 };
+ if (bc->hasVelocityBoundaryFlag(fdir))
+ {
+ const double forceTerm = f - fnbr;
+ correction[0] = forceTerm * boundaryVelocity[0];
+ correction[1] = forceTerm * boundaryVelocity[1];
+ correction[2] = forceTerm * boundaryVelocity[2];
+ }
+
+ forceX1 += (f + fnbr) * D3Q27System::DX1[invDir] - 2.0 * D3Q27System::WEIGTH[invDir] * rho - correction[0];
+ forceX2 += (f + fnbr) * D3Q27System::DX2[invDir] - 2.0 * D3Q27System::WEIGTH[invDir] * rho - correction[1];
+ forceX3 += (f + fnbr) * D3Q27System::DX3[invDir] - 2.0 * D3Q27System::WEIGTH[invDir] * rho - correction[2];
+ }
+ }
+ }
+
+ return force;
+}
+//////////////////////////////////////////////////////////////////////////
+void CalculateTorqueCoProcessor::addInteractor( SPtr<D3Q27Interactor> interactor )
+{
+ interactors.push_back(interactor);
+}
+//////////////////////////////////////////////////////////////////////////
+void CalculateTorqueCoProcessor::write(std::ofstream *fileObject, double value, char *separator)
+{
+ (*fileObject).width(12);
+ (*fileObject).precision(16);
+ (*fileObject).setf(std::ios::fixed);
+ (*fileObject) << value;
+ (*fileObject) << separator;
+}
+
+
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h
new file mode 100644
index 000000000..2eda4a063
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h
@@ -0,0 +1,54 @@
+/*
+ * D3Q27ForcesCoProcessor.h
+ *
+ * Created on: 29.09.2012
+ * Author: K. Kucher
+ */
+
+#ifndef CalculateTorqueCoProcessor_H
+#define CalculateTorqueCoProcessor_H
+
+#include <PointerDefinitions.h>
+#include <string>
+#include <vector>
+
+#include "CoProcessor.h"
+#include "UbTuple.h"
+#include "D3Q27System.h"
+
+class ForceCalculator;
+class Communicator;
+class Grid3D;
+class UbScheduler;
+class D3Q27Interactor;
+class DistributionArray3D;
+class BoundaryConditions;
+
+class CalculateTorqueCoProcessor: public CoProcessor
+{
+public:
+ //! Constructor
+ CalculateTorqueCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string &path, SPtr<Communicator> comm);
+ virtual ~CalculateTorqueCoProcessor();
+ void process(double step);
+ void addInteractor(SPtr<D3Q27Interactor> interactor);
+protected:
+ void collectData(double step);
+ void calculateForces();
+ UbTupleDouble3 getForces(int x1, int x2, int x3, SPtr<DistributionArray3D> distributions, SPtr<BoundaryConditions> bc);
+ void write(std::ofstream *fileObject, double value, char *separator);
+private:
+ std::string path;
+ SPtr<Communicator> comm;
+ std::vector<SPtr<D3Q27Interactor> > interactors;
+ double forceX1global;
+ double forceX2global;
+ double forceX3global;
+
+ typedef void(*CalcMacrosFct)(const LBMReal* const& /*f[27]*/, LBMReal& /*rho*/, LBMReal& /*vx1*/, LBMReal& /*vx2*/, LBMReal& /*vx3*/);
+ CalcMacrosFct calcMacrosFct;
+ LBMReal compressibleFactor;
+};
+
+
+#endif /* D3Q27ForcesCoProcessor_H */
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
index fe5f31788..72e7114e6 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
@@ -123,7 +123,7 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
datanames.push_back("viscosity");
//datanames.push_back("lambda");
//datanames.push_back("ShearRate");
- //datanames.push_back("collFactor");
+ datanames.push_back("omega");
//datanames.push_back("Fluxx");
//datanames.push_back("Fluxy");
//datanames.push_back("Fluxz");
@@ -200,11 +200,13 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
//distributionsF->getDistribution(f, ix1, ix2, ix3);
//LBMReal rho = D3Q27System::getDensity(f);
- //gammaDot = BinghamModel::getShearRate(f, collFactor);
+ //LBMReal gammaDot = D3Q27System::getShearRate(f, collFactor);
//LBMReal collFactorF = collFactor - 1e-6 / (gammaDot + one * 1e-9);
//collFactorF = (collFactorF < 0.5) ? 0.5 : collFactorF;
+ //LBMReal collFactorF = Thixotropy::getBinghamCollFactor(collFactor, gammaDot, rho);
+
//data[index++].push_back(lambda);
//data[index++].push_back(gammaDot);
//data[index++].push_back(collFactorF);
@@ -213,10 +215,13 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
LBMReal rho = D3Q27System::getDensity(f);
LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
//LBMReal omega = Thixotropy::getHerschelBulkleyCollFactor(collFactor, shearRate, rho);
- LBMReal omega = Thixotropy::getPowellEyringCollFactor(collFactor, shearRate, rho);
+ //LBMReal omega = Thixotropy::getPowellEyringCollFactor(collFactor, shearRate, rho);
+ LBMReal omega = Thixotropy::getBinghamCollFactor(collFactor, shearRate, rho);
LBMReal viscosity = (omega == 0) ? 0 : c1o3 * (c1/omega-c1o2);
+
data[index++].push_back(viscosity);
+ data[index++].push_back(omega);
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
index d94a7f37d..d240ed8ef 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
@@ -1171,6 +1171,8 @@ usage: ...
LBMReal Dxz = -three * collFactorF * mfbab;
LBMReal Dyz = -three * collFactorF * mfabb;
+
+ //TODO: may be factor 2
return sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
index b731337e6..50caabbe8 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
@@ -33,6 +33,7 @@ public:
void setMu0(LBMReal mu);
LBMReal getMu0() const;
+ static LBMReal getBinghamCollFactorOld(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
static LBMReal getBinghamCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
static LBMReal getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
static LBMReal getHerschelBulkleyCollFactorBackward(LBMReal shearRate, LBMReal drho);
@@ -56,9 +57,42 @@ inline LBMReal Thixotropy::getBinghamCollFactor(LBMReal omegaInf, LBMReal shearR
{
LBMReal cs2 = UbMath::one_over_sqrt3 * UbMath::one_over_sqrt3;
LBMReal rho = UbMath::one + drho;
- LBMReal omega = omegaInf * (UbMath::one - (omegaInf * tau0) / (shearRate * cs2 * rho + UbMath::Epsilon<LBMReal>::val()));
+ //LBMReal omega = omegaInf * (UbMath::one - (omegaInf * tau0) / (shearRate * cs2 * rho + UbMath::Epsilon<LBMReal>::val()));
+
+ //LBMReal omega = cs2 * cs2 * shearRate * shearRate * omegaInf * rho * rho / (cs2 * cs2 * shearRate * shearRate * rho * rho + cs2 * shearRate * omegaInf * rho * tau0+omegaInf*omegaInf*tau0*tau0);
+
+ LBMReal a = omegaInf * tau0 / (cs2 * shearRate * rho);
+ //10 iterations
+ //LBMReal omega = omegaInf / (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a))))))))));
+
+ //20 iterations
+ LBMReal omega = omegaInf / (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a))))))))))))))))))));
+
+ //LBMReal omega = omegaInf*cs2 * shearRate * rho / (cs2 * shearRate * rho + omegaInf * tau0);
+ //LBMReal shearRateNew = shearRate * (omega / omegaInf);
+
+ //for (int i = 0; i < 20; i++)
+ //{
+ // omega = omegaInf * cs2 * shearRateNew * rho / (cs2 * shearRateNew * rho + omegaInf * tau0);
+ // shearRateNew = shearRate * (omega / omegaInf);
+ //}
+ //omega = omegaInf * cs2 * shearRateNew * rho / (cs2 * shearRateNew * rho + omegaInf * tau0);
+
+ //if (omega < 0.2)
+ // omega = 0.2;
return omega;
}
+
+inline LBMReal Thixotropy::getBinghamCollFactorOld(LBMReal omegaInf, LBMReal shearRate, LBMReal drho)
+{
+ const LBMReal cs2 = UbMath::c1o3; // UbMath::one_over_sqrt3* UbMath::one_over_sqrt3;
+ LBMReal rho = UbMath::one + drho;
+
+ if (rho * cs2 * (UbMath::c1 / omegaInf - UbMath::c1o2) * shearRate < tau0)
+ return 0.0;
+ else
+ return omegaInf;
+}
//////////////////////////////////////////////////////////////////////////
inline LBMReal Thixotropy::getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho)
{
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
new file mode 100644
index 000000000..cafded555
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
@@ -0,0 +1,899 @@
+#include "ThixotropyModelLBMKernel2.h"
+#include "D3Q27System.h"
+#include "BCArray3D.h"
+#include "D3Q27EsoTwist3DSplittedVector.h"
+#include <math.h>
+#include "DataSet3D.h"
+#include "LBMKernel.h"
+
+#define PROOF_CORRECTNESS
+
+ThixotropyModelLBMKernel2::ThixotropyModelLBMKernel2() : forcingX1(0), forcingX2(0), forcingX3(0)
+{
+ compressible = false;
+ OxyyMxzz = 1.0;
+}
+
+ThixotropyModelLBMKernel2::~ThixotropyModelLBMKernel2()
+{
+}
+
+void ThixotropyModelLBMKernel2::calculate(int step)
+{
+ using namespace D3Q27System;
+
+ //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);
+ }
+ /////////////////////////////////////
+
+ localDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
+ nonLocalDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
+ zeroDistributionsF = dynamicPointerCast<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;
+
+
+ {
+ for (int x3 = minX3; x3 < maxX3; x3++)
+ {
+ for (int x2 = minX2; x2 < maxX2; x2++)
+ {
+ for (int x1 = minX1; x1 < maxX1; x1++)
+ {
+ if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3))
+ {
+ int x1p = x1 + 1;
+ int x2p = x2 + 1;
+ int x3p = x3 + 1;
+ //////////////////////////////////////////////////////////////////////////
+ //read distribution
+ // Cumulant (NSE part)
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+
+ //E N T
+ //c c c
+ //////////
+ //W S B
+ //a a a
+
+ //Rest ist b
+
+ //mfxyz
+ //a - negative
+ //b - null
+ //c - positive
+
+ // a b c
+ //-1 0 1
+
+ LBMReal mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3);
+ LBMReal mfbcb = (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3);
+ LBMReal mfbbc = (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3);
+ LBMReal mfccb = (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3);
+ LBMReal mfacb = (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3);
+ LBMReal mfcbc = (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3);
+ LBMReal mfabc = (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3);
+ LBMReal mfbcc = (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3);
+ LBMReal mfbac = (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3);
+ LBMReal mfccc = (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3);
+ LBMReal mfacc = (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3);
+ LBMReal mfcac = (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3);
+ LBMReal mfaac = (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3);
+
+ LBMReal mfabb = (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3);
+ LBMReal mfbab = (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3);
+ LBMReal mfbba = (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p);
+ LBMReal mfaab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3);
+ LBMReal mfcab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3);
+ LBMReal mfaba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p);
+ LBMReal mfcba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p);
+ LBMReal mfbaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p);
+ LBMReal mfbca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p);
+ LBMReal mfaaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p);
+ LBMReal mfcaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p);
+ LBMReal mfaca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p);
+ LBMReal mfcca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p);
+
+ LBMReal mfbbb = (*this->zeroDistributionsF)(x1, x2, x3);
+
+ LBMReal m0, m1, m2;
+
+ LBMReal rho = (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 vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ (mfcbb - mfabb));
+ LBMReal vvy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ (mfbcb - mfbab));
+ LBMReal vvz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ (mfbbc - mfbba));
+
+ LBMReal collFactorF = collFactor;
+
+ //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();
+
+ vvx += forcingX1 * deltaT * 0.5; // X
+ vvy += forcingX2 * deltaT * 0.5; // Y
+ vvz += forcingX3 * deltaT * 0.5; // Z
+ }
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ LBMReal oMdrho;
+
+ oMdrho = mfccc + mfaaa;
+ m0 = mfaca + mfcac;
+ m1 = mfacc + mfcaa;
+ m2 = mfaac + mfcca;
+ oMdrho += m0;
+ m1 += m2;
+ oMdrho += m1;
+ m0 = mfbac + mfbca;
+ m1 = mfbaa + mfbcc;
+ m0 += m1;
+ m1 = mfabc + mfcba;
+ m2 = mfaba + mfcbc;
+ m1 += m2;
+ m0 += m1;
+ m1 = mfacb + mfcab;
+ m2 = mfaab + mfccb;
+ m1 += m2;
+ m0 += m1;
+ oMdrho += m0;
+ m0 = mfabb + mfcbb;
+ m1 = mfbab + mfbcb;
+ m2 = mfbba + mfbbc;
+ m0 += m1 + m2;
+ m0 += mfbbb; //hat gefehlt
+ oMdrho = 1. - (oMdrho + m0);
+
+ LBMReal vx2;
+ LBMReal vy2;
+ LBMReal vz2;
+ vx2 = vvx * vvx;
+ vy2 = vvy * vvy;
+ vz2 = vvz * vvz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ LBMReal wadjust;
+ LBMReal qudricLimit = 0.01;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //Hin
+ ////////////////////////////////////////////////////////////////////////////////////
+ // mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ m2 = mfaaa + mfaac;
+ m1 = mfaac - mfaaa;
+ m0 = m2 + mfaab;
+ mfaaa = m0;
+ m0 += c1o36 * oMdrho;
+ mfaab = m1 - m0 * vvz;
+ mfaac = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaba + mfabc;
+ m1 = mfabc - mfaba;
+ m0 = m2 + mfabb;
+ mfaba = m0;
+ m0 += c1o9 * oMdrho;
+ mfabb = m1 - m0 * vvz;
+ mfabc = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaca + mfacc;
+ m1 = mfacc - mfaca;
+ m0 = m2 + mfacb;
+ mfaca = m0;
+ m0 += c1o36 * oMdrho;
+ mfacb = m1 - m0 * vvz;
+ mfacc = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbaa + mfbac;
+ m1 = mfbac - mfbaa;
+ m0 = m2 + mfbab;
+ mfbaa = m0;
+ m0 += c1o9 * oMdrho;
+ mfbab = m1 - m0 * vvz;
+ mfbac = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbba + mfbbc;
+ m1 = mfbbc - mfbba;
+ m0 = m2 + mfbbb;
+ mfbba = m0;
+ m0 += c4o9 * oMdrho;
+ mfbbb = m1 - m0 * vvz;
+ mfbbc = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbca + mfbcc;
+ m1 = mfbcc - mfbca;
+ m0 = m2 + mfbcb;
+ mfbca = m0;
+ m0 += c1o9 * oMdrho;
+ mfbcb = m1 - m0 * vvz;
+ mfbcc = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcaa + mfcac;
+ m1 = mfcac - mfcaa;
+ m0 = m2 + mfcab;
+ mfcaa = m0;
+ m0 += c1o36 * oMdrho;
+ mfcab = m1 - m0 * vvz;
+ mfcac = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcba + mfcbc;
+ m1 = mfcbc - mfcba;
+ m0 = m2 + mfcbb;
+ mfcba = m0;
+ m0 += c1o9 * oMdrho;
+ mfcbb = m1 - m0 * vvz;
+ mfcbc = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcca + mfccc;
+ m1 = mfccc - mfcca;
+ m0 = m2 + mfccb;
+ mfcca = m0;
+ m0 += c1o36 * oMdrho;
+ mfccb = m1 - m0 * vvz;
+ mfccc = m2 - 2. * m1 * vvz + vz2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ // mit 1/6, 0, 1/18, 2/3, 0, 2/9, 1/6, 0, 1/18 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ m2 = mfaaa + mfaca;
+ m1 = mfaca - mfaaa;
+ m0 = m2 + mfaba;
+ mfaaa = m0;
+ m0 += c1o6 * oMdrho;
+ mfaba = m1 - m0 * vvy;
+ mfaca = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaab + mfacb;
+ m1 = mfacb - mfaab;
+ m0 = m2 + mfabb;
+ mfaab = m0;
+ mfabb = m1 - m0 * vvy;
+ mfacb = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaac + mfacc;
+ m1 = mfacc - mfaac;
+ m0 = m2 + mfabc;
+ mfaac = m0;
+ m0 += c1o18 * oMdrho;
+ mfabc = m1 - m0 * vvy;
+ mfacc = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbaa + mfbca;
+ m1 = mfbca - mfbaa;
+ m0 = m2 + mfbba;
+ mfbaa = m0;
+ m0 += c2o3 * oMdrho;
+ mfbba = m1 - m0 * vvy;
+ mfbca = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbab + mfbcb;
+ m1 = mfbcb - mfbab;
+ m0 = m2 + mfbbb;
+ mfbab = m0;
+ mfbbb = m1 - m0 * vvy;
+ mfbcb = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfbac + mfbcc;
+ m1 = mfbcc - mfbac;
+ m0 = m2 + mfbbc;
+ mfbac = m0;
+ m0 += c2o9 * oMdrho;
+ mfbbc = m1 - m0 * vvy;
+ mfbcc = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcaa + mfcca;
+ m1 = mfcca - mfcaa;
+ m0 = m2 + mfcba;
+ mfcaa = m0;
+ m0 += c1o6 * oMdrho;
+ mfcba = m1 - m0 * vvy;
+ mfcca = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcab + mfccb;
+ m1 = mfccb - mfcab;
+ m0 = m2 + mfcbb;
+ mfcab = m0;
+ mfcbb = m1 - m0 * vvy;
+ mfccb = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfcac + mfccc;
+ m1 = mfccc - mfcac;
+ m0 = m2 + mfcbc;
+ mfcac = m0;
+ m0 += c1o18 * oMdrho;
+ mfcbc = m1 - m0 * vvy;
+ mfccc = m2 - 2. * m1 * vvy + vy2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ // mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ m2 = mfaaa + mfcaa;
+ m1 = mfcaa - mfaaa;
+ m0 = m2 + mfbaa;
+ mfaaa = m0;
+ m0 += 1. * oMdrho;
+ mfbaa = m1 - m0 * vvx;
+ mfcaa = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaba + mfcba;
+ m1 = mfcba - mfaba;
+ m0 = m2 + mfbba;
+ mfaba = m0;
+ mfbba = m1 - m0 * vvx;
+ mfcba = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaca + mfcca;
+ m1 = mfcca - mfaca;
+ m0 = m2 + mfbca;
+ mfaca = m0;
+ m0 += c1o3 * oMdrho;
+ mfbca = m1 - m0 * vvx;
+ mfcca = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaab + mfcab;
+ m1 = mfcab - mfaab;
+ m0 = m2 + mfbab;
+ mfaab = m0;
+ mfbab = m1 - m0 * vvx;
+ mfcab = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfabb + mfcbb;
+ m1 = mfcbb - mfabb;
+ m0 = m2 + mfbbb;
+ mfabb = m0;
+ mfbbb = m1 - m0 * vvx;
+ mfcbb = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfacb + mfccb;
+ m1 = mfccb - mfacb;
+ m0 = m2 + mfbcb;
+ mfacb = m0;
+ mfbcb = m1 - m0 * vvx;
+ mfccb = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfaac + mfcac;
+ m1 = mfcac - mfaac;
+ m0 = m2 + mfbac;
+ mfaac = m0;
+ m0 += c1o3 * oMdrho;
+ mfbac = m1 - m0 * vvx;
+ mfcac = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfabc + mfcbc;
+ m1 = mfcbc - mfabc;
+ m0 = m2 + mfbbc;
+ mfabc = m0;
+ mfbbc = m1 - m0 * vvx;
+ mfcbc = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m2 = mfacc + mfccc;
+ m1 = mfccc - mfacc;
+ m0 = m2 + mfbcc;
+ mfacc = m0;
+ m0 += c1o9 * oMdrho;
+ mfbcc = m1 - m0 * vvx;
+ mfccc = m2 - 2. * m1 * vvx + vx2 * m0;
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Cumulants
+ ////////////////////////////////////////////////////////////////////////////////////
+ LBMReal OxxPyyPzz = 1.; //omega2 or bulk viscosity
+ LBMReal OxyyPxzz = 1.;//-s9;//2+s9;//
+ //LBMReal OxyyMxzz = 1.;//2+s9;//
+ LBMReal O4 = 1.;
+ LBMReal O5 = 1.;
+ LBMReal O6 = 1.;
+
+ //Cum 4.
+ //LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
+ //LBMReal CUMbcb = mfbcb - ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
+ //LBMReal CUMbbc = mfbbc - ((mfaac + c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015
+
+ LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab);
+ LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb);
+ LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb);
+
+ LBMReal CUMcca = mfcca - ((mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
+ LBMReal CUMcac = mfcac - ((mfcaa * mfaac + 2. * mfbab * mfbab) + c1o3 * (mfcaa + mfaac) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
+ LBMReal CUMacc = mfacc - ((mfaac * mfaca + 2. * mfabb * mfabb) + c1o3 * (mfaac + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
+
+ //Cum 5.
+ LBMReal CUMbcc = mfbcc - (mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb + 2. * (mfbab * mfacb + mfbba * mfabc)) - c1o3 * (mfbca + mfbac) * oMdrho;
+ LBMReal CUMcbc = mfcbc - (mfaac * mfcba + mfcaa * mfabc + 4. * mfbab * mfbbb + 2. * (mfabb * mfcab + mfbba * mfbac)) - c1o3 * (mfcba + mfabc) * oMdrho;
+ LBMReal CUMccb = mfccb - (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) - c1o3 * (mfacb + mfcab) * oMdrho;
+
+ //Cum 6.
+ LBMReal CUMccc = mfccc + ((-4. * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - 4. * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+ - 2. * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb))
+ + (4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + 2. * (mfcaa * mfaca * mfaac)
+ + 16. * mfbba * mfbab * mfabb)
+ - c1o3 * (mfacc + mfcac + mfcca) * oMdrho - c1o9 * oMdrho * oMdrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * oMdrho * (1. - 2. * oMdrho) - c1o27 * oMdrho * oMdrho * (-2. * oMdrho)
+ + (2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3 * oMdrho) + c1o27 * oMdrho;
+
+ //2.
+ // linear combinations
+ LBMReal mxxPyyPzz = mfcaa + mfaca + mfaac;
+ LBMReal mxxMyy = mfcaa - mfaca;
+ LBMReal mxxMzz = mfcaa - mfaac;
+
+ LBMReal dxux = -c1o2 * collFactorF * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
+ LBMReal dyuy = dxux + collFactorF * c3o2 * mxxMyy;
+ LBMReal dzuz = dxux + collFactorF * c3o2 * mxxMzz;
+
+ LBMReal Dxy = -three * collFactorF * mfbba;
+ LBMReal Dxz = -three * collFactorF * mfbab;
+ LBMReal Dyz = -three * collFactorF * mfabb;
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ //non Newtonian fluid collision factor
+ LBMReal shearRate = sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+ //collFactorF = getThyxotropyCollFactor(collFactorF, shearRate, rho);
+ ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+
+ //relax
+ mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) - 3. * (1. - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+
+ LBMReal collFactorFyy = getThyxotropyCollFactor(collFactorF, std::sqrt(dxux*dxux + dyuy*dyuy) / (rho + one), rho);
+ mxxMyy += collFactorFyy * (-mxxMyy) - 3. * (1. - c1o2 * collFactorFyy) * (vx2 * dxux - vy2 * dyuy);
+
+ LBMReal collFactorFzz = getThyxotropyCollFactor(collFactorF, std::sqrt(dxux*dxux + dzuz*dzuz) / (rho + one), rho);
+ mxxMzz += collFactorFzz * (-mxxMzz) - 3. * (1. - c1o2 * collFactorFzz) * (vx2 * dxux - vz2 * dzuz);
+
+ mfabb += getThyxotropyCollFactor(collFactorF, std::abs(Dyz) / (rho + one), rho) * (-mfabb);
+ mfbab += getThyxotropyCollFactor(collFactorF, std::abs(Dxz) / (rho + one), rho) * (-mfbab);
+ mfbba += getThyxotropyCollFactor(collFactorF, std::abs(Dxy) / (rho + one), rho) * (-mfbba);
+
+ // linear combinations back
+ mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaca = c1o3 * (-2. * mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaac = c1o3 * (mxxMyy - 2. * mxxMzz + mxxPyyPzz);
+
+ //3.
+ // linear combinations
+ LBMReal mxxyPyzz = mfcba + mfabc;
+ LBMReal mxxyMyzz = mfcba - mfabc;
+
+ LBMReal mxxzPyyz = mfcab + mfacb;
+ LBMReal mxxzMyyz = mfcab - mfacb;
+
+ LBMReal mxyyPxzz = mfbca + mfbac;
+ LBMReal mxyyMxzz = mfbca - mfbac;
+
+ //relax
+ wadjust = OxyyMxzz + (1. - OxyyMxzz) * fabs(mfbbb) / (fabs(mfbbb) + qudricLimit);
+ mfbbb += wadjust * (-mfbbb);
+ wadjust = OxyyPxzz + (1. - OxyyPxzz) * fabs(mxxyPyzz) / (fabs(mxxyPyzz) + qudricLimit);
+ mxxyPyzz += wadjust * (-mxxyPyzz);
+ wadjust = OxyyMxzz + (1. - OxyyMxzz) * fabs(mxxyMyzz) / (fabs(mxxyMyzz) + qudricLimit);
+ mxxyMyzz += wadjust * (-mxxyMyzz);
+ wadjust = OxyyPxzz + (1. - OxyyPxzz) * fabs(mxxzPyyz) / (fabs(mxxzPyyz) + qudricLimit);
+ mxxzPyyz += wadjust * (-mxxzPyyz);
+ wadjust = OxyyMxzz + (1. - OxyyMxzz) * fabs(mxxzMyyz) / (fabs(mxxzMyyz) + qudricLimit);
+ mxxzMyyz += wadjust * (-mxxzMyyz);
+ wadjust = OxyyPxzz + (1. - OxyyPxzz) * fabs(mxyyPxzz) / (fabs(mxyyPxzz) + qudricLimit);
+ mxyyPxzz += wadjust * (-mxyyPxzz);
+ wadjust = OxyyMxzz + (1. - OxyyMxzz) * fabs(mxyyMxzz) / (fabs(mxyyMxzz) + qudricLimit);
+ mxyyMxzz += wadjust * (-mxyyMxzz);
+
+ // linear combinations back
+ 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.
+ CUMacc += O4 * (-CUMacc);
+ CUMcac += O4 * (-CUMcac);
+ CUMcca += O4 * (-CUMcca);
+
+ CUMbbc += O4 * (-CUMbbc);
+ CUMbcb += O4 * (-CUMbcb);
+ CUMcbb += O4 * (-CUMcbb);
+
+ //5.
+ CUMbcc += O5 * (-CUMbcc);
+ CUMcbc += O5 * (-CUMcbc);
+ CUMccb += O5 * (-CUMccb);
+
+ //6.
+ CUMccc += O6 * (-CUMccc);
+
+ //back cumulants to central moments
+ //4.
+ //mfcbb = CUMcbb + ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
+ //mfbcb = CUMbcb + ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
+ //mfbbc = CUMbbc + ((mfaac + c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015
+
+ mfcbb = CUMcbb + ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab);
+ mfbcb = CUMbcb + ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb);
+ mfbbc = CUMbbc + ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb);
+
+ mfcca = CUMcca + (mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho;
+ mfcac = CUMcac + (mfcaa * mfaac + 2. * mfbab * mfbab) + c1o3 * (mfcaa + mfaac) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho;
+ mfacc = CUMacc + (mfaac * mfaca + 2. * mfabb * mfabb) + c1o3 * (mfaac + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho;
+
+ //5.
+ mfbcc = CUMbcc + (mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb + 2. * (mfbab * mfacb + mfbba * mfabc)) + c1o3 * (mfbca + mfbac) * oMdrho;
+ mfcbc = CUMcbc + (mfaac * mfcba + mfcaa * mfabc + 4. * mfbab * mfbbb + 2. * (mfabb * mfcab + mfbba * mfbac)) + c1o3 * (mfcba + mfabc) * oMdrho;
+ mfccb = CUMccb + (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) + c1o3 * (mfacb + mfcab) * oMdrho;
+
+ //6.
+ mfccc = CUMccc - ((-4. * mfbbb * mfbbb
+ - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+ - 4. * (mfabb * mfcbb + mfbac * mfbca + mfbba * mfbbc)
+ - 2. * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb))
+ + (4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+ + 2. * (mfcaa * mfaca * mfaac)
+ + 16. * mfbba * mfbab * mfabb)
+ - c1o3 * (mfacc + mfcac + mfcca) * oMdrho - c1o9 * oMdrho * oMdrho
+ - c1o9 * (mfcaa + mfaca + mfaac) * oMdrho * (1. - 2. * oMdrho) - c1o27 * oMdrho * oMdrho * (-2. * oMdrho)
+ + (2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+ + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3 * oMdrho) - c1o27 * oMdrho;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //forcing
+ mfbaa = -mfbaa;
+ mfaba = -mfaba;
+ mfaab = -mfaab;
+ //////////////////////////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //back
+ ////////////////////////////////////////////////////////////////////////////////////
+ //mit 1, 0, 1/3, 0, 0, 0, 1/3, 0, 1/9 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ m0 = mfaac * c1o2 + mfaab * (vvz - c1o2) + (mfaaa + 1. * oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfaac - 2. * mfaab * vvz + mfaaa * (1. - vz2) - 1. * oMdrho * vz2;
+ m2 = mfaac * c1o2 + mfaab * (vvz + c1o2) + (mfaaa + 1. * oMdrho) * (vz2 + vvz) * c1o2;
+ mfaaa = m0;
+ mfaab = m1;
+ mfaac = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfabc * c1o2 + mfabb * (vvz - c1o2) + mfaba * (vz2 - vvz) * c1o2;
+ m1 = -mfabc - 2. * mfabb * vvz + mfaba * (1. - vz2);
+ m2 = mfabc * c1o2 + mfabb * (vvz + c1o2) + mfaba * (vz2 + vvz) * c1o2;
+ mfaba = m0;
+ mfabb = m1;
+ mfabc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfacc * c1o2 + mfacb * (vvz - c1o2) + (mfaca + c1o3 * oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfacc - 2. * mfacb * vvz + mfaca * (1. - vz2) - c1o3 * oMdrho * vz2;
+ m2 = mfacc * c1o2 + mfacb * (vvz + c1o2) + (mfaca + c1o3 * oMdrho) * (vz2 + vvz) * c1o2;
+ mfaca = m0;
+ mfacb = m1;
+ mfacc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfbac * c1o2 + mfbab * (vvz - c1o2) + mfbaa * (vz2 - vvz) * c1o2;
+ m1 = -mfbac - 2. * mfbab * vvz + mfbaa * (1. - vz2);
+ m2 = mfbac * c1o2 + mfbab * (vvz + c1o2) + mfbaa * (vz2 + vvz) * c1o2;
+ mfbaa = m0;
+ mfbab = m1;
+ mfbac = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbbc * c1o2 + mfbbb * (vvz - c1o2) + mfbba * (vz2 - vvz) * c1o2;
+ m1 = -mfbbc - 2. * mfbbb * vvz + mfbba * (1. - vz2);
+ m2 = mfbbc * c1o2 + mfbbb * (vvz + c1o2) + mfbba * (vz2 + vvz) * c1o2;
+ mfbba = m0;
+ mfbbb = m1;
+ mfbbc = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbcc * c1o2 + mfbcb * (vvz - c1o2) + mfbca * (vz2 - vvz) * c1o2;
+ m1 = -mfbcc - 2. * mfbcb * vvz + mfbca * (1. - vz2);
+ m2 = mfbcc * c1o2 + mfbcb * (vvz + c1o2) + mfbca * (vz2 + vvz) * c1o2;
+ mfbca = m0;
+ mfbcb = m1;
+ mfbcc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcac * c1o2 + mfcab * (vvz - c1o2) + (mfcaa + c1o3 * oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfcac - 2. * mfcab * vvz + mfcaa * (1. - vz2) - c1o3 * oMdrho * vz2;
+ m2 = mfcac * c1o2 + mfcab * (vvz + c1o2) + (mfcaa + c1o3 * oMdrho) * (vz2 + vvz) * c1o2;
+ mfcaa = m0;
+ mfcab = m1;
+ mfcac = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfcbc * c1o2 + mfcbb * (vvz - c1o2) + mfcba * (vz2 - vvz) * c1o2;
+ m1 = -mfcbc - 2. * mfcbb * vvz + mfcba * (1. - vz2);
+ m2 = mfcbc * c1o2 + mfcbb * (vvz + c1o2) + mfcba * (vz2 + vvz) * c1o2;
+ mfcba = m0;
+ mfcbb = m1;
+ mfcbc = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfccc * c1o2 + mfccb * (vvz - c1o2) + (mfcca + c1o9 * oMdrho) * (vz2 - vvz) * c1o2;
+ m1 = -mfccc - 2. * mfccb * vvz + mfcca * (1. - vz2) - c1o9 * oMdrho * vz2;
+ m2 = mfccc * c1o2 + mfccb * (vvz + c1o2) + (mfcca + c1o9 * oMdrho) * (vz2 + vvz) * c1o2;
+ mfcca = m0;
+ mfccb = m1;
+ mfccc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ //mit 1/6, 2/3, 1/6, 0, 0, 0, 1/18, 2/9, 1/18 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ m0 = mfaca * c1o2 + mfaba * (vvy - c1o2) + (mfaaa + c1o6 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfaca - 2. * mfaba * vvy + mfaaa * (1. - vy2) - c1o6 * oMdrho * vy2;
+ m2 = mfaca * c1o2 + mfaba * (vvy + c1o2) + (mfaaa + c1o6 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfaaa = m0;
+ mfaba = m1;
+ mfaca = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfacb * c1o2 + mfabb * (vvy - c1o2) + (mfaab + c2o3 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfacb - 2. * mfabb * vvy + mfaab * (1. - vy2) - c2o3 * oMdrho * vy2;
+ m2 = mfacb * c1o2 + mfabb * (vvy + c1o2) + (mfaab + c2o3 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfaab = m0;
+ mfabb = m1;
+ mfacb = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfacc * c1o2 + mfabc * (vvy - c1o2) + (mfaac + c1o6 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfacc - 2. * mfabc * vvy + mfaac * (1. - vy2) - c1o6 * oMdrho * vy2;
+ m2 = mfacc * c1o2 + mfabc * (vvy + c1o2) + (mfaac + c1o6 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfaac = m0;
+ mfabc = m1;
+ mfacc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfbca * c1o2 + mfbba * (vvy - c1o2) + mfbaa * (vy2 - vvy) * c1o2;
+ m1 = -mfbca - 2. * mfbba * vvy + mfbaa * (1. - vy2);
+ m2 = mfbca * c1o2 + mfbba * (vvy + c1o2) + mfbaa * (vy2 + vvy) * c1o2;
+ mfbaa = m0;
+ mfbba = m1;
+ mfbca = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbcb * c1o2 + mfbbb * (vvy - c1o2) + mfbab * (vy2 - vvy) * c1o2;
+ m1 = -mfbcb - 2. * mfbbb * vvy + mfbab * (1. - vy2);
+ m2 = mfbcb * c1o2 + mfbbb * (vvy + c1o2) + mfbab * (vy2 + vvy) * c1o2;
+ mfbab = m0;
+ mfbbb = m1;
+ mfbcb = m2;
+ /////////b//////////////////////////////////////////////////////////////////////////
+ m0 = mfbcc * c1o2 + mfbbc * (vvy - c1o2) + mfbac * (vy2 - vvy) * c1o2;
+ m1 = -mfbcc - 2. * mfbbc * vvy + mfbac * (1. - vy2);
+ m2 = mfbcc * c1o2 + mfbbc * (vvy + c1o2) + mfbac * (vy2 + vvy) * c1o2;
+ mfbac = m0;
+ mfbbc = m1;
+ mfbcc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcca * c1o2 + mfcba * (vvy - c1o2) + (mfcaa + c1o18 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfcca - 2. * mfcba * vvy + mfcaa * (1. - vy2) - c1o18 * oMdrho * vy2;
+ m2 = mfcca * c1o2 + mfcba * (vvy + c1o2) + (mfcaa + c1o18 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfcaa = m0;
+ mfcba = m1;
+ mfcca = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfccb * c1o2 + mfcbb * (vvy - c1o2) + (mfcab + c2o9 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfccb - 2. * mfcbb * vvy + mfcab * (1. - vy2) - c2o9 * oMdrho * vy2;
+ m2 = mfccb * c1o2 + mfcbb * (vvy + c1o2) + (mfcab + c2o9 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfcab = m0;
+ mfcbb = m1;
+ mfccb = m2;
+ /////////c//////////////////////////////////////////////////////////////////////////
+ m0 = mfccc * c1o2 + mfcbc * (vvy - c1o2) + (mfcac + c1o18 * oMdrho) * (vy2 - vvy) * c1o2;
+ m1 = -mfccc - 2. * mfcbc * vvy + mfcac * (1. - vy2) - c1o18 * oMdrho * vy2;
+ m2 = mfccc * c1o2 + mfcbc * (vvy + c1o2) + (mfcac + c1o18 * oMdrho) * (vy2 + vvy) * c1o2;
+ mfcac = m0;
+ mfcbc = m1;
+ mfccc = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ //mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ m0 = mfcaa * c1o2 + mfbaa * (vvx - c1o2) + (mfaaa + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcaa - 2. * mfbaa * vvx + mfaaa * (1. - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfcaa * c1o2 + mfbaa * (vvx + c1o2) + (mfaaa + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaaa = m0;
+ mfbaa = m1;
+ mfcaa = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcba * c1o2 + mfbba * (vvx - c1o2) + (mfaba + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcba - 2. * mfbba * vvx + mfaba * (1. - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfcba * c1o2 + mfbba * (vvx + c1o2) + (mfaba + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaba = m0;
+ mfbba = m1;
+ mfcba = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcca * c1o2 + mfbca * (vvx - c1o2) + (mfaca + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcca - 2. * mfbca * vvx + mfaca * (1. - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfcca * c1o2 + mfbca * (vvx + c1o2) + (mfaca + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaca = m0;
+ mfbca = m1;
+ mfcca = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcab * c1o2 + mfbab * (vvx - c1o2) + (mfaab + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcab - 2. * mfbab * vvx + mfaab * (1. - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfcab * c1o2 + mfbab * (vvx + c1o2) + (mfaab + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaab = m0;
+ mfbab = m1;
+ mfcab = m2;
+ ///////////b////////////////////////////////////////////////////////////////////////
+ m0 = mfcbb * c1o2 + mfbbb * (vvx - c1o2) + (mfabb + c4o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcbb - 2. * mfbbb * vvx + mfabb * (1. - vx2) - c4o9 * oMdrho * vx2;
+ m2 = mfcbb * c1o2 + mfbbb * (vvx + c1o2) + (mfabb + c4o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfabb = m0;
+ mfbbb = m1;
+ mfcbb = m2;
+ ///////////b////////////////////////////////////////////////////////////////////////
+ m0 = mfccb * c1o2 + mfbcb * (vvx - c1o2) + (mfacb + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfccb - 2. * mfbcb * vvx + mfacb * (1. - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfccb * c1o2 + mfbcb * (vvx + c1o2) + (mfacb + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfacb = m0;
+ mfbcb = m1;
+ mfccb = m2;
+ ////////////////////////////////////////////////////////////////////////////////////
+ ////////////////////////////////////////////////////////////////////////////////////
+ m0 = mfcac * c1o2 + mfbac * (vvx - c1o2) + (mfaac + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcac - 2. * mfbac * vvx + mfaac * (1. - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfcac * c1o2 + mfbac * (vvx + c1o2) + (mfaac + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfaac = m0;
+ mfbac = m1;
+ mfcac = m2;
+ ///////////c////////////////////////////////////////////////////////////////////////
+ m0 = mfcbc * c1o2 + mfbbc * (vvx - c1o2) + (mfabc + c1o9 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfcbc - 2. * mfbbc * vvx + mfabc * (1. - vx2) - c1o9 * oMdrho * vx2;
+ m2 = mfcbc * c1o2 + mfbbc * (vvx + c1o2) + (mfabc + c1o9 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfabc = m0;
+ mfbbc = m1;
+ mfcbc = m2;
+ ///////////c////////////////////////////////////////////////////////////////////////
+ m0 = mfccc * c1o2 + mfbcc * (vvx - c1o2) + (mfacc + c1o36 * oMdrho) * (vx2 - vvx) * c1o2;
+ m1 = -mfccc - 2. * mfbcc * vvx + mfacc * (1. - vx2) - c1o36 * oMdrho * vx2;
+ m2 = mfccc * c1o2 + mfbcc * (vvx + c1o2) + (mfacc + c1o36 * oMdrho) * (vx2 + vvx) * c1o2;
+ mfacc = m0;
+ mfbcc = m1;
+ mfccc = m2;
+
+ //////////////////////////////////////////////////////////////////////////
+ //proof correctness
+ //////////////////////////////////////////////////////////////////////////
+#ifdef PROOF_CORRECTNESS
+ LBMReal rho_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 = fabs(rho - rho_post);
+ LBMReal dif = rho - rho_post;
+#ifdef SINGLEPRECISION
+ if (dif > 10.0E-7 || dif < -10.0E-7)
+#else
+ if (dif > 10.0E-15 || dif < -10.0E-15)
+#endif
+ {
+ UB_THROW(UbException(UB_EXARGS, "rho=" + UbSystem::toString(rho) + ", rho_post=" + UbSystem::toString(rho_post)
+ + " dif=" + UbSystem::toString(dif)
+ + " rho is not correct for node " + UbSystem::toString(x1) + "," + UbSystem::toString(x2) + "," + UbSystem::toString(x3)));
+ //UBLOG(logERROR,"LBMKernelETD3Q27CCLB::collideAll(): rho is not correct for node "+UbSystem::toString(x1)+","+UbSystem::toString(x2)+","+UbSystem::toString(x3));
+ //exit(EXIT_FAILURE);
+ }
+#endif
+ //////////////////////////////////////////////////////////////////////////
+ //write distribution
+ //////////////////////////////////////////////////////////////////////////
+ (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) = mfabb;
+ (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) = mfbab;
+ (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) = mfbba;
+ (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) = mfaab;
+ (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab;
+ (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) = mfaba;
+ (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba;
+ (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa;
+ (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca;
+ (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa;
+ (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa;
+ (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca;
+ (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca;
+
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
+
+ (*this->zeroDistributionsF)(x1, x2, x3) = mfbbb;
+ }
+ }
+ }
+ }
+
+ }
+}
+
+//SPtr<LBMKernel> ThixotropyModelLBMKernel2::clone()
+//{
+// SPtr<LBMKernel> kernel(new ThixotropyModelLBMKernel2());
+// kernel->setNX(nx);
+// kernel->setCollisionFactor(collFactor);
+// collFactorF = collFactor;
+// dynamicPointerCast<ThixotropyModelLBMKernel2>(kernel)->initDataSet();
+// kernel->setBCProcessor(bcProcessor->clone(kernel));
+// kernel->setWithForcing(withForcing);
+// kernel->setForcingX1(muForcingX1);
+// kernel->setForcingX2(muForcingX2);
+// kernel->setForcingX3(muForcingX3);
+// kernel->setIndex(ix1, ix2, ix3);
+// kernel->setDeltaT(deltaT);
+//
+// return kernel;
+//}
+
+double ThixotropyModelLBMKernel2::getCalculationTime()
+{
+ return timer.getTotalTime();
+}
+
+void ThixotropyModelLBMKernel2::swapDistributions()
+{
+ LBMKernel::swapDistributions();
+}
+
+void ThixotropyModelLBMKernel2::initDataSet()
+{
+ SPtr<DistributionArray3D> df(new D3Q27EsoTwist3DSplittedVector(nx[0] + 2, nx[1] + 2, nx[2] + 2, -999.0));
+ dataSet->setFdistributions(df);
+}
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.h b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.h
new file mode 100644
index 000000000..17fa47797
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.h
@@ -0,0 +1,51 @@
+#ifndef ThixotropyModelLBMKernel2_H
+#define ThixotropyModelLBMKernel2_H
+
+#include "LBMKernel.h"
+#include "BCProcessor.h"
+#include "D3Q27System.h"
+#include "basics/utilities/UbTiming.h"
+#include "basics/container/CbArray4D.h"
+#include "basics/container/CbArray3D.h"
+
+class ThixotropyModelLBMKernel2;
+
+//! \brief Base class for model of thixotropy based on K16. Use Template Method design pattern for Implementation of different models.
+//! \author K. Kutscher, M. Geier
+class ThixotropyModelLBMKernel2 : public LBMKernel
+{
+public:
+ ThixotropyModelLBMKernel2();
+ virtual ~ThixotropyModelLBMKernel2();
+ void calculate(int step);
+ virtual SPtr<LBMKernel> clone() { UB_THROW(UbException("SPtr<LBMKernel> clone() - belongs in the derived class")); };
+ double getCalculationTime();
+
+ void swapDistributions();
+
+protected:
+ void initDataSet();
+
+ virtual LBMReal getThyxotropyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho) const { UB_THROW(UbException("LBMReal getThyxotropyCollFactor() - belongs in the derived class")); }
+
+ LBMReal f[D3Q27System::ENDF + 1];
+
+ UbTimer timer;
+
+ LBMReal OxyyMxzz;
+
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr localDistributionsF;
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributionsF;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr zeroDistributionsF;
+
+ mu::value_type muX1, muX2, muX3;
+ mu::value_type muDeltaT;
+ mu::value_type muNu;
+ LBMReal forcingX1;
+ LBMReal forcingX2;
+ LBMReal forcingX3;
+
+ bool test;
+};
+
+#endif
--
GitLab