diff --git a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
index 4f7951a75174ad7b78d6d3621652b16a58ac6628..28cf8de545ccb0d4a36a18a1046f600180847197 100644
--- a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
+++ b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
@@ -101,6 +101,7 @@ void bflow(string configname)
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()));
@@ -117,9 +118,12 @@ void bflow(string configname)
//fct.DefineConst("H", H);
SPtr<BCAdapter> velocityBCAdapter(new VelocityBCAdapter(true, false, false, fct, 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;
@@ -131,6 +135,8 @@ 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 CompressibleCumulant4thOrderViscosityLBMKernel());
SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
@@ -140,8 +146,8 @@ void bflow(string configname)
SPtr<Grid3D> grid(new Grid3D(comm));
grid->setPeriodicX1(false);
- grid->setPeriodicX2(false);
- grid->setPeriodicX3(false);
+ grid->setPeriodicX2(true);
+ grid->setPeriodicX3(true);
grid->setDeltaX(deltax);
grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
@@ -221,7 +227,7 @@ void bflow(string configname)
//wall interactors
SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, slipBCAdapter, Interactor3D::SOLID));
SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, slipBCAdapter, Interactor3D::SOLID));
-
+
SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, slipBCAdapter, Interactor3D::SOLID));
SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, slipBCAdapter, Interactor3D::SOLID));
@@ -237,10 +243,10 @@ void bflow(string configname)
InteractorsHelper intHelper(grid, metisVisitor);
intHelper.addInteractor(wallXminInt);
intHelper.addInteractor(wallXmaxInt);
- intHelper.addInteractor(wallZminInt);
- intHelper.addInteractor(wallZmaxInt);
intHelper.addInteractor(wallYminInt);
intHelper.addInteractor(wallYmaxInt);
+ intHelper.addInteractor(wallZminInt);
+ intHelper.addInteractor(wallZmaxInt);
intHelper.addInteractor(sphereInt);
intHelper.selectBlocks();
if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.cpp
index 55d9d2c005507ad021e4bc670ae10222c40b0c9f..380167005abc92d0137fd51fae81cf5fe6df8713 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.cpp
@@ -48,8 +48,6 @@ void SimpleSlipBCAlgorithm::applyBC()
{
//quadratic bounce back
const int invDir = D3Q27System::INVDIR[fdir];
- LBMReal q = bcPtr->getQ(invDir);// m+m q=0 stabiler
- //vx3=0;
LBMReal velocity = 0.0;
switch (invDir)
{
diff --git a/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp b/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp
index 362225638cf12f73c397493c9e38e839d681b754..b94840931d22d475b10e52cf22d6949c152f68e7 100644
--- a/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp
+++ b/src/cpu/VirtualFluidsCore/Grid/BasicCalculator.cpp
@@ -33,12 +33,12 @@
#include "BasicCalculator.h"
-#include "BCProcessor.h"
#include "Block3D.h"
-#include "Block3DConnector.h"
+#include "BCProcessor.h"
#include "LBMKernel.h"
-#include "UbLogger.h"
+#include "Block3DConnector.h"
#include "UbScheduler.h"
+#include "UbLogger.h"
#ifdef _OPENMP
#include <omp.h>
@@ -48,126 +48,141 @@
//#define TIMING
//#include "UbTiming.h"
-BasicCalculator::BasicCalculator(SPtr<Grid3D> grid, SPtr<UbScheduler> additionalGhostLayerUpdateScheduler,
- int numberOfTimeSteps)
- : Calculator(grid, additionalGhostLayerUpdateScheduler, numberOfTimeSteps)
+BasicCalculator::BasicCalculator(SPtr<Grid3D> grid, SPtr<UbScheduler> additionalGhostLayerUpdateScheduler, int numberOfTimeSteps) :
+ Calculator(grid, additionalGhostLayerUpdateScheduler, numberOfTimeSteps)
+{
+
+}
+//////////////////////////////////////////////////////////////////////////
+BasicCalculator::~BasicCalculator()
{
+
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::calculate()
{
- UBLOG(logDEBUG1, "OMPCalculator::calculate() - started");
- int calcStep = 0;
- try {
- int minInitLevel = minLevel;
- int maxInitLevel = maxLevel - minLevel;
- int straightStartLevel = minInitLevel;
- int internalIterations = 1 << (maxInitLevel - minInitLevel);
- // int forwardStartLevel;
- int threshold;
+ UBLOG(logDEBUG1, "OMPCalculator::calculate() - started");
+ int calcStep = 0;
+ try
+ {
+ int minInitLevel = minLevel;
+ int maxInitLevel = maxLevel - minLevel;
+ int straightStartLevel = minInitLevel;
+ int internalIterations = 1 << (maxInitLevel - minInitLevel);
+ int forwardStartLevel;
+ int threshold;
#ifdef TIMING
- UbTimer timer;
- double time[6];
+ UbTimer timer;
+ double time[6];
#endif
- for (calcStep = startTimeStep; calcStep <= numberOfTimeSteps; calcStep++) {
- //////////////////////////////////////////////////////////////////////////
+ for (calcStep = startTimeStep; calcStep <= numberOfTimeSteps; calcStep++)
+ {
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- UBLOG(logINFO, "calcStep = " << calcStep);
+ UBLOG(logINFO, "calcStep = " << calcStep);
#endif
- //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
- for (int staggeredStep = 1; staggeredStep <= internalIterations; staggeredStep++) {
- // forwardStartLevel = straightStartLevel;
- if (staggeredStep == internalIterations)
- straightStartLevel = minInitLevel;
- else {
- for (straightStartLevel = maxInitLevel, threshold = 1; (staggeredStep & threshold) != threshold;
- straightStartLevel--, threshold <<= 1)
- ;
- }
- //////////////////////////////////////////////////////////////////////////
+ for (int staggeredStep = 1; staggeredStep <= internalIterations; staggeredStep++)
+ {
+ forwardStartLevel = straightStartLevel;
+ if (staggeredStep == internalIterations) straightStartLevel = minInitLevel;
+ else
+ {
+ for (straightStartLevel = maxInitLevel, threshold = 1;
+ (staggeredStep & threshold) != threshold; straightStartLevel--, threshold <<= 1);
+ }
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- timer.resetAndStart();
+ timer.resetAndStart();
#endif
- //////////////////////////////////////////////////////////////////////////
- applyPreCollisionBC(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
+ applyPreCollisionBC(straightStartLevel, maxInitLevel);
- // do collision for all blocks
- calculateBlocks(straightStartLevel, maxInitLevel, calcStep);
- //////////////////////////////////////////////////////////////////////////
+ //do collision for all blocks
+ calculateBlocks(straightStartLevel, maxInitLevel, calcStep);
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- time[0] = timer.stop();
- UBLOG(logINFO, "calculateBlocks time = " << time[0]);
+ time[0] = timer.stop();
+ UBLOG(logINFO, "calculateBlocks time = " << time[0]);
#endif
- //////////////////////////////////////////////////////////////////////////
- //////////////////////////////////////////////////////////////////////////
- // exchange data between blocks
- exchangeBlockData(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //exchange data between blocks
+ exchangeBlockData(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- time[1] = timer.stop();
- UBLOG(logINFO, "exchangeBlockData time = " << time[1]);
+ time[1] = timer.stop();
+ UBLOG(logINFO, "exchangeBlockData time = " << time[1]);
#endif
- //////////////////////////////////////////////////////////////////////////
- applyPostCollisionBC(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ applyPostCollisionBC(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- time[2] = timer.stop();
- UBLOG(logINFO, "applyBCs time = " << time[2]);
+ time[2] = timer.stop();
+ UBLOG(logINFO, "applyBCs time = " << time[2]);
#endif
- //////////////////////////////////////////////////////////////////////////
- // swap distributions in kernel
- swapDistributions(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //swap distributions in kernel
+ swapDistributions(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- time[3] = timer.stop();
- UBLOG(logINFO, "swapDistributions time = " << time[3]);
+ time[3] = timer.stop();
+ UBLOG(logINFO, "swapDistributions time = " << time[3]);
#endif
- //////////////////////////////////////////////////////////////////////////
- if (refinement) {
- if (straightStartLevel < maxInitLevel)
- exchangeBlockData(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ if (refinement)
+ {
+ if (straightStartLevel < maxInitLevel)
+ exchangeBlockData(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- time[4] = timer.stop();
- UBLOG(logINFO, "refinement exchangeBlockData time = " << time[4]);
+ time[4] = timer.stop();
+ UBLOG(logINFO, "refinement exchangeBlockData time = " << time[4]);
#endif
- //////////////////////////////////////////////////////////////////////////
- // now ghost nodes have actual values
- // interpolation of interface nodes between grid levels
- interpolation(straightStartLevel, maxInitLevel);
- //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //now ghost nodes have actual values
+ //interpolation of interface nodes between grid levels
+ interpolation(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
#ifdef TIMING
- time[5] = timer.stop();
- UBLOG(logINFO, "refinement interpolation time = " << time[5]);
+ time[5] = timer.stop();
+ UBLOG(logINFO, "refinement interpolation time = " << time[5]);
#endif
- //////////////////////////////////////////////////////////////////////////
- }
- }
- // exchange data between blocks for visualization
- if (additionalGhostLayerUpdateScheduler->isDue(calcStep)) {
- exchangeBlockData(straightStartLevel, maxInitLevel);
+ //////////////////////////////////////////////////////////////////////////
}
- coProcess((double)(calcStep));
- // now ghost nodes have actual values
- }
- UBLOG(logDEBUG1, "OMPCalculator::calculate() - stoped");
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- UBLOG(logERROR, " step = " << calcStep);
- // throw e;
- // exit(EXIT_FAILURE);
- } catch (std::string &s) {
- UBLOG(logERROR, s);
- // exit(EXIT_FAILURE);
- // throw s;
- } catch (...) {
- UBLOG(logERROR, "unknown exception");
- // exit(EXIT_FAILURE);
- // throw;
- }
+ }
+ //exchange data between blocks for visualization
+ if (additionalGhostLayerUpdateScheduler->isDue(calcStep))
+ {
+ exchangeBlockData(straightStartLevel, maxInitLevel);
+ }
+ coProcess((double)(calcStep));
+ //now ghost nodes have actual values
+ }
+ UBLOG(logDEBUG1, "OMPCalculator::calculate() - stoped");
+ }
+ catch (std::exception & e)
+ {
+ UBLOG(logERROR, e.what());
+ UBLOG(logERROR, " step = " << calcStep);
+ //throw e;
+ //exit(EXIT_FAILURE);
+ }
+ catch (std::string & s)
+ {
+ UBLOG(logERROR, s);
+ //exit(EXIT_FAILURE);
+ //throw s;
+ }
+ catch (...)
+ {
+ UBLOG(logERROR, "unknown exception");
+ //exit(EXIT_FAILURE);
+ //throw;
+ }
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calcStep)
@@ -175,45 +190,50 @@ void BasicCalculator::calculateBlocks(int startLevel, int maxInitLevel, int calc
#ifdef _OPENMP
#pragma omp parallel
#endif
- {
- SPtr<Block3D> blockTemp;
- // startLevel bis maxInitLevel
- for (int level = startLevel; level <= maxInitLevel; level++) {
- // timer.resetAndStart();
- // call LBM kernel
- int size = (int)blocks[level].size();
+ {
+ SPtr<Block3D> blockTemp;
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level <= maxInitLevel; level++)
+ {
+ //timer.resetAndStart();
+ //call LBM kernel
+ int size = (int)blocks[level].size();
#ifdef _OPENMP
#pragma omp for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
- try {
- blockTemp = blocks[level][i];
- blockTemp->getKernel()->calculate(calcStep);
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- UBLOG(logERROR, blockTemp->toString() << " step = " << calcStep);
- std::exit(EXIT_FAILURE);
- }
+ for (int i = 0; i < size; i++)
+ {
+ try
+ {
+ blockTemp = blocks[level][i];
+ blockTemp->getKernel()->calculate(calcStep);
}
- // timer.stop();
- // UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " <<
- // timer.getTotalTime());
- }
- }
+ catch (std::exception & e)
+ {
+ UBLOG(logERROR, e.what());
+ UBLOG(logERROR, blockTemp->toString() << " step = " << calcStep);
+ std::exit(EXIT_FAILURE);
+ }
+ }
+ //timer.stop();
+ //UBLOG(logINFO, "level = " << level << " blocks = " << blocks[level].size() << " collision time = " << timer.getTotalTime());
+ }
+ }
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::exchangeBlockData(int startLevel, int maxInitLevel)
{
- // startLevel bis maxInitLevel
- for (int level = startLevel; level <= maxInitLevel; level++) {
- // connectorsPrepareLocal(localConns[level]);
- connectorsSendLocal(localConns[level]);
- // connectorsReceiveLocal(localConns[level]);
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level <= maxInitLevel; level++)
+ {
+ //connectorsPrepareLocal(localConns[level]);
+ connectorsSendLocal(localConns[level]);
+ //connectorsReceiveLocal(localConns[level]);
- connectorsPrepareRemote(remoteConns[level]);
- connectorsSendRemote(remoteConns[level]);
- connectorsReceiveRemote(remoteConns[level]);
- }
+ connectorsPrepareRemote(remoteConns[level]);
+ connectorsSendRemote(remoteConns[level]);
+ connectorsReceiveRemote(remoteConns[level]);
+ }
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::swapDistributions(int startLevel, int maxInitLevel)
@@ -221,156 +241,192 @@ void BasicCalculator::swapDistributions(int startLevel, int maxInitLevel)
#ifdef _OPENMP
#pragma omp parallel
#endif
- {
- // startLevel bis maxInitLevel
- for (int level = startLevel; level <= maxInitLevel; level++) {
- int size = (int)blocks[level].size();
+ {
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level <= maxInitLevel; level++)
+ {
+ int size = (int)blocks[level].size();
#ifdef _OPENMP
#pragma omp for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
- blocks[level][i]->getKernel()->swapDistributions();
- }
- }
- }
+ for (int i = 0; i < size; i++)
+ {
+ blocks[level][i]->getKernel()->swapDistributions();
+ }
+ }
+ }
}
//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsPrepareLocal(std::vector<SPtr<Block3DConnector>> &connectors)
+void BasicCalculator::connectorsPrepareLocal(std::vector< SPtr<Block3DConnector> >& connectors)
{
- int size = (int)connectors.size();
+ int size = (int)connectors.size();
#ifdef _OPENMP
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
- try {
- connectors[i]->prepareForReceive();
- connectors[i]->prepareForSend();
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- std::exit(EXIT_FAILURE);
- }
- }
+ for (int i = 0; i < size; i++)
+ {
+ try
+ {
+ connectors[i]->prepareForReceive();
+ connectors[i]->prepareForSend();
+ }
+ catch (std::exception & e)
+ {
+ UBLOG(logERROR, e.what());
+ std::exit(EXIT_FAILURE);
+ }
+ }
}
//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsSendLocal(std::vector<SPtr<Block3DConnector>> &connectors)
+void BasicCalculator::connectorsSendLocal(std::vector< SPtr<Block3DConnector> >& connectors)
{
- int size = (int)connectors.size();
+ int size = (int)connectors.size();
#ifdef _OPENMP
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
- try {
- connectors[i]->fillSendVectors();
- connectors[i]->sendVectors();
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- std::exit(EXIT_FAILURE);
- }
- }
+ for (int i = 0; i < size; i++)
+ {
+ try
+ {
+ connectors[i]->fillSendVectors();
+ connectors[i]->sendVectors();
+ }
+ catch (std::exception & e)
+ {
+ UBLOG(logERROR, e.what());
+ std::exit(EXIT_FAILURE);
+ }
+ }
}
//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsReceiveLocal(std::vector<SPtr<Block3DConnector>> &connectors)
+void BasicCalculator::connectorsReceiveLocal(std::vector< SPtr<Block3DConnector> >& connectors)
{
- int size = (int)connectors.size();
+ int size = (int)connectors.size();
#ifdef _OPENMP
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
- connectors[i]->receiveVectors();
- connectors[i]->distributeReceiveVectors();
- }
+ for (int i = 0; i < size; i++)
+ {
+ connectors[i]->receiveVectors();
+ connectors[i]->distributeReceiveVectors();
+ }
}
-void BasicCalculator::connectorsPrepareRemote(std::vector<SPtr<Block3DConnector>> &connectors)
+void BasicCalculator::connectorsPrepareRemote(std::vector< SPtr<Block3DConnector> >& connectors)
{
- int size = (int)connectors.size();
- for (int i = 0; i < size; i++) {
- connectors[i]->prepareForReceive();
- connectors[i]->prepareForSend();
- }
+ int size = (int)connectors.size();
+ for (int i = 0; i < size; i++)
+ {
+ connectors[i]->prepareForReceive();
+ connectors[i]->prepareForSend();
+ }
}
//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsSendRemote(std::vector<SPtr<Block3DConnector>> &connectors)
+void BasicCalculator::connectorsSendRemote(std::vector< SPtr<Block3DConnector> >& connectors)
{
- int size = (int)connectors.size();
- for (int i = 0; i < size; i++) {
- connectors[i]->fillSendVectors();
- connectors[i]->sendVectors();
- }
+ int size = (int)connectors.size();
+ for (int i = 0; i < size; i++)
+ {
+ connectors[i]->fillSendVectors();
+ connectors[i]->sendVectors();
+ }
}
//////////////////////////////////////////////////////////////////////////
-void BasicCalculator::connectorsReceiveRemote(std::vector<SPtr<Block3DConnector>> &connectors)
+void BasicCalculator::connectorsReceiveRemote(std::vector< SPtr<Block3DConnector> >& connectors)
{
- int size = (int)connectors.size();
- for (int i = 0; i < size; i++) {
- connectors[i]->receiveVectors();
- connectors[i]->distributeReceiveVectors();
- }
+ int size = (int)connectors.size();
+ for (int i = 0; i < size; i++)
+ {
+ connectors[i]->receiveVectors();
+ connectors[i]->distributeReceiveVectors();
+ }
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::interpolation(int startLevel, int maxInitLevel)
{
- for (int level = startLevel; level < maxInitLevel; level++) {
- connectorsPrepareLocal(localInterConns[level]);
- connectorsPrepareRemote(remoteInterConns[level]);
- }
+ for (int level = startLevel; level < maxInitLevel; level++)
+ {
+ connectorsPrepareLocal(localInterConns[level]);
+ connectorsPrepareRemote(remoteInterConns[level]);
+ }
- for (int level = startLevel; level < maxInitLevel; level++) {
- connectorsSendLocal(localInterConns[level]);
- connectorsSendRemote(remoteInterConns[level]);
- }
+ for (int level = startLevel; level < maxInitLevel; level++)
+ {
+ connectorsSendLocal(localInterConns[level]);
+ connectorsSendRemote(remoteInterConns[level]);
+ }
- for (int level = startLevel; level < maxInitLevel; level++) {
- connectorsReceiveLocal(localInterConns[level]);
- connectorsReceiveRemote(remoteInterConns[level]);
- }
+ for (int level = startLevel; level < maxInitLevel; level++)
+ {
+ connectorsReceiveLocal(localInterConns[level]);
+ connectorsReceiveRemote(remoteInterConns[level]);
+ }
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::applyPreCollisionBC(int startLevel, int maxInitLevel)
{
- // startLevel bis maxInitLevel
- for (int level = startLevel; level <= maxInitLevel; level++) {
- int size = (int)blocks[level].size();
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level <= maxInitLevel; level++)
+ {
+ int size = (int)blocks[level].size();
#ifdef _OPENMP
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
- try {
- blocks[level][i]->getKernel()->getBCProcessor()->applyPreCollisionBC();
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- exit(EXIT_FAILURE);
- } catch (std::string &s) {
- UBLOG(logERROR, s);
- exit(EXIT_FAILURE);
- } catch (...) {
- UBLOG(logERROR, "unknown exception");
- exit(EXIT_FAILURE);
- }
- }
- }
+ for (int i = 0; i < size; i++)
+ {
+ try
+ {
+ blocks[level][i]->getKernel()->getBCProcessor()->applyPreCollisionBC();
+ }
+ catch (std::exception & e)
+ {
+ UBLOG(logERROR, e.what());
+ exit(EXIT_FAILURE);
+ }
+ catch (std::string & s)
+ {
+ UBLOG(logERROR, s);
+ exit(EXIT_FAILURE);
+ }
+ catch (...)
+ {
+ UBLOG(logERROR, "unknown exception");
+ exit(EXIT_FAILURE);
+ }
+ }
+ }
}
//////////////////////////////////////////////////////////////////////////
void BasicCalculator::applyPostCollisionBC(int startLevel, int maxInitLevel)
{
- // startLevel bis maxInitLevel
- for (int level = startLevel; level <= maxInitLevel; level++) {
- int size = (int)blocks[level].size();
+ //startLevel bis maxInitLevel
+ for (int level = startLevel; level <= maxInitLevel; level++)
+ {
+ int size = (int)blocks[level].size();
#ifdef _OPENMP
#pragma omp parallel for schedule(OMP_SCHEDULE)
#endif
- for (int i = 0; i < size; i++) {
- try {
- blocks[level][i]->getKernel()->getBCProcessor()->applyPostCollisionBC();
- } catch (std::exception &e) {
- UBLOG(logERROR, e.what());
- exit(EXIT_FAILURE);
- } catch (std::string &s) {
- UBLOG(logERROR, s);
- exit(EXIT_FAILURE);
- } catch (...) {
- UBLOG(logERROR, "unknown exception");
- exit(EXIT_FAILURE);
- }
- }
- }
+ for (int i = 0; i < size; i++)
+ {
+ try
+ {
+ blocks[level][i]->getKernel()->getBCProcessor()->applyPostCollisionBC();
+ }
+ catch (std::exception & e)
+ {
+ UBLOG(logERROR, e.what());
+ exit(EXIT_FAILURE);
+ }
+ catch (std::string & s)
+ {
+ UBLOG(logERROR, s);
+ exit(EXIT_FAILURE);
+ }
+ catch (...)
+ {
+ UBLOG(logERROR, "unknown exception");
+ exit(EXIT_FAILURE);
+ }
+ }
+ }
}
+
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
index 9287b99c5be35a1ea608c672e77a733242c4b550..43713cbc676e2aacb6775bd8912ff0c1b59cfb77 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
@@ -107,24 +107,25 @@ void RheologyK17LBMKernel::calculate(int step)
int maxX2 = bcArrayMaxX2-ghostLayerWidth;
int maxX3 = bcArrayMaxX3-ghostLayerWidth;
- LBMReal omega = collFactor;
+ //LBMReal omega = collFactor;
- //magic parameter for rheology
- LBMReal a = 10;
- OxxPyyPzz = c1 / (a * ((c1 / omega) - c1o2) + c1o2);
+ ////magic parameter for rheology
+ //LBMReal a = 10;
+ //OxxPyyPzz = c1 / (a * ((c1 / omega) - c1o2) + c1o2);
+ //OxxPyyPzz = (OxxPyyPzz > c1) ? c1 : OxxPyyPzz;
- //LBMReal OxyyPxzz = eight*(-two+omega)*(one+two*omega)/(-eight-fourteen*omega+seven*omega*omega);//one;
- //LBMReal OxyyMxzz = eight*(-two+omega)*(-seven+four*omega)/(fiftysix-fifty*omega+nine*omega*omega);//one;
- //LBMReal Oxyz = twentyfour*(-two+omega)*(-two-seven*omega+three*omega*omega)/(fourtyeight+c152*omega-c130*omega*omega+twentynine*omega*omega*omega);
- LBMReal OxyyPxzz = 8.0*(omega-2.0)*(OxxPyyPzz*(3.0*omega-1.0)-5.0*omega)/(8.0*(5.0-2.0*omega)*omega+OxxPyyPzz*(8.0+omega*(9.0*omega-26.0)));
- LBMReal OxyyMxzz = 8.0*(omega-2.0)*(omega+OxxPyyPzz*(3.0*omega-7.0))/(OxxPyyPzz*(56.0-42.0*omega+9.0*omega*omega)-8.0*omega);
- LBMReal Oxyz = 24.0*(omega-2.0)*(4.0*omega*omega+omega*OxxPyyPzz*(18.0-13.0*omega)+OxxPyyPzz*OxxPyyPzz*(2.0+omega*(6.0*omega-11.0)))/(16.0*omega*omega*(omega-6.0)-2.0*omega*OxxPyyPzz*(216.0+5.0*omega*(9.0*omega-46.0))+OxxPyyPzz*OxxPyyPzz*(omega*(3.0*omega-10.0)*(15.0*omega-28.0)-48.0));
+ ////LBMReal OxyyPxzz = eight*(-two+omega)*(one+two*omega)/(-eight-fourteen*omega+seven*omega*omega);//one;
+ ////LBMReal OxyyMxzz = eight*(-two+omega)*(-seven+four*omega)/(fiftysix-fifty*omega+nine*omega*omega);//one;
+ ////LBMReal Oxyz = twentyfour*(-two+omega)*(-two-seven*omega+three*omega*omega)/(fourtyeight+c152*omega-c130*omega*omega+twentynine*omega*omega*omega);
+ //LBMReal OxyyPxzz = 8.0*(omega-2.0)*(OxxPyyPzz*(3.0*omega-1.0)-5.0*omega)/(8.0*(5.0-2.0*omega)*omega+OxxPyyPzz*(8.0+omega*(9.0*omega-26.0)));
+ //LBMReal OxyyMxzz = 8.0*(omega-2.0)*(omega+OxxPyyPzz*(3.0*omega-7.0))/(OxxPyyPzz*(56.0-42.0*omega+9.0*omega*omega)-8.0*omega);
+ //LBMReal Oxyz = 24.0*(omega-2.0)*(4.0*omega*omega+omega*OxxPyyPzz*(18.0-13.0*omega)+OxxPyyPzz*OxxPyyPzz*(2.0+omega*(6.0*omega-11.0)))/(16.0*omega*omega*(omega-6.0)-2.0*omega*OxxPyyPzz*(216.0+5.0*omega*(9.0*omega-46.0))+OxxPyyPzz*OxxPyyPzz*(omega*(3.0*omega-10.0)*(15.0*omega-28.0)-48.0));
- //LBMReal A = (four + two*omega - three*omega*omega) / (two - seven*omega + five*omega*omega);
- //LBMReal B = (four + twentyeight*omega - fourteen*omega*omega) / (six - twentyone*omega + fiveteen*omega*omega);
+ ////LBMReal A = (four + two*omega - three*omega*omega) / (two - seven*omega + five*omega*omega);
+ ////LBMReal B = (four + twentyeight*omega - fourteen*omega*omega) / (six - twentyone*omega + fiveteen*omega*omega);
- LBMReal A = (4.0*omega*omega+2.0*omega*OxxPyyPzz*(omega-6.0)+OxxPyyPzz*OxxPyyPzz*(omega*(10.0-3.0*omega)-4.0))/((omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
- LBMReal B = (4.0*omega*OxxPyyPzz*(9.0*omega-16.0)-4.0*omega*omega-2.0*OxxPyyPzz*OxxPyyPzz*(2.0+9.0*omega*(omega-2.0)))/(3.0*(omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
+ //LBMReal A = (4.0*omega*omega+2.0*omega*OxxPyyPzz*(omega-6.0)+OxxPyyPzz*OxxPyyPzz*(omega*(10.0-3.0*omega)-4.0))/((omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
+ //LBMReal B = (4.0*omega*OxxPyyPzz*(9.0*omega-16.0)-4.0*omega*omega-2.0*OxxPyyPzz*OxxPyyPzz*(2.0+9.0*omega*(omega-2.0)))/(3.0*(omega-OxxPyyPzz)*(OxxPyyPzz*(2.0+3.0*omega)-8.0*omega));
for (int x3 = minX3; x3 < maxX3; x3++)
{
@@ -206,7 +207,7 @@ void RheologyK17LBMKernel::calculate(int step)
(mfbbc-mfbba))/rho;
////////////////////////////////////////////////////////////////////////////////////
- LBMReal collFactorF = collFactor;
+ LBMReal omega = collFactor;
//forcing
///////////////////////////////////////////////////////////////////////////////////////////
@@ -605,7 +606,7 @@ void RheologyK17LBMKernel::calculate(int step)
///////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//incl. correction (hat noch nicht so gut funktioniert...Optimierungsbedarf??)
- LBMReal dxux = c1o2 * (-omega) *(mxxMyy+mxxMzz)+c1o2 * OxxPyyPzz * (mfaaa-mxxPyyPzz);
+ LBMReal dxux = c1o2 * (-omega) * (mxxMyy + mxxMzz);// +c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
LBMReal dyuy = dxux+omega * c3o2 * mxxMyy;
LBMReal dzuz = dxux+omega * c3o2 * mxxMzz;
@@ -615,19 +616,42 @@ void RheologyK17LBMKernel::calculate(int step)
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//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);
- collFactorF = Thixotropy::getHerschelBulkleyCollFactor(collFactorF, shearRate, drho);
+ LBMReal shearRate = sqrt(dxux * dxux + dyuy * dyuy + dzuz * dzuz + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (drho + c1);
+ //omega = getThyxotropyCollFactor(omega, shearRate, rho);
+ omega = Thixotropy::getHerschelBulkleyCollFactor(omega, shearRate, drho);
+ //omega = Thixotropy::getBinghamCollFactor(omega, shearRate, drho);
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ mxxMyy += omega * (-mxxMyy) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+ mxxMzz += omega * (-mxxMzz) - three * (one + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
+
+ mfabb += omega * (-mfabb);
+ mfbab += omega * (-mfbab);
+ mfbba += omega * (-mfbba);
+
+ if(omega < c1) { omega = c1; } //arbitrary limit (24.09.2020)
+
+ //magic parameter for rheology
+ LBMReal a = 10;
+ OxxPyyPzz = c1 / (a * ((c1 / omega) - c1o2) + c1o2);
+ OxxPyyPzz = (OxxPyyPzz > c1) ? c1 : OxxPyyPzz;
+
+ LBMReal OxyyPxzz = 8.0 * (omega - 2.0) * (OxxPyyPzz * (3.0 * omega - 1.0) - 5.0 * omega) / (8.0 * (5.0 - 2.0 * omega) * omega + OxxPyyPzz * (8.0 + omega * (9.0 * omega - 26.0)));
+ LBMReal OxyyMxzz = 8.0 * (omega - 2.0) * (omega + OxxPyyPzz * (3.0 * omega - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * omega + 9.0 * omega * omega) - 8.0 * omega);
+ LBMReal Oxyz = 24.0 * (omega - 2.0) * (4.0 * omega * omega + omega * OxxPyyPzz * (18.0 - 13.0 * omega) + OxxPyyPzz * OxxPyyPzz * (2.0 + omega * (6.0 * omega - 11.0))) / (16.0 * omega * omega * (omega - 6.0) - 2.0 * omega * OxxPyyPzz * (216.0 + 5.0 * omega * (9.0 * omega - 46.0)) + OxxPyyPzz * OxxPyyPzz * (omega * (3.0 * omega - 10.0) * (15.0 * omega - 28.0) - 48.0));
+
+ LBMReal A = (4.0 * omega * omega + 2.0 * omega * OxxPyyPzz * (omega - 6.0) + OxxPyyPzz * OxxPyyPzz * (omega * (10.0 - 3.0 * omega) - 4.0)) / ((omega - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * omega) - 8.0 * omega));
+ LBMReal B = (4.0 * omega * OxxPyyPzz * (9.0 * omega - 16.0) - 4.0 * omega * omega - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * omega * (omega - 2.0))) / (3.0 * (omega - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * omega) - 8.0 * omega));
+
//relax
- wadjust = OxxPyyPzz+(one-OxxPyyPzz)*fabs((mfaaa-mxxPyyPzz))/(fabs((mfaaa-mxxPyyPzz))+qudricLimitD);
- mxxPyyPzz += wadjust*(mfaaa-mxxPyyPzz)-three * (one-c1o2 * OxxPyyPzz) * (vx2 * dxux+vy2 * dyuy+vz2 * dzuz);
+ //wadjust = OxxPyyPzz+(one-OxxPyyPzz)*fabs((mfaaa-mxxPyyPzz))/(fabs((mfaaa-mxxPyyPzz))+qudricLimitD);
+ //mxxPyyPzz += wadjust*(mfaaa-mxxPyyPzz)-three * (one-c1o2 * OxxPyyPzz) * (vx2 * dxux+vy2 * dyuy+vz2 * dzuz);
+ mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) - three * (one - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
// mxxPyyPzz += OxxPyyPzz*(mfaaa-mxxPyyPzz)-three * (one-c1o2 * OxxPyyPzz) * (vx2 * dxux+vy2 * dyuy+vz2 * dzuz);//-magicBulk*OxxPyyPzz;
- mxxMyy += omega * (-mxxMyy)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vy2 * dyuy);
- mxxMzz += omega * (-mxxMzz)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vz2 * dzuz);
+ //mxxMyy += omega * (-mxxMyy)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vy2 * dyuy);
+ //mxxMzz += omega * (-mxxMzz)-three * (one+c1o2 * (-omega)) * (vx2 * dxux-vz2 * dzuz);
//////////////////////////////////////////////////////////////////////////
//limiter-Scheise Teil 2
@@ -644,9 +668,9 @@ void RheologyK17LBMKernel::calculate(int step)
//mxxMyy += -(-omega) * (-mxxMyy);
//mxxMzz += -(-omega) * (-mxxMzz);
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
- mfabb += omega * (-mfabb);
- mfbab += omega * (-mfbab);
- mfbba += omega * (-mfbba);
+ //mfabb += omega * (-mfabb);
+ //mfbab += omega * (-mfbab);
+ //mfbba += omega * (-mfbba);
//////////////////////////////////////////////////////////////////////////
//limiter-Scheise Teil 3
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h
index 03015cef8fe5d95725d6aedebc051d8ee4f46025..f2c211dffa5da0ca5e189759c75bdaad5a9a5006 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h
@@ -8,7 +8,7 @@
#include "basics/container/CbArray4D.h"
#include "basics/container/CbArray3D.h"
-//! \brief compressible cumulant LBM kernel.
+//! \brief compressible cumulant LBM kernel with rheological properties of shear and bulk viscosity for non-Newtonian fluids.
//! \details CFD solver that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model
//! \author K. Kutscher, M. Geier
class RheologyK17LBMKernel : public LBMKernel