diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 399ea713523aade56ebf3bff38cf4be31320d59c..f44fed21bd9abfd9657b716a93505aa3146b0d77 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -25,6 +25,7 @@ IF(${VFCPU_ENABLE_MultiphaseFlow})
add_subdirectory(${APPS_ROOT_CPU}/ShotcreteJet)
add_subdirectory(${APPS_ROOT_CPU}/ConcreteExtrusion)
add_subdirectory(${APPS_ROOT_CPU}/MultiphaseSymmetryTest)
+ add_subdirectory(${APPS_ROOT_CPU}/DamBreak)
ENDIF()
IF(${VFCPU_ENABLE_LiggghtsCoupling} AND ${VFCPU_ENABLE_MultiphaseFlow})
diff --git a/apps/cpu/DamBreak/CMakeLists.txt b/apps/cpu/DamBreak/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..7fd74ee6cc69cc9a31ac0fab99d0b01883a34f32
--- /dev/null
+++ b/apps/cpu/DamBreak/CMakeLists.txt
@@ -0,0 +1,3 @@
+PROJECT(DamBreak)
+
+vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} MultiphaseFlow NonNewtonianFluids FILES main.cpp )
diff --git a/apps/cpu/DamBreak/config.cfg b/apps/cpu/DamBreak/config.cfg
new file mode 100644
index 0000000000000000000000000000000000000000..016e34072c1ff69e284cc1743ad684cb1a382e82
--- /dev/null
+++ b/apps/cpu/DamBreak/config.cfg
@@ -0,0 +1,44 @@
+#pathname = d:/temp/MultiphaseDropletTest
+pathname = E:/Multiphase/DropletTest_Test
+
+numOfThreads = 4
+availMem = 10e9
+
+#Grid
+
+boundingBox = 0 256 512 768 0 3
+blocknx = 16 16 3
+
+dx = 1
+refineLevel = 0
+
+#Simulation
+uLB = 0 #0.001#0.005#0.005
+Re = 10
+nuL = 1e-2 #1e-5# 1.0e-5 #!1e-2
+nuG = 0.015811388300841892 #5e-2 #1e-4 # 1e-8 # 1.16e-4 #!1e-2
+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 = 25.6
+contactAngle = 110.0
+#gravity = 0.0
+gravity = -1.0348028606838648e-08 #-5.04e-6
+phi_L = 0.0
+phi_H = 1.0
+Phase-field Relaxation = 0.6
+Mobility = 0.056 # 0.01 ./. 0.08, fine correction of Phase-field Relaxation parameter, to activate it need to change in kernel tauH to tauH1
+
+
+logToFile = false
+
+newStart = true
+restartStep = 100000
+
+cpStart = 1000
+cpStep = 1000
+
+outTime = 100
+endTime = 10000
+
+rStep = 159990 #160000
\ No newline at end of file
diff --git a/apps/cpu/DamBreak/main.cpp b/apps/cpu/DamBreak/main.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..54db263d163ece870074be82f3b4e865a336a44f
--- /dev/null
+++ b/apps/cpu/DamBreak/main.cpp
@@ -0,0 +1,475 @@
+#include <iostream>
+#include <string>
+#include <memory>
+
+#if defined(__unix__)
+#include <stdio.h>
+#include <stdlib.h>
+#endif
+
+#include "VirtualFluids.h"
+#include "MultiphaseFlow/MultiphaseFlow.h"
+#include "NonNewtonianFluids/NonNewtonianFluids.h"
+
+using namespace std;
+
+void run(string configname)
+{
+ using namespace vf::lbm::dir;
+
+ try {
+ vf::basics::ConfigurationFile config;
+ config.load(configname);
+
+ string pathname = config.getValue<string>("pathname");
+ int numOfThreads = config.getValue<int>("numOfThreads");
+ vector<int> blocknx = config.getVector<int>("blocknx");
+ vector<real> boundingBox = config.getVector<real>("boundingBox");
+ real uLB = config.getValue<real>("uLB");
+ real nuL = config.getValue<real>("nuL");
+ real nuG = config.getValue<real>("nuG");
+ real densityRatio = config.getValue<real>("densityRatio");
+ real sigma = config.getValue<real>("sigma");
+ real interfaceThickness = config.getValue<real>("interfaceThickness");
+ real radius = config.getValue<real>("radius");
+ real theta = config.getValue<real>("contactAngle");
+ //double gr = config.getValue<double>("gravity");
+ real phiL = config.getValue<real>("phi_L");
+ real phiH = config.getValue<real>("phi_H");
+ real tauH = config.getValue<real>("Phase-field Relaxation");
+ real mob = config.getValue<real>("Mobility");
+
+ real endTime = config.getValue<real>("endTime");
+ real outTime = config.getValue<real>("outTime");
+ real availMem = config.getValue<real>("availMem");
+ int refineLevel = config.getValue<int>("refineLevel");
+ real Re = config.getValue<real>("Re");
+ real dx = config.getValue<real>("dx");
+ bool logToFile = config.getValue<bool>("logToFile");
+ real restartStep = config.getValue<real>("restartStep");
+ real cpStart = config.getValue<real>("cpStart");
+ real cpStep = config.getValue<real>("cpStep");
+ bool newStart = config.getValue<bool>("newStart");
+ //double rStep = config.getValue<double>("rStep");
+
+ SPtr<vf::parallel::Communicator> comm = vf::parallel::MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ if (myid == 0)
+ UBLOG(logINFO, "Droplet Test: 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());
+ }
+ }
+
+ 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
+
+ //Sleep(30000);
+
+ // LBMReal dLB = 0; // = length[1] / dx;
+ real rhoLB = 0.0;
+ real nuLB = nuL; //(uLB*dLB) / Re;
+
+ //diameter of circular droplet
+ real D = 2.0*radius;
+
+ //density retio
+ real r_rho = densityRatio;
+
+ //density of heavy fluid
+ real rho_h = 1.0;
+ //density of light fluid
+ real rho_l = rho_h / r_rho;
+
+ //kinimatic viscosity
+ real nu_h = nuL;
+ //LBMReal nu_l = nuG;
+ //#dynamic viscosity
+ real mu_h = rho_h * nu_h;
+
+ //gravity
+ real g_y = Re* Re* mu_h* mu_h / (rho_h * (rho_h - rho_l) * D * D * D);
+ //Eotvos number
+ real Eo = 10; // 0.100;
+ //surface tension
+ sigma =0*0.001;
+ //0 * rho_h *g_y *D *D / Eo;
+
+ //g_y = 0;
+
+ real beta = 12.0 * sigma / interfaceThickness;
+ real kappa = 1.5 * interfaceThickness * sigma;
+
+ if (myid == 0) {
+ //UBLOG(logINFO, "uLb = " << uLB);
+ //UBLOG(logINFO, "rho = " << rhoLB);
+ UBLOG(logINFO, "D = " << D);
+ UBLOG(logINFO, "nuL = " << nuL);
+ UBLOG(logINFO, "nuG = " << nuG);
+ UBLOG(logINFO, "Re = " << Re);
+ UBLOG(logINFO, "Eo = " << Eo);
+ UBLOG(logINFO, "g_y = " << g_y);
+ UBLOG(logINFO, "sigma = " << sigma);
+ UBLOG(logINFO, "dx = " << dx);
+ UBLOG(logINFO, "Preprocess - start");
+ }
+
+ SPtr<LBMUnitConverter> conv(new LBMUnitConverter());
+
+ //const int baseLevel = 0;
+
+
+ SPtr<Rheology> thix = Rheology::getInstance();
+ thix->setYieldStress(0*1e-1); //(1.0e-2);
+
+ SPtr<LBMKernel> kernel;
+
+ //kernel = SPtr<LBMKernel>(new MultiphaseScratchCumulantLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
+ //kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
+ //kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
+ kernel = make_shared< MultiphaseScaleDistributionLBMKernel>();
+ //kernel = SPtr<LBMKernel>(new MultiphaseSharpInterfaceLBMKernel());
+ // kernel = make_shared<MultiphaseSharpInterfaceLBMKernel>();
+
+ mu::Parser fgr;
+ fgr.SetExpr("-(rho-rho_l)*g_y");
+ fgr.DefineConst("rho_l", rho_l);
+ fgr.DefineConst("g_y", g_y);
+
+ kernel->setWithForcing(true);
+ kernel->setForcingX1(0.0);
+ kernel->setForcingX2(fgr);
+ kernel->setForcingX3(0.0);
+
+ kernel->setPhiL(phiL);
+ kernel->setPhiH(phiH);
+ kernel->setPhaseFieldRelaxation(tauH);
+ kernel->setMobility(mob);
+ kernel->setInterfaceWidth(interfaceThickness);
+
+
+ kernel->setCollisionFactorMultiphase(nuL, nuG);
+ kernel->setDensityRatio(densityRatio);
+ kernel->setMultiphaseModelParameters(beta, kappa);
+ kernel->setContactAngle(theta);
+ kernel->setSigma(sigma);
+
+ SPtr<BCSet> bcProc(new BCSet());
+ // BCSetPtr bcProc(new ThinWallBCSet());
+
+ kernel->setBCSet(bcProc);
+
+ SPtr<BC> noSlipBC(new NoSlipBC());
+ noSlipBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseNoSlipBCStrategy()));
+
+ SPtr<BC> outflowBC(new DensityBC(rhoLB));
+ outflowBC->setBCStrategy(SPtr<BCStrategy>(new MultiphasePressureBCStrategy()));
+
+ //////////////////////////////////////////////////////////////////////////////////
+ // BC visitor
+ MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
+ bcVisitor.addBC(noSlipBC);
+ bcVisitor.addBC(outflowBC);
+
+ SPtr<Grid3D> grid(new Grid3D(comm));
+ grid->setDeltaX(dx);
+ grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
+ grid->setPeriodicX1(false);//true);
+ grid->setPeriodicX2(false);
+ grid->setPeriodicX3(true);
+ grid->setGhostLayerWidth(2);
+
+ SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, DIR_MMM, MetisPartitioner::RECURSIVE));
+
+ //////////////////////////////////////////////////////////////////////////
+ // restart
+ SPtr<UbScheduler> rSch(new UbScheduler(cpStep, cpStart));
+ //SPtr<MPIIORestartSimulationObserver> rcp(new MPIIORestartSimulationObserver(grid, rSch, pathname, comm));
+ SPtr<MPIIOMigrationSimulationObserver> rcp(new MPIIOMigrationSimulationObserver(grid, rSch, metisVisitor, pathname, comm));
+ //SPtr<MPIIOMigrationBESimulationObserver> rcp(new MPIIOMigrationBESimulationObserver(grid, rSch, pathname, comm));
+ // rcp->setNu(nuLB);
+ // rcp->setNuLG(nuL, nuG);
+ // rcp->setDensityRatio(densityRatio);
+
+ rcp->setLBMKernel(kernel);
+ rcp->setBCSet(bcProc);
+ //////////////////////////////////////////////////////////////////////////
+
+ if (newStart) {
+
+ // bounding box
+ real g_minX1 = boundingBox[0];
+ real g_minX2 = boundingBox[2];
+ real g_minX3 = boundingBox[4];
+
+ real g_maxX1 = boundingBox[1];
+ real g_maxX2 = boundingBox[3];
+ real 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());
+
+
+
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ real dx2 = 2.0 * dx;
+ 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> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBC, Interactor3D::SOLID));
+
+ //walls in x direction
+
+ 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());
+
+ SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, noSlipBC, Interactor3D::SOLID));
+
+ SPtr<WriteBlocksSimulationObserver> ppblocks(new WriteBlocksSimulationObserver(
+ grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+
+ InteractorsHelper intHelper(grid, metisVisitor, true);
+ intHelper.addInteractor(wallYminInt);
+ intHelper.addInteractor(wallYmaxInt);
+ intHelper.addInteractor(wallXminInt);
+ intHelper.addInteractor(wallXmaxInt);
+ intHelper.selectBlocks();
+
+ ppblocks->update(0);
+ ppblocks.reset();
+
+ unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
+ int ghostLayer = 5;
+ 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);
+ real needMemAll =
+ real(numberOfNodesPerBlockWithGhostLayer * (27 * sizeof(real) + sizeof(int) + sizeof(float) * 4));
+ real needMem = needMemAll / real(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");
+ }
+
+ MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nuL, nuG, availMem, needMem);
+
+ grid->accept(kernelVisitor);
+
+ if (refineLevel > 0) {
+ SetUndefinedNodesBlockVisitor undefNodesVisitor;
+ grid->accept(undefNodesVisitor);
+ }
+
+
+ intHelper.setBC();
+
+ // initialization of distributions
+ real x1c = 2.5 * D; // (g_maxX1 - g_minX1-1)/2; //
+ real x2c = 12.5 * D; //(g_maxX2 - g_minX2-1)/2;
+ real x3c = 1.5; // 2.5 * D;
+ //1.5; // 2.5 * D; //(g_maxX3 - g_minX3-1)/2;
+ //LBMReal x3c = 2.5 * D;
+ mu::Parser fct1;
+ //fct1.SetExpr("0.5-0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+0*(x3-x3c)^2)-radius)/interfaceThickness)");
+ //fct1.SetExpr("(x1-x1c<-100 && x3>3 ? 1: 0)");
+ fct1.SetExpr("(x1-x1c<-100 ? 1: 0)");
+ // fct1.SetExpr("(((x1-x1c)*(x1-x1c+1)<180 ? 1: 0)*(x2>530? 1:0)*(sin(x2*0.3)>-0.1?1:0)-0.5*0)");
+ // *(x2 > 530 ? 1 : 0) ");
+ fct1.DefineConst("x1c", x1c);
+ fct1.DefineConst("x2c", x2c);
+ fct1.DefineConst("x3c", x3c);
+ fct1.DefineConst("radius", radius/2);
+ fct1.DefineConst("interfaceThickness", interfaceThickness);
+
+ mu::Parser fct2;
+ fct2.SetExpr("0.5*uLB-uLB*0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+0*(x3-x3c)^2)-radius)/interfaceThickness)");
+ //fct2.SetExpr("uLB");
+ fct2.DefineConst("uLB", uLB);
+ fct2.DefineConst("x1c", x1c);
+ fct2.DefineConst("x2c", x2c);
+ fct2.DefineConst("x3c", x3c);
+ fct2.DefineConst("radius", radius);
+ fct2.DefineConst("interfaceThickness", interfaceThickness);
+
+ //MultiphaseInitDistributionsBlockVisitor initVisitor(densityRatio);
+ MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
+ initVisitor.setPhi(fct1);
+ initVisitor.setVx2(fct2);
+ grid->accept(initVisitor);
+
+ // boundary conditions grid
+ {
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ SPtr<WriteBoundaryConditionsSimulationObserver> ppgeo(new WriteBoundaryConditionsSimulationObserver(
+ grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppgeo->update(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);
+
+ if (myid == 0)
+ UBLOG(logINFO, "Restart - end");
+ }
+
+ grid->accept(bcVisitor);
+
+ //TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
+ //grid->accept(setConnsVisitor);
+
+ //ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ //grid->accept(setConnsVisitor);
+
+ TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ grid->accept(setConnsVisitor);
+
+ SPtr<UbScheduler> visSch(new UbScheduler(outTime));
+ real t_ast, t;
+ t_ast = 2;
+ t = (int)(t_ast/std::sqrt(g_y/D));
+ visSch->addSchedule(t,t,t); //t=2
+ t_ast = 3;
+ t = (int)(t_ast/std::sqrt(g_y/D));
+ visSch->addSchedule(t,t,t); //t=3
+ t_ast = 4;
+ t = (int)(t_ast/std::sqrt(g_y/D));
+ visSch->addSchedule(t,t,t); //t=4
+ t_ast = 5;
+ t = (int)(t_ast/std::sqrt(g_y/D));
+ visSch->addSchedule(t,t,t); //t=5
+ t_ast = 6;
+ t = (int)(t_ast/std::sqrt(g_y/D));
+ visSch->addSchedule(t,t,t); //t=6
+ t_ast = 7;
+ t = (int)(t_ast/std::sqrt(g_y/D));
+ visSch->addSchedule(t,t,t); //t=7
+ t_ast = 9;
+ t = (int)(t_ast/std::sqrt(g_y/D));
+ visSch->addSchedule(t,t,t); //t=9
+
+ SPtr<WriteSharpInterfaceQuantitiesSimulationObserver> pp(new WriteSharpInterfaceQuantitiesSimulationObserver(
+ grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
+ //if(grid->getTimeStep() == 0)
+ pp->update(0);
+
+ SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
+ SPtr<NUPSCounterSimulationObserver> npr(new NUPSCounterSimulationObserver(grid, nupsSch, numOfThreads, comm));
+
+ omp_set_num_threads(numOfThreads);
+
+ SPtr<UbScheduler> stepGhostLayer(new UbScheduler(1));
+ SPtr<Simulation> simulation(new Simulation(grid, stepGhostLayer, endTime));
+ simulation->addSimulationObserver(npr);
+ simulation->addSimulationObserver(pp);
+ simulation->addSimulationObserver(rcp);
+
+
+ if (myid == 0)
+ UBLOG(logINFO, "Simulation-start");
+ simulation->run();
+ 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
+
+ } 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[])
+{
+ // 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/MultiphaseDropletTest/CMakeLists.txt b/apps/cpu/MultiphaseDropletTest/CMakeLists.txt
index c48e517aaad9ff75f6245101ae0a433a2bd49886..e21c6ccad6b9bbc87e1cf8971d3f6d344b466501 100644
--- a/apps/cpu/MultiphaseDropletTest/CMakeLists.txt
+++ b/apps/cpu/MultiphaseDropletTest/CMakeLists.txt
@@ -1,3 +1,3 @@
PROJECT(MultiphaseDropletTest)
-vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} MultiphaseFlow NonNewtonianFluids FILES droplet.cpp )
+vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} MultiphaseFlow NonNewtonianFluids FILES droplet.cpp )
diff --git a/apps/cpu/MultiphaseDropletTest/droplet.cpp b/apps/cpu/MultiphaseDropletTest/droplet.cpp
index 1ae1d18996dfd0e3ae7854e349631ed141cbbc58..6a8b41189ba46ed4edbf6fa4523a1ed1afcda5f2 100644
--- a/apps/cpu/MultiphaseDropletTest/droplet.cpp
+++ b/apps/cpu/MultiphaseDropletTest/droplet.cpp
@@ -30,7 +30,7 @@ void run(string configname)
real nuG = config.getValue<real>("nuG");
real densityRatio = config.getValue<real>("densityRatio");
real sigma = config.getValue<real>("sigma");
- int interfaceThickness = config.getValue<int>("interfaceThickness");
+ real interfaceThickness = config.getValue<real>("interfaceThickness");
real radius = config.getValue<real>("radius");
real theta = config.getValue<real>("contactAngle");
//double gr = config.getValue<double>("gravity");
@@ -115,9 +115,10 @@ void run(string configname)
//gravity
real g_y = Re* Re* mu_h* mu_h / (rho_h * (rho_h - rho_l) * D * D * D);
//Eotvos number
- real Eo = 100;
+ real Eo = 10; // 0.100;
//surface tension
- sigma = rho_h* g_y* D* D / Eo;
+ sigma =0*0.001;
+ //0 * rho_h *g_y *D *D / Eo;
//g_y = 0;
@@ -144,7 +145,7 @@ void run(string configname)
SPtr<Rheology> thix = Rheology::getInstance();
- thix->setYieldStress(0);
+ thix->setYieldStress(0*1e-1); //(1.0e-2);
SPtr<LBMKernel> kernel;
@@ -152,9 +153,9 @@ void run(string configname)
// kernel = SPtr<LBMKernel>(new MultiphaseCumulantLBMKernel());
//kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
//kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
- //kernel = make_shared< MultiphaseScaleDistributionLBMKernel>();
+ kernel = make_shared< MultiphaseScaleDistributionLBMKernel>();
//kernel = SPtr<LBMKernel>(new MultiphaseSharpInterfaceLBMKernel());
- kernel = make_shared<MultiphaseSharpInterfaceLBMKernel>();
+ // kernel = make_shared<MultiphaseSharpInterfaceLBMKernel>();
mu::Parser fgr;
fgr.SetExpr("-(rho-rho_l)*g_y");
@@ -177,6 +178,7 @@ void run(string configname)
kernel->setDensityRatio(densityRatio);
kernel->setMultiphaseModelParameters(beta, kappa);
kernel->setContactAngle(theta);
+ kernel->setSigma(sigma);
SPtr<BCSet> bcProc(new BCSet());
// BCSetPtr bcProc(new ThinWallBCSet());
@@ -198,7 +200,7 @@ void run(string configname)
SPtr<Grid3D> grid(new Grid3D(comm));
grid->setDeltaX(dx);
grid->setBlockNX(blocknx[0], blocknx[1], blocknx[2]);
- grid->setPeriodicX1(true);
+ grid->setPeriodicX1(false);//true);
grid->setPeriodicX2(false);
grid->setPeriodicX3(true);
grid->setGhostLayerWidth(2);
@@ -250,12 +252,24 @@ void run(string configname)
SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBC, Interactor3D::SOLID));
SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBC, Interactor3D::SOLID));
+ //walls in x direction
+
+ 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());
+
+ SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, noSlipBC, Interactor3D::SOLID));
+
SPtr<WriteBlocksSimulationObserver> ppblocks(new WriteBlocksSimulationObserver(
grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
InteractorsHelper intHelper(grid, metisVisitor, true);
intHelper.addInteractor(wallYminInt);
intHelper.addInteractor(wallYmaxInt);
+ intHelper.addInteractor(wallXminInt);
+ intHelper.addInteractor(wallXmaxInt);
intHelper.selectBlocks();
ppblocks->update(0);
@@ -300,21 +314,25 @@ void run(string configname)
intHelper.setBC();
// initialization of distributions
- real x1c = (g_maxX1 + g_minX1)/2; //
- real x2c = (g_maxX2 + g_minX2)/2;
- real x3c = (g_maxX3 + g_minX3) / 2;
+ real x1c = 2.5 * D; // (g_maxX1 - g_minX1-1)/2; //
+ real x2c = 12.5 * D; //(g_maxX2 - g_minX2-1)/2;
+ real x3c = 1.5; // 2.5 * D;
//1.5; // 2.5 * D; //(g_maxX3 - g_minX3-1)/2;
//LBMReal x3c = 2.5 * D;
mu::Parser fct1;
- fct1.SetExpr("0.5-0.5*tanh(2*(sqrt(0*(x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
+ //fct1.SetExpr("0.5-0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+0*(x3-x3c)^2)-radius)/interfaceThickness)");
+ //fct1.SetExpr("(x1-x1c<-100 && x3>3 ? 1: 0)");
+ fct1.SetExpr("(x1-x1c<-100 ? 1: 0)");
+ // fct1.SetExpr("(((x1-x1c)*(x1-x1c+1)<180 ? 1: 0)*(x2>530? 1:0)*(sin(x2*0.3)>-0.1?1:0)-0.5*0)");
+ // *(x2 > 530 ? 1 : 0) ");
fct1.DefineConst("x1c", x1c);
fct1.DefineConst("x2c", x2c);
fct1.DefineConst("x3c", x3c);
- fct1.DefineConst("radius", radius);
+ fct1.DefineConst("radius", radius/2);
fct1.DefineConst("interfaceThickness", interfaceThickness);
mu::Parser fct2;
- fct2.SetExpr("0.5*uLB-uLB*0.5*tanh(2*(sqrt(0*(x1-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
+ fct2.SetExpr("0.5*uLB-uLB*0.5*tanh(2*(sqrt((x1-x1c)^2+(x2-x2c)^2+0*(x3-x3c)^2)-radius)/interfaceThickness)");
//fct2.SetExpr("uLB");
fct2.DefineConst("uLB", uLB);
fct2.DefineConst("x1c", x1c);
diff --git a/src/cpu/MultiphaseFlow/BoundaryConditions/MultiphaseNoSlipBCStrategy.cpp b/src/cpu/MultiphaseFlow/BoundaryConditions/MultiphaseNoSlipBCStrategy.cpp
index d438f80f53fcefd553813b0c9517e0fe545af96f..ef2c7ac21d0ffb25172cba9730d2d1bc3ab527dc 100644
--- a/src/cpu/MultiphaseFlow/BoundaryConditions/MultiphaseNoSlipBCStrategy.cpp
+++ b/src/cpu/MultiphaseFlow/BoundaryConditions/MultiphaseNoSlipBCStrategy.cpp
@@ -88,11 +88,11 @@ void MultiphaseNoSlipBCStrategy::applyBC()
//quadratic bounce back
const int invDir = D3Q27System::INVDIR[fdir];
real fReturn = f[invDir];
- //distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
- distributions->setDistributionForDirection(fReturn, x1, x2, x3, invDir);//delay BB
+ distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
+ //distributions->setDistributionForDirection(fReturn, x1, x2, x3, invDir);//delay BB
real hReturn = h[invDir];
- // distributionsH->setDistributionForDirection(hReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
- distributionsH->setDistributionForDirection(hReturn, x1, x2, x3, invDir);//delay BB
+ distributionsH->setDistributionForDirection(hReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
+ //distributionsH->setDistributionForDirection(hReturn, x1, x2, x3, invDir);//delay BB
if (distributionsH2)
{
real h2Return = h2[invDir];
diff --git a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
index 84613be5c0932bdde23bcf35002ec3df9d79d6e4..e791f0e7e55a1b7c9cdb741017a47f331ba0a1fc 100644
--- a/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
+++ b/src/cpu/MultiphaseFlow/LBM/MultiphaseScaleDistributionLBMKernel.cpp
@@ -190,6 +190,61 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
int maxX1 = bcArrayMaxX1 - ghostLayerWidth;
int maxX2 = bcArrayMaxX2 - ghostLayerWidth;
int maxX3 = bcArrayMaxX3 - ghostLayerWidth;
+
+
+ //// 08.10.23 post collision BC from new liquid nodes
+ // this->swapDistributions();
+ // {
+ // real fff[27];
+ // SPtr<DistributionArray3D> distributionC = this->getDataSet()->getFdistributions();
+ // //for (int x3 = minX3 - 1; x3 < maxX3 + 1; x3++) {
+ // // for (int x2 = minX2 - 1; x2 < maxX2 + 1; x2++) {
+ // // for (int x1 = minX1 - 1; x1 < maxX1 + 1; x1++) {
+ // // if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3))
+ // for (int x3 = minX3 - ghostLayerWidth + 1; x3 < maxX3 + ghostLayerWidth - 1; x3++) {
+ // for (int x2 = minX2 - ghostLayerWidth + 1; x2 < maxX2 + ghostLayerWidth - 1; x2++) {
+ // for (int x1 = minX1 - ghostLayerWidth + 1; x1 < maxX1 + ghostLayerWidth - 1; x1++) {
+ // if (!bcArray->isSolid(x1, x2, x3) /* && !bcArray->isUndefined(x1, x2, x3)*/)
+
+ // {
+
+ // findNeighbors(phaseFieldOld, x1, x2, x3);
+ // findNeighbors2(phaseField, x1, x2, x3);
+ // if ((phi[DIR_000] <= phiLim) && (phi2[DIR_000] > phiLim)) {
+ // for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
+ // if ((phi2[fdir] > phiLim) &&
+ // ( !bcArray->isSolid(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir])
+ // &&
+ // !bcArray->isUndefined(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir]))) {
+ // distributionC->getDistributionInv(fff, x1 + D3Q27System::DX1[fdir],
+ // x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir]);
+ // real vx, vy, vz, rho;
+ // D3Q27System::calcIncompMacroscopicValues(fff, rho, vx, vy, vz);
+ // real feq = D3Q27System::getIncompFeqForDirection(fdir, rho, vx, vy, vz);
+ // // real feqZero=D3Q27System::getIncompFeqForDirection(fdir, 0.0,vx,vy,vz);
+ // // real feqZeroI = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], 0.0,
+ // // vx, vy, vz); distribution->setDistributionForDirection(
+ // // c1o2 * (feqZero+feqZeroI-feq+
+ // // distribution->getDistributionInvForDirection(
+ // // x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 +
+ // // D3Q27System::DX3[fdir],fdir)),
+ // // x1, x2, x3,
+ // // fdir);
+ // distributionC->setDistributionForDirection(feq, x1, x2, x3, fdir);
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ // }
+ // this->swapDistributions();
+ // //! 08.10.23
+
//real omegaDRho = 1.0;// 1.25;// 1.3;
for (int x3 = minX3 - ghostLayerWidth; x3 < maxX3 + ghostLayerWidth; x3++) {
for (int x2 = minX2 - ghostLayerWidth; x2 < maxX2 + ghostLayerWidth; x2++) {
@@ -252,6 +307,68 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
}
}
+// // 19.09.23 Velocity from Death Reckoning
+// for (int x3 = minX3 - ghostLayerWidth; x3 < maxX3 + ghostLayerWidth; x3++) {
+// for (int x2 = minX2 - ghostLayerWidth; x2 < maxX2 + ghostLayerWidth; x2++) {
+// for (int x1 = minX1 - ghostLayerWidth; x1 < maxX1 + ghostLayerWidth; x1++) {
+// if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3)) {
+// int x1p = x1 + 1;
+// int x2p = x2 + 1;
+// int x3p = x3 + 1;
+//
+// real mfcbb = (*this->localDistributionsH1)(D3Q27System::ET_E, x1, x2, x3);
+// real mfbcb = (*this->localDistributionsH1)(D3Q27System::ET_N, x1, x2, x3);
+// real mfbbc = (*this->localDistributionsH1)(D3Q27System::ET_T, x1, x2, x3);
+// real mfccb = (*this->localDistributionsH1)(D3Q27System::ET_NE, x1, x2, x3);
+// real mfacb = (*this->localDistributionsH1)(D3Q27System::ET_NW, x1p, x2, x3);
+// real mfcbc = (*this->localDistributionsH1)(D3Q27System::ET_TE, x1, x2, x3);
+// real mfabc = (*this->localDistributionsH1)(D3Q27System::ET_TW, x1p, x2, x3);
+// real mfbcc = (*this->localDistributionsH1)(D3Q27System::ET_TN, x1, x2, x3);
+// real mfbac = (*this->localDistributionsH1)(D3Q27System::ET_TS, x1, x2p, x3);
+// real mfccc = (*this->localDistributionsH1)(D3Q27System::ET_TNE, x1, x2, x3);
+// real mfacc = (*this->localDistributionsH1)(D3Q27System::ET_TNW, x1p, x2, x3);
+// real mfcac = (*this->localDistributionsH1)(D3Q27System::ET_TSE, x1, x2p, x3);
+// real mfaac = (*this->localDistributionsH1)(D3Q27System::ET_TSW, x1p, x2p, x3);
+// real mfabb = (*this->nonLocalDistributionsH1)(D3Q27System::ET_W, x1p, x2, x3);
+// real mfbab = (*this->nonLocalDistributionsH1)(D3Q27System::ET_S, x1, x2p, x3);
+// real mfbba = (*this->nonLocalDistributionsH1)(D3Q27System::ET_B, x1, x2, x3p);
+// real mfaab = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SW, x1p, x2p, x3);
+// real mfcab = (*this->nonLocalDistributionsH1)(D3Q27System::ET_SE, x1, x2p, x3);
+// real mfaba = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BW, x1p, x2, x3p);
+// real mfcba = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BE, x1, x2, x3p);
+// real mfbaa = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BS, x1, x2p, x3p);
+// real mfbca = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BN, x1, x2, x3p);
+// real mfaaa = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BSW, x1p, x2p, x3p);
+// real mfcaa = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BSE, x1, x2p, x3p);
+// real mfaca = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
+// real mfcca = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
+// real mfbbb = (*this->zeroDistributionsH1)(x1, x2, x3);
+//
+// findNeighbors(phaseField, x1, x2, x3);
+//
+// real dX1_phi = gradX1_phi();
+// real dX2_phi = gradX2_phi();
+// real dX3_phi = gradX3_phi();
+//
+// real denom = sqrt(dX1_phi * dX1_phi + dX2_phi * dX2_phi + dX3_phi * dX3_phi) + 1.0e-20; //+ 1e-9+1e-3;
+// real normX1 = dX1_phi / denom;
+// real normX2 = dX2_phi / denom;
+// real normX3 = dX3_phi / denom;
+// real cx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+// (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) + (mfcbb - mfabb));
+// real cy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+// (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) + (mfbcb - mfbab));
+// real cz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+// (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) + (mfbbc - mfbba));
+//
+// (*vxNode)(x1, x2, x3) =(c6o1 * cx -(normX1 * oneOverInterfaceScale * (phi[DIR_000] - phiH) * (phi[DIR_000] - phiL)) / (phiH - phiL)) /(c6o1 * phi[DIR_000]);
+// (*vyNode)(x1, x2, x3) =(c6o1 * cy -(normX2 * oneOverInterfaceScale * (phi[DIR_000] - phiH) * (phi[DIR_000] - phiL)) / (phiH - phiL)) /(c6o1 * phi[DIR_000]);
+// (*vzNode)(x1, x2, x3) =(c6o1 * cz -(normX3 * oneOverInterfaceScale * (phi[DIR_000] - phiH) * (phi[DIR_000] - phiL)) / (phiH - phiL)) /(c6o1 * phi[DIR_000]);
+// }
+// }
+// }
+// }
+////!19.09.23
this->swapDistributions();
for (int x3 = minX3 - ghostLayerWidth+1; x3 < maxX3 + ghostLayerWidth-1; x3++) {
for (int x2 = minX2 - ghostLayerWidth+1; x2 < maxX2 + ghostLayerWidth-1; x2++) {
@@ -612,16 +729,20 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
real feqNew = D3Q27System::getIncompFeqForDirection(DIR_000, rhoG,vx,vy,vz);
distribution->setDistributionForDirection(fL-feqOLD+feqNew, x1, x2, x3, DIR_000);
}
- D3Q27System::calcIncompMacroscopicValues(ff, rhoG, vx, vy, vz);
- ff[DIR_000] = vx * vx + vy * vy + vz * vz +
- (((ff[DIR_MM0] + ff[DIR_PP0]) + (ff[DIR_MP0] + ff[DIR_PM0])) + ((ff[DIR_0MM] + ff[DIR_0PP]) + (ff[DIR_0MP] + ff[DIR_0PM])) + ((ff[DIR_M0M] + ff[DIR_P0P]) + (ff[DIR_M0P] + ff[DIR_P0M])) +
- c2o1 * ((((ff[DIR_MMM] + ff[DIR_PPP]) + (ff[DIR_MMP] + ff[DIR_PPM]))) + (((ff[DIR_MPM] + ff[DIR_PMP]) + (ff[DIR_MPP] + ff[DIR_PMM])))));
- distribution->setDistributionForDirection(ff[DIR_000], x1, x2, x3, DIR_000);
+ //D3Q27System::calcIncompMacroscopicValues(ff, rhoG, vx, vy, vz);
+ //ff[DIR_000] = vx * vx + vy * vy + vz * vz +
+ // (((ff[DIR_MM0] + ff[DIR_PP0]) + (ff[DIR_MP0] + ff[DIR_PM0])) + ((ff[DIR_0MM] + ff[DIR_0PP]) + (ff[DIR_0MP] + ff[DIR_0PM])) + ((ff[DIR_M0M] + ff[DIR_P0P]) + (ff[DIR_M0P] + ff[DIR_P0M])) +
+ // c2o1 * ((((ff[DIR_MMM] + ff[DIR_PPP]) + (ff[DIR_MMP] + ff[DIR_PPM]))) + (((ff[DIR_MPM] + ff[DIR_PMP]) + (ff[DIR_MPP] + ff[DIR_PMM])))));
+ //distribution->setDistributionForDirection(ff[DIR_000], x1, x2, x3, DIR_000);
}
+
else {//no refill of gas required
rhoG = (*rhoNode)(x1, x2, x3);
if (phi2[DIR_000] <= phiLim) { // no refill liquid
+ // real ffRecover[27];
+ // distribution->getDistributionInv(ffRecover, x1, x2, x3);//only needed for BB (03.10.2023)
+
for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
if ((phi[fdir] > phiLim)) {
// real vxBC = ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
@@ -685,6 +806,28 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
real fBC = feqNew + (fL - feqOLD) * (c1o1 / collFactorG - c1o1) / (c1o1 / collFactorL - c1o1) ;
+ // 03.10.2023
+ //real qq = ((*phaseField)(x1, x2, x3) - c1o2) /
+ // ((*phaseField)(x1, x2, x3) - (*phaseField)(x1 + D3Q27System::DX1[fdir],
+ // x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir]));
+ //real fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
+ //real fBC = ((fG - collFactorG * fGEQ) / (c1o1 - collFactorG) * (c1o1 - qq) -
+ // c6o1 * WEIGTH[fdir] *
+ // ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir]) *
+ // D3Q27System::DX1[fdir] +
+ // (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir]) *
+ // D3Q27System::DX2[fdir] +
+ // (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir]) *
+ // D3Q27System::DX3[fdir]) +
+ // qq * (c2o1*fG -c6o1*WEIGTH[fdir]*(D3Q27System::DX1[fdir]*vx+D3Q27System::DX2[fdir]*vy+D3Q27System::DX3[fdir]*vz))) /
+ // (qq + c1o1);
+ //! 03.10.2023
+
+
//real fBC = (fGEQ - WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
//real qq = c1o1 - ((c1o1 - c2o1 * (*phaseField)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir])) /
@@ -766,8 +909,38 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
// eqBCN = eqBC;
// distribution->setDistributionForDirection(LaplacePressure* WEIGTH[fdir] + (fBC + fG - eqBC - eqG) / densityRatio + (eqBCN + eqGN) * (c1o1 - c1o1 / densityRatio*0) - fL - 0*(feqG - feqL - 2 * fL + 2 * feqL) * (c1o1 / densityRatio - c1o1) * vBC, x1, x2,
// x3, fdir);// (vxBC * D3Q27System::DX1[fdir] + vyBC * D3Q27System::DX2[fdir] + vzBC * D3Q27System::DX3[fdir]), x1, x2, x3, fdir);
- distribution->setDistributionForDirection(laplacePressureBC * WEIGTH[fdir] + (fBC + fG) / densityRatio + (eqBCN + eqGN) * (c1o1 - c1o1 / densityRatio) - fL, x1, x2, x3, fdir);
-
+ real vLink = -(D3Q27System::DX1[fdir] * (*vxNode)(x1 + D3Q27System::DX1[fdir],
+ x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir]) +
+ D3Q27System::DX2[fdir] * (*vyNode)(x1 + D3Q27System::DX1[fdir],
+ x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir]) +
+ D3Q27System::DX3[fdir] * (*vzNode)(x1 + D3Q27System::DX1[fdir],
+ x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir]));
+ distribution->setDistributionForDirection(
+ laplacePressureBC * WEIGTH[fdir] + (fBC + fG) / densityRatio +
+ c1o2* ( (eqBCN + eqGN) * (c1o1 -2* c1o1 / densityRatio) - fL) +
+ //alternative to antiBB is BB
+ c1o2*(fL- WEIGTH[fdir]*(c6o1*vLink
+ +c2o1*(*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir])))+
+ //(-fL*vLink/(c1o1-vLink))+
+ //!BB
+
+ 0*(fL - feqOLD)+
+ 0*(c2o1 - collFactorL) * (fL - feqOLD) / (c1o1 - collFactorL),
+ x1, x2, x3, fdir);
+ /* std::cout
+ << "eqBCN=" << eqBCN << " eqGN=" << eqGN
+ << (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir])
+ << " "
+ <<
+ (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir])
+ << " "<<(*vzNode)(x1 + D3Q27System::DX1[fdir],
+ x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir])<< "\n";*/
@@ -873,7 +1046,22 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
(*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
real fBC = feqNew + (fL - feqOLD) * (c1o1 / collFactorG - c1o1) / (c1o1 / collFactorL - c1o1) ;
-
+ //03.10.2023
+ //real qq = ((*phaseField)(x1, x2, x3)-c1o2)/((*phaseField)(x1, x2, x3)-(*phaseField)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+ //real fGEQ = D3Q27System::getIncompFeqForDirection(fdir, rhoG, vx, vy, vz);
+ //real fBC =
+ // ((fG - collFactorG * fGEQ) / (c1o1 - collFactorG) * (c1o1 - qq) -
+ // c6o1 * WEIGTH[fdir] *
+ // ((*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]) * D3Q27System::DX1[fdir] + (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]) * D3Q27System::DX2[fdir] +
+ // (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+ // x3 + D3Q27System::DX3[fdir]) *
+ // D3Q27System::DX3[fdir]) +
+ // qq * (c2o1 * fG -
+ // c6o1 * WEIGTH[fdir] *
+ // (D3Q27System::DX1[fdir] * vx + D3Q27System::DX2[fdir] * vy +
+ // D3Q27System::DX3[fdir] * vz))) /
+ // (qq + c1o1);
+ //!03.10.2023
//real fBC = (fGEQ - WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) - c6o1 * WEIGTH[fdir] * (vBC);
//real qq = c1o1-((c1o1 - c2o1 * (*phaseField)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir])) /
@@ -953,17 +1141,36 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
//distribution->setDistributionForDirection(LaplacePressure* WEIGTH[fdir] + (fBC + fG - eqBC - eqG) / densityRatio + (eqBCN + eqGN) * (c1o1 - c1o1 / densityRatio*0) - fL - 0*(feqG - feqL - 2 * fL + 2 * feqL) * (c1o1 / densityRatio - c1o1) * vBC, x1, x2, x3, fdir);// (vxBC * D3Q27System::DX1[fdir] + vyBC * D3Q27System::DX2[fdir] + vzBC * D3Q27System::DX3[fdir]), x1, x2, x3, fdir);
// fBC = (fG) / (densityRatio - c1o1) +
// ((densityRatio) / (densityRatio - c1o1)) * ((eqBCN + eqGN) * (c1o1 - c1o1 / densityRatio) - c2o1 * fL + (fGEQ - WEIGTH[fdir] * dvDir * (c1o1 / collFactorG - c1o1)) + laplacePressureBC * WEIGTH[fdir]);
- // 13.07.2023
- real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]),
- (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]),
- (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
- real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoG, (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]),
- (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
-
- real fBC = feqNew + (fL - feqOLD) * (c1o1 / collFactorG - c1o1) / (c1o1 / collFactorL - c1o1) ;
-
-
- distribution->setDistributionForDirection(laplacePressureBC* WEIGTH[fdir] + (fBC + fG) / densityRatio + (eqBCN + eqGN) * (c1o1 - c1o1 / densityRatio) - fL , x1, x2, x3, fdir);
+ //// 13.07.2023
+ //real feqOLD = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], (*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]),
+ // (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]),
+ // (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+ //real feqNew = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], rhoG, (*vxNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]),
+ // (*vyNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]), (*vzNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir]));
+
+ //real fBC = feqNew + (fL - feqOLD) * (c1o1 / collFactorG - c1o1) / (c1o1 / collFactorL - c1o1) ;
+
+ real vLink = -(D3Q27System::DX1[fdir] * (*vxNode)(x1 + D3Q27System::DX1[fdir],
+ x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir]) +
+ D3Q27System::DX2[fdir] * (*vyNode)(x1 + D3Q27System::DX1[fdir],
+ x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir]) +
+ D3Q27System::DX3[fdir] * (*vzNode)(x1 + D3Q27System::DX1[fdir],
+ x2 + D3Q27System::DX2[fdir],
+ x3 + D3Q27System::DX3[fdir]));
+ distribution->setDistributionForDirection(
+ laplacePressureBC * WEIGTH[fdir] + (fBC + fG) / densityRatio +
+ c1o2* ( (eqBCN + eqGN) * (c1o1 -2* c1o1 / densityRatio) - fL) +
+ //alternative to antiBB is BB
+ c1o2*(fL- WEIGTH[fdir]*(c6o1*vLink
+ +c2o1*(*rhoNode)(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 + D3Q27System::DX3[fdir])))+
+ //(-fL*vLink/(c1o1-vLink))+
+ //!BB
+
+ 0*(fL - feqOLD)+
+ 0*(c2o1 - collFactorL) * (fL - feqOLD) / (c1o1 - collFactorL),
+ x1, x2, x3, fdir);
// real number = 666;
@@ -3371,7 +3578,8 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
//real mu = 2 * beta * phi[DIR_000] * (phi[DIR_000] - 1) * (2 * phi[DIR_000] - 1) - kappa * nabla2_phi();
//----------- Calculating Macroscopic Values -------------
- real rho = phi[DIR_000] > phiLim ? rhoH : rhoL;//rhoH + rhoToPhi * (phi[DIR_000] - phiH); //Incompressible
+ //real rho = phi[DIR_000] > phiLim ? rhoH : rhoL;//rhoH + rhoToPhi * (phi[DIR_000] - phiH); //Incompressible
+ real rho =rhoL+ (rhoH-rhoL) * (phi[DIR_000] - phiL)/(phiH-phiL);
//real rho = rhoH + rhoToPhi * (tanh(pushInterface*(c2o1*phi[DIR_000]-c1o1))/tanh(pushInterface)*c1o2 +c1o2 - phiH); //Incompressible
///scaled phase field
@@ -3947,17 +4155,18 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
real Dxy = -c3o1 * collFactorM * mfbba;
real Dxz = -c3o1 * collFactorM * mfbab;
real Dyz = -c3o1 * collFactorM * mfabb;
- // if (phi[DIR_000] > phiLim)
+ //if ((phi[DIR_000] > phiLim)
+ //if (((phi[DIR_000] < 0.9) || (phi[DIR_000]>0.1))
if (((phi[DIR_000] > phiLim) && ((phi[DIR_P00] <= phiLim) || (phi[DIR_M00] <= phiLim) || (phi[DIR_00P] <= phiLim) || (phi[DIR_00M] <= phiLim) || (phi[DIR_0M0] <= phiLim) || (phi[DIR_0P0] <= phiLim) || (phi[DIR_PP0] <= phiLim) || (phi[DIR_PM0] <= phiLim) || (phi[DIR_P0P] <= phiLim) ||
(phi[DIR_P0M] <= phiLim) || (phi[DIR_MP0] <= phiLim) || (phi[DIR_MM0] <= phiLim) || (phi[DIR_M0P] <= phiLim) || (phi[DIR_M0M] <= phiLim) || (phi[DIR_0PM] <= phiLim) || (phi[DIR_0MM] <= phiLim) || (phi[DIR_0PP] <= phiLim) || (phi[DIR_0MP] <= phiLim) ||
(phi[DIR_PPP] <= phiLim) || (phi[DIR_PMP] <= phiLim) || (phi[DIR_MPP] <= phiLim) || (phi[DIR_MMP] <= phiLim) ||
(phi[DIR_PPM] <= phiLim) || (phi[DIR_PMM] <= phiLim) || (phi[DIR_MPM] <= phiLim) || (phi[DIR_MMM] <= phiLim)))
- /*
+
|| ((phi[DIR_000] <= phiLim) && ((phi[DIR_P00] > phiLim) || (phi[DIR_M00] > phiLim) || (phi[DIR_00P] > phiLim) || (phi[DIR_00M] > phiLim) || (phi[DIR_0M0] > phiLim) || (phi[DIR_0P0] > phiLim) || (phi[DIR_PP0] > phiLim) ||
(phi[DIR_PM0] > phiLim) || (phi[DIR_P0P] > phiLim) || (phi[DIR_P0M] > phiLim) || (phi[DIR_MP0] > phiLim) || (phi[DIR_MM0] > phiLim) || (phi[DIR_M0P] > phiLim) || (phi[DIR_M0M] > phiLim) ||
(phi[DIR_0PM] > phiLim) || (phi[DIR_0MM] > phiLim) || (phi[DIR_0PP] > phiLim) || (phi[DIR_0MP] > phiLim) || (phi[DIR_PPP] > phiLim) || (phi[DIR_PMP] > phiLim) || (phi[DIR_MPP] > phiLim) ||
(phi[DIR_MMP] > phiLim) || (phi[DIR_PPM] > phiLim) || (phi[DIR_PMM] > phiLim) || (phi[DIR_MPM] > phiLim) || (phi[DIR_MMM] > phiLim)))
- */){
+ ){
// {
/// QR eddyviscosity:
@@ -4581,9 +4790,13 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
//!
////////////////////////////////////////////////////////////////////////////////////
// second component
- real concentration =
- ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) + (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) + ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
- ////////////////////////////////////////////////////////////////////////////////////
+ // real concentration =
+ // ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) + (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) + ((mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc))) + mfbbb;
+ // 06.10.2023 filtered phi
+ real concentration = phi[DIR_000];
+ //for (int dir = 0; dir <= 26; dir++) {concentration += WEIGTH[dir] * phi[dir];}
+
+ ////////////////////////////////////////////////////////////////////////////////////
//bool inverseConcentration = false;
//if (true) //(concentration <= phiLim)
//{
@@ -4699,9 +4912,11 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
real MomXDenom = sqrt(MomX1 * MomX1 + MomX2 * MomX2 + MomX3 * MomX3) + 1.0e-100;
real scaleNorm = (normX1 * MomX1 + normX2 * MomX2 + normX3 * MomX3) / MomXDenom;
- scaleNorm = scaleNorm * scaleNorm; //(c1o2 + c1o2 * scaleNorm) * scaleNorm; // scaleNorm * (c1o1+(c2o1*concentration - c1o1) * (c2o1*concentration - c1o1) * (scaleNorm-c1o1)); // *scaleNorm* scaleNorm;
-
- if (phi[DIR_000] > phiLim)
+ //scaleNorm = scaleNorm * scaleNorm; //(c1o2 + c1o2 * scaleNorm) * scaleNorm; // scaleNorm * (c1o1+(c2o1*concentration - c1o1) * (c2o1*concentration - c1o1) * (scaleNorm-c1o1)); // *scaleNorm* scaleNorm;
+ //scaleNorm = scaleNorm * scaleNorm;
+ //scaleNorm = scaleNorm * scaleNorm;
+ //scaleNorm = scaleNorm * scaleNorm;
+ if (true)//(phi[DIR_000] > phiLim) //(true) // ((phi[DIR_000] > phiLim)||(normX1*vvx+normX2*vvy+normX3*vvz<0))
{
normX1 = (normX1 * (c1o1 - mixNormal) + mixNormal * MomX1 / MomXDenom) * scaleNorm;
@@ -4728,13 +4943,18 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
normX2 = (normX2 * (c1o1 - mixNormal) + mixNormal * MomX2 / MomXDenom) * scaleNorm;
normX3 = (normX3 * (c1o1 - mixNormal) + mixNormal * MomX3 / MomXDenom) * scaleNorm;
real scaleAdvectionContribution = (c1o1 - (concentration - phiL) / (phiH - phiL) * c2o1) ;
- scaleAdvectionContribution = scaleAdvectionContribution * scaleAdvectionContribution * scaleAdvectionContribution * scaleAdvectionContribution * scaleAdvectionContribution;
- cx = scaleAdvectionContribution * (cx * (c1o1 - omegaD) + omegaD * vvx * concentration + normX1 * (c1o1 - 0.5 * omegaD) * (phiH - concentration) * (concentration - phiL) * c1o3 * oneOverInterfaceScale / (phiH - phiL))
+ // scaleAdvectionContribution = scaleAdvectionContribution * scaleAdvectionContribution * scaleAdvectionContribution * scaleAdvectionContribution * scaleAdvectionContribution;
+ //scaleAdvectionContribution = scaleAdvectionContribution * scaleAdvectionContribution *
+ // scaleAdvectionContribution * scaleAdvectionContribution *
+ // scaleAdvectionContribution;
+ // scaleAdvectionContribution = 0;
+ // scaleAdvectionContribution = (concentration > 0.0001) ? 0 : 1;
+ cx = scaleAdvectionContribution * (cx * (c1o1 - omegaD) + omegaD * vvx * concentration + normX1 * (c1o1 - 0.5 * omegaD) * (phiH - concentration) * (concentration - phiL) * c1o3 * oneOverInterfaceScale / (phiH - phiL))
+(c1o1-scaleAdvectionContribution)*(cx - normX1FD* (oneOverInterfaceScale * (concentration - phiH) * (concentration - phiL)) / ((phiH - phiL)) * c1o3);
cy = scaleAdvectionContribution * (cy * (c1o1 - omegaD) + omegaD * vvy * concentration + normX2 * (c1o1 - 0.5 * omegaD) * (phiH - concentration) * (concentration - phiL) * c1o3 * oneOverInterfaceScale / (phiH - phiL))
+(c1o1 - scaleAdvectionContribution) *
(cy - normX2FD * (oneOverInterfaceScale * (concentration - phiH) * (concentration - phiL)) / ((phiH - phiL)) * c1o3);
- cz = scaleAdvectionContribution * (cz * (c1o1 - omegaD) + omegaD * vvx * concentration + normX3 * (c1o1 - 0.5 * omegaD) * (phiH - concentration) * (concentration - phiL) * c1o3 * oneOverInterfaceScale / (phiH - phiL))
+ cz = scaleAdvectionContribution * (cz * (c1o1 - omegaD) + omegaD * vvz * concentration + normX3 * (c1o1 - 0.5 * omegaD) * (phiH - concentration) * (concentration - phiL) * c1o3 * oneOverInterfaceScale / (phiH - phiL))
+(c1o1 - scaleAdvectionContribution) *
(cz - normX3FD * (oneOverInterfaceScale * (concentration - phiH) * (concentration - phiL)) / ((phiH - phiL)) * c1o3);
// cy = cy - normX2FD * (oneOverInterfaceScale * (concentration - phiH) * (concentration - phiL)) / ((phiH - phiL)) * c1o3;
@@ -5210,7 +5430,56 @@ void MultiphaseScaleDistributionLBMKernel::calculate(int step)
// }
// }
//}
-}
+
+//{
+// real fff[27];
+// SPtr<DistributionArray3D> distributionC = this->getDataSet()->getFdistributions();
+// // for (int x3 = minX3 - 1; x3 < maxX3 + 1; x3++) {
+// // for (int x2 = minX2 - 1; x2 < maxX2 + 1; x2++) {
+// // for (int x1 = minX1 - 1; x1 < maxX1 + 1; x1++) {
+// // if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3))
+// for (int x3 = minX3 - ghostLayerWidth + 1; x3 < maxX3 + ghostLayerWidth - 1; x3++) {
+// for (int x2 = minX2 - ghostLayerWidth + 1; x2 < maxX2 + ghostLayerWidth - 1; x2++) {
+// for (int x1 = minX1 - ghostLayerWidth + 1; x1 < maxX1 + ghostLayerWidth - 1; x1++) {
+// if (!bcArray->isSolid(x1, x2, x3) /* && !bcArray->isUndefined(x1, x2, x3)*/)
+//
+// {
+//
+// findNeighbors(phaseFieldOld, x1, x2, x3);
+// findNeighbors2(phaseField, x1, x2, x3);
+// if ((phi[DIR_000] <= phiLim) && (phi2[DIR_000] > phiLim)) {
+// for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
+// if ((phi2[fdir] > phiLim) &&
+// (!bcArray->isSolid(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+// x3 + D3Q27System::DX3[fdir]) &&
+// !bcArray->isUndefined(x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir],
+// x3 + D3Q27System::DX3[fdir]))) {
+// distributionC->getDistributionInv(fff, x1 + D3Q27System::DX1[fdir],
+// x2 + D3Q27System::DX2[fdir],
+// x3 + D3Q27System::DX3[fdir]);
+// real vx, vy, vz, rho;
+// D3Q27System::calcIncompMacroscopicValues(fff, rho, vx, vy, vz);
+// real feq = D3Q27System::getIncompFeqForDirection(fdir, rho, vx, vy, vz);
+// // real feqZero=D3Q27System::getIncompFeqForDirection(fdir, 0.0,vx,vy,vz);
+// // real feqZeroI = D3Q27System::getIncompFeqForDirection(D3Q27System::INVDIR[fdir], 0.0,
+// // vx, vy, vz); distribution->setDistributionForDirection(
+// // c1o2 * (feqZero+feqZeroI-feq+
+// // distribution->getDistributionInvForDirection(
+// // x1 + D3Q27System::DX1[fdir], x2 + D3Q27System::DX2[fdir], x3 +
+// // D3Q27System::DX3[fdir],fdir)),
+// // x1, x2, x3,
+// // fdir);
+// distributionC->setDistributionForDirection(feq, x1, x2, x3, fdir);
+// }
+// }
+// }
+// }
+// }
+// }
+// }
+// }
+
+ }
@@ -5448,18 +5717,34 @@ void MultiphaseScaleDistributionLBMKernel::findNeighbors2(CbArray3D<real, Indexe
for (int k = FSTARTDIR; k <= FENDDIR; k++) {
-
- if (!bcArray->isSolid(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k])) {
- phi2[k] = (*ph)(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k]);
- }
- else {
- //if (bcArray->getBC(x1, x2, x3)->hasVelocityBoundaryFlag(D3Q27System::INVDIR[k]))
- // phi2[k] = (*ph)(x1, x2, x3); // neutral wetting
- //else
- phi2[k] = 0.0; // unwetting
- // phi2[k] = (*ph)(x1, x2, x3) * 0.7;
- }
- }
+ if (!bcArray->isSolid(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k])) {
+ phi2[k] = (*ph)(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k]);
+ } else {
+ // if (bcArray->getBC(x1, x2, x3)->hasVelocityBoundaryFlag(D3Q27System::INVDIR[k]))
+ // phi[k] = (*ph)(x1, x2, x3); // neutral wetting
+ // else
+ // phi[k] = 0.0; // unwetting
+ // phi[k] = (*ph)(x1, x2, x3) * 0.7;
+ if (bcArray->getBC(x1, x2, x3)->hasNoSlipBoundaryFlag(D3Q27System::INVDIR[k])) {
+ if (bcArray->getBC(x1, x2, x3)->getNoSlipSecondaryOption(D3Q27System::INVDIR[k]) == 0)
+ phi2[k] = (*ph)(x1, x2, x3); // neutral wetting
+ else
+ phi2[k] = 0.0; // unwetting
+ }
+ }
+ }
+
+ // if (!bcArray->isSolid(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k])) {
+ // phi2[k] = (*ph)(x1 + DX1[k], x2 + DX2[k], x3 + DX3[k]);
+ // }
+ // else {
+ // //if (bcArray->getBC(x1, x2, x3)->hasVelocityBoundaryFlag(D3Q27System::INVDIR[k]))
+ // // phi2[k] = (*ph)(x1, x2, x3); // neutral wetting
+ // //else
+ // phi2[k] = 0.0; // unwetting
+ // // phi2[k] = (*ph)(x1, x2, x3) * 0.7;
+ // }
+ //}
}
void MultiphaseScaleDistributionLBMKernel::swapDistributions()