From ccb8cc57a5bdfba5bb931e9a458c82b0b4b5a064 Mon Sep 17 00:00:00 2001
From: kutscher <kutscher@irmb.tu-bs.de>
Date: Fri, 11 Sep 2020 14:38:33 +0200
Subject: [PATCH] add flow past sphere and flow past cylinder benchmarks
---
apps/cpu/Applications.cmake | 4 +-
apps/cpu/FlowAroundCylinder/CMakeLists.txt | 19 +-
apps/cpu/FlowAroundCylinder/cylinder.cfg | 8 +-
apps/cpu/FlowAroundCylinder/cylinder.cpp | 19 +-
apps/cpu/HerschelBulkleySphere/hbsphere.cfg | 18 +-
apps/cpu/HerschelBulkleySphere/hbsphere.cpp | 52 +++--
apps/cpu/sphere/sphere.cpp | 190 +++++++++----------
src/cpu/VirtualFluidsCore/Data/DataSet3D.h | 33 ++++
src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp | 8 +-
src/cpu/VirtualFluidsCore/LBM/Thixotropy.h | 9 +-
10 files changed, 188 insertions(+), 172 deletions(-)
diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 1e444b79b..2a40b73bd 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -11,8 +11,8 @@ add_subdirectory(${APPS_ROOT_CPU}/sphere)
#add_subdirectory(Applications/bananas2)
# add_subdirectory(Applications/plate)
# add_subdirectory(Applications/plate2)
-##add_subdirectory(Applications/FlowAroundCylinder)
-#add_subdirectory(${APPS_ROOT_CPU}/LaminarTubeFlow)
+add_subdirectory(${APPS_ROOT_CPU}/FlowAroundCylinder)
+add_subdirectory(${APPS_ROOT_CPU}/LaminarTubeFlow)
# add_subdirectory(Applications/LaminarTubeFlowConv)
#add_subdirectory(Applications/cylinderSt)
#add_subdirectory(Applications/mpichTest)
diff --git a/apps/cpu/FlowAroundCylinder/CMakeLists.txt b/apps/cpu/FlowAroundCylinder/CMakeLists.txt
index a4173b62b..54397d6b4 100644
--- a/apps/cpu/FlowAroundCylinder/CMakeLists.txt
+++ b/apps/cpu/FlowAroundCylinder/CMakeLists.txt
@@ -5,21 +5,6 @@ CMAKE_MINIMUM_REQUIRED(VERSION 2.8)
########################################################
PROJECT(cylinder)
-INCLUDE(${APPS_ROOT}/IncludsList.cmake)
+INCLUDE(${APPS_ROOT_CPU}/IncludsList.cmake)
-#################################################################
-### LOCAL FILES ###
-#################################################################
-FILE(GLOB SPECIFIC_FILES ${CMAKE_CURRENT_SOURCE_DIR}/*.h
- ${CMAKE_CURRENT_SOURCE_DIR}/*.cpp
- ${CMAKE_CURRENT_SOURCE_DIR}/*.hpp )
-
-SET(ALL_SOURCES ${ALL_SOURCES} ${SPECIFIC_FILES})
-SOURCE_GROUP(src FILES ${SPECIFIC_FILES})
-
-SET(CAB_ADDITIONAL_LINK_LIBRARIES VirtualFluids)
-
-#################################################################
-### CREATE PROJECT ###
-#################################################################
-CREATE_CAB_PROJECT(cylinder BINARY)
+vf_add_library(BUILDTYPE binary DEPENDS VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} FILES cylinder.cpp )
diff --git a/apps/cpu/FlowAroundCylinder/cylinder.cfg b/apps/cpu/FlowAroundCylinder/cylinder.cfg
index b895c7a49..0a7066ed9 100644
--- a/apps/cpu/FlowAroundCylinder/cylinder.cfg
+++ b/apps/cpu/FlowAroundCylinder/cylinder.cfg
@@ -2,11 +2,11 @@ pathOut = d:/temp/cylinder_test
numOfThreads = 4
availMem = 15e9
-refineLevel = 1
+refineLevel = 0
blockNx = 25 41 41
uLB = 0.001
-dx = 0.005
-#dx = 0.01
+#dx = 0.005
+dx = 0.01
logToFile = false
@@ -16,7 +16,7 @@ restartStep = 1000
cpStart = 1000
cpStep = 1000
-outTime = 100
+outTime = 10000
endTime = 100000
nupsStep = 100 100 10000000
\ No newline at end of file
diff --git a/apps/cpu/FlowAroundCylinder/cylinder.cpp b/apps/cpu/FlowAroundCylinder/cylinder.cpp
index da443cc37..8e0073b77 100644
--- a/apps/cpu/FlowAroundCylinder/cylinder.cpp
+++ b/apps/cpu/FlowAroundCylinder/cylinder.cpp
@@ -259,9 +259,9 @@ void run(string configname)
intHelper.setBC();
- //domain decomposition
- PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
- grid->accept(pqPartVisitor);
+ ////domain decomposition
+ //PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
+ //grid->accept(pqPartVisitor);
//initialization of distributions
InitDistributionsBlockVisitor initVisitor;
@@ -297,10 +297,11 @@ void run(string configname)
SPtr<CoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, stepSch, pathOut, WbWriterVtkXmlBinary::getInstance(), conv, comm));
- //double area = (2.0*radius*H)/(dx*dx);
- //double v = 4.0*uLB/9.0;
- //CalculateForcesCoProcessor fp(grid, stepSch, pathOut + "/results/forces.txt", comm, v, area);
- //fp.addInteractor(cylinderInt);
+ double area = (2.0*radius*H)/(dx*dx);
+ double v = 4.0*uLB/9.0;
+ SPtr<UbScheduler> forceSch(new UbScheduler(100));
+ SPtr<CalculateForcesCoProcessor> fp = make_shared<CalculateForcesCoProcessor>(grid, forceSch, pathOut + "/results/forces.txt", comm, v, area);
+ fp->addInteractor(cylinderInt);
SPtr<UbScheduler> nupsSch(new UbScheduler(nupsStep[0], nupsStep[1], nupsStep[2]));
std::shared_ptr<CoProcessor> nupsCoProcessor(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
@@ -309,7 +310,9 @@ void run(string configname)
SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
calculator->addCoProcessor(nupsCoProcessor);
- calculator->addCoProcessor(writeMQCoProcessor);
+ calculator->addCoProcessor(fp);
+ calculator->addCoProcessor(writeMQCoProcessor);
+
if(myid == 0) UBLOG(logINFO,"Simulation-start");
calculator->calculate();
if(myid == 0) UBLOG(logINFO,"Simulation-end");
diff --git a/apps/cpu/HerschelBulkleySphere/hbsphere.cfg b/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
index ebcc5d1f8..639737953 100644
--- a/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
+++ b/apps/cpu/HerschelBulkleySphere/hbsphere.cfg
@@ -1,18 +1,18 @@
outputPath = d:/temp/Herschel-BulkleySphere
-numOfThreads = 1
+numOfThreads = 4
availMem = 8e9
logToFile = false
-blocknx = 5 5 5
-boundingBox = 40 40 40 #30*20=600**3=216000000
+blocknx = 50 50 50
+boundingBox = 150 50 50 #30*20=600**3=216000000
deltax = 1
-refineLevel = 1
+refineLevel = 0
radius = 5
velocity = 1e-3
-n = 0.3
+n = 0.9
Re = 1
Bn = 0.01
@@ -20,8 +20,8 @@ Bn = 0.01
newStart = true
restartStep = 100000
-cpStart = 10000
-cpStep = 10000
+cpStart = 100000
+cpStep = 100000
-outTime = 1
-endTime = 1000
\ No newline at end of file
+outTime = 10000
+endTime = 1000000
\ No newline at end of file
diff --git a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
index b437a2736..19ffd60b5 100644
--- a/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
+++ b/apps/cpu/HerschelBulkleySphere/hbsphere.cpp
@@ -61,21 +61,21 @@ void bflow(string configname)
//bounding box
- //double g_minX1 = 0;
- //double g_minX2 = 0;
- //double g_minX3 = 0;
+ 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_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_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 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;
@@ -84,13 +84,20 @@ void bflow(string configname)
double Gamma = U / d;
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 HerschelBulkleyModelNoSlipBCAlgorithm()));
@@ -117,6 +124,8 @@ void bflow(string configname)
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);
@@ -129,7 +138,7 @@ void bflow(string configname)
if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), outputPath + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
//sphere
- SPtr<GbObject3D> sphere(new GbSphere3D(0, 0, 0, radius));
+ SPtr<GbObject3D> sphere(new GbSphere3D(75, 25, 25, radius));
GbSystem3D::writeGeoObject(sphere.get(), outputPath + "/geo/sphere", WbWriterVtkXmlBinary::getInstance());
SPtr<D3Q27Interactor> sphereInt(new D3Q27Interactor(sphere, grid, noSlipBCAdapter, Interactor3D::SOLID));
@@ -147,13 +156,13 @@ void bflow(string configname)
//UBLOG(logINFO, "forcing = " << forcing);
UBLOG(logINFO, "rho = " << rhoLB);
UBLOG(logINFO, "U = " << U);
- UBLOG(logINFO, "Re = " << Re);
- UBLOG(logINFO, "Bn = " << Bn);
+ 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 levels = " << refineLevel + 1);
UBLOG(logINFO, "number of threads = " << numOfThreads);
UBLOG(logINFO, "number of processes = " << comm->getNumberOfProcesses());
UBLOG(logINFO, "Preprozess - start");
@@ -219,7 +228,7 @@ void bflow(string configname)
intHelper.addInteractor(wallYmaxInt);
intHelper.addInteractor(wallXminInt);
intHelper.addInteractor(wallXmaxInt);
- //intHelper.addInteractor(sphereInt);
+ intHelper.addInteractor(sphereInt);
intHelper.selectBlocks();
if (myid == 0) UBLOG(logINFO, "deleteSolidBlocks - end");
//////////////////////////////////////
@@ -281,6 +290,7 @@ void bflow(string configname)
//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);
@@ -297,10 +307,10 @@ void bflow(string configname)
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
//SPtr<UbScheduler> visSch(new UbScheduler(10,1));
SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlASCII::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
- writeMQCoProcessor->process(0);
+ //writeMQCoProcessor->process(0);
- double area = 4*UbMath::PI*radius*radius;
- SPtr<UbScheduler> forceSch(new UbScheduler(1000));
+ 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);
diff --git a/apps/cpu/sphere/sphere.cpp b/apps/cpu/sphere/sphere.cpp
index c33a3d696..400e0e28a 100644
--- a/apps/cpu/sphere/sphere.cpp
+++ b/apps/cpu/sphere/sphere.cpp
@@ -9,37 +9,36 @@ void run(string configname)
{
try
{
-
- //Sleep(30000);
-
- string machine = QUOTEME(CAB_MACHINE);
-
SPtr<Communicator> comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
int mybundle = comm->getBundleID();
int root = comm->getRoot();
- ConfigurationFile config;
- config.load(configname);
-
- string pathname = config.getValue<string>("path");
- double availMem = config.getValue<double>("memory");
- string metafile = config.getValue<string>("metafile");
- double outstep = config.getValue<double>("outstep");
- double endstep = config.getValue<double>("endstep");
- int numOfThreads = config.getValue<int>("threads");
- const int refineLevel = config.getValue<int>("level");
-
- bool test = config.getValue<bool>("test");
-
- LBMReal radius = 10;
- LBMReal uLB = 0.01;
+ //ConfigurationFile config;
+ //config.load(configname);
+
+ //string pathname = config.getValue<string>("path");
+ //double availMem = config.getValue<double>("memory");
+ //string metafile = config.getValue<string>("metafile");
+ //double outstep = config.getValue<double>("outstep");
+ //double endstep = config.getValue<double>("endstep");
+ //int numOfThreads = config.getValue<int>("threads");
+ //const int refineLevel = config.getValue<int>("level");
+
+ string outputPath = "d:/temp/sphereParab4";
+ double availMem = 8e9;
+ double outstep = 10000;
+ double endstep = 1e6;
+ int numOfThreads = 4;
+ omp_set_num_threads(numOfThreads);
+ int refineLevel = 0;
+
+ LBMReal radius = 5;
+ LBMReal uLB = 1e-3;
LBMReal Re = 1;
LBMReal rhoLB = 0.0;
- //LBMReal nuLB = (uLB*2.0*radius)/Re;
- //LBMReal nuLB = (uLB*L2)/Re;
- LBMReal nuLB = 0.168666666667/100;
+ LBMReal nuLB = (uLB*2.0*radius)/Re;
double dp_LB = 1e-6;
double rhoLBinflow = dp_LB*3.0;
@@ -47,11 +46,17 @@ void run(string configname)
SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
+ double H = 50;
+ //mu::Parser fct;
+ //fct.SetExpr("U");
+ //fct.DefineConst("U", uLB);
mu::Parser fct;
- fct.SetExpr("U");
+ fct.SetExpr("16*U*x2*x3*(H-x2)*(H-x3)/H^4");
fct.DefineConst("U", uLB);
+ fct.DefineConst("H", H);
SPtr<BCAdapter> velBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
- velBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
+ //velBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
+ velBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new SimpleVelocityBCAlgorithm()));
SPtr<BCAdapter> denBCAdapter(new DensityBCAdapter(rhoLB));
denBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
@@ -61,39 +66,34 @@ void run(string configname)
bcVisitor.addBC(velBCAdapter);
bcVisitor.addBC(denBCAdapter);
- SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
-
double dx = 1;
- const int blocknx1 = 5;
- const int blocknx2 = 5;
- const int blocknx3 = 5;
-
- const int gridNx1 = 8;//18;
- const int gridNx2 = 8;// 11;
- const int gridNx3 = 8;// 11;
-
- //const int blocknx1 = 40;
- //const int blocknx2 = 40;
- //const int blocknx3 = 40;
-
- //const int gridNx1 = 2;
- //const int gridNx2 = 2;
- //const int gridNx3 = 2;
+ const int blocknx1 = 50;
+ const int blocknx2 = 50;
+ const int blocknx3 = 50;
- double L1 = gridNx1*blocknx1;
- double L2, L3;
- L2 = gridNx2*blocknx1;
- L3 = gridNx3*blocknx1;
+ const int gridNx1 = 150;
+ const int gridNx2 = H;
+ const int gridNx3 = H;
+ double L1, L2, L3;
+ L1 = gridNx1;
+ L2 = gridNx2;
+ L3 = gridNx3;
+ SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
SPtr<Grid3D> grid(new Grid3D(comm));
grid->setDeltaX(dx);
grid->setBlockNX(blocknx1, blocknx2, blocknx3);
+ //sphere
+ //SPtr<GbObject3D> sphere(new GbSphere3D(L1 * 0.5, L2 * 0.5, L3 * 0.5, radius));
+ SPtr<GbObject3D> sphere(new GbSphere3D(75, 25, 25, radius));
+ GbSystem3D::writeGeoObject(sphere.get(), outputPath + "/geo/sphere", WbWriterVtkXmlBinary::getInstance());
+ SPtr<D3Q27Interactor> sphereInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(sphere, grid, noSlipBCAdapter, Interactor3D::SOLID));
- if (grid->getTimeStep() == 0)
+ if (true)
{
const int baseLevel = 0;
@@ -123,64 +123,48 @@ void run(string configname)
}
SPtr<GbObject3D> gridCube(new GbCuboid3D(d_minX1, d_minX2, d_minX3, d_maxX1, d_maxX2, d_maxX3));
- if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), outputPath + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
GenBlocksGridVisitor genBlocks(gridCube);
grid->accept(genBlocks);
- //sphere
- //SPtr<GbObject3D> sphereRef(new GbSphere3D(L1/4.0, L2*0.5, L3*0.5, radius+1.0));
- //GbSystem3D::writeGeoObject(sphereRef.get(),pathname + "/geo/sphereRef", WbWriterVtkXmlBinary::getInstance());
-
-
- //sphere
- SPtr<GbObject3D> sphere(new GbSphere3D(L1*0.5, L2*0.5, L3*0.5, radius));
- //SPtr<GbObject3D> sphere(new GbSphere3D(L1/2.0-4.0, L2*0.5+4.0, L3*0.5+4.0, radius));
- //SPtr<GbObject3D> sphere(new GbCuboid3D(L1/4.0-radius, L2/2.0-radius, L3/2.0-radius, L1/4.0+radius, L2/2.0+radius, L3/2.0+radius));
- GbSystem3D::writeGeoObject(sphere.get(), pathname + "/geo/sphere", WbWriterVtkXmlBinary::getInstance());
-
double off = 0.0;
SPtr<GbObject3D> refCube(new GbCuboid3D(sphere->getX1Minimum() - off, sphere->getX2Minimum() - off, sphere->getX3Minimum(),
sphere->getX1Maximum() + off, sphere->getX2Maximum() + off, sphere->getX3Maximum()));
- if (myid == 0) GbSystem3D::writeGeoObject(refCube.get(), pathname + "/geo/refCube", WbWriterVtkXmlBinary::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(refCube.get(), outputPath + "/geo/refCube", WbWriterVtkXmlBinary::getInstance());
if (refineLevel > 0)
{
if (myid == 0) UBLOG(logINFO, "Refinement - start");
RefineCrossAndInsideGbObjectHelper refineHelper(grid, refineLevel, comm);
- refineHelper.addGbObject(sphere, refineLevel);
- //refineHelper.addGbObject(refCube, refineLevel);
+ //refineHelper.addGbObject(sphere, refineLevel);
+ refineHelper.addGbObject(refCube, refineLevel);
refineHelper.refine();
if (myid == 0) UBLOG(logINFO, "Refinement - end");
}
//walls
GbCuboid3DPtr addWallYmin(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_minX2, d_maxX3 + 4.0*blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(addWallYmin.get(), pathname + "/geo/addWallYmin", WbWriterVtkXmlASCII::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallYmin.get(), outputPath + "/geo/addWallYmin", WbWriterVtkXmlASCII::getInstance());
GbCuboid3DPtr addWallZmin(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_minX3));
- if (myid == 0) GbSystem3D::writeGeoObject(addWallZmin.get(), pathname + "/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallZmin.get(), outputPath + "/geo/addWallZmin", WbWriterVtkXmlASCII::getInstance());
GbCuboid3DPtr addWallYmax(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_maxX2, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(addWallYmax.get(), pathname + "/geo/addWallYmax", WbWriterVtkXmlASCII::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallYmax.get(), outputPath + "/geo/addWallYmax", WbWriterVtkXmlASCII::getInstance());
GbCuboid3DPtr addWallZmax(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_maxX3, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(addWallZmax.get(), pathname + "/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(addWallZmax.get(), outputPath + "/geo/addWallZmax", WbWriterVtkXmlASCII::getInstance());
//inflow
GbCuboid3DPtr geoInflow(new GbCuboid3D(d_minX1 - 4.0*blockLength, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_minX1, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), pathname + "/geo/geoInflow", WbWriterVtkXmlASCII::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), outputPath + "/geo/geoInflow", WbWriterVtkXmlASCII::getInstance());
//outflow
GbCuboid3DPtr geoOutflow(new GbCuboid3D(d_maxX1, d_minX2 - 4.0*blockLength, d_minX3 - 4.0*blockLength, d_maxX1 + 4.0*blockLength, d_maxX2 + 4.0*blockLength, d_maxX3 + 4.0*blockLength));
- if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), outputPath + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
- SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
-
-
-
- //sphere
- SPtr<D3Q27Interactor> sphereInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(sphere, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), outputPath, WbWriterVtkXmlBinary::getInstance(), comm));
//walls
SPtr<D3Q27Interactor> addWallYminInt(new D3Q27Interactor(addWallYmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
@@ -204,7 +188,7 @@ void run(string configname)
SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::B));
InteractorsHelper intHelper(grid, metisVisitor);
- //intHelper.addInteractor(sphereInt);
+ intHelper.addInteractor(sphereInt);
intHelper.addInteractor(addWallYminInt);
intHelper.addInteractor(addWallZminInt);
intHelper.addInteractor(addWallYmaxInt);
@@ -217,11 +201,6 @@ void run(string configname)
//PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
//grid->accept(pqPartVisitor);
- ////set connectors
- //InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- //SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- //grid->accept(setConnsVisitor);
-
ppblocks->process(0);
ppblocks.reset();
@@ -241,8 +220,8 @@ void run(string configname)
UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
}
- //SPtr<LBMKernel> kernel(new IncompressibleCumulantLBMKernel());
- SPtr<LBMKernel> kernel(new CompressibleCumulantLBMKernel());
+ SPtr<LBMKernel> kernel(new IncompressibleCumulantLBMKernel());
+ //SPtr<LBMKernel> kernel(new CompressibleCumulantLBMKernel());
SPtr<BCProcessor> bcProcessor(new BCProcessor());
@@ -270,14 +249,14 @@ void run(string configname)
//initialization of distributions
InitDistributionsBlockVisitor initVisitor;
- //initVisitor.setVx1(fct);
+ initVisitor.setVx1(fct);
//initVisitor.setRho(fctRoh);
grid->accept(initVisitor);
//Postrozess
SPtr<UbScheduler> geoSch(new UbScheduler(1));
SPtr<CoProcessor> ppgeo(
- new WriteBoundaryConditionsCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ new WriteBoundaryConditionsCoProcessor(grid, geoSch, outputPath, WbWriterVtkXmlBinary::getInstance(), comm));
ppgeo->process(0);
ppgeo.reset();
@@ -286,28 +265,34 @@ void run(string configname)
else
{
UBLOG(logINFO, "SetConnectors - start, id=" << myid);
+ }
- //set connectors
- //SPtr<InterpolationProcessor> iProcessor(new IncompressibleOffsetInterpolationProcessor());
- SPtr<CompressibleOffsetMomentsInterpolationProcessor> iProcessor(new CompressibleOffsetMomentsInterpolationProcessor());
- SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- //SPtr<ConnectorFactory> cFactory(new Block3DConnectorFactory());
- //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, cFactory);
- grid->accept(setConnsVisitor);
+
- UBLOG(logINFO, "SetConnectors - end, id=" << myid);
- }
+ UBLOG(logINFO, "SetConnectors - start, id=" << myid);
+ //set connectors
+ SPtr<InterpolationProcessor> iProcessor(new IncompressibleOffsetInterpolationProcessor());
+ //SPtr<CompressibleOffsetMomentsInterpolationProcessor> iProcessor(new CompressibleOffsetMomentsInterpolationProcessor());
+ SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
+ grid->accept(setConnsVisitor);
+ UBLOG(logINFO, "SetConnectors - end, id=" << myid);
SPtr<UbScheduler> stepSch(new UbScheduler(outstep));
//stepSch->addSchedule(10000, 0, 1000000);
- SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, stepSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv,comm));
+ SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, stepSch, outputPath, WbWriterVtkXmlBinary::getInstance(), conv,comm));
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
SPtr<CoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
+ double area = UbMath::PI * radius * radius;
+ SPtr<UbScheduler> forceSch(new UbScheduler(1));
+ SPtr<CalculateForcesCoProcessor> fp = make_shared<CalculateForcesCoProcessor>(grid, forceSch, outputPath + "/forces/forces.txt", comm, uLB, area);
+ fp->addInteractor(sphereInt);
+
SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endstep));
calculator->addCoProcessor(npr);
+ calculator->addCoProcessor(fp);
calculator->addCoProcessor(writeMQCoProcessor);
@@ -336,14 +321,15 @@ int main(int argc, char* argv[])
{
if (argv != NULL)
{
- if (argv[1] != NULL)
- {
- run(string(argv[1]));
- }
- else
- {
- cout << "Configuration file is missing!" << endl;
- }
+ //if (argv[1] != NULL)
+ //{
+ // run(string(argv[1]));
+ //}
+ //else
+ //{
+ // cout << "Configuration file is missing!" << endl;
+ //}
+ run("sphere");
}
}
diff --git a/src/cpu/VirtualFluidsCore/Data/DataSet3D.h b/src/cpu/VirtualFluidsCore/Data/DataSet3D.h
index d0ed4c60f..e869c39cd 100644
--- a/src/cpu/VirtualFluidsCore/Data/DataSet3D.h
+++ b/src/cpu/VirtualFluidsCore/Data/DataSet3D.h
@@ -1,3 +1,36 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file DataSet3D.h
+//! \ingroup Data
+//! \author Konstantin Kutscher
+//=======================================================================================
+
#ifndef DataSet3D_h
#define DataSet3D_h
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp
index 1c8784df0..a90f00afc 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.cpp
@@ -18,7 +18,7 @@ void Thixotropy::setYieldStress(LBMReal yieldStress)
{
tau0 = yieldStress;
}
-LBMReal Thixotropy::getYieldStress()
+LBMReal Thixotropy::getYieldStress() const
{
return tau0;
}
@@ -26,7 +26,7 @@ void Thixotropy::setViscosityParameter(LBMReal kParameter)
{
k = kParameter;
}
-LBMReal Thixotropy::getViscosityParameter()
+LBMReal Thixotropy::getViscosityParameter() const
{
return k;
}
@@ -34,7 +34,7 @@ void Thixotropy::setPowerIndex(LBMReal index)
{
n = index;
}
-LBMReal Thixotropy::getPowerIndex()
+LBMReal Thixotropy::getPowerIndex() const
{
return n;
}
@@ -43,7 +43,7 @@ void Thixotropy::setOmegaMin(LBMReal omega)
{
omegaMin = omega;
}
-LBMReal Thixotropy::getOmegaMin()
+LBMReal Thixotropy::getOmegaMin() const
{
return omegaMin;
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
index 54517d283..677e488e6 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
@@ -12,16 +12,16 @@ public:
Thixotropy& operator=(Thixotropy const&) = delete;
static SPtr<Thixotropy> getInstance();
void setYieldStress(LBMReal tau0);
- LBMReal getYieldStress();
+ LBMReal getYieldStress() const;
void setViscosityParameter(LBMReal k);
- LBMReal getViscosityParameter();
+ LBMReal getViscosityParameter() const;
void setPowerIndex(LBMReal n);
- LBMReal getPowerIndex();
+ LBMReal getPowerIndex() const;
void setOmegaMin(LBMReal omegaMin);
- LBMReal getOmegaMin();
+ LBMReal getOmegaMin() const;
static LBMReal getBinghamCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
static LBMReal getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho);
@@ -53,7 +53,6 @@ inline LBMReal Thixotropy::getHerschelBulkleyCollFactor(LBMReal omegaInf, LBMRea
LBMReal gammaDot = shearRate;
LBMReal omega = omegaInf;
LBMReal epsilon = 1;
- LBMReal omegaMin = 0.7;
while (epsilon > 1e-10)
{
--
GitLab