diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 2a40b73bd14ef909048725fd35332a5a3cb24825..d097031970a174f146b25126d80d14aea84752ad 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -65,4 +65,6 @@ add_subdirectory(${APPS_ROOT_CPU}/LaminarTubeFlow)
#add_subdirectory(Applications/OrganPipe)
#add_subdirectory(Applications/LidDrivenCavity)
add_subdirectory(${APPS_ROOT_CPU}/HerschelBulkleySphere)
+add_subdirectory(${APPS_ROOT_CPU}/HerschelBulkleyModel)
+add_subdirectory(${APPS_ROOT_CPU}/rheometer)
diff --git a/apps/cpu/HerschelBulkleyModel/CMakeLists.txt b/apps/cpu/HerschelBulkleyModel/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..323d1e24ab308f0090f2849dcc9f216b5658c567
--- /dev/null
+++ b/apps/cpu/HerschelBulkleyModel/CMakeLists.txt
@@ -0,0 +1,10 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
+
+########################################################
+## C++ PROJECT ###
+########################################################
+PROJECT(HerschelBulkleyModel)
+
+INCLUDE(${APPS_ROOT_CPU}/IncludsList.cmake)
+
+vf_add_library(BUILDTYPE binary DEPENDS VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} FILES hbflow.cpp )
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleyModel/hbf1.cfg b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..830fe2ae98a124a7a81752aa183eeb496294b7e7
--- /dev/null
+++ b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
@@ -0,0 +1,32 @@
+pathname = d:/temp/Herschel-BulkleyConvergenceStudy/HerschelBulkley32_n_0.4_Re_1_Bn_0.01_K15
+
+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 = 100000
+endTime = 50000000
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleyModel/hbflow.cpp b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..0e5897fa6d94aae0342cd2e03ba66dc7c2f0231d
--- /dev/null
+++ b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
@@ -0,0 +1,305 @@
+#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 = (U * d) / (Re * std::pow(Gamma, n - 1));
+ double tau0 = 0;// 1e-12;// 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()));
+
+ //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());
+ 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/HerschelBulkleySphere/hbsphere.cfg b/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
index 953f2c2872e2dff565438bc5d42aad24b0359e31..ee4b49ef79ab0a589f355725a5b5bd61a376ee4c 100644
--- a/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
+++ b/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
@@ -1,11 +1,11 @@
-outputPath = d:/temp/Herschel-Bulkley_sphere_MP_test
+outputPath = d:/temp/Herschel-Bulkley_sphere_Re_200_n_0.3_Bn_0.01_slip
numOfThreads = 4
availMem = 8e9
logToFile = false
blocknx = 25 25 25
-boundingBox = 250 50 50 #30*20=600**3=216000000
+boundingBox = 125 50 50
deltax = 1
refineLevel = 0
@@ -14,8 +14,8 @@ radius = 5
sphereCenter = 25 25 25
velocity = 1e-3
-n = 0.9
-Re = 1
+n = 0.3
+Re = 200
Bn = 0.01
@@ -25,5 +25,5 @@ restartStep = 100000
cpStart = 100000
cpStep = 100000
-outTime = 10000
-endTime = 1000000
\ No newline at end of file
+outTime = 10000
+endTime = 100000
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
index a31267160369ca0947ed78dd4780e141e37e2941..bfbf4d9ab58f168d9383c42b36aec614e8b85831 100644
--- a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
+++ b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
@@ -162,10 +162,10 @@ void bflow(string configname)
//////////////////////////////////////////////////////////////////////////
//restart
SPtr<UbScheduler> mSch(new UbScheduler(cpStep, cpStart));
- SPtr<MPIIOMigrationBECoProcessor> restartCoProcessor(new MPIIOMigrationBECoProcessor(grid, mSch, outputPath, comm));
+ SPtr<MPIIOMigrationCoProcessor> restartCoProcessor(new MPIIOMigrationCoProcessor(grid, mSch, outputPath, comm));
restartCoProcessor->setLBMKernel(kernel);
restartCoProcessor->setBCProcessor(bcProc);
- restartCoProcessor->setNu(k);
+ //restartCoProcessor->setNu(k);
//////////////////////////////////////////////////////////////////////////
if (myid == 0)
diff --git a/apps/cpu/rheometer/CMakeLists.txt b/apps/cpu/rheometer/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..66f007ecc6a0c33aa9bf17bed0f6fd171ac1e045
--- /dev/null
+++ b/apps/cpu/rheometer/CMakeLists.txt
@@ -0,0 +1,10 @@
+CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
+
+########################################################
+## C++ PROJECT ###
+########################################################
+PROJECT(rheometer)
+
+INCLUDE(${APPS_ROOT_CPU}/IncludsList.cmake)
+
+vf_add_library(BUILDTYPE binary DEPENDS VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} FILES rheometer.cpp )
\ No newline at end of file
diff --git a/apps/cpu/rheometer/rheometer.cfg b/apps/cpu/rheometer/rheometer.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..5f95613e5a25f27881a1e261eb821c6edb5cd6fb
--- /dev/null
+++ b/apps/cpu/rheometer/rheometer.cfg
@@ -0,0 +1,29 @@
+outputPath = d:/temp/rheometer
+
+numOfThreads = 4
+availMem = 8e9
+logToFile = false
+
+blocknx = 8 8 2
+boundingBox = 32 32 2
+deltax = 1
+
+refineLevel = 0
+
+radius = 5
+sphereCenter = 25 25 25
+
+velocity = 1e-3
+n = 0.3
+Re = 10
+Bn = 0.01
+
+
+newStart = true
+restartStep = 100000
+
+cpStart = 100000
+cpStep = 100000
+
+outTime = 1000
+endTime = 100000
\ No newline at end of file
diff --git a/apps/cpu/rheometer/rheometer.cpp b/apps/cpu/rheometer/rheometer.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..69de586fbb901883209cd772ef0cf45fa4031c1d
--- /dev/null
+++ b/apps/cpu/rheometer/rheometer.cpp
@@ -0,0 +1,408 @@
+#include <iostream>
+#include <string>
+
+#include <VirtualFluids.h>
+
+using namespace std;
+
+
+void bflow(string configname)
+{
+ try
+ {
+ ConfigurationFile config;
+ config.load(configname);
+
+ string outputPath = config.getValue<string>("outputPath");
+ 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 radius = config.getValue<double>("radius");
+ 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 n = config.getValue<double>("n");
+ double Re = config.getValue<double>("Re");
+ double Bn = config.getValue<double>("Bn");
+ vector<double> sphereCenter = config.getVector<double>("sphereCenter");
+
+ SPtr<Communicator> comm = MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ if (logToFile)
+ {
+#if defined(__unix__)
+ if (myid == 0)
+ {
+ const char* str = outputPath.c_str();
+ mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ }
+#endif
+
+ if (myid == 0)
+ {
+ stringstream logFilename;
+ logFilename << outputPath + "/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 = -boundingBox[0]/2.0;
+ //double g_minX2 = -boundingBox[1] / 2.0;
+ //double g_minX3 = -boundingBox[2]/2.0;
+
+ //double g_maxX1 = boundingBox[0]/2.0;
+ //double g_maxX2 = boundingBox[1]/2.0;
+ //double g_maxX3 = boundingBox[2]/2.0;
+
+ double blockLength = 3.0 * deltax;
+
+ double d = 2.0 * radius;
+ double U = velocity;
+ double Gamma = U / d;
+
+ double k = (U * d) / (Re);
+
+ //double k = (U * d) / (Re * std::pow(Gamma, n - 1));
+ double tau0 = 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 forcing = 8e-7;
+
+ double omegaMin = 1.0e-8;
+
+ SPtr<Thixotropy> thix = Thixotropy::getInstance();
+ thix->setPowerIndex(n);
+ thix->setViscosityParameter(k);
+ thix->setYieldStress(tau0);
+ thix->setOmegaMin(omegaMin);
+
+ SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
+ //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new HerschelBulkleyModelNoSlipBCAlgorithm()));
+ //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.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);
+
+ 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 SimpleVelocityBCAlgorithm()));
+ //velocityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityWithDensityBCAlgorithm()));
+
+ //SPtr<BCAdapter> densityBCAdapter(new DensityBCAdapter());
+ //densityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
+ ////densityBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonReflectingOutflowBCAlgorithm()));
+
+
+ //BS visitor
+ BoundaryConditionsBlockVisitor bcVisitor;
+ bcVisitor.addBC(noSlipBCAdapter);
+ //bcVisitor.addBC(slipBCAdapter);
+ bcVisitor.addBC(velocityBCAdapter);
+ //bcVisitor.addBC(densityBCAdapter);
+
+ SPtr<BCProcessor> bcProc;
+ bcProc = SPtr<BCProcessor>(new BCProcessor());
+
+ 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());
+ kernel->setBCProcessor(bcProc);
+ //kernel->setForcingX1(forcing);
+ //kernel->setWithForcing(true);
+
+ SPtr<Grid3D> grid(new Grid3D(comm));
+ grid->setPeriodicX1(false);
+ 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(), outputPath + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+
+ //sphere
+ //SPtr<GbObject3D> sphere(new GbSphere3D(sphereCenter[0], sphereCenter[1], sphereCenter[2], radius));
+ //GbSystem3D::writeGeoObject(sphere.get(), outputPath + "/geo/sphere", WbWriterVtkXmlBinary::getInstance());
+ //SPtr<D3Q27Interactor> sphereInt(new D3Q27Interactor(sphere, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ //////////////////////////////////////////////////////////////////////////
+ //restart
+ SPtr<UbScheduler> mSch(new UbScheduler(cpStep, cpStart));
+ SPtr<MPIIOMigrationCoProcessor> restartCoProcessor(new MPIIOMigrationCoProcessor(grid, mSch, outputPath, comm));
+ restartCoProcessor->setLBMKernel(kernel);
+ restartCoProcessor->setBCProcessor(bcProc);
+ //restartCoProcessor->setNu(k);
+ //////////////////////////////////////////////////////////////////////////
+
+ if (myid == 0)
+ {
+ UBLOG(logINFO, "Parameters:");
+ //UBLOG(logINFO, "forcing = " << forcing);
+ UBLOG(logINFO, "rho = " << rhoLB);
+ 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, "number of threads = " << numOfThreads);
+ UBLOG(logINFO, "number of processes = " << comm->getNumberOfProcesses());
+ UBLOG(logINFO, "blocknx = " << blocknx[0] << " " << blocknx[1] << " " << blocknx[2]);
+ UBLOG(logINFO, "boundingBox = " << boundingBox[0] << " " << boundingBox[1] << " " << boundingBox[2]);
+ UBLOG(logINFO, "sphereCenter = " << sphereCenter[0] << " " << sphereCenter[1] << " " << sphereCenter[2]);
+ UBLOG(logINFO, "output path = " << outputPath);
+ UBLOG(logINFO, "Preprozess - start");
+ }
+
+ if (newStart)
+ {
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ if (refineLevel > 0)
+ {
+ GbCuboid3DPtr refCube(new GbCuboid3D(-10, -10, -10, 0, 0, 0));
+ if (myid == 0) GbSystem3D::writeGeoObject(refCube.get(), outputPath + "/geo/refCube", WbWriterVtkXmlASCII::getInstance());
+
+ if (myid == 0) UBLOG(logINFO, "Refinement - start");
+ RefineCrossAndInsideGbObjectHelper refineHelper(grid, refineLevel, comm);
+ //refineHelper.addGbObject(sphere, refineLevel);
+ refineHelper.addGbObject(refCube, refineLevel);
+ refineHelper.refine();
+ 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));
+ if (myid == 0) GbSystem3D::writeGeoObject(wallZmin.get(), outputPath + "/geo/wallZmin", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr wallZmax(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_maxX3, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ if (myid == 0) GbSystem3D::writeGeoObject(wallZmax.get(), outputPath + "/geo/wallZmax", WbWriterVtkXmlASCII::getInstance());
+
+ //GbCuboid3DPtr wallYmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_minX2, g_maxX3 + blockLength));
+ //if (myid == 0) GbSystem3D::writeGeoObject(wallYmin.get(), outputPath + "/geo/wallYmin", WbWriterVtkXmlASCII::getInstance());
+
+ //GbCuboid3DPtr wallYmax(new GbCuboid3D(g_minX1 - blockLength, g_maxX2, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ //if (myid == 0) GbSystem3D::writeGeoObject(wallYmax.get(), outputPath + "/geo/wallYmax", WbWriterVtkXmlASCII::getInstance());
+
+ //GbCuboid3DPtr wallXmin(new GbCuboid3D(g_minX1 - blockLength, g_minX2 - blockLength, g_minX3 - blockLength, g_minX1, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ //if (myid == 0) GbSystem3D::writeGeoObject(wallXmin.get(), outputPath + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
+
+ //GbCuboid3DPtr wallXmax(new GbCuboid3D(g_maxX1, g_minX2 - blockLength, g_minX3 - blockLength, g_maxX1 + blockLength, g_maxX2 + blockLength, g_maxX3 + blockLength));
+ //if (myid == 0) GbSystem3D::writeGeoObject(wallXmax.get(), outputPath + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
+
+ //wall interactors
+ SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ //SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, slipBCAdapter, Interactor3D::SOLID));
+ //SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, slipBCAdapter, Interactor3D::SOLID));
+
+ //SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, velocityBCAdapter, Interactor3D::SOLID));
+ //SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, densityBCAdapter, 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(wallXminInt);
+ //intHelper.addInteractor(wallXmaxInt);
+ //intHelper.addInteractor(wallYminInt);
+ //intHelper.addInteractor(wallYmaxInt);
+ //intHelper.addInteractor(wallZminInt);
+ //intHelper.addInteractor(wallZmaxInt);
+ intHelper.addInteractor(statorInt);
+ intHelper.addInteractor(rotorInt);
+ intHelper.selectBlocks();
+ if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
+ //////////////////////////////////////
+
+ SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), outputPath, 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);
+
+ if (refineLevel > 0)
+ {
+ SetUndefinedNodesBlockVisitor undefNodesVisitor;
+ grid->accept(undefNodesVisitor);
+ }
+
+ //BC
+ intHelper.setBC();
+
+ //initialization of distributions
+ InitDistributionsBlockVisitor initVisitor;
+ grid->accept(initVisitor);
+
+
+ if (myid == 0) UBLOG(logINFO, "Preprozess - end");
+ }
+ else
+ {
+ restartCoProcessor->restart((int)restartStep);
+ grid->setTimeStep(restartStep);
+ //SetBcBlocksBlockVisitor v(sphereInt);
+ //grid->accept(v);
+ //sphereInt->initInteractor();
+ }
+
+ omp_set_num_threads(numOfThreads);
+
+ //set connectors
+ InterpolationProcessorPtr iProcessor(new ThixotropyInterpolationProcessor());
+ static_pointer_cast<ThixotropyInterpolationProcessor>(iProcessor)->setOmegaMin(thix->getOmegaMin());
+ 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, outputPath, WbWriterVtkXmlASCII::getInstance(), comm);
+ ppgeo.process(0);
+
+ 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<UbScheduler> visSch(new UbScheduler(10,1));
+ 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<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(writeMQCoProcessor);
+ //calculator->addCoProcessor(writeThixotropicMQCoProcessor);
+ 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/src/cpu/VirtualFluids.h b/src/cpu/VirtualFluids.h
index 2b4672355e8fce7dcd81a00ab3b98fb14ce74aa0..77a8a551e33b52b03fdf1c6d90d619495387194d 100644
--- a/src/cpu/VirtualFluids.h
+++ b/src/cpu/VirtualFluids.h
@@ -131,6 +131,7 @@
#include <BoundaryConditions/BinghamModelNoSlipBCAlgorithm.h>
#include <BoundaryConditions/HerschelBulkleyModelNoSlipBCAlgorithm.h>
#include <BoundaryConditions/SimpleSlipBCAlgorithm.h>
+#include <BoundaryConditions/PowellEyringModelNoSlipBCAlgorithm.h>
#include <Connectors/Block3DConnector.h>
#include <Connectors/Block3DConnectorFactory.h>
@@ -224,6 +225,8 @@
#include <LBM/HerschelBulkleyModelLBMKernel.h>
#include <LBM/ThixotropyInterpolationProcessor.h>
#include <LBM/RheologyK17LBMKernel.h>
+#include <LBM/PowellEyringModelLBMKernel.h>
+
#include <geometry3d/CoordinateTransformation3D.h>
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
index 4ff5745026f87aa428b64fd084954bfff4f4a4de..f9427d41a140a275ac45bf35c5ae5a0c049706f7 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCAlgorithm.h
@@ -65,7 +65,7 @@ public:
static const char HerschelBulkleyModelNoSlipBCAlgorithm = 15;
static const char SimpleVelocityBCAlgorithm = 16;
static const char SimpleSlipBCAlgorithm = 17;
-
+ static const char PowellEyringModelNoSlipBCAlgorithm = 18;
public:
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/PowellEyringModelNoSlipBCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/PowellEyringModelNoSlipBCAlgorithm.h
new file mode 100644
index 0000000000000000000000000000000000000000..f94266be4d87b33fc5858e1f48263d3cdbc082c5
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/PowellEyringModelNoSlipBCAlgorithm.h
@@ -0,0 +1,27 @@
+#ifndef PowellEyringModelNoSlipBCAlgorithm_h__
+#define PowellEyringModelNoSlipBCAlgorithm_h__
+
+#include "ThixotropyNoSlipBCAlgorithm.h"
+#include "Thixotropy.h"
+
+class PowellEyringModelNoSlipBCAlgorithm : public ThixotropyNoSlipBCAlgorithm
+{
+public:
+ PowellEyringModelNoSlipBCAlgorithm()
+ {
+ BCAlgorithm::type = BCAlgorithm::PowellEyringModelNoSlipBCAlgorithm;
+ BCAlgorithm::preCollision = true;
+ }
+ ~PowellEyringModelNoSlipBCAlgorithm() {}
+ SPtr<BCAlgorithm> clone() override
+ {
+ SPtr<BCAlgorithm> bc(new PowellEyringModelNoSlipBCAlgorithm());
+ return bc;
+ }
+protected:
+ LBMReal getThyxotropyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho) const override
+ {
+ return Thixotropy::getHerschelBulkleyCollFactor(omegaInf, shearRate, drho);
+ }
+};
+#endif // PowellEyringModelNoSlipBCAlgorithm_h__
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
index 4ff7be254a6964f6ece23d209105a32c58938741..fe5f317886271365aaf5d7bcb8304de1711ee95b 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
@@ -212,7 +212,8 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
distributionsF->getDistribution(f, ix1, ix2, ix3);
LBMReal rho = D3Q27System::getDensity(f);
LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
- LBMReal omega = Thixotropy::getHerschelBulkleyCollFactor(collFactor, shearRate, rho);
+ //LBMReal omega = Thixotropy::getHerschelBulkleyCollFactor(collFactor, shearRate, rho);
+ LBMReal omega = Thixotropy::getPowellEyringCollFactor(collFactor, shearRate, rho);
LBMReal viscosity = (omega == 0) ? 0 : c1o3 * (c1/omega-c1o2);
data[index++].push_back(viscosity);
diff --git a/src/cpu/VirtualFluidsCore/LBM/PowellEyringModelLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/PowellEyringModelLBMKernel.h
new file mode 100644
index 0000000000000000000000000000000000000000..02932d808e2cbc242617405979435c78fd6386d9
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/PowellEyringModelLBMKernel.h
@@ -0,0 +1,38 @@
+#ifndef PowellEyringModelLBMKernel_H
+#define PowellEyringModelLBMKernel_H
+
+#include "ThixotropyModelLBMKernel.h"
+#include "Thixotropy.h"
+
+//! \brief Cumulant LBM kernel + Herschel-Bulkley plastic model
+//! \author K. Kutscher, M. Geier
+class PowellEyringModelLBMKernel : public ThixotropyModelLBMKernel
+{
+public:
+ PowellEyringModelLBMKernel() {};
+ ~PowellEyringModelLBMKernel() {};
+ SPtr<LBMKernel> clone() override
+ {
+ SPtr<LBMKernel> kernel(new PowellEyringModelLBMKernel());
+ kernel->setNX(nx);
+ kernel->setCollisionFactor(collFactor);
+ dynamicPointerCast<PowellEyringModelLBMKernel>(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;
+ }
+protected:
+ LBMReal getThyxotropyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho) const override
+ {
+ return Thixotropy::getPowellEyringCollFactor(omegaInf, shearRate, drho);
+ }
+};
+
+
+#endif
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp
index a90f00afc8bd33f92f9be9ebe4bb592e83781f67..f28c34f0ddc8be17807a32118074a0c07f36bcfa 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp
@@ -5,6 +5,9 @@ LBMReal Thixotropy::tau0 = 0;
LBMReal Thixotropy::k = 0;
LBMReal Thixotropy::n = 1;
LBMReal Thixotropy::omegaMin = 0;
+LBMReal Thixotropy::beta = 0;
+LBMReal Thixotropy::c = 0;
+LBMReal Thixotropy::mu0 = 0;
//////////////////////////////////////////////////////////////////////////
SPtr<Thixotropy> Thixotropy::getInstance()
@@ -48,6 +51,36 @@ LBMReal Thixotropy::getOmegaMin() const
return omegaMin;
}
+void Thixotropy::setBeta(LBMReal PowellEyringBeta)
+{
+ beta = PowellEyringBeta;
+}
+
+LBMReal Thixotropy::getBeta() const
+{
+ return beta;
+}
+
+void Thixotropy::setC(LBMReal PowellEyringC)
+{
+ c = PowellEyringC;
+}
+
+LBMReal Thixotropy::getC() const
+{
+ return c;
+}
+
+void Thixotropy::setMu0(LBMReal mu)
+{
+ mu0 = mu;
+}
+
+LBMReal Thixotropy::getMu0() const
+{
+ return mu0;
+}
+
Thixotropy::Thixotropy()
{
}
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
index 677e488e68359eef7774b21384ca2a14af93e1b5..b731337e696451a365db73ebc200d7fc841c7514 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
@@ -4,6 +4,7 @@
#include <PointerDefinitions.h>
#include <LBMSystem.h>
#include <UbMath.h>
+#include <math.h>
class Thixotropy
{
@@ -23,9 +24,19 @@ public:
void setOmegaMin(LBMReal omegaMin);
LBMReal getOmegaMin() const;
+ void setBeta(LBMReal PowellEyringBeta);
+ LBMReal getBeta() const;
+
+ void setC(LBMReal PowellEyringC);
+ LBMReal getC() const;
+
+ void setMu0(LBMReal mu);
+ LBMReal getMu0() const;
+
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);
+ static LBMReal getPowellEyringCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
private:
Thixotropy();
@@ -35,6 +46,9 @@ private:
static LBMReal k;
static LBMReal n;
static LBMReal omegaMin;
+ static LBMReal beta;
+ static LBMReal c;
+ static LBMReal mu0;
};
//////////////////////////////////////////////////////////////////////////
@@ -53,12 +67,19 @@ inline LBMReal Thixotropy::getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMRea
LBMReal gammaDot = shearRate;
LBMReal omega = omegaInf;
LBMReal epsilon = 1;
+ LBMReal gammaDotPowN = std::pow(gammaDot, n);
while (epsilon > 1e-10)
{
LBMReal omegaOld = omega;
- LBMReal gammaDotPowN = std::pow(gammaDot, n);
- LBMReal omegaByOmegaInfPowN = std::pow(omega / omegaInf, n);
+ LBMReal omegaByOmegaInfPowN = std::pow(omega / omegaInf, n);/*
+ LBMReal gammaDotPowOneMinusN = std::pow(gammaDot,1- n);
+ LBMReal omegaByOmegaInfPowOneMinusN = std::pow(omega / omegaInf, 1-n);
+ LBMReal numeratorA = (2.0* k * omegaInf + cs2 * gammaDotPowOneMinusN * omegaByOmegaInfPowOneMinusN *omegaInf* rho );
+ LBMReal numeratorB = ( cs2 * gammaDot * ( - 2.0) * rho + 2.0 * omegaInf * tau0);
+ LBMReal denominatorA = (2.0 * k * n * omegaInf + cs2 * gammaDot * rho * omegaInf* gammaDotPowOneMinusN * omegaByOmegaInfPowOneMinusN) + UbMath::Epsilon<LBMReal>::val();
+ LBMReal denominatorB = (2.0 * k * n * gammaDotPowN * omegaByOmegaInfPowN * omegaInf + cs2 * gammaDot * rho * omega) + UbMath::Epsilon<LBMReal>::val();
+ omega = omega - omega *( numeratorA / denominatorA+ numeratorB / denominatorB);*/
LBMReal numerator = (2.0 * gammaDotPowN * k * omegaByOmegaInfPowN * omegaInf + cs2 * gammaDot * (omega - 2.0) * rho + 2.0 * omegaInf * tau0);
LBMReal denominator = (2.0 * k * n * gammaDotPowN * omegaByOmegaInfPowN * omegaInf + cs2 * gammaDot * rho * omega) + UbMath::Epsilon<LBMReal>::val();
omega = omega - omega * numerator / denominator;
@@ -78,4 +99,30 @@ inline LBMReal Thixotropy::getHerschelBulkleyCollFactorBackward(LBMReal shearRat
return 1.0 / ((tau0 + k * std::pow(gamma, n)) / (cs2 * rho * gamma) + UbMath::c1o2);
}
+//////////////////////////////////////////////////////////////////////////
+inline LBMReal Thixotropy::getPowellEyringCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho)
+{
+ using namespace UbMath;
+ LBMReal cs2 = c1o3; // UbMath::one_over_sqrt3* UbMath::one_over_sqrt3;
+ LBMReal rho = c1 + drho;
+ LBMReal gammaDot = shearRate;
+ LBMReal omega = omegaInf;
+ LBMReal epsilon = 1;
+
+ while (epsilon > 1e-10)
+ {
+ LBMReal omegaOld = omega;
+ epsilon = std::abs(omega - omegaOld);
+
+ LBMReal numerator = c*sqrt(c1+(gammaDot*gammaDot*omega*omega)/(c*c*omegaInf*omegaInf))*(beta*(c2*gammaDot*mu0*omega+cs2*gammaDot*(omega-c2)*rho+c2*omegaInf*tau0)+c2*omegaInf*(asinh((gammaDot*omega)/(c*omegaInf))));
+
+ LBMReal denominator = gammaDot*(c2+beta*c*sqrt(c1+(gammaDot*gammaDot*omega*omega)/(c*c*omegaInf*omegaInf))*(c2*mu0+cs2*rho)) + UbMath::Epsilon<LBMReal>::val();
+
+ omega = omega - numerator / denominator;
+
+ omega = (omega < UbMath::zeroReal) ? UbMath::c1o2 * omegaOld : omega;
+ }
+
+ return omega;
+}
#endif