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..426891d2cd7fee2dde95796a0e3e3f7a7a0b5483 100644
--- a/apps/cpu/JetBreakup/JetBreakup.cfg
+++ b/apps/cpu/JetBreakup/JetBreakup.cfg
@@ -1,39 +1,34 @@
-pathname = d:/temp/Multiphase
-pathGeo = d:/Projects/VirtualFluids-Multiphase/source/Applications/Multiphase/backup
-geoFile = JetBreakup2.ASCII.stl
+pathname = d:/temp/JetBreakupBenchmark_D50
+#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
+uLB = 0.01 #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)
+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
+deltaXfactor = 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 +39,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..262eafca939fab9d95272758c1bea7abe7c43e5f 100644
--- a/apps/cpu/JetBreakup/JetBreakup.cpp
+++ b/apps/cpu/JetBreakup/JetBreakup.cpp
@@ -1,516 +1,582 @@
#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 uLB = config.getValue<double>("uLB");
+ // 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 deltaXfactor = config.getValue<double>("deltaXfactor");
+ 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");
+ double dx = D/deltaXfactor;
+ 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");
+
+ double beta = 12 * sigma / interfaceWidth;
+ double kappa = 1.5 * interfaceWidth * sigma;
+
+ 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();
-
+ if (myid == 0) {
+ const char *str = pathname.c_str();
+ mkdir(str, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+ }
+#endif
- ppblocks->process(0);
- ppblocks.reset();
+ 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;
+ }
+
+ // LBMReal dLB = 0; // = length[1] / dx;
+ LBMReal rhoLB = 0.0;
+ LBMReal nuLB = nuL; //(uLB*dLB) / Re;
+
+ 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(nuL, nuG);
+ kernel->setDensityRatio(densityRatio);
+ 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", uLB);
+
+ 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, "uLb = " << uLB);
+ UBLOG(logINFO, "rho = " << rhoLB);
+ UBLOG(logINFO, "nuLb = " << nuLB);
+ UBLOG(logINFO, "Re = " << Re);
+ UBLOG(logINFO, "dx = " << dx);
+ UBLOG(logINFO, "Preprocess - start");
+ }
- 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());
+ 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");
+ }
- 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++)
+ MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, 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
{
- int nobl = grid->getNumberOfBlocks(level);
- UBLOG(logINFO, "Number of blocks for level " << level << " = " << nobl);
- UBLOG(logINFO, "Number of nodes for level " << level << " = " << nobl * numberOfNodesPerBlock);
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ SPtr<WriteBoundaryConditionsCoProcessor> ppgeo(new WriteBoundaryConditionsCoProcessor(
+ grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppgeo->process(0);
+ ppgeo.reset();
}
- 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();
-
-
- 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)
- {
- UBLOG(logINFO, "Parameters:");
- UBLOG(logINFO, "uLb = " << uLB);
- UBLOG(logINFO, "rho = " << rhoLB);
- UBLOG(logINFO, "nuLb = " << nuLB);
- UBLOG(logINFO, "Re = " << Re);
- UBLOG(logINFO, "dx = " << dx);
- 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);
-
-
-
-
-
-
- //UbSchedulerPtr bcSch(new UbScheduler(1, 12000, 12000));
- //TimeDependentBCCoProcessorPtr inflowF2 (new TimeDependentBCCoProcessor(grid,bcSch));
- //inflowF2->addInteractor(inflowF2_1Int);
- //inflowF2->addInteractor(inflowF2_2Int);
-
- //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;
- }
+ if (myid == 0)
+ UBLOG(logINFO, "Preprocess - end");
+ } else {
+ if (myid == 0) {
+ UBLOG(logINFO, "Parameters:");
+ UBLOG(logINFO, "uLb = " << uLB);
+ UBLOG(logINFO, "rho = " << rhoLB);
+ UBLOG(logINFO, "nuLb = " << nuLB);
+ UBLOG(logINFO, "Re = " << Re);
+ UBLOG(logINFO, "dx = " << dx);
+ //UBLOG(logINFO, "number of levels = " << refineLevel + 1);
+ UBLOG(logINFO, "numOfThreads = " << numOfThreads);
+ UBLOG(logINFO, "path = " << pathname);
+ }
+ 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/(uLB/(deltaXfactor)));
+ 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/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 6eb283ae25524714a8b4680b937038ecab01fffb..0182f569cec34dc2f993803e274dd853e194cfc2 100644
--- a/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
+++ b/apps/cpu/RisingBubble2D/RisingBubble2D.cpp
@@ -174,7 +174,7 @@ void run(string configname)
SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseNoSlipBCAlgorithm()));
SPtr<BCAdapter> slipBCAdapter(new SlipBCAdapter());
- noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseSlipBCAlgorithm()));
+ slipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new MultiphaseSlipBCAlgorithm()));
//////////////////////////////////////////////////////////////////////////////////
// BC visitor
MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
@@ -194,55 +194,58 @@ void run(string configname)
//////////////////////////////////////////////////////////////////////////
// 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);
+ 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, metisVisitor, pathname, comm));
+ //rcp->setNu(nuLB);
+ //rcp->setNuLG(nuL, nuG);
+ //rcp->setDensityRatio(densityRatio);
rcp->setLBMKernel(kernel);
rcp->setBCProcessor(bcProc);
//////////////////////////////////////////////////////////////////////////
+ // bounding box
+ 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];
+
+ double dx2 = 2.0 * dx;
+ GbCuboid3DPtr wallXmin(
+ new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_minX1, g_maxX2 + dx2, g_maxX3 + dx2));
+ GbSystem3D::writeGeoObject(wallXmin.get(), pathname + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallXmax(
+ new GbCuboid3D(g_maxX1, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
+ GbSystem3D::writeGeoObject(wallXmax.get(), pathname + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
+
+ GbCuboid3DPtr wallYmin(
+ new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_minX2, g_maxX3 + dx2));
+ GbSystem3D::writeGeoObject(wallYmin.get(), pathname + "/geo/wallYmin", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallYmax(
+ new GbCuboid3D(g_minX1 - dx2, g_maxX2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
+ GbSystem3D::writeGeoObject(wallYmax.get(), pathname + "/geo/wallYmax", WbWriterVtkXmlASCII::getInstance());
+
+ SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, slipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, slipBCAdapter, Interactor3D::SOLID));
+
+ SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
if (newStart) {
- // bounding box
- 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
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());
-
-
+ WbWriterVtkXmlBinary::getInstance());
GenBlocksGridVisitor genBlocks(gridCube);
grid->accept(genBlocks);
- double dx2 = 2.0 * dx;
- GbCuboid3DPtr wallXmin(new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_minX1, g_maxX2 + dx2, g_maxX3 + dx2));
- GbSystem3D::writeGeoObject(wallXmin.get(), pathname + "/geo/wallXmin", WbWriterVtkXmlASCII::getInstance());
- GbCuboid3DPtr wallXmax(new GbCuboid3D(g_maxX1, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
- GbSystem3D::writeGeoObject(wallXmax.get(), pathname + "/geo/wallXmax", WbWriterVtkXmlASCII::getInstance());
-
- GbCuboid3DPtr wallYmin(new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_minX2, g_maxX3 + dx2));
- GbSystem3D::writeGeoObject(wallYmin.get(), pathname + "/geo/wallYmin", WbWriterVtkXmlASCII::getInstance());
- GbCuboid3DPtr wallYmax(new GbCuboid3D(g_minX1 - dx2, g_maxX2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
- GbSystem3D::writeGeoObject(wallYmax.get(), pathname + "/geo/wallYmax", WbWriterVtkXmlASCII::getInstance());
-
- SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
SPtr<WriteBlocksCoProcessor> ppblocks(new WriteBlocksCoProcessor(
grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
@@ -349,7 +352,28 @@ void run(string configname)
}
rcp->restart((int)restartStep);
- grid->setTimeStep(restartStep);
+ //grid->setTimeStep(restartStep);
+
+ //rcp->readBlocks((int)restartStep);
+ //rcp->readDataSet((int)restartStep);
+ //rcp->readBoundaryConds((int)restartStep);
+ //grid->setTimeStep((int)restartStep);
+
+ //SetBcBlocksBlockVisitor v2(wallXminInt);
+ //grid->accept(v2);
+ //wallXminInt->initInteractor();
+
+ //SetBcBlocksBlockVisitor v3(wallXmaxInt);
+ //grid->accept(v3);
+ //wallXmaxInt->initInteractor();
+
+ //SetBcBlocksBlockVisitor v4(wallYminInt);
+ //grid->accept(v4);
+ //wallYminInt->initInteractor();
+
+ //SetBcBlocksBlockVisitor v1(wallYmaxInt);
+ //grid->accept(v1);
+ //wallYmaxInt->initInteractor();
if (myid == 0)
UBLOG(logINFO, "Restart - end");
@@ -392,8 +416,9 @@ void run(string configname)
SPtr<WriteMultiphaseQuantitiesCoProcessor> pp(new WriteMultiphaseQuantitiesCoProcessor(
grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
- if(grid->getTimeStep() == 0)
- pp->process(0);
+ //if(grid->getTimeStep() == 0)
+ //pp->process(0);
+ pp->process(grid->getTimeStep());
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
SPtr<NUPSCounterCoProcessor> npr(new NUPSCounterCoProcessor(grid, nupsSch, numOfThreads, comm));
diff --git a/apps/cpu/ViskomatXL/viskomat.cfg b/apps/cpu/ViskomatXL/viskomat.cfg
index c46760239722b590e6f08fc1dd6505168780c84c..3e84b68cf857a8fd43c41d0687717a24be6d7b90 100644
--- a/apps/cpu/ViskomatXL/viskomat.cfg
+++ b/apps/cpu/ViskomatXL/viskomat.cfg
@@ -1,24 +1,24 @@
-outputPath = d:/temp/viskomatXL_dx_0.5
-#geoPath = d:/Projects/TRR277/Project/WP1/Rheometer/Aileen
-geoPath = d:/Projects/TRR277/Project/WP1/Rheometer
-#geoFile = fishboneT.stl
-geoFile = cylinder.stl
+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 = 1
availMem = 15e9
logToFile = false
-#blocknx = 14 14 14
+blocknx = 14 14 14
#blocknx = 14 15 15
#blocknx = 35 83 83
-#boundingBox = 0 140 -82.5 82.5 -82.5 82.5
+boundingBox = 0 140 -82.5 82.5 -82.5 82.5
-blocknx = 32 12 12
-boundingBox = 0 32 -12 12 -12 12
+#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 = 0.5
+deltax = 1
refineLevel = 0
@@ -27,11 +27,11 @@ mu = 5 # Pa s
N = 80 # rpm
tau0 = 20e-7
-newStart = true
-restartStep = 1.3e6
+newStart = false
+restartStep = 10
-cpStart = 100000
-cpStep = 100000
+cpStart = 10
+cpStep = 10
-outTime = 1
-endTime = 200
\ No newline at end of file
+outTime = 10000
+endTime = 20
\ No newline at end of file
diff --git a/apps/cpu/ViskomatXL/viskomat.cpp b/apps/cpu/ViskomatXL/viskomat.cpp
index bef8fa53a92ef3aac3ac7a055ab1778b03c12ef4..45694a4e453245d96ca9a57970c591e2cec66283 100644
--- a/apps/cpu/ViskomatXL/viskomat.cpp
+++ b/apps/cpu/ViskomatXL/viskomat.cpp
@@ -363,8 +363,12 @@ void bflow(string configname)
}
else
{
- restartCoProcessor->restart((int)restartStep);
- grid->setTimeStep(restartStep);
+ //restartCoProcessor->restart((int)restartStep);
+
+ restartCoProcessor->readBlocks((int)restartStep);
+ restartCoProcessor->readDataSet((int)restartStep);
+ //restartCoProcessor->readBoundaryConds((int)restartStep);
+ grid->setTimeStep((int)restartStep);
SetBcBlocksBlockVisitor v2(wallXmaxInt);
grid->accept(v2);
@@ -381,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);
@@ -405,7 +413,7 @@ void bflow(string configname)
SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
//writeMQCoProcessor->process(100);
- SPtr<UbScheduler> forceSch(new UbScheduler(1000));
+ SPtr<UbScheduler> forceSch(new UbScheduler(100));
SPtr<CalculateTorqueCoProcessor> fp = make_shared<CalculateTorqueCoProcessor>(grid, forceSch, outputPath + "/torque/TorqueRotor.csv", comm);
fp->addInteractor(rotorInt);
SPtr<CalculateTorqueCoProcessor> fp2 = make_shared<CalculateTorqueCoProcessor>(grid, forceSch, outputPath + "/torque/TorqueStator.csv", comm);
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/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;