add Powell-Eyring model / add rheometer setup

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)
 

apps/cpu/HerschelBulkleyModel/CMakeLists.txt

+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

apps/cpu/HerschelBulkleyModel/hbf1.cfg

+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

apps/cpu/HerschelBulkleyModel/hbflow.cpp

+#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;
+}

apps/cpu/HerschelBulkleySphere/hbsphere.cfg

-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

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)

apps/cpu/rheometer/CMakeLists.txt

+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

apps/cpu/rheometer/rheometer.cfg

+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

apps/cpu/rheometer/rheometer.cpp

+#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;
+}

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>

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:

src/cpu/VirtualFluidsCore/BoundaryConditions/PowellEyringModelNoSlipBCAlgorithm.h

+#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

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);

src/cpu/VirtualFluidsCore/LBM/PowellEyringModelLBMKernel.h

+#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

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

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