From 440e465229e72e53c0b8759e00c26648f2b4f60b Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Tue, 2 Feb 2021 11:35:21 +0100
Subject: [PATCH] fix RheologyK17LBMKernel
---
.gitignore | 1 +
apps/cpu/HerschelBulkleyModel/hbf1.cfg | 17 ++++---
apps/cpu/HerschelBulkleyModel/hbflow.cpp | 44 +++++++++++++------
apps/cpu/rheometer/rheometer.cfg | 10 ++---
apps/cpu/rheometer/rheometer.cpp | 2 +-
.../LBM/RheologyK17LBMKernel.cpp | 4 +-
src/cpu/VirtualFluidsCore/LBM/Thixotropy.h | 2 +-
7 files changed, 52 insertions(+), 28 deletions(-)
diff --git a/.gitignore b/.gitignore
index a5b8d0a8b..4879d82bb 100644
--- a/.gitignore
+++ b/.gitignore
@@ -1,4 +1,5 @@
# build directories
+buildDev/
build/
buildDev/
bin/
diff --git a/apps/cpu/HerschelBulkleyModel/hbf1.cfg b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
index b5f8ef94d..2eed172d5 100644
--- a/apps/cpu/HerschelBulkleyModel/hbf1.cfg
+++ b/apps/cpu/HerschelBulkleyModel/hbf1.cfg
@@ -1,4 +1,6 @@
-pathname = d:/temp/BinghamNewConvergenceAcousticScaling_it1/128
+pathname = d:/temp/BinghamNewNewConvergenceAcousticScaling_20it/256
+
+#pathname = d:/temp/Newton45/32
numOfThreads = 4
availMem = 8e9
@@ -18,15 +20,18 @@ blocknx = 2 32 2
boundingBox = 2 32 2
deltax = 1
-#nuLB = 0.05
-forcing = 8e-6
+scaleFactor = 8
+
+nuLB = 0.005
+forcing = 6e-9
+tau0 = 3e-8
+
#n = 0.8
-#tau0 = 3e-6
velocity = 1e-3
n = 0.4
Re = 1
Bn = 0.01
-outTime = 10000
-endTime = 300000
\ No newline at end of file
+outTime = 1000
+endTime = 1000000
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleyModel/hbflow.cpp b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
index e0a24aba9..d83f0f6d0 100644
--- a/apps/cpu/HerschelBulkleyModel/hbflow.cpp
+++ b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
@@ -17,7 +17,7 @@ void bflow(string configname)
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 nuLB = config.getValue<double>("nuLB");
double endTime = config.getValue<double>("endTime");
double outTime = config.getValue<double>("outTime");
double availMem = config.getValue<double>("availMem");
@@ -31,11 +31,12 @@ void bflow(string configname)
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 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");
+ double scaleFactor = config.getValue<double>("scaleFactor");
SPtr<Communicator> comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
@@ -63,13 +64,13 @@ void bflow(string configname)
SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
//bounding box
- //double g_minX1 = 0;
- //double g_minX2 = 0;
- //double g_minX3 = 0;
+ //double g_minX1 = -0.707107+1.0;
+ //double g_minX2 = -1.41421+1.0;
+ //double g_minX3 = 1.0;
//double g_maxX1 = boundingBox[0];
//double g_maxX2 = boundingBox[1];
- //double g_maxX3 = boundingBox[2];
+ //double g_maxX3 = boundingBox[2]+1.0;
double g_minX1 = 0.0;
double g_minX2 = -boundingBox[1]/2.0;
@@ -79,6 +80,8 @@ void bflow(string configname)
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;
@@ -93,7 +96,7 @@ void bflow(string configname)
double U = velocity;
double Gamma = U / d;
- double scaleFactor = 4.0;
+ //double scaleFactor = 2.0;
// Diffusive Scaling
@@ -105,13 +108,13 @@ void bflow(string configname)
// Acoustic Scaling
- double k = 0.005 * scaleFactor;
- double tau0 = 3e-5;
+ double k = nuLB * scaleFactor;
+ //double tau0 = 3e-5;
forcing /= scaleFactor;
endTime *= scaleFactor;
deltax /= scaleFactor;
- outTime = endTime;
+ //outTime = endTime;
double beta = 14;
double c = 10; // 1.0 / 6.0;
@@ -130,11 +133,11 @@ void bflow(string configname)
SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
//noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new HerschelBulkleyModelNoSlipBCAlgorithm()));
//noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new PowellEyringModelNoSlipBCAlgorithm()));
- noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelNoSlipBCAlgorithm()));
+ //noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new BinghamModelNoSlipBCAlgorithm()));
//BS visitor
BoundaryConditionsBlockVisitor bcVisitor;
- bcVisitor.addBC(noSlipBCAdapter);
+ //bcVisitor.addBC(noSlipBCAdapter);
SPtr<BCProcessor> bcProc;
bcProc = SPtr<BCProcessor>(new BCProcessor());
@@ -142,6 +145,11 @@ void bflow(string configname)
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new HerschelBulkleyModelLBMKernel());
SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new RheologyK17LBMKernel());
//SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new BinghamModelLBMKernel());
+ //SPtr<LBMKernel> kernel = SPtr<LBMKernel>(new CompressibleCumulant4thOrderViscosityLBMKernel());
+
+ //double forcingXY = forcing / sqrt(2.0);
+ //kernel->setForcingX1(forcingXY);
+ //kernel->setForcingX2(forcingXY);
kernel->setForcingX1(forcing);
kernel->setWithForcing(true);
kernel->setBCProcessor(bcProc);
@@ -192,6 +200,13 @@ void bflow(string configname)
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());
+ //GbCuboid3DPtr addWallXY(new GbCuboid3D(g_minX1 - 0.5*deltax, g_minX2 - 0.5 * deltax, g_minX3 - 0.5 * deltax, g_maxX1 + 0.5 * deltax, g_maxX2 + 0.5 * deltax, g_maxX3 + 0.5 * deltax));
+ //if (myid == 0) GbSystem3D::writeGeoObject(addWallXY.get(), pathname + "/geo/addWallXY", WbWriterVtkXmlASCII::getInstance());
+
+ //SPtr<GbTriFaceMesh3D> wall45(new GbTriFaceMesh3D());
+ //wall45->readMeshFromSTLFile("d:/Projects/TRR277/Project/WP1/PF45/Geo/wall45_ASCII.stl", true);
+
+
//wall interactors
//SPtr<D3Q27Interactor> addWallZminInt(new D3Q27Interactor(addWallZmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
//SPtr<D3Q27Interactor> addWallZmaxInt(new D3Q27Interactor(addWallZmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
@@ -199,9 +214,11 @@ void bflow(string configname)
SPtr<D3Q27Interactor> addWallYminInt(new D3Q27Interactor(addWallYmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
SPtr<D3Q27Interactor> addWallYmaxInt(new D3Q27Interactor(addWallYmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ //SPtr<D3Q27TriFaceMeshInteractor> wall45Int(new D3Q27TriFaceMeshInteractor(wall45, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
////////////////////////////////////////////
//METIS
- SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::KWAY));
+ SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::RECURSIVE));
////////////////////////////////////////////
/////delete solid blocks
if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - start");
@@ -210,6 +227,7 @@ void bflow(string configname)
//intHelper.addInteractor(addWallZmaxInt);
intHelper.addInteractor(addWallYminInt);
intHelper.addInteractor(addWallYmaxInt);
+ //intHelper.addInteractor(wall45Int);
intHelper.selectBlocks();
if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
//////////////////////////////////////
diff --git a/apps/cpu/rheometer/rheometer.cfg b/apps/cpu/rheometer/rheometer.cfg
index c65d95096..f5d9fef26 100644
--- a/apps/cpu/rheometer/rheometer.cfg
+++ b/apps/cpu/rheometer/rheometer.cfg
@@ -1,4 +1,4 @@
-outputPath = d:/temp/rheometer/rheometerBinghamqQBB/rheometerBingham_nu_1.5e-3_newBC_K17
+outputPath = d:/temp/rheometer/rheometerBinghamqQBB/rheometerBingham_tau_20e-7_nu_1.5e-3_new_lim
viscosityPath = d:/Projects/VirtualFluidsCombined/apps/cpu/rheometer
@@ -15,11 +15,11 @@ deltax = 1
refineLevel = 0
-OmegaLB = 1e-4
-tau0 = 5e-7
+OmegaLB = 4e-5
+tau0 = 20e-7
-resolution = 128
-scaleFactor = 4
+resolution = 32
+scaleFactor = 1
newStart = true
restartStep = 100000
diff --git a/apps/cpu/rheometer/rheometer.cpp b/apps/cpu/rheometer/rheometer.cpp
index 47851a944..96c50e5c3 100644
--- a/apps/cpu/rheometer/rheometer.cpp
+++ b/apps/cpu/rheometer/rheometer.cpp
@@ -67,7 +67,7 @@ void bflow(string configname)
OmegaLB /= scaleFactor;
nuLB *=scaleFactor;
endTime *= scaleFactor;
- outTime = endTime;
+ //outTime = endTime;
cpStart = endTime;
cpStep = endTime;
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
index 9df968e12..42c324007 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
@@ -616,7 +616,7 @@ 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) / (drho + c1);
+ LBMReal shearRate = sqrt(c2 * (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);
@@ -639,7 +639,7 @@ void RheologyK17LBMKernel::calculate(int step)
if(omega < c1) { omega = c1; } //arbitrary limit (24.09.2020)
- omega = collFactor;
+ //omega = collFactor;
//magic parameter for rheology
LBMReal a = 10;
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
index 0b2aee946..a5ad23dab 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
@@ -71,7 +71,7 @@ inline LBMReal Thixotropy::getBinghamCollFactor(LBMReal omegaInf, LBMReal shearR
LBMReal omega = omegaInf*cs2 * shearRate * rho / (cs2 * shearRate * rho + omegaInf * tau0);
LBMReal shearRateNew = shearRate * (omega / omegaInf);
- for (int i = 0; i < 3; i++)
+ for (int i = 0; i < 20; i++)
{
omega = omegaInf * cs2 * shearRateNew * rho / (cs2 * shearRateNew * rho + omegaInf * tau0);
shearRateNew = shearRate * (omega / omegaInf);
--
GitLab