diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 210268968abe6fc0abf3eac0ba97e4310a9b4b32..436db3c5d3cc2f9c3bc867ce633bb6d9585ef0a8 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -10,6 +10,7 @@ add_subdirectory(${APPS_ROOT_CPU}/FlowAroundCylinder)
add_subdirectory(${APPS_ROOT_CPU}/LaminarTubeFlow)
add_subdirectory(${APPS_ROOT_CPU}/MultiphaseDropletTest)
add_subdirectory(${APPS_ROOT_CPU}/RisingBubble2D)
+add_subdirectory(${APPS_ROOT_CPU}/JetBreakup)
#add_subdirectory(tests)
#add_subdirectory(Applications/gridRf)
#add_subdirectory(Applications/greenvortex)
diff --git a/apps/cpu/JetBreakup/JetBreakup.cfg b/apps/cpu/JetBreakup/JetBreakup.cfg
index 22d20f7d5667ae30c2f3405334566c4d109e7d9f..8d39e8a730434f4afc831a2796deaee6225d3de1 100644
--- a/apps/cpu/JetBreakup/JetBreakup.cfg
+++ b/apps/cpu/JetBreakup/JetBreakup.cfg
@@ -1,39 +1,35 @@
-pathname = d:/temp/Multiphase
-pathGeo = d:/Projects/VirtualFluids-Multiphase/source/Applications/Multiphase/backup
-geoFile = JetBreakup2.ASCII.stl
+pathname = d:/temp/JetBreakupBenchmark_D50_2
+#pathGeo = d:/Projects/VirtualFluids-Multiphase/source/Applications/Multiphase/backup
+pathGeo = d:/Projects/VirtualFluidsCombined/apps/cpu/Multiphase/backup
+#geoFile = JetBreakupR.ASCII.stl
+#geoFile = inlet1.stl
+geoFile = tubeTransformed.stl
+
numOfThreads = 4
availMem = 10e9
#Grid
-
-#boundingBox = -1.0 121.0 0.5 629.0 -1.0 121.0 #(Jet Breakup) (Original with inlet length)
-#boundingBox = -60.5 60.5 -1.0 -201.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
-#blocknx = 22 20 22
-
-boundingBox = -60.5 60.5 -1.0 -21.0 -60.5 60.5 #(Jet Breakup2) (Original without inlet length)
-blocknx = 22 20 22
-
-
-dx = 0.5
-refineLevel = 0
+blocknx = 10 10 10
#Simulation
-uLB = 0.05 #inlet velocity
-uF2 = 0.0001
-Re = 10
-nuL = 1.0e-5 #!1e-2
-nuG = 1.16e-4 #!1e-2
-densityRatio = 10 #30
-sigma = 4.66e-3 #surface tension 1e-4 ./. 1e-5
-interfaceThickness = 5
-radius = 615.0 (Jet Breakup)
+case = 2
+U_LB = 0.01 #inlet velocity
+#uF2 = 0.0001
+#Re = 10
+#nuL =0.00016922169811320757# 1.0e-5 #!1e-2
+#nuG =0.00016922169811320757# 1.16e-4 #!1e-2
+#densityRatio = 24.579710144927535
+#sigma = 1.7688679245283022e-07
+interfaceWidth = 5
+
+D = 0.0001 # m
+D_LB = 50
+
contactAngle = 110.0
-gravity = 0.0
-#gravity = -5.04e-6
phi_L = 0.0
phi_H = 1.0
Phase-field Relaxation = 0.6
-Mobility = 0.02 # 0.01 ./. 0.08, fine correction of Phase-field Relaxation parameter, to activate it need to change in kernel tauH to tauH1
+Mobility = 0.02 # 0.01 ./. 0.08, fine correction of Phase-field Relaxation parameter, to activate it need to change in kernel tauH to tauH1
logToFile = false
@@ -44,5 +40,5 @@ restartStep = 100000
cpStart = 100000
cpStep = 100000
-outTime = 1
-endTime = 200000000
\ No newline at end of file
+outTime = 1000
+endTime = 36000
\ No newline at end of file
diff --git a/apps/cpu/JetBreakup/JetBreakup.cpp b/apps/cpu/JetBreakup/JetBreakup.cpp
index eb7d705537e4307e4ca1066ac9d06dafb72449f4..549748f23c0e14afd09b1aa93ccda762d237e37b 100644
--- a/apps/cpu/JetBreakup/JetBreakup.cpp
+++ b/apps/cpu/JetBreakup/JetBreakup.cpp
@@ -1,516 +1,594 @@
#include <iostream>
+#include <memory>
#include <string>
#include "VirtualFluids.h"
using namespace std;
+void setInflowBC(double x1, double x2, double x3, double radius, int dir)
+{
+
+}
void run(string configname)
{
- try
- {
- vf::basics::ConfigurationFile config;
- config.load(configname);
-
- string pathname = config.getString("pathname");
- string pathGeo = config.getString("pathGeo");
- string geoFile = config.getString("geoFile");
- int numOfThreads = config.getInt("numOfThreads");
- vector<int> blocknx = config.getVector<int>("blocknx");
- vector<double> boundingBox = config.getVector<double>("boundingBox");
- //vector<double> length = config.getVector<double>("length");
- double uLB = config.getDouble("uLB");
- double uF2 = config.getDouble("uF2");
- double nuL = config.getDouble("nuL");
- double nuG = config.getDouble("nuG");
- double densityRatio = config.getDouble("densityRatio");
- double sigma = config.getDouble("sigma");
- int interfaceThickness = config.getInt("interfaceThickness");
- double radius = config.getDouble("radius");
- double theta = config.getDouble("contactAngle");
- double gr = config.getDouble("gravity");
- double phiL = config.getDouble("phi_L");
- double phiH = config.getDouble("phi_H");
- double tauH = config.getDouble("Phase-field Relaxation");
- double mob = config.getDouble("Mobility");
-
-
- double endTime = config.getDouble("endTime");
- double outTime = config.getDouble("outTime");
- double availMem = config.getDouble("availMem");
- int refineLevel = config.getInt("refineLevel");
- double Re = config.getDouble("Re");
- double dx = config.getDouble("dx");
- bool logToFile = config.getBool("logToFile");
- double restartStep = config.getDouble("restartStep");
- double cpStart = config.getValue<double>("cpStart");
- double cpStep = config.getValue<double>("cpStep");
- bool newStart = config.getValue<bool>("newStart");
-
- double beta = 12 * sigma / interfaceThickness;
- double kappa = 1.5 * interfaceThickness * sigma;
-
- CommunicatorPtr comm = vf::mpi::MPICommunicator::getInstance();
- int myid = comm->getProcessID();
-
- if (logToFile)
- {
+ try {
+
+ // Sleep(30000);
+
+ vf::basics::ConfigurationFile config;
+ config.load(configname);
+
+ string pathname = config.getValue<string>("pathname");
+ string pathGeo = config.getValue<string>("pathGeo");
+ string geoFile = config.getValue<string>("geoFile");
+ int numOfThreads = config.getValue<int>("numOfThreads");
+ vector<int> blocknx = config.getVector<int>("blocknx");
+ //vector<double> boundingBox = config.getVector<double>("boundingBox");
+ // vector<double> length = config.getVector<double>("length");
+ double U_LB = config.getValue<double>("U_LB");
+ // double uF2 = config.getValue<double>("uF2");
+ //double nuL = config.getValue<double>("nuL");
+ //double nuG = config.getValue<double>("nuG");
+ //double densityRatio = config.getValue<double>("densityRatio");
+ //double sigma = config.getValue<double>("sigma");
+ int interfaceWidth = config.getValue<int>("interfaceWidth");
+ //double D = config.getValue<double>("D");
+ double theta = config.getValue<double>("contactAngle");
+ double D_LB = config.getValue<double>("D_LB");
+ double phiL = config.getValue<double>("phi_L");
+ double phiH = config.getValue<double>("phi_H");
+ double tauH = config.getValue<double>("Phase-field Relaxation");
+ double mob = config.getValue<double>("Mobility");
+
+ double endTime = config.getValue<double>("endTime");
+ double outTime = config.getValue<double>("outTime");
+ double availMem = config.getValue<double>("availMem");
+ //int refineLevel = config.getValue<int>("refineLevel");
+ //double Re = config.getValue<double>("Re");
+
+ bool logToFile = config.getValue<bool>("logToFile");
+ double restartStep = config.getValue<double>("restartStep");
+ double cpStart = config.getValue<double>("cpStart");
+ double cpStep = config.getValue<double>("cpStep");
+ bool newStart = config.getValue<bool>("newStart");
+
+
+
+ int caseN = config.getValue<int>("case");
+
+ SPtr<vf::mpi::Communicator> comm = vf::mpi::MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ if (myid == 0)
+ UBLOG(logINFO, "Jet Breakup: Start!");
+
+ 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());
- }
- }
-
- //Sleep(30000);
-
- LBMReal dLB; // = length[1] / dx;
- LBMReal rhoLB = 0.0;
- LBMReal nuLB = nuL; //(uLB*dLB) / Re;
-
- LBMUnitConverterPtr conv = LBMUnitConverterPtr(new LBMUnitConverter());
-
- const int baseLevel = 0;
-
-
-
- Grid3DPtr grid(new Grid3D(comm));
- //grid->setPeriodicX1(true);
- //grid->setPeriodicX2(true);
- //grid->setPeriodicX3(true);
- //////////////////////////////////////////////////////////////////////////
- //restart
- UbSchedulerPtr rSch(new UbScheduler(cpStep, cpStart));
- //RestartCoProcessor rp(grid, rSch, comm, pathname, RestartCoProcessor::TXT);
- MPIIORestart1CoProcessor rcp(grid, rSch, pathname, comm);
- //////////////////////////////////////////////////////////////////////////
-
-
-
-
-
- mu::Parser fctF1;
- //fctF1.SetExpr("vy1*(1-((x1-x0)^2+(x3-z0)^2)/(R^2))");
- //fctF1.SetExpr("vy1*(1-(sqrt((x1-x0)^2+(x3-z0)^2)/R))^0.1");
- fctF1.SetExpr("vy1");
- fctF1.DefineConst("vy1", -uLB);
- fctF1.DefineConst("R", 8.0);
- fctF1.DefineConst("x0", 0.0);
- fctF1.DefineConst("z0", 0.0);
-
-
- if (newStart)
- {
-
- //bounding box
- /*double g_minX1 = 0.0;
- double g_minX2 = -length[1] / 2.0;
- double g_minX3 = -length[2] / 2.0;
-
- double g_maxX1 = length[0];
- double g_maxX2 = length[1] / 2.0;
- double g_maxX3 = length[2] / 2.0;*/
-
- double g_minX1 = boundingBox[0];
- double g_minX2 = boundingBox[2];
- double g_minX3 = boundingBox[4];
-
- double g_maxX1 = boundingBox[1];
- double g_maxX2 = boundingBox[3];
- double g_maxX3 = boundingBox[5];
-
- //geometry
-
- //GbObject3DPtr innerCube(new GbCuboid3D(g_minX1+2, g_minX2+2, g_minX3+2, g_maxX1-2, g_maxX2-2, g_maxX3-2));
-
- //GbObject3DPtr cylinder1(new GbCylinder3D(g_minX1 - 2.0*dx, g_maxX2/2, g_maxX3/2, g_minX1 + 12.0*dx, g_maxX2/2, g_maxX3/2, radius));
- //GbObject3DPtr cylinder2(new GbCylinder3D(g_minX1 + 12.0*dx, g_maxX2/2, g_maxX3/2, g_maxX1 + 2.0*dx, g_maxX2/2, g_maxX3/2, dLB / 2.0));
-
- //GbObject3DPtr cylinder(new GbCylinder3D(g_minX1 - 2.0*dx, g_maxX2/2, g_maxX3/2, g_maxX1 + 2.0*dx, g_maxX2/2, g_maxX3/2, dLB / 2.0));
- //GbObject3DPtr cylinders(new GbObject3DManager());
- //GbObject3DPtr cylinders1(new GbObjectGroup3D());
-
-
-
-
- GbObject3DPtr 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());
-
- GbTriFaceMesh3DPtr cylinder;
- if (myid == 0) UBLOG(logINFO, "Read geoFile:start");
- //cylinder = GbTriFaceMesh3DPtr(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile2(pathGeo+"/"+geoFile, "geoCylinders", GbTriFaceMesh3D::KDTREE_SAHPLIT, false));
- cylinder = GbTriFaceMesh3DPtr(GbTriFaceMesh3DCreator::getInstance()->readMeshFromSTLFile(pathGeo + "/" + geoFile, "geoCylinders", GbTriFaceMesh3D::KDTREE_SAHPLIT));
- GbSystem3D::writeGeoObject(cylinder.get(), pathname + "/geo/Stlgeo", WbWriterVtkXmlBinary::getInstance());
-
-
-
- //inflow
- //GbCuboid3DPtr geoInflowF1(new GbCuboid3D(40.0, 628.0, 40.0, 80, 631.0, 80.0)); // For JetBreakup (Original)
- //GbCuboid3DPtr geoInflowF1(new GbCuboid3D(g_minX1-2.0*dx, g_minX2-2.0*dx, g_minX3-2.0*dx, g_maxX1+2.0*dx, g_minX2+2.0*dx, g_maxX3+2.0*dx));
- //if (myid == 0) GbSystem3D::writeGeoObject(geoInflowF1.get(), pathname + "/geo/geoInflowF1", WbWriterVtkXmlASCII::getInstance());
-
-
- ////outflow
- ////GbCuboid3DPtr geoOutflow(new GbCuboid3D(-1.0, -1, -1.0, 121.0, 1.0, 121.0)); // For JetBreakup (Original)
- //GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_minX1-2.0*dx, g_maxX2, g_minX3-2.0*dx, g_maxX1+2.0*dx, g_maxX2+2.0*dx, g_maxX3+2.0*dx));
- //if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr geoInflowF1(new GbCuboid3D(g_minX1, g_minX2-0.5*dx, g_minX3, g_maxX1, g_minX2 - 1.0*dx, g_maxX3));
- if (myid==0) GbSystem3D::writeGeoObject(geoInflowF1.get(), pathname+"/geo/geoInflowF1", WbWriterVtkXmlASCII::getInstance());
-
-
- //outflow
- //GbCuboid3DPtr geoOutflow(new GbCuboid3D(-1.0, -1, -1.0, 121.0, 1.0, 121.0)); // For JetBreakup (Original)
- GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_minX1, g_maxX2-1*dx, g_minX3, g_maxX1, g_maxX2, g_maxX3));
- if (myid==0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname+"/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
-
- double blockLength = blocknx[0] * dx;
-
-
-
- if (myid == 0)
- {
- UBLOG(logINFO, "uLb = " << uLB);
- UBLOG(logINFO, "rho = " << rhoLB);
- UBLOG(logINFO, "nuLb = " << nuLB);
- UBLOG(logINFO, "Re = " << Re);
- UBLOG(logINFO, "dx = " << dx);
- UBLOG(logINFO, "Preprocess - start");
- }
-
- grid->setDeltaX(dx);
- grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
-
- grid->setPeriodicX1(false);
- grid->setPeriodicX2(false);
- grid->setPeriodicX3(false);
-
-
-
- GenBlocksGridVisitor genBlocks(gridCube);
- grid->accept(genBlocks);
-
-
-
-
- //BC Adapter
- //////////////////////////////////////////////////////////////////////////////
- BCAdapterPtr noSlipBCAdapter(new NoSlipBCAdapter());
- noSlipBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new NoSlipBCAlgorithmMultiphase()));
-
-
- BCAdapterPtr denBCAdapter(new DensityBCAdapter(rhoLB));
- denBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new NonReflectingOutflowBCAlgorithmMultiphase()));
-
- double r = 5.0; //boost::dynamic_pointer_cast<GbCylinder3D>(cylinder)->getRadius();
- double cx1 = g_minX1;
- double cx2 = 0.0; //cylinder->getX2Centroid();
- double cx3 = 0.0; //cylinder->getX3Centroid();
-
-
-
- mu::Parser fctPhi_F1;
- fctPhi_F1.SetExpr("phiH");
- fctPhi_F1.DefineConst("phiH", phiH);
-
- mu::Parser fctPhi_F2;
- fctPhi_F2.SetExpr("phiL");
- fctPhi_F2.DefineConst("phiL", phiL);
-
- mu::Parser fctvel_F2_init;
- fctvel_F2_init.SetExpr("U");
- fctvel_F2_init.DefineConst("U", 0);
-
- //fct.SetExpr("U");
- //fct.DefineConst("U", uLB);
- //BCAdapterPtr velBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
-
- BCAdapterPtr velBCAdapterF1(new VelocityBCAdapterMultiphase(false, true, false, fctF1, phiH, 0.0, endTime));
-
- //BCAdapterPtr velBCAdapterF2_1_init(new VelocityBCAdapterMultiphase(false, false, true, fctF2_1, phiH, 0.0, endTime));
- //BCAdapterPtr velBCAdapterF2_2_init(new VelocityBCAdapterMultiphase(false, false, true, fctF2_2, phiH, 0.0, endTime));
-
- //BCAdapterPtr velBCAdapterF2_1_init(new VelocityBCAdapterMultiphase(false, false, true, fctvel_F2_init, phiL, 0.0, endTime));
- //BCAdapterPtr velBCAdapterF2_2_init(new VelocityBCAdapterMultiphase(false, false, true, fctvel_F2_init, phiL, 0.0, endTime));
-
- velBCAdapterF1->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithmMultiphase()));
- //velBCAdapterF2_1_init->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithmMultiphase()));
- //velBCAdapterF2_2_init->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithmMultiphase()));
-
-
- //velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityWithDensityBCAlgorithm()));
- //mu::Parser fct;
- //fct.SetExpr("U");
- //fct.DefineConst("U", uLB);
- //BCAdapterPtr velBCAdapter(new VelocityBCAdapter(true, false, false, fct, 0, BCFunction::INFCONST));
- //velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new NonReflectingVelocityBCAlgorithm()));
-
-
- //////////////////////////////////////////////////////////////////////////////////
- //BC visitor
- BoundaryConditionsBlockVisitorMultiphase bcVisitor;
- bcVisitor.addBC(noSlipBCAdapter);
- bcVisitor.addBC(denBCAdapter);
- bcVisitor.addBC(velBCAdapterF1);
- //bcVisitor.addBC(velBCAdapterF2_1_init);
- //bcVisitor.addBC(velBCAdapterF2_2_init);
-
-
-
- WriteBlocksCoProcessorPtr ppblocks(new WriteBlocksCoProcessor(grid, UbSchedulerPtr(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
-
- ppblocks->process(0);
-
- Interactor3DPtr tubes(new D3Q27TriFaceMeshInteractor(cylinder, grid, noSlipBCAdapter, Interactor3D::SOLID));
-
- D3Q27InteractorPtr inflowF1Int = D3Q27InteractorPtr(new D3Q27Interactor(geoInflowF1, grid, velBCAdapterF1, Interactor3D::SOLID));
-
- //D3Q27InteractorPtr inflowF2_1Int_init = D3Q27InteractorPtr(new D3Q27Interactor(geoInflowF2_1, grid, velBCAdapterF2_1_init, Interactor3D::SOLID));
-
- //D3Q27InteractorPtr inflowF2_2Int_init = D3Q27InteractorPtr(new D3Q27Interactor(geoInflowF2_2, grid, velBCAdapterF2_2_init, Interactor3D::SOLID));
-
- D3Q27InteractorPtr outflowInt = D3Q27InteractorPtr(new D3Q27Interactor(geoOutflow, grid, denBCAdapter, Interactor3D::SOLID));
-
- //SetSolidBlockVisitor visitor1(inflowF2_1Int, SetSolidBlockVisitor::BC);
- //grid->accept(visitor1);
- //SetSolidBlockVisitor visitor2(inflowF2_2Int, SetSolidBlockVisitor::BC);
- //grid->accept(visitor2);
-
-
- Grid3DVisitorPtr metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW));
- InteractorsHelper intHelper(grid, metisVisitor);
- intHelper.addInteractor(tubes);
- intHelper.addInteractor(inflowF1Int);
- intHelper.addInteractor(outflowInt);
- intHelper.selectBlocks();
-
-
- ppblocks->process(0);
- ppblocks.reset();
-
- unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
- int ghostLayer = 3;
- unsigned long long numberOfNodesPerBlock = (unsigned long long)(blocknx[0]) * (unsigned long long)(blocknx[1]) * (unsigned long long)(blocknx[2]);
- unsigned long long numberOfNodes = numberOfBlocks * numberOfNodesPerBlock;
- unsigned long long numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (blocknx[0] + ghostLayer) * (blocknx[1] + ghostLayer) * (blocknx[2] + ghostLayer);
- double needMemAll = double(numberOfNodesPerBlockWithGhostLayer * (27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
- double needMem = needMemAll / double(comm->getNumberOfProcesses());
-
- if (myid == 0)
- {
- UBLOG(logINFO, "Number of blocks = " << numberOfBlocks);
- UBLOG(logINFO, "Number of nodes = " << numberOfNodes);
- 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 * numberOfNodesPerBlock);
+ if (myid == 0) {
+ const char *str = pathname.c_str();
+ mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
}
- UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
- UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
- UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
- }
-
- LBMKernelPtr kernel;
-
- kernel = LBMKernelPtr(new MultiphaseCumulantLBMKernel(blocknx[0], blocknx[1], blocknx[2], MultiphaseCumulantLBMKernel::NORMAL));
-
- kernel->setWithForcing(true);
- kernel->setForcingX1(0.0);
- kernel->setForcingX2(gr);
- kernel->setForcingX3(0.0);
-
- kernel->setPhiL(phiL);
- kernel->setPhiH(phiH);
- kernel->setPhaseFieldRelaxation(tauH);
- kernel->setMobility(mob);
-
- BCProcessorPtr bcProc(new BCProcessor());
- //BCProcessorPtr bcProc(new ThinWallBCProcessor());
-
- kernel->setBCProcessor(bcProc);
-
- SetKernelBlockVisitorMultiphase kernelVisitor(kernel, nuL, nuG, densityRatio, beta, kappa, theta, availMem, needMem);
-
- grid->accept(kernelVisitor);
-
- if (refineLevel > 0)
- {
- SetUndefinedNodesBlockVisitor undefNodesVisitor;
- grid->accept(undefNodesVisitor);
- }
-
- //inflowF2_1Int->initInteractor();
- //inflowF2_2Int->initInteractor();
-
- intHelper.setBC();
+#endif
-
- grid->accept(bcVisitor);
-
- //initialization of distributions
- LBMReal x1c = radius; //g_minX1; //radius; //19; //(g_maxX1+g_minX1)/2;
- LBMReal x2c = (g_maxX2 + g_minX2) / 2; //g_minX2 + 2;
- LBMReal x3c = (g_maxX3 + g_minX3) / 2;
- mu::Parser fct1;
-
- //fct1.SetExpr("0.5-0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
- //fct1.SetExpr("phiM-phiM*tanh((sqrt((x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/(interfaceThickness*phiM))");
-
- //fct1.SetExpr("0.5*(phiH + phiL)-0.5*(phiH - phiL)*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
-
-
- //fct1.SetExpr("0.5*(phiH + phiL) + 0.5*(phiH - phiL)*tanh(2*((x2-radius))/interfaceThickness)");
- fct1.SetExpr("phiL");
- fct1.DefineConst("x1c", x1c);
- fct1.DefineConst("x2c", x2c);
- fct1.DefineConst("x3c", x3c);
- fct1.DefineConst("phiL", phiL);
- fct1.DefineConst("phiH", phiH);
- fct1.DefineConst("radius", radius);
- fct1.DefineConst("interfaceThickness", interfaceThickness);
-
- mu::Parser fct2;
- //fct2.SetExpr("vx1*(1-((x2-y0)^2+(x3-z0)^2)/(R^2))");
- fct2.SetExpr("vx1");
- fct2.DefineConst("R", 10.0);
- fct2.DefineConst("vx1", uLB);
- fct2.DefineConst("y0", 1.0);
- fct2.DefineConst("z0", 31.0);
- /*fct2.SetExpr("0.5*uLB-uLB*0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
- fct2.DefineConst("uLB", uLB);
- fct2.DefineConst("x1c", x1c);
- fct2.DefineConst("x2c", x2c);
- fct2.DefineConst("x3c", x3c);
- fct2.DefineConst("radius", radius);
- fct2.DefineConst("interfaceThickness", interfaceThickness);*/
-
-
- InitDistributionsBlockVisitorMultiphase initVisitor(densityRatio, interfaceThickness, radius);
- initVisitor.setPhi(fct1);
- //initVisitor.setVx1(fct2);
- grid->accept(initVisitor);
-
- //set connectors
- InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- //InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor());
- SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- //ConnectorFactoryPtr factory(new Block3DConnectorFactory());
- //ConnectorBlockVisitor setConnsVisitor(comm, nuLB, iProcessor, factory);
- grid->accept(setConnsVisitor);
-
- //domain decomposition for threads
- //PQueuePartitioningGridVisitor pqPartVisitor(numOfThreads);
- //grid->accept(pqPartVisitor);
-
-
-
-
- //boundary conditions grid
- {
- UbSchedulerPtr geoSch(new UbScheduler(1));
- WriteBoundaryConditionsCoProcessorPtr ppgeo(
- new WriteBoundaryConditionsCoProcessor(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
- ppgeo->process(0);
- ppgeo.reset();
- }
-
- if (myid == 0) UBLOG(logINFO, "Preprocess - end");
- }
- else
- {
- if (myid == 0)
- {
+ if (myid == 0) {
+ stringstream logFilename;
+ logFilename << pathname + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt";
+ UbLog::output_policy::setStream(logFilename.str());
+ }
+ }
+
+ // Sleep(30000);
+
+ double rho_h, rho_l, r_rho, mu_h, mu_l, Uo, D, sigma;
+
+ switch (caseN) {
+ case 1:
+ //density of heavy fluid (kg/m^3)
+ rho_h = 848;
+ //density of light fluid (kg/m^3)
+ rho_l = 34.5;
+ //density retio
+ r_rho = rho_h / rho_l;
+ //dynamic viscosity of heavy fluid (Pa · s)
+ mu_h = 2.87e-3;
+ //dynamic viscosity of light fluid (Pa · s)
+ mu_l = 1.97e-5;
+ //velocity (m/s)
+ Uo = 100;
+ //diameter of jet (m)
+ D = 0.0001;
+ //surface tension (N/m)
+ sigma = 0.03;
+ break;
+ case 2:
+ // density of heavy fluid (kg/m^3)
+ rho_h = 848;
+ // density of light fluid (kg/m^3)
+ rho_l = 1.205;
+ // density retio
+ r_rho = rho_h / rho_l;
+ // dynamic viscosity of heavy fluid (Pa · s)
+ mu_h = 2.87e-3;
+ // dynamic viscosity of light fluid (Pa · s)
+ mu_l = 1.84e-5;
+ // velocity (m/s)
+ Uo = 200;
+ // diameter of jet (m)
+ D = 0.0001;
+ // surface tension (N/m)
+ sigma = 0.03;
+ break;
+ }
+
+ double Re = rho_h * Uo * D / mu_h;
+ double We = rho_h * Uo * Uo * D / sigma;
+
+ double dx = D / D_LB;
+ double nu_h = U_LB * D_LB / Re;
+ double nu_l = nu_h;
+
+ double rho_h_LB = 1;
+ //surface tension
+ double sigma_LB = rho_h_LB * U_LB * U_LB * D_LB / We;
+
+ // LBMReal dLB = 0; // = length[1] / dx;
+ LBMReal rhoLB = 0.0;
+ LBMReal nuLB = nu_l; //(uLB*dLB) / Re;
+
+ double beta = 12.0 * sigma_LB / interfaceWidth;
+ double kappa = 1.5 * interfaceWidth * sigma_LB;
+
+ if (myid == 0) {
UBLOG(logINFO, "Parameters:");
- UBLOG(logINFO, "uLb = " << uLB);
+ UBLOG(logINFO, "U_LB = " << U_LB);
UBLOG(logINFO, "rho = " << rhoLB);
- UBLOG(logINFO, "nuLb = " << nuLB);
+ UBLOG(logINFO, "nu_l = " << nu_l);
+ UBLOG(logINFO, "nu_h = " << nu_h);
UBLOG(logINFO, "Re = " << Re);
+ UBLOG(logINFO, "We = " << We);
UBLOG(logINFO, "dx = " << dx);
- UBLOG(logINFO, "number of levels = " << refineLevel + 1);
+ UBLOG(logINFO, "sigma = " << sigma);
+ UBLOG(logINFO, "density ratio = " << r_rho);
+ // UBLOG(logINFO, "number of levels = " << refineLevel + 1);
UBLOG(logINFO, "numOfThreads = " << numOfThreads);
UBLOG(logINFO, "path = " << pathname);
- }
-
- rcp.restart((int)restartStep);
- grid->setTimeStep(restartStep);
-
- //BCAdapterPtr velBCAdapter(new VelocityBCAdapter());
- //velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityBCAlgorithm()));
- //velBCAdapter->setBcAlgorithm(BCAlgorithmPtr(new VelocityWithDensityBCAlgorithm()));
- //bcVisitor.addBC(velBCAdapter);
- //grid->accept(bcVisitor);
-
- //set connectors
- //InterpolationProcessorPtr iProcessor(new IncompressibleOffsetInterpolationProcessor());
- InterpolationProcessorPtr iProcessor(new CompressibleOffsetInterpolationProcessor());
- SetConnectorsBlockVisitor setConnsVisitor(comm, true, D3Q27System::ENDDIR, nuLB, iProcessor);
- grid->accept(setConnsVisitor);
-
- if (myid == 0) UBLOG(logINFO, "Restart - end");
- }
- UbSchedulerPtr visSch(new UbScheduler(outTime));
- WriteMacroscopicQuantitiesCoProcessor pp(grid, visSch, pathname, WbWriterVtkXmlASCII::getInstance(), conv, comm);
-
- UbSchedulerPtr nupsSch(new UbScheduler(10, 30, 100));
- NUPSCounterCoProcessor npr(grid, nupsSch, numOfThreads, comm);
-
-
-
-
-
+ }
+
+ SPtr<LBMUnitConverter> conv(new LBMUnitConverter());
+
+ // const int baseLevel = 0;
+
+ SPtr<LBMKernel> kernel;
+
+ // kernel = SPtr<LBMKernel>(new MultiphaseScratchCumulantLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsCumulantLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
+ kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
+
+ kernel->setWithForcing(true);
+ kernel->setForcingX1(0.0);
+ kernel->setForcingX2(0.0);
+ kernel->setForcingX3(0.0);
+
+ kernel->setPhiL(phiL);
+ kernel->setPhiH(phiH);
+ kernel->setPhaseFieldRelaxation(tauH);
+ kernel->setMobility(mob);
+
+ // nuL, nuG, densityRatio, beta, kappa, theta,
+
+ kernel->setCollisionFactorMultiphase(nu_h, nu_l);
+ kernel->setDensityRatio(r_rho);
+ kernel->setMultiphaseModelParameters(beta, kappa);
+ kernel->setContactAngle(theta);
+ kernel->setInterfaceWidth(interfaceWidth);
+ dynamicPointerCast<MultiphasePressureFilterLBMKernel>(kernel)->setPhaseFieldBC(0.0);
+
+ SPtr<BCProcessor> bcProc(new BCProcessor());
+ // BCProcessorPtr bcProc(new ThinWallBCProcessor());
+
+ kernel->setBCProcessor(bcProc);
+
+ SPtr<Grid3D> grid(new Grid3D(comm));
+ // grid->setPeriodicX1(true);
+ // grid->setPeriodicX2(true);
+ // grid->setPeriodicX3(true);
+ grid->setGhostLayerWidth(2);
+
+ SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(
+ comm, MetisPartitioningGridVisitor::LevelBased, D3Q27System::BSW, MetisPartitioner::RECURSIVE));
+
+ //////////////////////////////////////////////////////////////////////////
+ // restart
+ SPtr<UbScheduler> rSch(new UbScheduler(cpStep, cpStart));
+ // SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(grid, rSch, pathname, comm));
+ SPtr<MPIIOMigrationCoProcessor> rcp(new MPIIOMigrationCoProcessor(grid, rSch, metisVisitor, pathname, comm));
+ // SPtr<MPIIOMigrationBECoProcessor> rcp(new MPIIOMigrationBECoProcessor(grid, rSch, pathname, comm));
+ // rcp->setNu(nuLB);
+ // rcp->setNuLG(nuL, nuG);
+ // rcp->setDensityRatio(densityRatio);
+
+ rcp->setLBMKernel(kernel);
+ rcp->setBCProcessor(bcProc);
+ //////////////////////////////////////////////////////////////////////////
+ // BC Adapter
+ //////////////////////////////////////////////////////////////////////////////
+ mu::Parser fctF1;
+ // fctF1.SetExpr("vy1*(1-((x1-x0)^2+(x3-z0)^2)/(R^2))");
+ // fctF1.SetExpr("vy1*(1-(sqrt((x1-x0)^2+(x3-z0)^2)/R))^0.1");
+ fctF1.SetExpr("vy1");
+ fctF1.DefineConst("vy1", 0.0);
+ fctF1.DefineConst("R", 8.0);
+ fctF1.DefineConst("x0", 0.0);
+ fctF1.DefineConst("z0", 0.0);
+ // SPtr<BCAdapter> velBCAdapterF1(
+ // new MultiphaseVelocityBCAdapter(false, true, false, fctF1, phiH, 0.0, BCFunction::INFCONST));
+
+ mu::Parser fctF2;
+ fctF2.SetExpr("vy1");
+ fctF2.DefineConst("vy1", U_LB);
+
+ double startTime = 30;
+ SPtr<BCAdapter> velBCAdapterF1(
+ new MultiphaseVelocityBCAdapter(true, false, false, fctF1, phiH, 0.0, startTime));
+ SPtr<BCAdapter> velBCAdapterF2(
+ new MultiphaseVelocityBCAdapter(true, false, false, fctF2, phiH, startTime, endTime));
+
+ SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseNoSlipBCAlgorithm()));
+
+ SPtr<BCAdapter> denBCAdapter(new DensityBCAdapter(rhoLB));
+ denBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseNonReflectingOutflowBCAlgorithm()));
+
+ mu::Parser fctPhi_F1;
+ fctPhi_F1.SetExpr("phiH");
+ fctPhi_F1.DefineConst("phiH", phiH);
+
+ mu::Parser fctPhi_F2;
+ fctPhi_F2.SetExpr("phiL");
+ fctPhi_F2.DefineConst("phiL", phiL);
+
+ mu::Parser fctvel_F2_init;
+ fctvel_F2_init.SetExpr("U");
+ fctvel_F2_init.DefineConst("U", 0);
+
+ velBCAdapterF1->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseVelocityBCAlgorithm()));
+ //////////////////////////////////////////////////////////////////////////////////
+ // BC visitor
+ MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
+ bcVisitor.addBC(noSlipBCAdapter);
+ bcVisitor.addBC(denBCAdapter); // Ohne das BB?
+ bcVisitor.addBC(velBCAdapterF1);
+
+ //SPtr<D3Q27Interactor> inflowF1Int;
+ //SPtr<D3Q27Interactor> cylInt;
+
+ SPtr<D3Q27Interactor> inflowInt;
+
+ if (newStart) {
+
+ // if (newStart) {
+
+ // bounding box
+ double g_minX1 = 0;
+ double g_minX2 = 0;
+ double g_minX3 = 0;
+
+ double g_maxX1 = 8.0*D;
+ double g_maxX2 = 2.5*D;
+ double g_maxX3 = 2.5*D;
+
+ // geometry
+ 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());
+
+ //if (myid == 0)
+ // UBLOG(logINFO, "Read geoFile:start");
+ //SPtr<GbTriFaceMesh3D> cylinder = make_shared<GbTriFaceMesh3D>();
+ //cylinder->readMeshFromSTLFileBinary(pathGeo + "/" + geoFile, false);
+ //GbSystem3D::writeGeoObject(cylinder.get(), pathname + "/geo/Stlgeo", WbWriterVtkXmlBinary::getInstance());
+ //if (myid == 0)
+ // UBLOG(logINFO, "Read geoFile:stop");
+ // inflow
+ // GbCuboid3DPtr geoInflowF1(new GbCuboid3D(g_minX1, g_minX2 - 0.5 * dx, g_minX3, g_maxX1, g_minX2 - 1.0 *
+ // dx, g_maxX3));
+ //GbCuboid3DPtr geoInflowF1(new GbCuboid3D(g_minX1 * 0.5 - dx, g_minX2 - dx, g_minX3 * 0.5 - dx,
+ // g_maxX1 * 0.5 + dx, g_minX2, g_maxX3 * 0.5 + dx));
+ //if (myid == 0)
+ // GbSystem3D::writeGeoObject(geoInflowF1.get(), pathname + "/geo/geoInflowF1",
+ // WbWriterVtkXmlASCII::getInstance());
+
+ GbCylinder3DPtr geoInflow(new GbCylinder3D(g_minX1 - 2.0*dx, g_maxX2 / 2.0, g_maxX3 / 2.0, g_minX1,
+ g_maxX2 / 2.0,
+ g_maxX3 / 2.0, D / 2.0));
+ if (myid == 0)
+ GbSystem3D::writeGeoObject(geoInflow.get(), pathname + "/geo/geoInflow",
+ WbWriterVtkXmlASCII::getInstance());
+
+ GbCylinder3DPtr geoSolid(new GbCylinder3D(g_minX1 - 2.0 * dx, g_maxX2 / 2.0, g_maxX3 / 2.0, g_minX1-dx,
+ g_maxX2 / 2.0, g_maxX3 / 2.0, 1.5*D / 2.0));
+ if (myid == 0)
+ GbSystem3D::writeGeoObject(geoSolid.get(), pathname + "/geo/geoSolid",
+ WbWriterVtkXmlASCII::getInstance());
+
+
+ // GbCylinder3DPtr cylinder2(
+ // new GbCylinder3D(0.0, g_minX2 - 2.0 * dx / 2.0, 0.0, 0.0, g_minX2 + 4.0 * dx, 0.0, 8.0+2.0*dx));
+ // if (myid == 0)
+ // GbSystem3D::writeGeoObject(cylinder2.get(), pathname + "/geo/cylinder2",
+ // WbWriterVtkXmlASCII::getInstance());
+ // outflow
+ // GbCuboid3DPtr geoOutflow(new GbCuboid3D(-1.0, -1, -1.0, 121.0, 1.0, 121.0)); // For JetBreakup (Original)
+ // GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_minX1, g_maxX2 - 40 * dx, g_minX3, g_maxX1, g_maxX2, g_maxX3));
+ GbCuboid3DPtr geoOutflow(new GbCuboid3D(g_minX1, g_maxX2, g_minX3, g_maxX1, g_maxX2 + dx, g_maxX3));
+ if (myid == 0)
+ GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow",
+ WbWriterVtkXmlASCII::getInstance());
+
+ // double blockLength = blocknx[0] * dx;
+
+ if (myid == 0) {
+ UBLOG(logINFO, "Preprocess - start");
+ }
- //UbSchedulerPtr bcSch(new UbScheduler(1, 12000, 12000));
- //TimeDependentBCCoProcessorPtr inflowF2 (new TimeDependentBCCoProcessor(grid,bcSch));
- //inflowF2->addInteractor(inflowF2_1Int);
- //inflowF2->addInteractor(inflowF2_2Int);
+ grid->setDeltaX(dx);
+ grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
+
+ grid->setPeriodicX1(false);
+ grid->setPeriodicX2(false);
+ grid->setPeriodicX3(false);
+
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ SPtr<WriteBlocksCoProcessor> ppblocks(new WriteBlocksCoProcessor(
+ grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+
+ //SPtr<Interactor3D> tubes(new D3Q27TriFaceMeshInteractor(cylinder, grid, noSlipBCAdapter,
+ // Interactor3D::SOLID, Interactor3D::POINTS));
+
+ // inflowF1Int =
+ // SPtr<D3Q27Interactor>(new D3Q27Interactor(cylinder1, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ // inflowF1Int->addBCAdapter(velBCAdapterF2);
+
+ SPtr<D3Q27Interactor> outflowInt(new D3Q27Interactor(geoOutflow, grid, denBCAdapter, Interactor3D::SOLID));
+
+ // Create boundary conditions geometry
+ GbCuboid3DPtr wallXmin(
+ new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_minX1, g_maxX2 + 2.0*dx, g_maxX3));
+ GbSystem3D::writeGeoObject(wallXmin.get(), pathname + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallXmax(
+ new GbCuboid3D(g_maxX1, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3));
+ GbSystem3D::writeGeoObject(wallXmax.get(), pathname + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallZmin(
+ new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_minX3));
+ GbSystem3D::writeGeoObject(wallZmin.get(), pathname + "/geo/wallZmin", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallZmax(
+ new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_maxX3, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3 + 2.0*dx));
+ GbSystem3D::writeGeoObject(wallZmax.get(), pathname + "/geo/wallZmax", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallYmin(
+ new GbCuboid3D(g_minX1 - 2.0*dx, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_minX2, g_maxX3));
+ GbSystem3D::writeGeoObject(wallYmin.get(), pathname + "/geo/wallYmin", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallYmax(
+ new GbCuboid3D(g_minX1 - 2.0*dx, g_maxX2, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3));
+ GbSystem3D::writeGeoObject(wallYmax.get(), pathname + "/geo/wallYmax", WbWriterVtkXmlASCII::getInstance());
+
+ // Add boundary conditions to grid generator
+ SPtr<D3Q27Interactor> wallXminInt(
+ new D3Q27Interactor(wallXmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXmaxInt(
+ new D3Q27Interactor(wallXmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ 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, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallYmaxInt(
+ new D3Q27Interactor(wallYmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ //cylInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(cylinder1, grid, velBCAdapterF1, Interactor3D::SOLID));
+ //cylInt->addBCAdapter(velBCAdapterF2);
+ // SPtr<D3Q27Interactor> cyl2Int(new D3Q27Interactor(cylinder2, grid, noSlipBCAdapter,
+ // Interactor3D::SOLID));
+
+ inflowInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, velBCAdapterF1, Interactor3D::SOLID));
+ inflowInt->addBCAdapter(velBCAdapterF2);
+
+ SPtr<D3Q27Interactor> solidInt =
+ SPtr<D3Q27Interactor>(new D3Q27Interactor(geoSolid, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ InteractorsHelper intHelper(grid, metisVisitor, true);
+ //intHelper.addInteractor(cylInt);
+ //intHelper.addInteractor(tubes);
+ // intHelper.addInteractor(outflowInt);
+ // intHelper.addInteractor(cyl2Int);
+
+ intHelper.addInteractor(wallXminInt);
+ intHelper.addInteractor(wallXmaxInt);
+ intHelper.addInteractor(wallZminInt);
+ intHelper.addInteractor(wallZmaxInt);
+ intHelper.addInteractor(wallYminInt);
+ intHelper.addInteractor(wallYmaxInt);
+ intHelper.addInteractor(inflowInt);
+ //intHelper.addInteractor(solidInt);
+
+ intHelper.selectBlocks();
+
+ ppblocks->process(0);
+ ppblocks.reset();
+
+ unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
+ int ghostLayer = 3;
+ unsigned long long numberOfNodesPerBlock =
+ (unsigned long long)(blocknx[0]) * (unsigned long long)(blocknx[1]) * (unsigned long long)(blocknx[2]);
+ unsigned long long numberOfNodes = numberOfBlocks * numberOfNodesPerBlock;
+ unsigned long long numberOfNodesPerBlockWithGhostLayer =
+ numberOfBlocks * (blocknx[0] + ghostLayer) * (blocknx[1] + ghostLayer) * (blocknx[2] + ghostLayer);
+ double needMemAll =
+ double(numberOfNodesPerBlockWithGhostLayer * (27 * sizeof(double) + sizeof(int) + sizeof(float) * 4));
+ double needMem = needMemAll / double(comm->getNumberOfProcesses());
+
+ if (myid == 0) {
+ UBLOG(logINFO, "Number of blocks = " << numberOfBlocks);
+ UBLOG(logINFO, "Number of nodes = " << numberOfNodes);
+ 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 * numberOfNodesPerBlock);
+ }
+ UBLOG(logINFO, "Necessary memory = " << needMemAll << " bytes");
+ UBLOG(logINFO, "Necessary memory per process = " << needMem << " bytes");
+ UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
+ }
- //CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, visSch,CalculationManager::MPI));
- CalculationManagerPtr calculation(new CalculationManager(grid, numOfThreads, endTime, visSch));
- if (myid == 0) UBLOG(logINFO, "Simulation-start");
- calculation->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;
- }
+ MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nu_h, nu_l, availMem, needMem);
+
+ grid->accept(kernelVisitor);
+
+ //if (refineLevel > 0) {
+ // SetUndefinedNodesBlockVisitor undefNodesVisitor;
+ // grid->accept(undefNodesVisitor);
+ //}
+
+ intHelper.setBC();
+
+ // initialization of distributions
+ mu::Parser fct1;
+ fct1.SetExpr("phiL");
+ fct1.DefineConst("phiL", phiL);
+ // MultiphaseInitDistributionsBlockVisitor initVisitor(interfaceThickness);
+ MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
+ initVisitor.setPhi(fct1);
+ grid->accept(initVisitor);
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ //{
+ // std::vector<std::vector<SPtr<Block3D>>> blockVector;
+ // int gridRank = comm->getProcessID();
+ // int minInitLevel = grid->getCoarsestInitializedLevel();
+ // int maxInitLevel = grid->getFinestInitializedLevel();
+ // blockVector.resize(maxInitLevel + 1);
+ // for (int level = minInitLevel; level <= maxInitLevel; level++) {
+ // grid->getBlocks(level, gridRank, true, blockVector[level]);
+ //}
+ // for (int level = minInitLevel; level <= maxInitLevel; level++) {
+ // for (SPtr<Block3D> block : blockVector[level]) {
+ // if (block) {
+ // int ix1 = block->getX1();
+ // int ix2 = block->getX2();
+ // int ix3 = block->getX3();
+ // int level = block->getLevel();
+
+ // for (int dir = 0; dir < D3Q27System::ENDDIR; dir++) {
+ // SPtr<Block3D> neighBlock = grid->getNeighborBlock(dir, ix1, ix2, ix3, level);
+
+ // if (!neighBlock) {
+
+ // }
+ // }
+ // }
+ // }
+ //}
+ // SPtr<Block3D> block = grid->getBlock(0, 0, 0, 0);
+ // SPtr<LBMKernel> kernel = dynamicPointerCast<LBMKernel>(block->getKernel());
+ // SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
+
+ // for (int ix3 = 0; ix3 <= 13; ix3++) {
+ // for (int ix2 = 0; ix2 <= 13; ix2++) {
+ // for (int ix1 = 0; ix1 <= 13; ix1++) {
+ // if (ix1 == 0 || ix2 == 0 || ix3 == 0 || ix1 == 13 || ix2 == 13 || ix3 == 13)
+ // bcArray->setUndefined(ix1, ix2, ix3);
+ // }
+ // }
+ // }
+ //}
+ ////////////////////////////////////////////////////////////////////////////////////////////
+ // boundary conditions grid
+ {
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ SPtr<WriteBoundaryConditionsCoProcessor> ppgeo(new WriteBoundaryConditionsCoProcessor(
+ grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppgeo->process(0);
+ ppgeo.reset();
+ }
+ if (myid == 0)
+ UBLOG(logINFO, "Preprocess - end");
+ } else {
+ rcp->restart((int)restartStep);
+ grid->setTimeStep(restartStep);
+
+ if (myid == 0)
+ UBLOG(logINFO, "Restart - end");
+ }
+
+ // TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
+ // grid->accept(setConnsVisitor);
+
+ // ThreeDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
+
+ grid->accept(bcVisitor);
+
+ // ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ grid->accept(setConnsVisitor);
+
+ SPtr<UbScheduler> visSch(new UbScheduler(outTime));
+ double t_ast, t;
+ t_ast = 7.19;
+ t = (int)(t_ast/(U_LB/(D_LB)));
+ visSch->addSchedule(t,t,t); //t=7.19
+ SPtr<WriteMultiphaseQuantitiesCoProcessor> pp(new WriteMultiphaseQuantitiesCoProcessor(
+ grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
+ pp->process(0);
+
+ SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
+ SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
+
+ SPtr<UbScheduler> timeBCSch(new UbScheduler(1, startTime, startTime));
+ auto timeDepBC = make_shared<TimeDependentBCCoProcessor>(TimeDependentBCCoProcessor(grid, timeBCSch));
+ timeDepBC->addInteractor(inflowInt);
+
+#ifdef _OPENMP
+ omp_set_num_threads(numOfThreads);
+#endif
+
+ SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
+ SPtr<Calculator> calculator(new BasicCalculator(grid, stepGhostLayer, endTime));
+ calculator->addCoProcessor(npr);
+ calculator->addCoProcessor(pp);
+ calculator->addCoProcessor(timeDepBC);
+ calculator->addCoProcessor(rcp);
+
+ 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[])
+int main(int argc, char *argv[])
{
- //Sleep(30000);
- if (argv != NULL)
- {
- if (argv[1] != NULL)
- {
- run(string(argv[1]));
- }
- else
- {
- cout << "Configuration file is missing!" << endl;
- }
- }
-
+ // Sleep(30000);
+ if (argv != NULL) {
+ if (argv[1] != NULL) {
+ run(string(argv[1]));
+ } else {
+ cout << "Configuration file is missing!" << endl;
+ }
+ }
}
-
diff --git a/apps/cpu/Multiphase/Multiphase.cfg b/apps/cpu/Multiphase/Multiphase.cfg
index c294ea68ce96c751030380d52d16eb35d06f9faa..b2f435db04ce51f915c3994b8418ba97b49c4843 100644
--- a/apps/cpu/Multiphase/Multiphase.cfg
+++ b/apps/cpu/Multiphase/Multiphase.cfg
@@ -1,11 +1,11 @@
-pathname = d:/temp/MultiphaseNew5
+pathname = d:/temp/JetBreakup
#pathGeo = d:/Projects/VirtualFluids-Multiphase/source/Applications/Multiphase/backup
pathGeo = d:/Projects/VirtualFluidsCombined/apps/cpu/Multiphase/backup
#geoFile = JetBreakupR.ASCII.stl
#geoFile = inlet1.stl
geoFile = tubeTransformed.stl
-numOfThreads = 4
+numOfThreads = 1
availMem = 10e9
#Grid
@@ -22,21 +22,23 @@ availMem = 10e9
#boundingBox = -40e-3 40e-3 1.0e-3 11.0e-3 -403-3 40e-3 #(Jet Breakup2) (Original without inlet length)
#blocknx = 20 20 20
-boundingBox = 6.0e-3 46.0e-3 -5e-3 5e-3 -5e-3 5e-3
-blocknx = 20 20 20
+#boundingBox = 6.0e-3 46.0e-3 -5e-3 5e-3 -5e-3 5e-3
+#blocknx = 20 20 20
+boundingBox = 0 9 0 9 0 9
+blocknx = 10 10 10
-dx = 1.66666666667e-4
+dx = 1 #1.66666666667e-4
refineLevel = 0
#Simulation
uLB = 0.005 #inlet velocity
#uF2 = 0.0001
Re = 10
-nuL =1e-2# 1.0e-5 #!1e-2
-nuG =1e-2# 1.16e-4 #!1e-2
+nuL =1e-3# 1.0e-5 #!1e-2
+nuG =1e-6# 1.16e-4 #!1e-2
densityRatio = 1000
-sigma = 1e-5 #4.66e-3 #surface tension 1e-4 ./. 1e-5
-interfaceThickness = 5
+sigma = 0 #1e-5 #4.66e-3 #surface tension 1e-4 ./. 1e-5
+interfaceWidth = 5
radius = 615.0 (Jet Breakup)
contactAngle = 110.0
gravity = 0.0
diff --git a/apps/cpu/Multiphase/Multiphase.cpp b/apps/cpu/Multiphase/Multiphase.cpp
index 77eb317e98ea019c983b3f150ca7f4c98d5a034a..9623cd730b7a80548a5ca14c149ae3583aba8c7b 100644
--- a/apps/cpu/Multiphase/Multiphase.cpp
+++ b/apps/cpu/Multiphase/Multiphase.cpp
@@ -28,7 +28,7 @@ void run(string configname)
double nuG = config.getValue<double>("nuG");
double densityRatio = config.getValue<double>("densityRatio");
double sigma = config.getValue<double>("sigma");
- int interfaceThickness = config.getValue<int>("interfaceThickness");
+ int interfaceWidth = config.getValue<int>("interfaceWidth");
//double radius = config.getValue<double>("radius");
double theta = config.getValue<double>("contactAngle");
double gr = config.getValue<double>("gravity");
@@ -49,8 +49,8 @@ void run(string configname)
double cpStep = config.getValue<double>("cpStep");
bool newStart = config.getValue<bool>("newStart");
- double beta = 12 * sigma / interfaceThickness;
- double kappa = 1.5 * interfaceThickness * sigma;
+ double beta = 12 * sigma / interfaceWidth;
+ double kappa = 1.5 * interfaceWidth * sigma;
SPtr<vf::mpi::Communicator> comm = vf::mpi::MPICommunicator::getInstance();
int myid = comm->getProcessID();
@@ -89,7 +89,8 @@ void run(string configname)
//kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
//kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsCumulantLBMKernel());
//kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel());
- kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
+ kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
kernel->setWithForcing(true);
kernel->setForcingX1(0.0);
@@ -107,6 +108,7 @@ void run(string configname)
kernel->setDensityRatio(densityRatio);
kernel->setMultiphaseModelParameters(beta, kappa);
kernel->setContactAngle(theta);
+ kernel->setInterfaceWidth(interfaceWidth);
SPtr<BCProcessor> bcProc(new BCProcessor());
// BCProcessorPtr bcProc(new ThinWallBCProcessor());
@@ -114,9 +116,9 @@ void run(string configname)
kernel->setBCProcessor(bcProc);
SPtr<Grid3D> grid(new Grid3D(comm));
- // grid->setPeriodicX1(true);
- // grid->setPeriodicX2(true);
- // grid->setPeriodicX3(true);
+ //grid->setPeriodicX1(true);
+ //grid->setPeriodicX2(true);
+ //grid->setPeriodicX3(true);
grid->setGhostLayerWidth(2);
@@ -347,8 +349,7 @@ void run(string configname)
UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
}
- MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, densityRatio, beta, kappa, theta, availMem,
- needMem);
+ MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, availMem, needMem);
grid->accept(kernelVisitor);
@@ -363,10 +364,52 @@ void run(string configname)
mu::Parser fct1;
fct1.SetExpr("phiL");
fct1.DefineConst("phiL", phiL);
- MultiphaseInitDistributionsBlockVisitor initVisitor(interfaceThickness);
+ //MultiphaseInitDistributionsBlockVisitor initVisitor(interfaceThickness);
+ MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
initVisitor.setPhi(fct1);
grid->accept(initVisitor);
-
+///////////////////////////////////////////////////////////////////////////////////////////
+ //{
+ // std::vector<std::vector<SPtr<Block3D>>> blockVector;
+ // int gridRank = comm->getProcessID();
+ // int minInitLevel = grid->getCoarsestInitializedLevel();
+ // int maxInitLevel = grid->getFinestInitializedLevel();
+ // blockVector.resize(maxInitLevel + 1);
+ // for (int level = minInitLevel; level <= maxInitLevel; level++) {
+ // grid->getBlocks(level, gridRank, true, blockVector[level]);
+ //}
+ // for (int level = minInitLevel; level <= maxInitLevel; level++) {
+ // for (SPtr<Block3D> block : blockVector[level]) {
+ // if (block) {
+ // int ix1 = block->getX1();
+ // int ix2 = block->getX2();
+ // int ix3 = block->getX3();
+ // int level = block->getLevel();
+
+ // for (int dir = 0; dir < D3Q27System::ENDDIR; dir++) {
+ // SPtr<Block3D> neighBlock = grid->getNeighborBlock(dir, ix1, ix2, ix3, level);
+
+ // if (!neighBlock) {
+
+ // }
+ // }
+ // }
+ // }
+ //}
+ // SPtr<Block3D> block = grid->getBlock(0, 0, 0, 0);
+ // SPtr<LBMKernel> kernel = dynamicPointerCast<LBMKernel>(block->getKernel());
+ // SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
+
+ // for (int ix3 = 0; ix3 <= 13; ix3++) {
+ // for (int ix2 = 0; ix2 <= 13; ix2++) {
+ // for (int ix1 = 0; ix1 <= 13; ix1++) {
+ // if (ix1 == 0 || ix2 == 0 || ix3 == 0 || ix1 == 13 || ix2 == 13 || ix3 == 13)
+ // bcArray->setUndefined(ix1, ix2, ix3);
+ // }
+ // }
+ // }
+ //}
+ ////////////////////////////////////////////////////////////////////////////////////////////
// boundary conditions grid
{
SPtr<UbScheduler> geoSch(new UbScheduler(1));
@@ -403,9 +446,10 @@ void run(string configname)
//ThreeDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
- grid->accept(bcVisitor);
+ grid->accept(bcVisitor);
- ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ //ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
grid->accept(setConnsVisitor);
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
diff --git a/apps/cpu/MultiphaseDropletTest/droplet.cpp b/apps/cpu/MultiphaseDropletTest/droplet.cpp
index 49a887e98d4c7f2b1a33ee474756e8de7e62a09b..29a2228bb7a6b7a6808e6e76162a9d3a20c23126 100644
--- a/apps/cpu/MultiphaseDropletTest/droplet.cpp
+++ b/apps/cpu/MultiphaseDropletTest/droplet.cpp
@@ -46,7 +46,7 @@ void run(string configname)
double cpStart = config.getValue<double>("cpStart");
double cpStep = config.getValue<double>("cpStep");
bool newStart = config.getValue<bool>("newStart");
- double rStep = config.getValue<double>("rStep");
+ //double rStep = config.getValue<double>("rStep");
SPtr<vf::mpi::Communicator> comm = vf::mpi::MPICommunicator::getInstance();
int myid = comm->getProcessID();
@@ -71,19 +71,19 @@ void run(string configname)
std::string fileName = "./LastTimeStep" + std::to_string((int)boundingBox[1]) + ".txt";
-#if defined(__unix__)
- double lastTimeStep = 0;
- //if (!newStart)
- {
- std::ifstream ifstr(fileName);
- ifstr >> lastTimeStep;
- restartStep = lastTimeStep;
- if(endTime >= lastTimeStep)
- endTime = lastTimeStep + rStep;
- else
- return;
- }
-#endif
+//#if defined(__unix__)
+// double lastTimeStep = 0;
+// //if (!newStart)
+// {
+// std::ifstream ifstr(fileName);
+// ifstr >> lastTimeStep;
+// restartStep = lastTimeStep;
+// if(endTime >= lastTimeStep)
+// endTime = lastTimeStep + rStep;
+// else
+// return;
+// }
+//#endif
//Sleep(30000);
@@ -271,8 +271,7 @@ void run(string configname)
UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
}
- MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, densityRatio, beta, kappa, theta, availMem,
- needMem);
+ MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, availMem, needMem);
grid->accept(kernelVisitor);
@@ -402,23 +401,23 @@ void run(string configname)
if (myid == 0)
UBLOG(logINFO, "Simulation-end");
-#if defined(__unix__)
- //if (!newStart)
- //{
- if (myid == 0)
- {
- std::ofstream ostr(fileName);
- ostr << endTime;
- cout << "start sbatch\n";
- //system("./start.sh");
- //system("echo test!");
- std::string str = "sbatch startJob" + std::to_string((int)boundingBox[1]) + ".sh";
- //system("sbatch startJob512.sh");
- system(str.c_str());
- }
- //MPI_Barrier((MPI_Comm)comm->getNativeCommunicator());
- //}
-#endif
+//#if defined(__unix__)
+// //if (!newStart)
+// //{
+// if (myid == 0)
+// {
+// std::ofstream ostr(fileName);
+// ostr << endTime;
+// cout << "start sbatch\n";
+// //system("./start.sh");
+// //system("echo test!");
+// std::string str = "sbatch startJob" + std::to_string((int)boundingBox[1]) + ".sh";
+// //system("sbatch startJob512.sh");
+// system(str.c_str());
+// }
+// //MPI_Barrier((MPI_Comm)comm->getNativeCommunicator());
+// //}
+//#endif
} catch (std::exception &e) {
cerr << e.what() << endl << flush;
diff --git a/apps/cpu/RisingBubble2D/RisingBubble2D.cfg b/apps/cpu/RisingBubble2D/RisingBubble2D.cfg
index 659a3a2d2e137f8bdb773271eb1e26fc21924816..d0635ea272199311a6e09a68e77ae7ca59a239f0 100644
--- a/apps/cpu/RisingBubble2D/RisingBubble2D.cfg
+++ b/apps/cpu/RisingBubble2D/RisingBubble2D.cfg
@@ -5,8 +5,12 @@ availMem = 10e9
#Grid
-boundingBox = 0 160 0 320 0 3
-blocknx = 16 16 3
+#boundingBox = 0 160 0 320 0 3
+#blocknx = 16 16 3
+#blocknx = 80 80 3
+
+boundingBox = 0 20 0 20 0 3
+blocknx = 20 20 3
dx = 1
refineLevel = 0
@@ -21,7 +25,7 @@ nuG = 1e-3
densityRatio = 10
sigma = 1.0850694444444444e-06 #1e-10 #1e-6 # 1e-5 #4.66e-3 #surface tension 1e-4 ./. 1e-5
interfaceThickness = 4.096
-radius = 40
+radius = 5 #40
contactAngle = 110.0
phi_L = 0.0
phi_H = 1.0
@@ -37,7 +41,7 @@ restartStep = 10
cpStart = 10
cpStep = 10
-outTime = 10000
-endTime = 20
+outTime = 100000
+endTime = 13
rStep = 159990 #160000
\ No newline at end of file
diff --git a/apps/cpu/RisingBubble2D/RisingBubble2D.cpp b/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
index 1c2b551ee5c9472bdef1438138ba98687fa9ba71..9001b64f00a30ee405797a5022ae68597e7c97ad 100644
--- a/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
+++ b/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
@@ -85,7 +85,7 @@ void run(string configname)
// }
//#endif
- //Sleep(30000);
+ //Sleep(20000);
// LBMReal dLB = 0; // = length[1] / dx;
LBMReal rhoLB = 0.0;
@@ -196,10 +196,14 @@ void run(string configname)
SPtr<UbScheduler> rSch(new UbScheduler(cpStep, cpStart));
SPtr<MPIIORestartCoProcessor> rcp(new MPIIORestartCoProcessor(grid, rSch, pathname, comm));
//SPtr<MPIIOMigrationCoProcessor> rcp(new MPIIOMigrationCoProcessor(grid, rSch, metisVisitor, pathname, comm));
+<<<<<<< HEAD
//SPtr<MPIIOMigrationBECoProcessor> rcp(new MPIIOMigrationBECoProcessor(grid, rSch, pathname, comm));
+=======
+ //SPtr<MPIIOMigrationBECoProcessor> rcp(new MPIIOMigrationBECoProcessor(grid, rSch, metisVisitor, pathname, comm));
+>>>>>>> 2d5b026e0450341bc450d1ab0085bd1940db01c7
// rcp->setNu(nuLB);
- // rcp->setNuLG(nuL, nuG);
- // rcp->setDensityRatio(densityRatio);
+ // rcp->setNuLG(nuL, nuG);
+ //rcp->setDensityRatio(densityRatio);
rcp->setLBMKernel(kernel);
rcp->setBCProcessor(bcProc);
@@ -283,8 +287,7 @@ void run(string configname)
UBLOG(logINFO, "Available memory per process = " << availMem << " bytes");
}
- MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, densityRatio, beta, kappa, theta, availMem,
- needMem);
+ MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, availMem, needMem);
grid->accept(kernelVisitor);
diff --git a/apps/cpu/ViskomatXL/viskomat.cfg b/apps/cpu/ViskomatXL/viskomat.cfg
index 75833bef60efc726a7d2e46fd11cdb386aa40104..3e84b68cf857a8fd43c41d0687717a24be6d7b90 100644
--- a/apps/cpu/ViskomatXL/viskomat.cfg
+++ b/apps/cpu/ViskomatXL/viskomat.cfg
@@ -1,33 +1,37 @@
-outputPath = d:/temp/viskomatXL
+outputPath = d:/temp/viskomatXL_restart_test
geoPath = d:/Projects/TRR277/Project/WP1/Rheometer/Aileen
#geoPath = d:/Projects/TRR277/Project/WP1/Rheometer
geoFile = fishboneT.stl
#geoFile = cylinder.stl
-numOfThreads = 4
+numOfThreads = 1
availMem = 15e9
logToFile = false
blocknx = 14 14 14
+#blocknx = 14 15 15
+#blocknx = 35 83 83
boundingBox = 0 140 -82.5 82.5 -82.5 82.5
+#blocknx = 32 12 12
+#boundingBox = 0 32 -12 12 -12 12
+#boundingBox = 0 64 -24 24 -24 24
+#boundingBox = 0 64 -24 24 -24 24
+
deltax = 1
+
refineLevel = 0
-#nuLB = 1.5e-4
OmegaLB = 1e-5
-N = 80 #rpm
-mu = 1 # Pa s
+mu = 5 # Pa s
+N = 80 # rpm
tau0 = 20e-7
-resolution = 32
-scaleFactor = 1
-
newStart = false
-restartStep = 50
+restartStep = 10
-cpStart = 50
-cpStep = 50
+cpStart = 10
+cpStep = 10
outTime = 10000
-endTime = 100 #0000
\ No newline at end of file
+endTime = 20
\ No newline at end of file
diff --git a/apps/cpu/ViskomatXL/viskomat.cpp b/apps/cpu/ViskomatXL/viskomat.cpp
index c29e8699d04cafaa654026ecba6a6e65c068240c..ad07752eda626eadab395f3dd2f890f4a109cbc8 100644
--- a/apps/cpu/ViskomatXL/viskomat.cpp
+++ b/apps/cpu/ViskomatXL/viskomat.cpp
@@ -111,16 +111,18 @@ void bflow(string configname)
//// rotation around X-axis
mu::Parser fctVy;
- fctVy.SetExpr("-Omega*(x3-z0-r)");
+ fctVy.SetExpr("-Omega*(x3-z0-r)/deltax");
fctVy.DefineConst("Omega", OmegaLB);
fctVy.DefineConst("r", 0.5 * (g_maxX3 - g_minX3));
fctVy.DefineConst("z0", g_minX3);
+ fctVy.DefineConst("deltax", deltax);
mu::Parser fctVz;
- fctVz.SetExpr("Omega*(x2-y0-r)");
+ fctVz.SetExpr("Omega*(x2-y0-r)/deltax");
fctVz.DefineConst("Omega", OmegaLB);
fctVz.DefineConst("r", 0.5 * (g_maxX2 - g_minX2));
fctVz.DefineConst("y0", g_minX2);
+ fctVz.DefineConst("deltax", deltax);
mu::Parser fctVx;
fctVx.SetExpr("0.0");
@@ -218,6 +220,8 @@ void bflow(string configname)
SPtr<GbTriFaceMesh3D> stator = make_shared<GbTriFaceMesh3D>();
stator->readMeshFromSTLFileBinary(geoPath + "/" + geoFile, false);
+ //stator->scale(2.0, 2.0, 2.0);
+ //stator->translate(16.0, 0.0, 0.0);
//stator->translate(4.0, -73.0, -6.0);
SPtr<D3Q27Interactor> statorInt = SPtr<D3Q27TriFaceMeshInteractor>(
@@ -360,7 +364,11 @@ void bflow(string configname)
else
{
restartCoProcessor->restart((int)restartStep);
- grid->setTimeStep(restartStep);
+
+ //restartCoProcessor->readBlocks((int)restartStep);
+ //restartCoProcessor->readDataSet((int)restartStep);
+ ////restartCoProcessor->readBoundaryConds((int)restartStep);
+ //grid->setTimeStep((int)restartStep);
SetBcBlocksBlockVisitor v2(wallXmaxInt);
grid->accept(v2);
@@ -377,6 +385,10 @@ void bflow(string configname)
SetBcBlocksBlockVisitor v1(wallXminInt);
grid->accept(v1);
wallXminInt->initInteractor();
+
+ SPtr<CoProcessor> ppblocks(new WriteBlocksCoProcessor(grid, SPtr<UbScheduler>(new UbScheduler(1)), outputPath,
+ WbWriterVtkXmlBinary::getInstance(), comm));
+ ppblocks->process(1);
}
omp_set_num_threads(numOfThreads);
diff --git a/apps/cpu/rheometer/rheometer.cfg b/apps/cpu/rheometer/rheometer.cfg
index 9eec8c6ded9b7a5ab8d1e6177c43354a4514ccc3..9b739bc67ed42d46c89adaefab1b020ad67da660 100644
--- a/apps/cpu/rheometer/rheometer.cfg
+++ b/apps/cpu/rheometer/rheometer.cfg
@@ -1,4 +1,4 @@
-#outputPath = d:/temp/rheometer/rheometerBinghamqQBB/rheometerBingham_tau_20e-7_nu_1.5e-3_new_lim_test
+#outputPath = d:/temp/rheometerTest
outputPath = d:/temp/Taylor-CouetteFlowIncompCum
viscosityPath = d:/Projects/VirtualFluidsCombined/apps/cpu/rheometer
@@ -8,7 +8,7 @@ logToFile = false
blocknx = 16 16 1 #8 8 1
#boundingBox = 32 32 1
-deltax = 1
+deltax = 0.5
#boundingBox = 0.02 0.02 0.00125
#deltax = 0.000625
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNoSlipBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNoSlipBCAlgorithm.cpp
index 51fc2d5abdfe7cfa5bafb7ae21571c684cd26b02..76128e0e72336ce6a056178e38708da5e888da92 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNoSlipBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseNoSlipBCAlgorithm.cpp
@@ -85,8 +85,12 @@ void MultiphaseNoSlipBCAlgorithm::applyBC()
//quadratic bounce back
const int invDir = D3Q27System::INVDIR[fdir];
LBMReal fReturn = f[invDir];
+ //if (UbMath::isNaN(fReturn))
+ //UBLOG(logINFO, "fReturn: " << fReturn);
distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
LBMReal hReturn = h[invDir];
+ //if (UbMath::isNaN(hReturn))
+ //UBLOG(logINFO, "hReturn: " << hReturn);
distributionsH->setDistributionForDirection(hReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
}
}
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseSlipBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseSlipBCAlgorithm.cpp
index beba9a256b869a37828efe44e886bb988bf9fa71..483cf3db8878a4b2b78e29327c2066bb8e79ec83 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseSlipBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseSlipBCAlgorithm.cpp
@@ -133,8 +133,8 @@ void MultiphaseSlipBCAlgorithm::applyBC()
LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity*rho)/(1.0+q));
distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
- LBMReal hReturn = ((1.0-q)/(1.0+q))*((h[invDir]-heq[invDir])/(1.0-collFactorPh)+heq[invDir])+((q/(1.0+q))*(h[invDir]+h[fdir]));
- //LBMReal hReturn = h[invDir];
+ //LBMReal hReturn = ((1.0-q)/(1.0+q))*((h[invDir]-heq[invDir])/(1.0-collFactorPh)+heq[invDir])+((q/(1.0+q))*(h[invDir]+h[fdir]));
+ LBMReal hReturn = h[invDir];
distributionsH->setDistributionForDirection(hReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
}
}
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
index 58c359887fe8b4263d8140038dc03754aeab74bb..70b0cceffbf3b0567cb3080c5861957a6ea2d2cb 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/MultiphaseVelocityBCAlgorithm.cpp
@@ -100,7 +100,7 @@ void MultiphaseVelocityBCAlgorithm::applyBC()
else if (bcPtr->hasVelocityBoundaryFlag(D3Q27System::S)) { nx2 += 1; }
else if (bcPtr->hasVelocityBoundaryFlag(D3Q27System::T)) { nx3 -= 1; }
else if (bcPtr->hasVelocityBoundaryFlag(D3Q27System::B)) { nx3 += 1; }
- else UB_THROW(UbException(UB_EXARGS, "Danger...no orthogonal BC-Flag on velocity boundary..."));
+ //else UB_THROW(UbException(UB_EXARGS, "Danger...no orthogonal BC-Flag on velocity boundary..."));
phiBC = bcPtr->getBoundaryPhaseField();
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
index a3ccfd51f966a4e89b5aba63c8e637f97a013eca..b2c7466f7cd6e7d5dd0aeb0baa152bfb6ced93ae 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
@@ -82,8 +82,6 @@ void CalculateTorqueCoProcessor::calculateForces()
torqueX2global = 0.0;
torqueX3global = 0.0;
- int counter = 0;
-
for(SPtr<D3Q27Interactor> interactor : interactors)
{
double x1Centre = interactor->getGbObject3D()->getX1Centroid();
@@ -103,17 +101,6 @@ void CalculateTorqueCoProcessor::calculateForces()
SPtr<ILBMKernel> kernel = block->getKernel();
- //if (kernel->getCompressible())
- //{
- // calcMacrosFct = &D3Q27System::calcCompMacroscopicValues;
- // compressibleFactor = 1.0;
- //}
- //else
- //{
- // calcMacrosFct = &D3Q27System::calcIncompMacroscopicValues;
- // compressibleFactor = 0.0;
- //}
-
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
distributions->swap();
@@ -151,29 +138,14 @@ void CalculateTorqueCoProcessor::calculateForces()
torqueX1 += ry * Fz - rz * Fy;
torqueX2 += rz * Fx - rx * Fz;
torqueX3 += rx * Fy - ry * Fx;
-
-
- counter++;
- //UBLOG(logINFO, "x1="<<(worldCoordinates[1] - x2Centre)<<",x2=" << (worldCoordinates[2] - x3Centre)<< ",x3=" << (worldCoordinates[0] - x1Centre) <<" forceX3 = " << forceX3);
}
}
- //if we have got discretization with more level
- // deltaX is LBM deltaX and equal LBM deltaT
- //double deltaX = 0.5; // LBMSystem::getDeltaT(block->getLevel()); //grid->getDeltaT(block);
- //double deltaXquadrat = deltaX*deltaX;
- //torqueX1 *= deltaXquadrat;
- //torqueX2 *= deltaXquadrat;
- //torqueX3 *= deltaXquadrat;
distributions->swap();
torqueX1global += torqueX1;
torqueX2global += torqueX2;
torqueX3global += torqueX3;
-
- //UBLOG(logINFO, "torqueX1global = " << torqueX1global);
-
- //UBLOG(logINFO, "counter = " << counter);
}
}
std::vector<double> values;
@@ -182,8 +154,6 @@ void CalculateTorqueCoProcessor::calculateForces()
values.push_back(torqueX2global);
values.push_back(torqueX3global);
- //UBLOG(logINFO, "counter = " << counter);
-
rvalues = comm->gather(values);
if (comm->getProcessID() == comm->getRoot())
{
@@ -206,9 +176,6 @@ UbTupleDouble3 CalculateTorqueCoProcessor::getForces(int x1, int x2, int x3, SP
LBMReal fs[D3Q27System::ENDF + 1];
distributions->getDistributionInv(fs, x1, x2, x3);
- //LBMReal /*rho = 0.0,*/ vx1 = 0.0, vx2 = 0.0, vx3 = 0.0, drho = 0.0;
- //calcMacrosFct(fs, drho, vx1, vx2, vx3);
- //rho = 1.0 + drho * compressibleFactor;
if(bc)
{
@@ -226,22 +193,9 @@ UbTupleDouble3 CalculateTorqueCoProcessor::getForces(int x1, int x2, int x3, SP
f = dynamicPointerCast<EsoTwist3D>(distributions)->getDistributionInvForDirection(x1, x2, x3, invDir);
fnbr = dynamicPointerCast<EsoTwist3D>(distributions)->getDistributionInvForDirection(x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
- // Vector3D boundaryVelocity;
- // boundaryVelocity[0] = bc->getBoundaryVelocityX1();
- // boundaryVelocity[1] = bc->getBoundaryVelocityX2();
- // boundaryVelocity[2] = bc->getBoundaryVelocityX3();
- //double correction[3] = { 0.0, 0.0, 0.0 };
- // if (bc->hasVelocityBoundaryFlag(fdir))
- // {
- // const double forceTerm = f - fnbr;
- // correction[0] = forceTerm * boundaryVelocity[0];
- // correction[1] = forceTerm * boundaryVelocity[1];
- // correction[2] = forceTerm * boundaryVelocity[2];
- // }
-
- forceX1 += (f + fnbr) * D3Q27System::DX1[invDir];// - 2.0 * D3Q27System::WEIGTH[invDir] * rho - correction[0];
- forceX2 += (f + fnbr) * D3Q27System::DX2[invDir];// - 2.0 * D3Q27System::WEIGTH[invDir] * rho - correction[1];
- forceX3 += (f + fnbr) * D3Q27System::DX3[invDir];// - 2.0 * D3Q27System::WEIGTH[invDir] * rho - correction[2];
+ forceX1 += (f + fnbr) * D3Q27System::DX1[invDir];
+ forceX2 += (f + fnbr) * D3Q27System::DX2[invDir];
+ forceX3 += (f + fnbr) * D3Q27System::DX3[invDir];
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h
index 872bccdc4c1098bb3d33c0efe6c996e96152439e..e488b442b60b2f726747a521e51cad9d4bacdbe9 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.h
@@ -43,10 +43,6 @@ private:
double torqueX1global;
double torqueX2global;
double torqueX3global;
-
- typedef void(*CalcMacrosFct)(const LBMReal* const& /*f[27]*/, LBMReal& /*rho*/, LBMReal& /*vx1*/, LBMReal& /*vx2*/, LBMReal& /*vx3*/);
- CalcMacrosFct calcMacrosFct;
- LBMReal compressibleFactor;
};
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/MPIIOCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/MPIIOCoProcessor.cpp
index 70bde9e97c2928f7c91a347653dbcfaff9c5e501..a16f32c7d9e0d83dff90a55bb139d4115285a196 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/MPIIOCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/MPIIOCoProcessor.cpp
@@ -55,7 +55,7 @@ MPIIOCoProcessor::MPIIOCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const
//-----------------------------------------------------------------------
MPI_Datatype typesBC[3] = { MPI_LONG_LONG_INT, MPI_FLOAT, MPI_CHAR };
- int blocksBC[3] = { 5, 33, 1 };
+ int blocksBC[3] = { 5, 34, 1 };
MPI_Aint offsetsBC[3], lbBC, extentBC;
offsetsBC[0] = 0;
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
index 1798bc68d26797833d371ae86435351c123b236a..b7e4d4b9a6cc3276728b5d07d89ac3c723ab027c 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/MPIIORestartCoProcessor.cpp
@@ -2500,10 +2500,12 @@ void MPIIORestartCoProcessor::readArray(int step, Arrays arrType, std::string fn
}
}
- MPI_File_read_at(file_handler, (MPI_Offset)(read_offset + sizeof(dataSetParam)), dataSetSmallArray, blocksCount, dataSetSmallType, MPI_STATUS_IGNORE);
+ MPI_File_read_at(file_handler, (MPI_Offset)(read_offset + sizeof(dataSetParam)), dataSetSmallArray, (int)blocksCount, dataSetSmallType, MPI_STATUS_IGNORE);
if (doubleCountInBlock > 0)
- MPI_File_read_at(file_handler, (MPI_Offset)(read_offset + sizeof(dataSetParam) + blocksCount * sizeof(DataSetSmallRestart)),
- &doubleValuesArray[0], blocksCount, dataSetDoubleType, MPI_STATUS_IGNORE);
+ MPI_File_read_at(
+ file_handler,
+ (MPI_Offset)(read_offset + sizeof(dataSetParam) + blocksCount * sizeof(DataSetSmallRestart)),
+ &doubleValuesArray[0], (int)blocksCount, dataSetDoubleType, MPI_STATUS_IGNORE);
MPI_File_close(&file_handler);
MPI_Type_free(&dataSetDoubleType);
diff --git a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
index 0dc57aabaacb6864d6b01e228716921bc5d10393..7372beed5261a6b933a274002a9cb173c718dff2 100644
--- a/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/MultiphasePressureFilterLBMKernel.cpp
@@ -47,8 +47,8 @@ MultiphasePressureFilterLBMKernel::MultiphasePressureFilterLBMKernel() { this->c
//////////////////////////////////////////////////////////////////////////
void MultiphasePressureFilterLBMKernel::initDataSet()
{
- SPtr<DistributionArray3D> f(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.9));
- SPtr<DistributionArray3D> h(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.9)); // For phase-field
+ SPtr<DistributionArray3D> f(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
+ SPtr<DistributionArray3D> h(new D3Q27EsoTwist3DSplittedVector( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0)); // For phase-field
//SPtr<PhaseFieldArray3D> divU1(new PhaseFieldArray3D( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr pressure(new CbArray3D<LBMReal, IndexerX3X2X1>( nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
@@ -58,7 +58,7 @@ void MultiphasePressureFilterLBMKernel::initDataSet()
//dataSet->setPhaseField(divU1);
dataSet->setPressureField(pressure);
- phaseField = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, -999.0));
+ phaseField = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
divU = CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr(new CbArray3D<LBMReal, IndexerX3X2X1>(nx[0] + 4, nx[1] + 4, nx[2] + 4, 0.0));
}
@@ -242,7 +242,7 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
mfcca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p);
mfbbb = (*this->zeroDistributionsF)(x1, x2, x3);
-
+
LBMReal rhoH = 1.0;
LBMReal rhoL = 1.0 / densityRatio;
@@ -300,7 +300,7 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
}
////!filter
- //bool firstTime = true;
+
for (int x3 = minX3; x3 < maxX3; x3++) {
for (int x2 = minX2; x2 < maxX2; x2++) {
for (int x1 = minX1; x1 < maxX1; x1++) {
@@ -476,6 +476,124 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
vvy += mu * dX2_phi * c1o2 / rho ;
vvz += mu * dX3_phi * c1o2 / rho;
+ //Abbas
+ LBMReal pStar = ((((((mfaaa + mfccc) + (mfaac + mfcca)) + ((mfcac + mfaca) + (mfcaa + mfacc)))
+ + (((mfaab + mfccb) + (mfacb + mfcab)) + ((mfaba + mfcbc) + (mfabc + mfcba)) + ((mfbaa + mfbcc) + (mfbac + mfbca))))
+ + ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb) * c1o3;
+
+ LBMReal M200 = ((((((mfaaa + mfccc) + (mfaac + mfcca)) + ((mfcac + mfaca) + (mfcaa + mfacc)))
+ + (((mfaab + mfccb) + (mfacb + mfcab)) + ((mfaba + mfcbc) + (mfabc + mfcba))))
+ + ((mfabb + mfcbb))));
+ LBMReal M020 = ((((((mfaaa + mfccc) + (mfaac + mfcca)) + ((mfcac + mfaca) + (mfcaa + mfacc)))
+ + (((mfaab + mfccb) + (mfacb + mfcab)) + ((mfbaa + mfbcc) + (mfbac + mfbca))))
+ + ((mfbab + mfbcb))));
+ LBMReal M002 = ((((((mfaaa + mfccc) + (mfaac + mfcca)) + ((mfcac + mfaca) + (mfcaa + mfacc)))
+ + (+((mfaba + mfcbc) + (mfabc + mfcba)) + ((mfbaa + mfbcc) + (mfbac + mfbca))))
+ + ((mfbba + mfbbc))));
+
+ LBMReal M110 = ((((((mfaaa + mfccc) + (-mfcac - mfaca)) + ((mfaac + mfcca) + (-mfcaa - mfacc)))
+ + (((mfaab + mfccb) + (-mfacb - mfcab))))
+ ));
+ LBMReal M101 = ((((((mfaaa + mfccc) - (mfaac + mfcca)) + ((mfcac + mfaca) - (mfcaa + mfacc)))
+ + (((mfaba + mfcbc) + (-mfabc - mfcba))))
+ ));
+ LBMReal M011 = ((((((mfaaa + mfccc) - (mfaac + mfcca)) + ((mfcaa + mfacc) - (mfcac + mfaca)))
+ + (((mfbaa + mfbcc) + (-mfbac - mfbca))))
+ ));
+ LBMReal vvxI = vvx;
+ LBMReal vvyI = vvy;
+ LBMReal vvzI = vvz;
+
+ LBMReal collFactorStore = collFactorM;
+ LBMReal stress;
+ for (int iter = 0; iter < 1; iter++) {
+ LBMReal OxxPyyPzz = 1.0;
+ LBMReal mxxPyyPzz = (M200 - vvxI * vvxI) + (M020 - vvyI * vvyI) + (M002 - vvzI * vvzI);
+ mxxPyyPzz -= c3 * pStar;
+
+ LBMReal mxxMyy = (M200 - vvxI * vvxI) - (M020 - vvyI * vvyI);
+ LBMReal mxxMzz = (M200 - vvxI * vvxI) - (M002 - vvzI * vvzI);
+ LBMReal mxy = M110 - vvxI * vvyI;
+ LBMReal mxz = M101 - vvxI * vvzI;
+ LBMReal myz = M011 - vvyI * vvzI;
+
+ ///////Bingham
+ //LBMReal dxux = -c1o2 * collFactorM * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (/*mfaaa*/ -mxxPyyPzz);
+ //LBMReal dyuy = dxux + collFactorM * c3o2 * mxxMyy;
+ //LBMReal dzuz = dxux + collFactorM * c3o2 * mxxMzz;
+ //LBMReal Dxy = -three * collFactorM * mxy;
+ //LBMReal Dxz = -three * collFactorM * mxz;
+ //LBMReal Dyz = -three * collFactorM * myz;
+
+ //LBMReal tau0 = phi[REST] * 1.0e-7;//(phi[REST]>0.01)?1.0e-6: 0;
+ //LBMReal shearRate =fabs(pStar)*0.0e-2+ sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho);
+ //collFactorM = collFactorM * (UbMath::one - (collFactorM * tau0) / (shearRate * c1o3 /* *rho*/ + 1.0e-15));
+ //collFactorM = (collFactorM < -1000000) ? -1000000 : collFactorM;
+ ////if(collFactorM < 0.1) {
+ //// int test = 1;
+ ////}
+ //////!Bingham
+
+
+ mxxMyy *= c1 - collFactorM * c1o2;
+ mxxMzz *= c1 - collFactorM * c1o2;
+ mxy *= c1 - collFactorM * c1o2;
+ mxz *= c1 - collFactorM * c1o2;
+ myz *= c1 - collFactorM * c1o2;
+ mxxPyyPzz *= c1 - OxxPyyPzz * c1o2;
+ //mxxPyyPzz += c3 * pStar;
+ LBMReal mxx = (mxxMyy + mxxMzz + mxxPyyPzz) * c1o3;
+ LBMReal myy = (-c2 * mxxMyy + mxxMzz + mxxPyyPzz) * c1o3;
+ LBMReal mzz = (mxxMyy - c2 * mxxMzz + mxxPyyPzz) * c1o3;
+ vvxI = vvx - (mxx * dX1_phi + mxy * dX2_phi + mxz * dX3_phi) * rhoToPhi / (rho);
+ vvyI = vvy - (mxy * dX1_phi + myy * dX2_phi + myz * dX3_phi) * rhoToPhi / (rho);
+ vvzI = vvz - (mxz * dX1_phi + myz * dX2_phi + mzz * dX3_phi) * rhoToPhi / (rho);
+
+
+
+ }
+
+
+ forcingX1 += c2 * (vvxI - vvx);
+ forcingX2 += c2 * (vvyI - vvy);
+ forcingX3 += c2 * (vvzI - vvz);
+
+ mfabb += c1o2 * (-forcingX1) * c2o9;
+ mfbab += c1o2 * (-forcingX2) * c2o9;
+ mfbba += c1o2 * (-forcingX3) * c2o9;
+ mfaab += c1o2 * (-forcingX1 - forcingX2) * c1o18;
+ mfcab += c1o2 * (forcingX1 - forcingX2) * c1o18;
+ mfaba += c1o2 * (-forcingX1 - forcingX3) * c1o18;
+ mfcba += c1o2 * (forcingX1 - forcingX3) * c1o18;
+ mfbaa += c1o2 * (-forcingX2 - forcingX3) * c1o18;
+ mfbca += c1o2 * (forcingX2 - forcingX3) * c1o18;
+ mfaaa += c1o2 * (-forcingX1 - forcingX2 - forcingX3) * c1o72;
+ mfcaa += c1o2 * (forcingX1 - forcingX2 - forcingX3) * c1o72;
+ mfaca += c1o2 * (-forcingX1 + forcingX2 - forcingX3) * c1o72;
+ mfcca += c1o2 * (forcingX1 + forcingX2 - forcingX3) * c1o72;
+ mfcbb += c1o2 * (forcingX1)*c2o9;
+ mfbcb += c1o2 * (forcingX2)*c2o9;
+ mfbbc += c1o2 * (forcingX3)*c2o9;
+ mfccb += c1o2 * (forcingX1 + forcingX2) * c1o18;
+ mfacb += c1o2 * (-forcingX1 + forcingX2) * c1o18;
+ mfcbc += c1o2 * (forcingX1 + forcingX3) * c1o18;
+ mfabc += c1o2 * (-forcingX1 + forcingX3) * c1o18;
+ mfbcc += c1o2 * (forcingX2 + forcingX3) * c1o18;
+ mfbac += c1o2 * (-forcingX2 + forcingX3) * c1o18;
+ mfccc += c1o2 * (forcingX1 + forcingX2 + forcingX3) * c1o72;
+ mfacc += c1o2 * (-forcingX1 + forcingX2 + forcingX3) * c1o72;
+ mfcac += c1o2 * (forcingX1 - forcingX2 + forcingX3) * c1o72;
+ mfaac += c1o2 * (-forcingX1 - forcingX2 + forcingX3) * c1o72;
+
+
+
+ vvx = vvxI;
+ vvy = vvyI;
+ vvz = vvzI;
+
+ //!Abbas
+
+
LBMReal vx2;
LBMReal vy2;
LBMReal vz2;
@@ -758,7 +876,7 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
LBMReal OxyyPxzz = 8.0 * (collFactorM - 2.0) * (OxxPyyPzz * (3.0 * collFactorM - 1.0) - 5.0 * collFactorM) / (8.0 * (5.0 - 2.0 * collFactorM) * collFactorM + OxxPyyPzz * (8.0 + collFactorM * (9.0 * collFactorM - 26.0)));
LBMReal OxyyMxzz = 8.0 * (collFactorM - 2.0) * (collFactorM + OxxPyyPzz * (3.0 * collFactorM - 7.0)) / (OxxPyyPzz * (56.0 - 42.0 * collFactorM + 9.0 * collFactorM * collFactorM) - 8.0 * collFactorM);
- // LBMReal Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
+ LBMReal Oxyz = 24.0 * (collFactorM - 2.0) * (4.0 * collFactorM * collFactorM + collFactorM * OxxPyyPzz * (18.0 - 13.0 * collFactorM) + OxxPyyPzz * OxxPyyPzz * (2.0 + collFactorM * (6.0 * collFactorM - 11.0))) / (16.0 * collFactorM * collFactorM * (collFactorM - 6.0) - 2.0 * collFactorM * OxxPyyPzz * (216.0 + 5.0 * collFactorM * (9.0 * collFactorM - 46.0)) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (3.0 * collFactorM - 10.0) * (15.0 * collFactorM - 28.0) - 48.0));
LBMReal A = (4.0 * collFactorM * collFactorM + 2.0 * collFactorM * OxxPyyPzz * (collFactorM - 6.0) + OxxPyyPzz * OxxPyyPzz * (collFactorM * (10.0 - 3.0 * collFactorM) - 4.0)) / ((collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
//FIXME: warning C4459: declaration of 'B' hides global declaration (message : see declaration of 'D3Q27System::B' )
LBMReal BB = (4.0 * collFactorM * OxxPyyPzz * (9.0 * collFactorM - 16.0) - 4.0 * collFactorM * collFactorM - 2.0 * OxxPyyPzz * OxxPyyPzz * (2.0 + 9.0 * collFactorM * (collFactorM - 2.0))) / (3.0 * (collFactorM - OxxPyyPzz) * (OxxPyyPzz * (2.0 + 3.0 * collFactorM) - 8.0 * collFactorM));
@@ -848,7 +966,7 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
LBMReal mxyyMxzz = mfbca - mfbac;
//relax
- wadjust = OxyyMxzz + (1. - OxyyMxzz) * fabs(mfbbb) / (fabs(mfbbb) + qudricLimit);
+ wadjust = Oxyz + (1. - Oxyz) * fabs(mfbbb) / (fabs(mfbbb) + qudricLimit);
mfbbb += wadjust * (-mfbbb);
wadjust = OxyyPxzz + (1. - OxyyPxzz) * fabs(mxxyPyzz) / (fabs(mxxyPyzz) + qudricLimit);
mxxyPyzz += wadjust * (-mxxyPyzz);
@@ -925,13 +1043,13 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
////////////////////////////////////////////////////////////////////////////////////
//forcing
- mfbaa = -mfbaa;
- mfaba = -mfaba;
- mfaab = -mfaab;
+ //mfbaa = -mfbaa;
+ //mfaba = -mfaba;
+ //mfaab = -mfaab;
//////////////////////////////////////////////////////////////////////////////////////
- mfbaa += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (2 * dxux * dX1_phi + Dxy * dX2_phi + Dxz * dX3_phi) / (rho);
- mfaba += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (Dxy * dX1_phi + 2 * dyuy * dX2_phi + Dyz * dX3_phi) / (rho);
- mfaab += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (Dxz * dX1_phi + Dyz * dX2_phi + 2 * dyuy * dX3_phi) / (rho);
+ //mfbaa += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (2 * dxux * dX1_phi + Dxy * dX2_phi + Dxz * dX3_phi) / (rho);
+ //mfaba += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (Dxy * dX1_phi + 2 * dyuy * dX2_phi + Dyz * dX3_phi) / (rho);
+ //mfaab += c1o3 * (c1 / collFactorM - c1o2) * rhoToPhi * (Dxz * dX1_phi + Dyz * dX2_phi + 2 * dyuy * dX3_phi) / (rho);
////////////////////////////////////////////////////////////////////////////////////
//back
////////////////////////////////////////////////////////////////////////////////////
@@ -1141,6 +1259,38 @@ void MultiphasePressureFilterLBMKernel::calculate(int step)
mfbcc = m1;
mfccc = m2;
+ ////forcing
+
+ mfabb += c1o2 * (-forcingX1) * c2o9;
+ mfbab += c1o2 * (-forcingX2) * c2o9;
+ mfbba += c1o2 * (-forcingX3) * c2o9;
+ mfaab += c1o2 * (-forcingX1 - forcingX2) * c1o18;
+ mfcab += c1o2 * (forcingX1 - forcingX2) * c1o18;
+ mfaba += c1o2 * (-forcingX1 - forcingX3) * c1o18;
+ mfcba += c1o2 * (forcingX1 - forcingX3) * c1o18;
+ mfbaa += c1o2 * (-forcingX2 - forcingX3) * c1o18;
+ mfbca += c1o2 * (forcingX2 - forcingX3) * c1o18;
+ mfaaa += c1o2 * (-forcingX1 - forcingX2 - forcingX3) * c1o72;
+ mfcaa += c1o2 * (forcingX1 - forcingX2 - forcingX3) * c1o72;
+ mfaca += c1o2 * (-forcingX1 + forcingX2 - forcingX3) * c1o72;
+ mfcca += c1o2 * (forcingX1 + forcingX2 - forcingX3) * c1o72;
+ mfcbb += c1o2 * (forcingX1)*c2o9;
+ mfbcb += c1o2 * (forcingX2)*c2o9;
+ mfbbc += c1o2 * (forcingX3)*c2o9;
+ mfccb += c1o2 * (forcingX1 + forcingX2) * c1o18;
+ mfacb += c1o2 * (-forcingX1 + forcingX2) * c1o18;
+ mfcbc += c1o2 * (forcingX1 + forcingX3) * c1o18;
+ mfabc += c1o2 * (-forcingX1 + forcingX3) * c1o18;
+ mfbcc += c1o2 * (forcingX2 + forcingX3) * c1o18;
+ mfbac += c1o2 * (-forcingX2 + forcingX3) * c1o18;
+ mfccc += c1o2 * (forcingX1 + forcingX2 + forcingX3) * c1o72;
+ mfacc += c1o2 * (-forcingX1 + forcingX2 + forcingX3) * c1o72;
+ mfcac += c1o2 * (forcingX1 - forcingX2 + forcingX3) * c1o72;
+ mfaac += c1o2 * (-forcingX1 - forcingX2 + forcingX3) * c1o72;
+
+
+
+
//////////////////////////////////////////////////////////////////////////
//proof correctness
//////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/Utilities/CheckpointConverter.h b/src/cpu/VirtualFluidsCore/Utilities/CheckpointConverter.h
index 6d5b3c686dbbc2d688df1b59dedba9c696c2959d..6fe24772d574a6db67428a820027971b4c7fd230 100644
--- a/src/cpu/VirtualFluidsCore/Utilities/CheckpointConverter.h
+++ b/src/cpu/VirtualFluidsCore/Utilities/CheckpointConverter.h
@@ -32,7 +32,7 @@ protected:
private:
MPI_Datatype gridParamType, block3dType;
MPI_Datatype dataSetParamType, dataSetTypeRead, dataSetTypeWrite;
- MPI_Datatype boundCondType, boundCondType1000, dataSetDoubleType;
+ MPI_Datatype boundCondType, boundCondType1000;
MPIIODataStructures::boundCondParam boundCondParamStr;
};
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
index 11863b6912dee88a161db75ba81db44c5ec586de..afe288e8c70a7ab3b0ceadae721b68db6a0102f4 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.cpp
@@ -19,9 +19,8 @@
// }
//}
//////////////////////////////////////////////////////////////////////////
-MultiphaseSetKernelBlockVisitor::MultiphaseSetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nuL, LBMReal nuG, LBMReal densityRatio, LBMReal beta, LBMReal kappa,
- LBMReal contactAngle, double availMem, double needMem, MultiphaseSetKernelBlockVisitor::Action action /*= SetKernelBlockVisitor::New*/) :
- Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(kernel), nuL(nuL), nuG(nuG), densityRatio(densityRatio), beta(beta), kappa(kappa), contactAngle(contactAngle), action(action), dataSetFlag(true)
+MultiphaseSetKernelBlockVisitor::MultiphaseSetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nuL, LBMReal nuG, double availMem, double needMem, MultiphaseSetKernelBlockVisitor::Action action /*= SetKernelBlockVisitor::New*/) :
+ Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(kernel), nuL(nuL), nuG(nuG), action(action), dataSetFlag(true)
{
if (needMem > availMem)
{
@@ -36,10 +35,6 @@ void MultiphaseSetKernelBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> blo
LBMReal collFactorL = LBMSystem::calcCollisionFactor(nuL, block->getLevel());
LBMReal collFactorG = LBMSystem::calcCollisionFactor(nuG, block->getLevel());
kernel->setCollisionFactorMultiphase(collFactorL, collFactorG);
-
- //kernel->setDensityRatio(densityRatio);
- //kernel->setMultiphaseModelParameters(beta, kappa);
- //kernel->setContactAngle(contactAngle);
kernel->setIndex(block->getX1(), block->getX2(), block->getX3());
kernel->setDeltaT(LBMSystem::getDeltaT(block->getLevel()));
diff --git a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.h b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.h
index 9f952e2dddaf67fed5fa10cd8874938e584b5114..24d2b35c3a85b80e793b94d61feceb58b607ff19 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/MultiphaseSetKernelBlockVisitor.h
@@ -42,11 +42,7 @@ class MultiphaseSetKernelBlockVisitor : public Block3DVisitor
public:
enum Action { NewKernel, ChangeKernel, ChangeKernelWithData};
public:
- //SetKernelBlockVisitor(LBMKernel3DPtr kernel, LBMReal nue);
-
- //SetKernelBlockVisitor(LBMKernel3DPtr kernel, LBMReal nue, double availMem, double needMem);
-
- MultiphaseSetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nuL, LBMReal nuG, LBMReal densityRatio, LBMReal beta, LBMReal kappa, LBMReal contactAngle, double availMem, double needMem,
+ MultiphaseSetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nuL, LBMReal nuG, double availMem, double needMem,
MultiphaseSetKernelBlockVisitor::Action action = MultiphaseSetKernelBlockVisitor::NewKernel);
virtual ~MultiphaseSetKernelBlockVisitor() {}
@@ -59,10 +55,6 @@ private:
SPtr<LBMKernel> kernel;
LBMReal nuL;
LBMReal nuG;
- LBMReal densityRatio;
- LBMReal beta;
- LBMReal kappa;
- LBMReal contactAngle;
Action action;
bool dataSetFlag;
};
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
index 42ed287c6d9ddfe30e38dce3a05814c385cdefd0..f39a95afc23e5351e0331f789a0e22ae5ff5ccf2 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.h
@@ -70,8 +70,7 @@ protected:
void setSameLevelConnectors(SPtr<Grid3D> grid, SPtr<Block3D> block);
void setRemoteConnectors(SPtr<Block3D> sblock, SPtr<Block3D> tblock, int dir);
std::shared_ptr<vf::mpi::Communicator> comm;
- int dirs;
- int gridRank;
+ int gridRank{0};
};
template <class T1, class T2>