diff --git a/apps/cpu/FallingSphere/in.lbdem b/apps/cpu/FallingSphere/in.lbdem
index 4a526a9e3a0296938b0128d928897f003fe019c3..5d0cfb511774c5270918f4e7007224d36080232c 100644
--- a/apps/cpu/FallingSphere/in.lbdem
+++ b/apps/cpu/FallingSphere/in.lbdem
@@ -42,7 +42,7 @@ fix xwalls2 all wall/gran model hertz tangential history primitive type 1 xplane
fix ywalls1 all wall/gran model hertz tangential history primitive type 1 yplane 0.
fix ywalls2 all wall/gran model hertz tangential history primitive type 1 yplane 1.
fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.
-fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane 2.
+fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane 10.
create_atoms 1 single 0.5 0.5 9.75
#create_atoms 1 single 0.38 0.05 0.05
diff --git a/apps/cpu/JetBreakup/JetBreakup.cpp b/apps/cpu/JetBreakup/JetBreakup.cpp
index f4b74ca379edaf8840cb6875ca0ff9fc7f296509..9a1331995657a7a7e43eeb948134c154061c431b 100644
--- a/apps/cpu/JetBreakup/JetBreakup.cpp
+++ b/apps/cpu/JetBreakup/JetBreakup.cpp
@@ -2,14 +2,13 @@
#include <memory>
#include <string>
-#include "VirtualFluids.h"
#include "MultiphaseFlow/MultiphaseFlow.h"
+#include "VirtualFluids.h"
using namespace std;
void setInflowBC(real x1, real x2, real x3, real radius, int dir)
{
-
}
void run(string configname)
@@ -24,20 +23,20 @@ void run(string configname)
config.load(configname);
string pathname = config.getValue<string>("pathname");
- //string pathGeo = config.getValue<string>("pathGeo");
- //string geoFile = config.getValue<string>("geoFile");
+ // 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");
+ // vector<double> boundingBox = config.getVector<double>("boundingBox");
+ // vector<double> length = config.getVector<double>("length");
real U_LB = config.getValue<real>("U_LB");
// double uF2 = config.getValue<double>("uF2");
- //double nuL = config.getValue<double>("nuL");
- //double nuG = config.getValue<double>("nuG");
- //double densityRatio = config.getValue<double>("densityRatio");
- //double sigma = config.getValue<double>("sigma");
+ // 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 D = config.getValue<double>("D");
real theta = config.getValue<real>("contactAngle");
real D_LB = config.getValue<real>("D_LB");
real phiL = config.getValue<real>("phi_L");
@@ -48,24 +47,21 @@ void run(string configname)
real endTime = config.getValue<real>("endTime");
real outTime = config.getValue<real>("outTime");
real availMem = config.getValue<real>("availMem");
- //int refineLevel = config.getValue<int>("refineLevel");
- //double Re = config.getValue<double>("Re");
-
+ // int refineLevel = config.getValue<int>("refineLevel");
+ // double Re = config.getValue<double>("Re");
+
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");
-
-
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 (myid == 0) UBLOG(logINFO, "Jet Breakup: Start!");
if (logToFile) {
#if defined(__unix__)
@@ -84,25 +80,25 @@ void run(string configname)
// Sleep(30000);
- real rho_h=0, rho_l=0, r_rho=0, mu_h=0, /*mu_l,*/ Uo=0, D=0, sigma=0;
+ real rho_h = 0, rho_l = 0, r_rho = 0, mu_h = 0, /*mu_l,*/ Uo = 0, D = 0, sigma = 0;
switch (caseN) {
- case 1:
- //density of heavy fluid (kg/m^3)
- rho_h = 848;
- //density of light fluid (kg/m^3)
+ case 1:
+ // density of heavy fluid (kg/m^3)
+ rho_h = 848;
+ // density of light fluid (kg/m^3)
rho_l = 34.5;
- //density ratio
+ // density ratio
r_rho = rho_h / rho_l;
- //dynamic viscosity of heavy fluid (Pa � s)
+ // 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)
+ // dynamic viscosity of light fluid (Pa � s)
+ // mu_l = 1.97e-5;
+ // velocity (m/s)
Uo = 100;
- //diameter of jet (m)
+ // diameter of jet (m)
D = 0.0001;
- //surface tension (N/m)
+ // surface tension (N/m)
sigma = 0.03;
break;
case 2:
@@ -115,7 +111,7 @@ void run(string configname)
// dynamic viscosity of heavy fluid (Pa � s)
mu_h = 2.87e-3;
// dynamic viscosity of light fluid (Pa � s)
- //mu_l = 1.84e-5;
+ // mu_l = 1.84e-5;
// velocity (m/s)
Uo = 200;
// diameter of jet (m)
@@ -133,14 +129,14 @@ void run(string configname)
// dynamic viscosity of heavy fluid (Pa � s)
mu_h = 2.87e-3;
// dynamic viscosity of light fluid (Pa � s)
- //mu_l = 1.84e-5;
+ // 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;
+ break;
}
real Re = rho_h * Uo * D / mu_h;
@@ -149,14 +145,15 @@ void run(string configname)
real dx = D / D_LB;
real nu_h = U_LB * D_LB / Re;
real nu_l = nu_h;
+ nu_h *= 100;
real rho_h_LB = 1;
- //surface tension
+ // surface tension
real sigma_LB = rho_h_LB * U_LB * U_LB * D_LB / We;
// LBMReal dLB = 0; // = length[1] / dx;
real rhoLB = 0.0;
- //LBMReal nuLB = nu_l; //(uLB*dLB) / Re;
+ // LBMReal nuLB = nu_l; //(uLB*dLB) / Re;
real beta = 12.0 * sigma_LB / interfaceWidth;
real kappa = 1.5 * interfaceWidth * sigma_LB;
@@ -177,6 +174,19 @@ void run(string configname)
UBLOG(logINFO, "path = " << pathname);
}
+ // bounding box
+ real g_minX1 = 0;
+ real g_minX2 = 0;
+ real g_minX3 = 0;
+
+ // double g_maxX1 = 8.0*D;
+ // double g_maxX2 = 2.5*D;
+ // double g_maxX3 = 2.5*D;
+
+ real g_maxX1 = 2.0 * D;
+ real g_maxX2 = 2.5 * D;
+ real g_maxX3 = 2.5 * D;
+
SPtr<LBMUnitConverter> conv(new LBMUnitConverter());
// const int baseLevel = 0;
@@ -188,8 +198,10 @@ void run(string configname)
// kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsCumulantLBMKernel());
// kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsVelocityCumulantLBMKernel());
// kernel = SPtr<LBMKernel>(new MultiphaseTwoPhaseFieldsPressureFilterLBMKernel());
- //kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
- kernel = SPtr<LBMKernel>(new MultiphaseSimpleVelocityBaseExternalPressureLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphasePressureFilterLBMKernel());
+ // kernel = SPtr<LBMKernel>(new MultiphaseSimpleVelocityBaseExternalPressureLBMKernel());
+ kernel = make_shared<MultiphaseScaleDistributionLBMKernel>();
+ //kernel = make_shared<MultiphaseSharpInterfaceLBMKernel>();
kernel->setWithForcing(true);
kernel->setForcingX1(0.0);
@@ -208,7 +220,8 @@ void run(string configname)
kernel->setMultiphaseModelParameters(beta, kappa);
kernel->setContactAngle(theta);
kernel->setInterfaceWidth(interfaceWidth);
- //dynamicPointerCast<MultiphasePressureFilterLBMKernel>(kernel)->setPhaseFieldBC(0.0);
+ // dynamicPointerCast<MultiphasePressureFilterLBMKernel>(kernel)->setPhaseFieldBC(0.0);
+ kernel->setSigma(sigma);
SPtr<BCSet> bcProc(new BCSet());
// BCSetPtr bcProc(new ThinWallBCSet());
@@ -221,8 +234,7 @@ void run(string configname)
// grid->setPeriodicX3(true);
grid->setGhostLayerWidth(2);
- SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(
- comm, MetisPartitioningGridVisitor::LevelBased, DIR_MMM, MetisPartitioner::RECURSIVE));
+ SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, DIR_MMM, MetisPartitioner::RECURSIVE));
//////////////////////////////////////////////////////////////////////////
// restart
@@ -240,25 +252,26 @@ void run(string configname)
// 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<BC> velBCF1(
- // new MultiphaseVelocityBC(false, true, false, fctF1, phiH, 0.0, BCFunction::INFCONST));
+ // u_actual(r) = U_max * (4 / (D^2) * (R^2 - r^2))
+ fctF1.SetExpr("vx1*(1-((x2-y0)^2+(x3-z0)^2)/(R^2))");
+ // fctF1.SetExpr("vx1*(1-(sqrt((x2-y0)^2+(x3-z0)^2))/(R))");
+ // fctF1.SetExpr("vy1*(1-(sqrt((x2-x0)^2+(x3-z0)^2)/R))^0.1");
+ // fctF1.SetExpr("vx1");
+ //fctF1.DefineConst("vx1", 0);
+ fctF1.DefineConst("vx1", U_LB);
+ fctF1.DefineConst("R", 0.5 * D);
+ fctF1.DefineConst("y0", (g_minX2 + g_maxX2) / 2);
+ fctF1.DefineConst("z0", (g_minX3 + g_maxX3) / 2);
+ SPtr<BC> velBCF1(new MultiphaseVelocityBC(true, false, false, fctF1, phiH, 0.0, BCFunction::INFCONST));
mu::Parser fctF2;
- fctF2.SetExpr("vy1");
- fctF2.DefineConst("vy1", U_LB);
+ fctF2.SetExpr("vx1");
+ fctF2.DefineConst("vx1", U_LB);
- real startTime = 1;
- SPtr<BC> velBCF1(
- new MultiphaseVelocityBC(true, false, false, fctF1, phiH, 0.0, startTime));
- SPtr<BC> velBCF2(
- new MultiphaseVelocityBC(true, false, false, fctF2, phiH, startTime, endTime));
+ // real startTime = 1;
+ // SPtr<BC> velBCF1(new MultiphaseVelocityBC(true, false, false, fctF1, phiH, 0.0, startTime));
+ // SPtr<BC> velBCF2(new MultiphaseVelocityBC(true, false, false, fctF2, phiH, startTime, endTime));
+ SPtr<BC> velBCF2(new MultiphaseVelocityBC(true, false, false, fctF2, phiH*0.2, 0.0, BCFunction::INFCONST));
SPtr<BC> noSlipBC(new NoSlipBC());
noSlipBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseNoSlipBCStrategy()));
@@ -279,15 +292,17 @@ void run(string configname)
fctvel_F2_init.DefineConst("U", 0);
velBCF1->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+ velBCF2->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
//////////////////////////////////////////////////////////////////////////////////
// BC visitor
MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
bcVisitor.addBC(noSlipBC);
- bcVisitor.addBC(denBC); // Ohne das BB?
+ // bcVisitor.addBC(denBC); // Ohne das BB?
bcVisitor.addBC(velBCF1);
+ bcVisitor.addBC(velBCF2);
- //SPtr<D3Q27Interactor> inflowF1Int;
- //SPtr<D3Q27Interactor> cylInt;
+ // SPtr<D3Q27Interactor> inflowF1Int;
+ // SPtr<D3Q27Interactor> cylInt;
SPtr<D3Q27Interactor> inflowInt;
@@ -295,54 +310,31 @@ void run(string configname)
// if (newStart) {
- // bounding box
- real g_minX1 = 0;
- real g_minX2 = 0;
- real g_minX3 = 0;
-
- //double g_maxX1 = 8.0*D;
- //double g_maxX2 = 2.5*D;
- //double g_maxX3 = 2.5*D;
-
- real g_maxX1 = 1.0 * D; // 8.0 * D;
- real g_maxX2 = 2.0 * D;
- real g_maxX3 = 2.0 * 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());
+ if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), pathname + "/geo/gridCube", WbWriterVtkXmlBinary::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());
+ // 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 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 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));
@@ -352,8 +344,8 @@ void run(string configname)
// 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_maxX1, g_minX2 - 2.0*dx, g_minX3 - 2.0*dx, g_maxX1 + 2.0*dx, g_maxX2 + 2.0*dx, g_maxX3));
- if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr geoOutflow(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));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), pathname + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
// double blockLength = blocknx[0] * dx;
@@ -371,11 +363,10 @@ void run(string configname)
GenBlocksGridVisitor genBlocks(gridCube);
grid->accept(genBlocks);
- SPtr<WriteBlocksSimulationObserver> ppblocks(new WriteBlocksSimulationObserver(
- grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ SPtr<WriteBlocksSimulationObserver> ppblocks(new WriteBlocksSimulationObserver(grid, SPtr<UbScheduler>(new UbScheduler(1)), pathname, WbWriterVtkXmlBinary::getInstance(), comm));
- //SPtr<Interactor3D> tubes(new D3Q27TriFaceMeshInteractor(cylinder, grid, noSlipBC,
- // Interactor3D::SOLID, Interactor3D::POINTS));
+ // SPtr<Interactor3D> tubes(new D3Q27TriFaceMeshInteractor(cylinder, grid, noSlipBC,
+ // Interactor3D::SOLID, Interactor3D::POINTS));
// inflowF1Int =
// SPtr<D3Q27Interactor>(new D3Q27Interactor(cylinder1, grid, noSlipBC, Interactor3D::SOLID));
@@ -384,64 +375,55 @@ void run(string configname)
SPtr<D3Q27Interactor> outflowInt(new D3Q27Interactor(geoOutflow, grid, denBC, 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));
+ 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));
+ 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));
+ 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));
+ 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));
+ 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));
+ 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, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallXmaxInt(
- new D3Q27Interactor(wallXmax, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallZminInt(
- new D3Q27Interactor(wallZmin, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallZmaxInt(
- new D3Q27Interactor(wallZmax, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallYminInt(
- new D3Q27Interactor(wallYmin, grid, noSlipBC, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> wallYmaxInt(
- new D3Q27Interactor(wallYmax, grid, noSlipBC, Interactor3D::SOLID));
-
- //cylInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(cylinder1, grid, velBCF1, Interactor3D::SOLID));
- //cylInt->addBC(velBCF2);
- // SPtr<D3Q27Interactor> cyl2Int(new D3Q27Interactor(cylinder2, grid, noSlipBC,
- // Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXminInt(new D3Q27Interactor(wallXmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallXmaxInt(new D3Q27Interactor(wallXmax, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallYminInt(new D3Q27Interactor(wallYmin, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallYmaxInt(new D3Q27Interactor(wallYmax, grid, noSlipBC, Interactor3D::SOLID));
+
+ // cylInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(cylinder1, grid, velBCF1, Interactor3D::SOLID));
+ // cylInt->addBC(velBCF2);
+ // SPtr<D3Q27Interactor> cyl2Int(new D3Q27Interactor(cylinder2, grid, noSlipBC,
+ // Interactor3D::SOLID));
inflowInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, velBCF1, Interactor3D::SOLID));
- inflowInt->addBC(velBCF2);
+ // inflowInt->addBC(velBCF2);
+
+ GbCylinder3DPtr geoAirInflow(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 * 0.9 / 2.0));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoAirInflow.get(), pathname + "/geo/geoAirInflow", WbWriterVtkXmlASCII::getInstance());
+ SPtr<D3Q27Interactor> inflowAirInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoAirInflow, grid, velBCF2, Interactor3D::SOLID));
- SPtr<D3Q27Interactor> solidInt =
- SPtr<D3Q27Interactor>(new D3Q27Interactor(geoSolid, grid, noSlipBC, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> solidInt = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoSolid, grid, noSlipBC, Interactor3D::SOLID));
InteractorsHelper intHelper(grid, metisVisitor, true);
- //intHelper.addInteractor(cylInt);
- //intHelper.addInteractor(tubes);
- intHelper.addInteractor(outflowInt);
- // intHelper.addInteractor(cyl2Int);
+ // intHelper.addInteractor(cylInt);
+ // intHelper.addInteractor(tubes);
+ // intHelper.addInteractor(outflowInt);
+ // intHelper.addInteractor(cyl2Int);
intHelper.addInteractor(wallXminInt);
- //intHelper.addInteractor(wallXmaxInt);
+ intHelper.addInteractor(wallXmaxInt);
intHelper.addInteractor(wallZminInt);
intHelper.addInteractor(wallZmaxInt);
intHelper.addInteractor(wallYminInt);
intHelper.addInteractor(wallYmaxInt);
- intHelper.addInteractor(inflowInt);
- //intHelper.addInteractor(solidInt);
+ //intHelper.addInteractor(inflowInt);
+ intHelper.addInteractor(inflowAirInt);
intHelper.selectBlocks();
@@ -450,13 +432,10 @@ void run(string configname)
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 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));
+ 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) {
@@ -478,31 +457,44 @@ void run(string configname)
grid->accept(kernelVisitor);
- //if (refineLevel > 0) {
- // SetUndefinedNodesBlockVisitor undefNodesVisitor;
- // grid->accept(undefNodesVisitor);
- //}
+ // if (refineLevel > 0) {
+ // SetUndefinedNodesBlockVisitor undefNodesVisitor;
+ // grid->accept(undefNodesVisitor);
+ // }
intHelper.setBC();
// initialization of distributions
- //mu::Parser fct1;
- //fct1.SetExpr("phiL");
- //fct1.DefineConst("phiL", phiL);
- real x1c = 0; // (g_maxX1 - g_minX1-1)/2; //
- real x2c = (g_maxX2 - g_minX2)/2;
- real x3c = (g_maxX3 - g_minX3)/2;
-
+ // mu::Parser fct1;
+ // fct1.SetExpr("phiL");
+ // fct1.DefineConst("phiL", phiL);
+ real x1c = 0; // (g_maxX1 - g_minX1-1)/2; //
+ real x2c = (g_minX2 + g_maxX2) / 2;
+ real x3c = (g_minX3 + g_maxX3) / 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("(0.5-0.5*tanh(2*(sqrt((x1/radius-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness))-(0.5-0.5*tanh(2*(sqrt((x1/radius-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius2)/interfaceThickness))");
+ //fct1.SetExpr("x1 < 4*dx ? (0.5-0.5*tanh(2*(1-4*dx)*(sqrt((x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)) : 0");
fct1.DefineConst("x1c", x1c);
fct1.DefineConst("x2c", x2c);
fct1.DefineConst("x3c", x3c);
- fct1.DefineConst("radius", 0.5*D);
- fct1.DefineConst("interfaceThickness", interfaceWidth*dx);
+ fct1.DefineConst("dx", dx);
+ fct1.DefineConst("radius", 0.5 * D /* + 2. * dx*/);
+ fct1.DefineConst("radius2", D/4.);
+ fct1.DefineConst("interfaceThickness", interfaceWidth * dx);
+
+ mu::Parser fct2;
+ // fct1.SetExpr("0.5-0.5*tanh(2*(sqrt((x1/radius-x1c)^2+(x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)");
+ fct2.SetExpr("x1 < 4*dx ? (0.5-0.5*tanh(2*(1-4*dx)*(sqrt((x2-x2c)^2+(x3-x3c)^2)-radius)/interfaceThickness)) : 0");
+ fct2.DefineConst("x1c", x1c);
+ fct2.DefineConst("x2c", x2c);
+ fct2.DefineConst("x3c", x3c);
+ fct2.DefineConst("dx", dx);
+ fct2.DefineConst("radius", 0.5 * D /* + 2. * dx*/);
+ fct2.DefineConst("interfaceThickness", interfaceWidth * dx);
MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
- initVisitor.setPhi(fct1);
+ //initVisitor.setPhi(fct1);
grid->accept(initVisitor);
///////////////////////////////////////////////////////////////////////////////////////////
//{
@@ -549,22 +541,19 @@ void run(string configname)
// boundary conditions grid
{
SPtr<UbScheduler> geoSch(new UbScheduler(1));
- SPtr<WriteBoundaryConditionsSimulationObserver> ppgeo(new WriteBoundaryConditionsSimulationObserver(
- grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
+ SPtr<WriteBoundaryConditionsSimulationObserver> ppgeo(new WriteBoundaryConditionsSimulationObserver(grid, geoSch, pathname, WbWriterVtkXmlBinary::getInstance(), comm));
ppgeo->update(0);
ppgeo.reset();
}
- if (myid == 0)
- UBLOG(logINFO, "Preprocess - end");
+ if (myid == 0) UBLOG(logINFO, "Preprocess - end");
} else {
rcp->restart((int)restartStep);
grid->setTimeStep(restartStep);
- if (myid == 0)
- UBLOG(logINFO, "Restart - end");
+ if (myid == 0) UBLOG(logINFO, "Restart - end");
}
-
+
// TwoDistributionsSetConnectorsBlockVisitor setConnsVisitor(comm);
// grid->accept(setConnsVisitor);
@@ -572,25 +561,26 @@ void run(string configname)
grid->accept(bcVisitor);
- ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
- //TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ // ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
grid->accept(setConnsVisitor);
SPtr<UbScheduler> visSch(new UbScheduler(outTime));
real t_ast, t;
t_ast = 7.19;
- t = (int)(t_ast/(U_LB/(D_LB)));
- visSch->addSchedule(t,t,t); //t=7.19
- SPtr<WriteMultiphaseQuantitiesSimulationObserver> pp(new WriteMultiphaseQuantitiesSimulationObserver(
- grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
+ t = (int)(t_ast / (U_LB / (D_LB)));
+ visSch->addSchedule(t, t, t); // t=7.19
+ // SPtr<WriteMultiphaseQuantitiesSimulationObserver> pp(new WriteMultiphaseQuantitiesSimulationObserver(
+ // grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
+ SPtr<WriteSharpInterfaceQuantitiesSimulationObserver> pp(new WriteSharpInterfaceQuantitiesSimulationObserver(grid, visSch, pathname, WbWriterVtkXmlBinary::getInstance(), conv, comm));
pp->update(0);
SPtr<UbScheduler> nupsSch(new UbScheduler(10, 30, 100));
SPtr<NUPSCounterSimulationObserver> npr(new NUPSCounterSimulationObserver(grid, nupsSch, numOfThreads, comm));
- SPtr<UbScheduler> timeBCSch(new UbScheduler(1, startTime, startTime));
- auto timeDepBC = make_shared<TimeDependentBCSimulationObserver>(TimeDependentBCSimulationObserver(grid, timeBCSch));
- timeDepBC->addInteractor(inflowInt);
+ // SPtr<UbScheduler> timeBCSch(new UbScheduler(1, startTime, startTime));
+ // auto timeDepBC = make_shared<TimeDependentBCSimulationObserver>(TimeDependentBCSimulationObserver(grid, timeBCSch));
+ // timeDepBC->addInteractor(inflowInt);
#ifdef _OPENMP
omp_set_num_threads(numOfThreads);
@@ -600,14 +590,12 @@ void run(string configname)
SPtr<Simulation> simulation(new Simulation(grid, stepGhostLayer, endTime));
simulation->addSimulationObserver(npr);
simulation->addSimulationObserver(pp);
- simulation->addSimulationObserver(timeDepBC);
+ // simulation->addSimulationObserver(timeDepBC);
simulation->addSimulationObserver(rcp);
- if (myid == 0)
- UBLOG(logINFO, "Simulation-start");
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
simulation->run();
- if (myid == 0)
- UBLOG(logINFO, "Simulation-end");
+ if (myid == 0) UBLOG(logINFO, "Simulation-end");
} catch (std::exception &e) {
cerr << e.what() << endl << flush;
} catch (std::string &s) {
diff --git a/apps/cpu/Nozzle/CMakeLists.txt b/apps/cpu/Nozzle/CMakeLists.txt
index f9e90c7e95115d4b4dfc436721a192100ad95e98..9764934d46b39abca2ceda133e3434a78980f294 100644
--- a/apps/cpu/Nozzle/CMakeLists.txt
+++ b/apps/cpu/Nozzle/CMakeLists.txt
@@ -1,3 +1,3 @@
PROJECT(Nozzle)
-vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} LiggghtsCoupling MultiphaseFlow NonNewtonianFluids FILES nozzle.cpp )
+vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} LiggghtsCoupling MultiphaseFlow NonNewtonianFluids FILES nozzleSinglePhase.cpp )
diff --git a/apps/cpu/Nozzle/nozzle.cpp b/apps/cpu/Nozzle/nozzle.cpp.1
similarity index 98%
rename from apps/cpu/Nozzle/nozzle.cpp
rename to apps/cpu/Nozzle/nozzle.cpp.1
index dc4e9e40ac63332fa0b54aa06aa3c7b32ec4381a..0feddc57b709be57c05000da0ef8075b6941b8c0 100644
--- a/apps/cpu/Nozzle/nozzle.cpp
+++ b/apps/cpu/Nozzle/nozzle.cpp.1
@@ -126,8 +126,8 @@ int main(int argc, char *argv[])
if (myid == 0) VF_LOG_INFO("Yield stress = {} Pa", tau0);
if (myid == 0) VF_LOG_INFO("Yield stress LB = {} ", tau0_LB);
- // SPtr<BC> noSlipBC(new NoSlipBC());
- // noSlipBC->setBCStrategy(SPtr<BCStrategy>(new NoSlipBCStrategy()));
+ //SPtr<BC> noSlipBC(new NoSlipBC());
+ //noSlipBC->setBCStrategy(SPtr<BCStrategy>(new NoSlipBCStrategy()));
SPtr<BC> noSlipBC(new NoSlipBC());
noSlipBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseNoSlipBCStrategy()));
@@ -161,8 +161,8 @@ int main(int argc, char *argv[])
// fct.DefineConst("U", uLB_ref * ((N + 3) / (N + 1)));
// fct.DefineConst("NplusOne", N + 1.0);
- // SPtr<BC> inflowConcreteBC(new VelocityBC(false, false, true, fct, 0, BCFunction::INFCONST));
- // inflowConcreteBC->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowConcreteBC(new VelocityBC(false, false, true, fct, 0, BCFunction::INFCONST));
+ //inflowConcreteBC->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
SPtr<BC> inflowConcreteBC(new MultiphaseVelocityBC(false, false, true, fct, phiH, 0, BCFunction::INFCONST));
inflowConcreteBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
@@ -444,7 +444,7 @@ int main(int argc, char *argv[])
bcVisitor.addBC(noSlipBC);
bcVisitor.addBC(inflowConcreteBC);
bcVisitor.addBC(inflowAirBC1);
- //bcVisitor.addBC(outflowBC);
+ bcVisitor.addBC(outflowBC);
// SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
// SPtr<LBMKernel> kernel = make_shared<CumulantK17LBMKernel>();
@@ -484,7 +484,7 @@ int main(int argc, char *argv[])
string geoPath = "d:/Projects/TRR277/Project/WP4/NozzleGeo";
- string outputPath = "d:/temp/NozzleFlowTest_SharpInterface_With_Part_without_gravity_0.6_2";
+ string outputPath = "d:/temp/NozzleFlowTest_SharpInterface_With_Part_fraction_0.23";
UbSystem::makeDirectory(outputPath);
UbSystem::makeDirectory(outputPath + "/liggghts");
@@ -525,7 +525,7 @@ int main(int argc, char *argv[])
//LatticeDecomposition lDec((g_maxX1 - g_minX1) / dx, (g_maxX2 - g_minX2) / dx, (g_maxX3 - g_minX3) / dx, wrapper.lmp, grid);
SPtr<UbScheduler> lScheduler = make_shared<UbScheduler>(1);
- SPtr<LiggghtsCouplingSimulationObserver> lcSimulationObserver = make_shared<LiggghtsCouplingSimulationObserver>(grid, lScheduler, comm, wrapper, demSubsteps, units);
+ SPtr<LiggghtsCouplingSimulationObserver> lcSimulationObserver = make_shared<LiggghtsCouplingSimulationObserver>(grid, lScheduler, comm, wrapper, demSubsteps, unitsAir);
//SPtr<Grid3DVisitor> partVisitor = make_shared<LiggghtsPartitioningGridVisitor>(std::ceil((g_maxX1 - g_minX1) / dx), std::ceil((g_maxX2 - g_minX2) / dx), std::ceil((g_maxX3 - g_minX3) / dx), wrapper.lmp);
SPtr<Grid3DVisitor> partVisitor = make_shared<LiggghtsPartitioningGridVisitor>(26,26,420, wrapper.lmp);
@@ -598,7 +598,7 @@ int main(int argc, char *argv[])
///////////////////////////////////////////////////////////
// inflow
//GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, 0.20105, -1.30181 + 0.0005, 0.390872 - 0.00229, 0.23, 0.013));
- GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, 0.175, -1.30181 + 0.0005, 0.390872 - 0.00229, 0.23, 0.013));
+ GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, 0.175, -1.30181 + 0.0005, 0.390872 - 0.00229, 0.21, 0.013));
if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), outputPath + "/geo/geoInflow", WbWriterVtkXmlBinary::getInstance());
SPtr<D3Q27Interactor> intrInflow = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, inflowConcreteBC, Interactor3D::SOLID));
///////////////////////////////////////////////////////////
@@ -674,7 +674,7 @@ int main(int argc, char *argv[])
intHelper.addInteractor(intrFluidArea);
intHelper.addInteractor(intrNozzleVolcanNozzle2);
// intHelper.addInteractor(intrBox);
- //intHelper.addInteractor(intrInflow);
+ intHelper.addInteractor(intrInflow);
// intHelper.addInteractor(intrAirInflow);
intHelper.addInteractor(intAirInlet1);
intHelper.addInteractor(intAirInlet2);
@@ -711,8 +711,8 @@ int main(int argc, char *argv[])
double x1c = -1.31431 + R;
double x2c = 0.375582 + R;
- double Ri = 5; //0.05;
- double x3c = 0.1225 + Ri;
+ double Ri = 5;
+ double x3c = 0.136 + Ri;
//R = 0.2 - 0.145; // 0.078-0.04; // 0.2;
@@ -726,7 +726,7 @@ int main(int argc, char *argv[])
fct1.DefineConst("interfaceThickness", interfaceThickness * dx);
MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
- //initVisitor.setPhi(fct1);
+ initVisitor.setPhi(fct1);
grid->accept(initVisitor);
//InitDistributionsBlockVisitor initVisitor;
diff --git a/apps/cpu/Nozzle/nozzleSinglePhase.cpp b/apps/cpu/Nozzle/nozzleSinglePhase.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ed89cb57374e543630134296ae4ba24d02481db5
--- /dev/null
+++ b/apps/cpu/Nozzle/nozzleSinglePhase.cpp
@@ -0,0 +1,798 @@
+#include <iostream>
+#include <memory>
+#include <string>
+
+#include "VirtualFluids.h"
+
+#include "LiggghtsCoupling/LiggghtsCoupling.h"
+
+#include "MultiphaseFlow/MultiphaseFlow.h"
+
+#include "NonNewtonianFluids/NonNewtonianFluids.h"
+
+using namespace std;
+
+int main(int argc, char *argv[])
+{
+ //Sleep(30000);
+
+ try {
+
+ std::shared_ptr<vf::mpi::Communicator> comm = vf::mpi::MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+ // bounding box
+ // double g_minX1 = -1341.81e-3;
+ // double g_minX2 = 348.087e-3;
+ // double g_minX3 = -210e-3;
+
+ // double g_maxX1 = -1260.81e-3;
+ // double g_maxX2 = 429.087e-3;
+ // double g_maxX3 = 214.5e-3;
+
+ // double g_minX1 = -1341.81e-3 + 10e-3;
+ // double g_minX2 = 0.360872;
+ // double g_minX3 = -210e-3;
+
+ // double g_maxX1 = -1260.81e-3 - 10e-3;
+ // double g_maxX2 = 0.416302;
+ // double g_maxX3 = 210e-3;
+
+ // int blockNX[3] = { 10, 10, 10 };
+
+ double g_minX1 = -1.31431;
+ double g_minX2 = 0.375582;
+ double g_minX3 = -0.21; //-210e-3 - 0.2 - 6e-3; //- 1e-3;
+
+ double g_maxX1 = -1.28831;
+ double g_maxX2 = 0.401582;
+ double g_maxX3 = 0.175;//0.21;
+
+ int blockNX[3] = { 26, 26, 35 };
+
+ double dx = 1e-3;
+
+ double uLB_ref = 0.0001;
+ // double rhoLB = 0.0;
+
+ // concrete
+ double d_part = 1e-3;
+ double V = 0.4*10.; // flow rate [m^3/h]
+ double D = 0.026; // shotcrete inlet diameter [m]
+ double R = D / 2.0; // radius [m]
+ double A = UbMath::PI * R * R;
+ double u = V / 3600 / A;
+ double muConcrete = 2.1133054011798826; // [Pa s]
+ double rhoAir = 1.2041; // [kg/m^3]
+ double tau0 = 715.218181094648; // Pa
+ double rhoConcrete = 2400; // [kg/m^3]
+ double nu = muConcrete / rhoConcrete;
+
+ // double Re_D = d_part * u / nu;
+ // if (myid == 0) UBLOG(logINFO, "Re_D = " << Re_D);
+ //
+ SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(d_part, 1., 2400, d_part / dx, uLB_ref);
+ if (myid == 0) std::cout << units->toString() << std::endl;
+
+ double interfaceThickness = 3; // 4.096;
+ double sigma = 0.3; //0.03;
+ double Re = rhoConcrete * u * d_part / muConcrete;
+ double We = rhoConcrete * u * u * d_part / sigma;
+
+ double u_LB_con = u * units->getFactorVelocityWToLb();
+ double nu_h_LB = nu * units->getFactorViscosityWToLb(); // uLB_ref * d_part * units->getFactorLentghWToLb() / Re;
+ double nu_l_LB = 0; // = nu_h_LB;
+
+ double rho_h_LB = 1;
+
+ // surface tension
+ double sigma_LB = rho_h_LB *u_LB_con *u_LB_con *d_part * units->getFactorLentghWToLb() / We;
+
+ // LBMReal dLB = 0; // = length[1] / dx;
+ LBMReal rhoLB = 0.0;
+ // LBMReal nuLB = nu_l; //(uLB_ref*dLB) / Re;
+
+ double beta = 12.0 * sigma_LB / interfaceThickness;
+ double kappa = 1.5 * interfaceThickness * sigma_LB;
+
+ double phiL = 0.0;
+ double phiH = 1.0;
+ double tauH = 0.6; // Phase - field Relaxation
+ double mob = 0.02; // Mobility
+ // double nuL = 1e-2;
+ // double nuG = 0.015811388300841892;
+ double densityRatio = rhoConcrete / rhoAir;
+ // double sigma_old = 1.0850694444444444e-06;
+ //
+ // double beta_old = 12.0 * sigma / interfaceThickness;
+ // double kappa_old = 1.5 * interfaceThickness * sigma;
+
+ double theta = 110; // contact angle
+
+ // https://civilsir.com/density-of-cement-sand-and-aggregate-in-kg-m3-list-of-material-density/
+
+ // SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(r_p, 1.480, 2060, r_p/dx);
+ // SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(r_p, LBMUnitConverter::AIR_20C, r_p / dx);
+ // SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(d_part, 1., 1000, d_part / dx, std::abs(uLB_ref));
+ // SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(d_part, 1., 1000, d_part / dx, std::abs(uLB_ref));
+ // SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(d_part, 1., 2400, d_part / dx, uRef);
+
+ double Bm = (tau0 * d_part) / (muConcrete * u);
+ double tau0_LB = Bm * nu_h_LB * u_LB_con / (d_part * units->getFactorLentghWToLb());
+
+ SPtr<Rheology> thix = Rheology::getInstance();
+ thix->setYieldStress(tau0_LB);
+
+ if (myid == 0) VF_LOG_INFO("Yield stress = {} Pa", tau0);
+ if (myid == 0) VF_LOG_INFO("Yield stress LB = {} ", tau0_LB);
+
+ SPtr<BC> noSlipBC(new NoSlipBC());
+ noSlipBC->setBCStrategy(SPtr<BCStrategy>(new NoSlipBCStrategy()));
+ //SPtr<BC> noSlipBC(new NoSlipBC());
+ //noSlipBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseNoSlipBCStrategy()));
+
+ // concrete inflow boundary condition
+ mu::Parser fct;
+ fct.SetExpr("U");
+ fct.DefineConst("U", -u_LB_con);
+ if (myid == 0) VF_LOG_INFO("Concrete inflow velocity = {} m/s", u);
+ if (myid == 0) VF_LOG_INFO("Concrete inflow velocity = {} dx/dt", u_LB_con);
+ if (myid == 0) VF_LOG_INFO("Concrete Re = {}", Re);
+
+ // // Å tigler, J. (2014). Analytical velocity profile in tube for laminar and turbulent flow. Engineering
+ // // Mechanics, 21(6), 371-379.
+ // double cx1 = -1.31431 + R;
+ // double cx2 = 0.375582 + R;
+ // //double cx3 = 0.20105 + R;
+ // double L = g_maxX1 - g_minX1;
+ // double p_concrete = 1e5; // Pa = 1 Bar
+ // double p1 = p_concrete * units->getFactorPressureWToLb();
+ // double p2 = 0.0;
+ // double drhoLB = 1.0 + rhoLB;
+ // double muLB = drhoLB * nuLB;
+ // double N = R * R / 2 * muLB * uLB_ref * (p1 - p2) / L - 3;
+
+ // // mu::Parser fct;
+ // fct.SetExpr("U*(1-(((((x2-y0)^2+(x1-x0)^2)^0.5)/R)^NplusOne))");
+ // fct.DefineConst("x0", cx1);
+ // fct.DefineConst("y0", cx2);
+ // //fct.DefineConst("z0", cx3);
+ // fct.DefineConst("R", R);
+ // fct.DefineConst("U", uLB_ref * ((N + 3) / (N + 1)));
+ // fct.DefineConst("NplusOne", N + 1.0);
+
+ SPtr<BC> inflowConcreteBC(new VelocityBC(false, false, true, fct, 0, BCFunction::INFCONST));
+ inflowConcreteBC->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowConcreteBC(new MultiphaseVelocityBC(false, false, true, fct, phiH, 0, BCFunction::INFCONST));
+ //inflowConcreteBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+
+ // air inflow boundary condition
+ // Å tigler, J. (2014). Analytical velocity profile in tube for laminar and turbulent flow. Engineering
+ // Mechanics, 21(6), 371-379.
+ // SPtr<LBMUnitConverter> unitsAir = std::make_shared<LBMUnitConverter>(d_part, LBMUnitConverter::AIR_20C, d_part / dx);
+ // SPtr<LBMUnitConverter> unitsAir = std::make_shared<LBMUnitConverter>(d_part, 1., 1.2041, d_part / dx, uLB_ref);
+ // double V = 40; // flow rate [m^3/h]
+ // double D = 0.0166; // air inlet diameter [m]
+ // double R = D / 2.0; // radius [m]
+ // double A = UbMath::PI * R * R;
+ // double u = V / 3600 / A;
+ // double uLB_ref = u * unitsAir->getFactorVelocityWToLb();
+ //// double cx1 = -1.2788 + R;
+ // double cx2 = 0.3803 + R;
+ // double cx3 = 0.1517 + R;
+ // double L = g_maxX1 - g_minX1;
+ // double p_air = 7e5; // Pa = 7 Bar
+ // double p1 = p_air;
+ // double p2 = 0.0;
+ // double mu = 17.2e-6; // Pa s, air 20° C
+ // double N = R * R / 2 * mu * u * (p1 - p2) / L - 3;
+ // if (myid == 0) VF_LOG_INFO("Air inflow velocity = {} m/s", u);
+ // if (myid == 0) VF_LOG_INFO("Air inflow velocity = {} dx/dt", uLB_ref);
+ //
+
+ // double nu = mu / rhoConcrete;
+ // double Re = d_part * u / nu;
+ // if (myid == 0) VF_LOG_INFO("Re_air = {}", Re);
+
+ // double nuLB = d_part * unitsAir->getFactorLentghWToLb() * uLB_ref / Re;
+ // if (myid == 0) VF_LOG_INFO("nuLB_air = {}", nuLB);
+ // nu_l_LB = nuLB;
+
+ SPtr<LBMUnitConverter> unitsAir = std::make_shared<LBMUnitConverter>(d_part, 1., 1.2041, d_part / dx, uLB_ref);
+ //SPtr<LBMUnitConverter> unitsAir = std::make_shared<LBMUnitConverter>(d_part, LBMUnitConverter::AIR_20C, d_part / dx);
+ double V_air = 40./6.; // flow rate [m^3/h] //10.
+ double D_air = 0.00553; // air inlet diameter [m]
+ double R_air = D_air / 2.0; // radius [m]
+ double A_air = UbMath::PI * (R_air * R_air);
+ double u_air = V_air / 3600 / A_air;
+ double uLB_air = u_air * unitsAir->getFactorVelocityWToLb();
+ // double cx1 = -1.2788 + R;
+ double cx2 = 0.385822 + R_air;
+ double cx3 = 0.135562 + R_air;
+ double L_air = 0.00747;
+ double p_air = 7e5; // Pa = 7 Bar
+ double p1 = p_air;
+ double p2 = 1e5;
+ double mu_air = 17.2e-6; // Pa s, air 20° C
+ double rho_air = 1.2041; // [kg/m^3]
+ double Re_inlet = D_air * u_air * rho_air / mu_air;
+ double lambda = 0.3164 / pow(Re_inlet, 0.25);
+ double deltaP = (lambda / (2. * R_air)) * (rho_air * pow(u_air, 2) / 2.); // Darcy friction factor (Rohrreibungszahl)
+ double N = pow(R_air, 2) / (2. * mu_air * u_air) * deltaP - 3.;
+ // double N = R_air * R_air / 2 * mu_air * u_air * (p1 - p2) / L_air - 3;
+ if (myid == 0) VF_LOG_INFO("Air inflow velocity = {} m/s", u_air);
+ if (myid == 0) VF_LOG_INFO("Air inflow velocity = {} dx/dt", uLB_air);
+
+ double nu_air = mu_air / rho_air;
+ double Re_air = d_part * u_air / nu_air;
+ if (myid == 0) VF_LOG_INFO("Air Re = {}", Re_air);
+
+ double nuLB_air = nu_air * unitsAir->getFactorViscosityWToLb(); // d_part * unitsAir->getFactorLentghWToLb() * uLB_air / Re_air;
+ if (myid == 0) VF_LOG_INFO("nuLB_air = {}", nuLB_air);
+ nu_l_LB = nuLB_air;
+
+ if (myid == 0) VF_LOG_INFO("nu_h = {}", nu_h_LB);
+ if (myid == 0) VF_LOG_INFO("nu_l = {}", nu_l_LB);
+ if (myid == 0) VF_LOG_INFO("sigma_LB = {}", sigma_LB);
+
+ double p_air_LB = p_air * unitsAir->getFactorPressureWToLb();
+ if (myid == 0) VF_LOG_INFO("p_air_LB = {}", p_air_LB);
+
+ // mu::Parser fctVx1;
+ ////fctVx1.SetExpr("U");
+ ////fctVx1.DefineConst("U", uLB_air);
+ // mu::Parser fctVx2;
+ // fctVx2.SetExpr("U");
+ // fctVx2.DefineConst("U", 0);
+ // mu::Parser fctVx3;
+ ////fctVx3.SetExpr("U");
+ ////fctVx3.DefineConst("U", -uLB_air);
+
+ double cx1 = 0;
+ double alpha = 0;
+ double gamma = 0;
+ double U = uLB_air;// * ((N + 3.) / (N + 1.));
+
+ mu::Parser fctVx1;
+ //fctVx1.SetExpr("U*cos(alpha*_pi/180)*(1-(((((x1-x0)^2+(x2-y0)^2+(x3-z0)^2)^0.5)/R)^NplusOne))");
+ //fctVx1.SetExpr("U*cos(alpha*_pi/180)*(1-(((((x1-x0)^2+(x2-y0)^2+(x3-z0)^2)^0.5)/R)))");
+ // fctVx1.SetExpr("(((x1-x0)^2+(x2-y0)^2+(x3-z0)^2)^0.5/R)^NplusOne");
+ fctVx1.SetExpr("U*cos(alpha*_pi/180)");
+ fctVx1.DefineConst("x0", cx1);
+ fctVx1.DefineConst("y0", cx2);
+ fctVx1.DefineConst("z0", cx3);
+ fctVx1.DefineConst("R", R_air);
+ fctVx1.DefineConst("U", U); //* ((N + 3.) / (N + 1.)));
+ fctVx1.DefineConst("NplusOne", N + 1.0);
+ fctVx1.DefineConst("alpha", alpha);
+ fctVx1.DefineConst("gamma", gamma);
+
+ mu::Parser fctVx2;
+ //fctVx2.SetExpr("U*sin(alpha*_pi/180)*(1-(((((x1-x0)^2+(x2-y0)^2+(x3-z0)^2)^0.5)/R)^NplusOne))");
+ //fctVx2.SetExpr("U*sin(alpha*_pi/180)*(1-(((((x1-x0)^2+(x2-y0)^2+(x3-z0)^2)^0.5)/R)))");
+ fctVx2.SetExpr("U*sin(alpha*_pi/180)");
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("y0", cx2);
+ fctVx2.DefineConst("z0", cx3);
+ fctVx2.DefineConst("R", R_air);
+ fctVx2.DefineConst("U", U); //* ((N + 3.) / (N + 1.)));
+ fctVx2.DefineConst("NplusOne", N + 1.0);
+ fctVx2.DefineConst("alpha", alpha);
+
+ mu::Parser fctVx3;
+ //fctVx3.SetExpr("U*(1-(((((x1-x0)^2+(x2-y0)^2+(x3-z0)^2)^0.5)/R)^NplusOne))");
+ //fctVx3.SetExpr("U*cos(alpha*_pi/180)*(1-(((((x1-x0)^2+(x2-y0)^2+(x3-z0)^2)^0.5)/R)))");
+ fctVx3.SetExpr("U");
+ fctVx3.DefineConst("x0", cx1);
+ fctVx3.DefineConst("y0", cx2);
+ fctVx3.DefineConst("z0", cx3);
+ fctVx3.DefineConst("R", R_air);
+ fctVx3.DefineConst("U", -U); //* ((N + 3.) / (N + 1.)));
+ fctVx3.DefineConst("NplusOne", N + 1.0);
+ fctVx3.DefineConst("alpha", alpha);
+ fctVx3.DefineConst("gamma", gamma);
+
+ // SPtr<BC> inflowAirBC1(new VelocityBC(true, false, false, fct, 0, BCFunction::INFCONST));
+ // inflowAirBC1->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ // t = U * sin(alpha * _pi / 180) * (1 - (((((x1 - x0) ^ 2 + (x2 - y0) ^ 2 + (x3 - z0) ^ 2) ^ 0.5) / R) ^ NplusOne));
+ cx1 = -1.31416;
+ cx2 = 0.388684;
+ cx3 = 0.138177;
+ alpha = 0;
+ gamma = 225;
+ fctVx1.DefineConst("x0", cx1);
+ fctVx1.DefineConst("y0", cx2);
+ fctVx1.DefineConst("z0", cx3);
+ fctVx1.DefineConst("alpha", alpha);
+ fctVx1.DefineConst("gamma", gamma);
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("y0", cx2);
+ fctVx2.DefineConst("z0", cx3);
+ fctVx2.DefineConst("alpha", alpha);
+ fctVx2.DefineConst("gamma", gamma);
+ fctVx3.DefineConst("x0", cx1);
+ fctVx3.DefineConst("y0", cx2);
+ fctVx3.DefineConst("z0", cx3);
+ fctVx3.DefineConst("alpha", alpha);
+ fctVx3.DefineConst("gamma", gamma);
+
+ SPtr<BC> inflowAirBC1(new VelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, 0, BCFunction::INFCONST));
+ inflowAirBC1->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowAirBC1(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
+ //inflowAirBC1->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+
+ fctVx1.DefineVar("x1", &cx1);
+ fctVx1.DefineVar("x2", &cx2);
+ fctVx1.DefineVar("x3", &cx3);
+ fctVx2.DefineVar("x1", &cx1);
+ fctVx2.DefineVar("x2", &cx2);
+ fctVx2.DefineVar("x3", &cx3);
+ fctVx3.DefineVar("x1", &cx1);
+ fctVx3.DefineVar("x2", &cx2);
+ fctVx3.DefineVar("x3", &cx3);
+
+ if (myid == 0) {
+
+ VF_LOG_INFO("fctVx1 = {}", fctVx1.Eval());
+ VF_LOG_INFO("fctVx2 = {}", fctVx2.Eval());
+ VF_LOG_INFO("fctVx3 = {}", fctVx3.Eval());
+ VF_LOG_INFO("N = {}", N);
+ VF_LOG_INFO("NplusOne = {}", N + 1.0);
+ // return 0;
+ }
+ cx1 = -1.31303;
+ cx2 = 0.377234;
+ cx3 = 0.138174;
+ alpha = 60;
+ fctVx1.DefineConst("x0", cx1);
+ fctVx1.DefineConst("y0", cx2);
+ fctVx1.DefineConst("z0", cx3);
+ fctVx1.DefineConst("alpha", alpha);
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("y0", cx2);
+ fctVx2.DefineConst("z0", cx3);
+ fctVx2.DefineConst("alpha", alpha);
+ fctVx3.DefineConst("x0", cx1);
+ fctVx3.DefineConst("y0", cx2);
+ fctVx3.DefineConst("z0", cx3);
+ SPtr<BC> inflowAirBC2(new VelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, 0, BCFunction::INFCONST));
+ inflowAirBC2->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowAirBC2(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
+ //inflowAirBC2->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+
+ cx1 = -1.2948374155694822;
+ cx2 = 0.37733728717266285;
+ cx3 = 0.13840460401111598;
+ alpha = 120;
+ fctVx1.DefineConst("x0", cx1);
+ fctVx1.DefineConst("y0", cx2);
+ fctVx1.DefineConst("z0", cx3);
+ fctVx1.DefineConst("alpha", alpha);
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("y0", cx2);
+ fctVx2.DefineConst("z0", cx3);
+ fctVx2.DefineConst("alpha", alpha);
+ fctVx3.DefineConst("x0", cx1);
+ fctVx3.DefineConst("y0", cx2);
+ fctVx3.DefineConst("z0", cx3);
+ SPtr<BC> inflowAirBC3(new VelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, 0, BCFunction::INFCONST));
+ inflowAirBC3->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowAirBC3(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
+ //inflowAirBC3->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+
+ cx1 = -1.28847;
+ cx2 = 0.3885;
+ cx3 = 0.1385;
+ alpha = 180;
+ fctVx1.DefineConst("x0", cx1);
+ fctVx1.DefineConst("y0", cx2);
+ fctVx1.DefineConst("z0", cx3);
+ fctVx1.DefineConst("alpha", alpha);
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("y0", cx2);
+ fctVx2.DefineConst("z0", cx3);
+ fctVx2.DefineConst("alpha", alpha);
+ fctVx3.DefineConst("x0", cx1);
+ fctVx3.DefineConst("y0", cx2);
+ fctVx3.DefineConst("z0", cx3);
+ SPtr<BC> inflowAirBC4(new VelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, 0, BCFunction::INFCONST));
+ inflowAirBC4->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowAirBC4(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
+ //inflowAirBC4->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+
+ cx1 = -1.294771417778694;
+ cx2 = 0.399787947463142;
+ cx3 = 0.1383429692754194;
+ alpha = 240;
+ fctVx1.DefineConst("x0", cx1);
+ fctVx1.DefineConst("y0", cx2);
+ fctVx1.DefineConst("z0", cx3);
+ fctVx1.DefineConst("alpha", alpha);
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("y0", cx2);
+ fctVx2.DefineConst("z0", cx3);
+ fctVx2.DefineConst("alpha", alpha);
+ fctVx3.DefineConst("x0", cx1);
+ fctVx3.DefineConst("y0", cx2);
+ fctVx3.DefineConst("z0", cx3);
+ SPtr<BC> inflowAirBC5(new VelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, 0, BCFunction::INFCONST));
+ inflowAirBC5->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowAirBC5(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
+ //inflowAirBC5->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+
+ cx1 = -1.3077338898450492;
+ cx2 = 0.3998516560596088;
+ cx3 = 0.13843501416896437;
+ alpha = 300;
+ fctVx1.DefineConst("x0", cx1);
+ fctVx1.DefineConst("y0", cx2);
+ fctVx1.DefineConst("z0", cx3);
+ fctVx1.DefineConst("alpha", alpha);
+ fctVx2.DefineConst("x0", cx1);
+ fctVx2.DefineConst("y0", cx2);
+ fctVx2.DefineConst("z0", cx3);
+ fctVx2.DefineConst("alpha", alpha);
+ fctVx3.DefineConst("x0", cx1);
+ fctVx3.DefineConst("y0", cx2);
+ fctVx3.DefineConst("z0", cx3);
+ SPtr<BC> inflowAirBC6(new VelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, 0, BCFunction::INFCONST));
+ inflowAirBC6->setBCStrategy(SPtr<BCStrategy>(new VelocityBCStrategy()));
+ //SPtr<BC> inflowAirBC6(new MultiphaseVelocityBC(true, true, true, fctVx1, fctVx2, fctVx3, phiL, 0, BCFunction::INFCONST));
+ //inflowAirBC6->setBCStrategy(SPtr<BCStrategy>(new MultiphaseVelocityBCStrategy()));
+
+ // Pressure BC for air inlet
+ // SPtr<BC> inflowAirBC1(new DensityBC(p_air_LB));
+ // inflowAirBC1->setBCStrategy(SPtr<BCStrategy>(new MultiphasePressureBCStrategy()));
+
+ SPtr<BC> outflowBC(new DensityBC(rhoLB));
+ outflowBC->setBCStrategy(SPtr<BCStrategy>(new NonEqDensityBCStrategy()));
+ //outflowBC->setBCStrategy(SPtr<BCStrategy>(new MultiphasePressureBCStrategy()));
+
+ // SPtr<BC> outflowBC(new DensityBC(rhoLB));
+ //outflowBC->setBCStrategy(SPtr<BCStrategy>(new MultiphaseNonReflectingOutflowBCStrategy()));
+ //////////////////////////////////////////////////////////////////////////////////
+ // BC visitor
+ BoundaryConditionsBlockVisitor bcVisitor;
+ //MultiphaseBoundaryConditionsBlockVisitor bcVisitor;
+ bcVisitor.addBC(noSlipBC);
+ //bcVisitor.addBC(inflowConcreteBC);
+ bcVisitor.addBC(inflowAirBC1);
+ bcVisitor.addBC(outflowBC);
+
+ // SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
+ //SPtr<LBMKernel> kernel = make_shared<CumulantK17LBMKernel>();
+ // SPtr<LBMKernel> kernel = make_shared<MultiphaseTwoPhaseFieldsPressureFilterLBMKernel>();
+ // SPtr<LBMKernel> kernel = make_shared<MultiphaseSimpleVelocityBaseExternalPressureLBMKernel>();
+ //SPtr<LBMKernel> kernel = make_shared<MultiphaseSharpInterfaceLBMKernel>();
+ SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
+ //SPtr<LBMKernel> kernel = make_shared<IBsharpInterfaceLBMKernel>();
+
+ kernel->setWithForcing(false);
+ kernel->setForcingX1(0.0);
+ kernel->setForcingX2(0.0);
+ kernel->setForcingX3(0.0);
+
+ kernel->setPhiL(phiL);
+ kernel->setPhiH(phiH);
+ kernel->setPhaseFieldRelaxation(tauH);
+ kernel->setMobility(mob);
+ kernel->setInterfaceWidth(interfaceThickness);
+
+ kernel->setCollisionFactorMultiphase(nu_h_LB, nu_l_LB);
+ kernel->setDensityRatio(densityRatio);
+ kernel->setMultiphaseModelParameters(beta, kappa);
+ kernel->setContactAngle(theta);
+ kernel->setSigma(sigma_LB);
+
+ SPtr<BCSet> bcProc = make_shared<BCSet>();
+ kernel->setBCSet(bcProc);
+
+ SPtr<Grid3D> grid = make_shared<Grid3D>(comm);
+ grid->setPeriodicX1(false);
+ grid->setPeriodicX2(false);
+ grid->setPeriodicX3(false);
+ grid->setDeltaX(dx);
+ grid->setBlockNX(blockNX[0], blockNX[1], blockNX[2]);
+ grid->setGhostLayerWidth(1);
+
+ string geoPath = "d:/Projects/TRR277/Project/WP4/NozzleGeo";
+
+ string outputPath = "d:/temp/NozzleFlowTest_SinglePhase_With_Part_fraction_0.13_vel_0.2_ogl";
+ UbSystem::makeDirectory(outputPath);
+ UbSystem::makeDirectory(outputPath + "/liggghts");
+
+ // if (myid == 0) {
+ // stringstream logFilename;
+ // logFilename << outputPath + "/logfile" + UbSystem::toString(UbSystem::getTimeStamp()) + ".txt";
+ // UbLog::output_policy::setStream(logFilename.str());
+ // }
+
+ /////////////////////////////////////////////////////////////////////
+ //LIGGGHTS things
+ /////////////////////////////////////////////////////////////////////
+ string inFile1 = "d:/Projects/TRR277/Project/WP4/Config/in.nozzle";
+ // string inFile2 = "d:/Projects/VirtualFluids_LIGGGHTS_coupling/apps/cpu/LiggghtsApp/in2.lbdem";
+ MPI_Comm mpi_comm = *(MPI_Comm *)(comm->getNativeCommunicator());
+ LiggghtsCouplingWrapper wrapper(argv, mpi_comm);
+
+ double v_frac = 0.1;
+ double dt_phys = units->getFactorTimeLbToW();
+ int demSubsteps = 10;
+ double dt_dem = dt_phys / (double)demSubsteps;
+ int vtkSteps = 1000;
+ string demOutDir = outputPath + "/liggghts";
+
+ // wrapper.execCommand("echo none");
+
+ // wrapper.execFile((char*)inFile1.c_str());
+
+ //// set timestep and output directory
+ wrapper.setVariable("t_step", dt_dem);
+ wrapper.setVariable("dmp_stp", vtkSteps * demSubsteps);
+ wrapper.setVariable("dmp_dir", demOutDir);
+
+ wrapper.execFile((char *)inFile1.c_str());
+ //wrapper.runUpto(demSubsteps - 1);
+ // wrapper.runUpto(1000);
+
+ //LatticeDecomposition lDec((g_maxX1 - g_minX1) / dx, (g_maxX2 - g_minX2) / dx, (g_maxX3 - g_minX3) / dx, wrapper.lmp, grid);
+
+ SPtr<UbScheduler> lScheduler = make_shared<UbScheduler>(1);
+ SPtr<LiggghtsCouplingSimulationObserver> lcSimulationObserver = make_shared<LiggghtsCouplingSimulationObserver>(grid, lScheduler, comm, wrapper, demSubsteps, unitsAir);
+ //SPtr<Grid3DVisitor> partVisitor = make_shared<LiggghtsPartitioningGridVisitor>(std::ceil((g_maxX1 - g_minX1) / dx), std::ceil((g_maxX2 - g_minX2) / dx), std::ceil((g_maxX3 - g_minX3) / dx), wrapper.lmp);
+ SPtr<Grid3DVisitor> partVisitor = make_shared<LiggghtsPartitioningGridVisitor>(26,26,420, wrapper.lmp);
+
+ /////////////////////////////////////////////////////////////////////
+ /////////////////////////////////////////////////////////////////////
+
+ SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased, vf::lbm::dir::DIR_MMM, MetisPartitioner::RECURSIVE));
+
+ SPtr<GbObject3D> gridCube = make_shared<GbCuboid3D>(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3);
+ if (myid == 0) GbSystem3D::writeGeoObject(gridCube.get(), outputPath + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ // geo
+ //////////////////////////////////////////////////////////
+ int accuracy = Interactor3D::EDGES;
+ ///////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshNozzleAirDistributor = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAirDistributor:start");
+ meshNozzleAirDistributor->readMeshFromSTLFileASCII(geoPath + "/01_Nozzle_Air_Distributor.stl", false);
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAirDistributor:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleAirDistributor.get(), outputPath + "/geo/meshNozzleAirDistributor", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrNozzleAirDistributor = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleAirDistributor, grid, noSlipBC, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshNozzleAirInlet = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAirInlet:start");
+ meshNozzleAirInlet->readMeshFromSTLFileASCII(geoPath + "/02_Nozzle_Air_Inlet.stl", false);
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAirInlet:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleAirInlet.get(), outputPath + "/geo/meshNozzleAirInlet", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrNozzleAirInlet = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleAirInlet, grid, noSlipBC, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshNozzleSpacer = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleSpacer:start");
+ meshNozzleSpacer->readMeshFromSTLFileASCII(geoPath + "/03_Nozzle_Spacer.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleSpacer:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleSpacer.get(), outputPath + "/geo/meshNozzleSpacer", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrNozzleSpacer = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleSpacer, grid, noSlipBC, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshNozzleAccDistributor = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAccDistributor:start");
+ meshNozzleAccDistributor->readMeshFromSTLFileASCII(geoPath + "/04_Nozzle_Acc_Distributor.stl", false);
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAccDistributor:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleAccDistributor.get(), outputPath + "/geo/meshNozzleAccDistributor", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrNozzleAccDistributor = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleAccDistributor, grid, noSlipBC, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshNozzleAccInlet = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAccInlet:start");
+ meshNozzleAccInlet->readMeshFromSTLFileASCII(geoPath + "/05_Nozzle_Acc_Inlet.stl", false);
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleAccInlet:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleAccInlet.get(), outputPath + "/geo/meshNozzleAccInlet", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrNozzleAccInlet = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleAccInlet, grid, noSlipBC, Interactor3D::SOLID, (Interactor3D::Accuracy)accuracy);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshNozzleVolcanNozzle1 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleVolcanNozzle1:start");
+ meshNozzleVolcanNozzle1->readMeshFromSTLFileBinary(geoPath + "/06_1_Nozzle_Volcan_Nozzle.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleVolcanNozzle1:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleVolcanNozzle1.get(), outputPath + "/geo/meshNozzleVolcanNozzle1", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrNozzleVolcanNozzle1 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleVolcanNozzle1, grid, noSlipBC, Interactor3D::SOLID, Interactor3D::EDGES);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshNozzleVolcanNozzle2 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleVolcanNozzle2:start");
+ meshNozzleVolcanNozzle2->readMeshFromSTLFileBinary(geoPath + "/06_2_Nozzle_Volcan_Nozzle.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshNozzleVolcanNozzle2:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshNozzleVolcanNozzle2.get(), outputPath + "/geo/meshNozzleVolcanNozzle2", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrNozzleVolcanNozzle2 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshNozzleVolcanNozzle2, grid, noSlipBC, Interactor3D::SOLID, Interactor3D::POINTS);
+ ///////////////////////////////////////////////////////////
+ // box
+ SPtr<D3Q27Interactor> intrBox = SPtr<D3Q27Interactor>(new D3Q27Interactor(gridCube, grid, noSlipBC, Interactor3D::INVERSESOLID));
+ ///////////////////////////////////////////////////////////
+ // inflow
+ //GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, 0.20105, -1.30181 + 0.0005, 0.390872 - 0.00229, 0.23, 0.013));
+ GbCylinder3DPtr geoInflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, 0.175, -1.30181 + 0.0005, 0.390872 - 0.00229, 0.21, 0.013));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), outputPath + "/geo/geoInflow", WbWriterVtkXmlBinary::getInstance());
+ SPtr<D3Q27Interactor> intrInflow = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, inflowConcreteBC, Interactor3D::SOLID));
+ ///////////////////////////////////////////////////////////
+ // outflow
+ //GbCylinder3DPtr geoOutflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, -0.22, -1.30181 + 0.0005, 0.390872 - 0.00229, -0.21, 0.013));
+ //GbCylinder3DPtr geoOutflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, -0.426, -1.30181 + 0.0005, 0.390872 - 0.00229, -0.415, 0.013));
+ GbCylinder3DPtr geoOutflow(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, g_minX3, -1.30181 + 0.0005, 0.390872 - 0.00229, g_minX3+2.*dx, 0.013));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), outputPath + "/geo/geoOutflow", WbWriterVtkXmlBinary::getInstance());
+ SPtr<D3Q27Interactor> intrOutflow = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflow, grid, outflowBC, Interactor3D::SOLID));
+ ///////////////////////////////////////////////////////////
+ // SPtr<GbTriFaceMesh3D> geoAirInlet = std::make_shared<GbTriFaceMesh3D>();
+ // if (myid == 0) UBLOG(logINFO, "Read Air_Inlet:start");
+ // geoAirInlet->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet.stl", true);
+ // if (myid == 0) UBLOG(logINFO, "Read Air_Inlet:end");
+ // if (myid == 0) GbSystem3D::writeGeoObject(geoAirInlet.get(), outputPath + "/geo/geoAirInlet", WbWriterVtkXmlBinary::getInstance());
+ // SPtr<Interactor3D> intrAirInlet = std::make_shared<D3Q27TriFaceMeshInteractor>(geoAirInlet, grid, inflowAirBC1, Interactor3D::SOLID, Interactor3D::EDGES);
+ /////////////////////////////////////////////////////////////
+ // Fluid area
+ GbCylinder3DPtr geoFluidArea(new GbCylinder3D(-1.30181 + 0.0005, 0.390872 - 0.00229, g_minX3, -1.30181 + 0.0005, 0.390872 - 0.00229, g_maxX3, 0.013));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoFluidArea.get(), outputPath + "/geo/geoFluidArea", WbWriterVtkXmlBinary::getInstance());
+ SPtr<D3Q27Interactor> intrFluidArea = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoFluidArea, grid, noSlipBC, Interactor3D::INVERSESOLID));
+ ///////////////////////////////////////////////////////////
+ ///////////////////////////////////////////////////////////
+ GbCylinder3DPtr geoAirInflow(new GbCylinder3D(-1.31431 - 0.0005, 0.388587, 0.1383275, -1.31431, 0.388587, 0.1383275, 0.002765));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoAirInflow.get(), outputPath + "/geo/geoAirInlet", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intrAirInflow = std::make_shared<D3Q27Interactor>(geoAirInflow, grid, inflowAirBC1, Interactor3D::SOLID, Interactor3D::EDGES);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshAirInlet1 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet1:start");
+ meshAirInlet1->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_1.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet1:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet1.get(), outputPath + "/geo/meshAirInlet1", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intAirInlet1 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet1, grid, inflowAirBC1, Interactor3D::SOLID, Interactor3D::POINTS);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshAirInlet2 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet2:start");
+ meshAirInlet2->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_2.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet2:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet2.get(), outputPath + "/geo/meshAirInlet2", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intAirInlet2 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet2, grid, inflowAirBC2, Interactor3D::SOLID, Interactor3D::POINTS);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshAirInlet3 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet3:start");
+ meshAirInlet3->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_3.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet3:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet3.get(), outputPath + "/geo/meshAirInlet3", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intAirInlet3 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet3, grid, inflowAirBC3, Interactor3D::SOLID, Interactor3D::POINTS);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshAirInlet4 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet4:start");
+ meshAirInlet4->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_4.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet4:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet4.get(), outputPath + "/geo/meshAirInlet4", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intAirInlet4 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet4, grid, inflowAirBC4, Interactor3D::SOLID, Interactor3D::POINTS);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshAirInlet5 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet5:start");
+ meshAirInlet5->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_5.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet5:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet5.get(), outputPath + "/geo/meshAirInlet5", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intAirInlet5 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet5, grid, inflowAirBC5, Interactor3D::SOLID, Interactor3D::POINTS);
+ ///////////////////////////////////////////////////////////
+ SPtr<GbTriFaceMesh3D> meshAirInlet6 = std::make_shared<GbTriFaceMesh3D>();
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet6:start");
+ meshAirInlet6->readMeshFromSTLFileASCII(geoPath + "/Air_Inlet_6.stl", true);
+ if (myid == 0) UBLOG(logINFO, "Read meshAirInlet6:end");
+ if (myid == 0) GbSystem3D::writeGeoObject(meshAirInlet6.get(), outputPath + "/geo/meshAirInlet6", WbWriterVtkXmlBinary::getInstance());
+ SPtr<Interactor3D> intAirInlet6 = std::make_shared<D3Q27TriFaceMeshInteractor>(meshAirInlet6, grid, inflowAirBC6, Interactor3D::SOLID, Interactor3D::POINTS);
+ ///////////////////////////////////////////////////////////
+
+ InteractorsHelper intHelper(grid, partVisitor, false);
+
+ intHelper.addInteractor(intrFluidArea);
+ intHelper.addInteractor(intrNozzleVolcanNozzle2);
+ // intHelper.addInteractor(intrBox);
+ //intHelper.addInteractor(intrInflow);
+ // intHelper.addInteractor(intrAirInflow);
+ intHelper.addInteractor(intAirInlet1);
+ intHelper.addInteractor(intAirInlet2);
+ intHelper.addInteractor(intAirInlet3);
+ intHelper.addInteractor(intAirInlet4);
+ intHelper.addInteractor(intAirInlet5);
+ intHelper.addInteractor(intAirInlet6);
+ intHelper.addInteractor(intrOutflow);
+
+ // intHelper.addInteractor(intrNozzleAirDistributor);
+ // intHelper.addInteractor(intrNozzleAirInlet);
+ // intHelper.addInteractor(intrNozzleSpacer);
+ // intHelper.addInteractor(intrNozzleAccDistributor);
+ // intHelper.addInteractor(intrNozzleAccInlet);
+ // intHelper.addInteractor(intrNozzleVolcanNozzle1);
+
+ intHelper.selectBlocks();
+
+ SPtr<SimulationObserver> ppblocks = make_shared<WriteBlocksSimulationObserver>(grid, SPtr<UbScheduler>(new UbScheduler(1)), outputPath, WbWriterVtkXmlBinary::getInstance(), comm);
+ ppblocks->update(0);
+ ppblocks.reset();
+
+ // if (myid == 0) UBLOG(logINFO, Utilities::toString(grid, comm->getNumberOfProcesses()));
+
+ SetKernelBlockVisitor kernelVisitor(kernel, nu_l_LB, comm->getNumberOfProcesses());
+ //MultiphaseSetKernelBlockVisitor kernelVisitor(kernel, nu_h_LB, nu_l_LB, 1e9, 1);
+ grid->accept(kernelVisitor);
+
+
+ intHelper.setBC();
+
+ // InitDistributionsBlockVisitor initVisitor;
+ // grid->accept(initVisitor);
+
+ double x1c = -1.31431 + R;
+ double x2c = 0.375582 + R;
+ double Ri = 5;
+ double x3c = 0.136 + Ri;
+
+ //R = 0.2 - 0.145; // 0.078-0.04; // 0.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(" 0.5 - 0.5 * tanh(2 * (sqrt((x1 - x1c) ^ 2 + (x2 - x2c) ^ 2 + (x3 - x3c) ^ 2) - radius) / interfaceThickness)");
+ fct1.DefineConst("x1c", x1c);
+ fct1.DefineConst("x2c", x2c);
+ fct1.DefineConst("x3c", x3c);
+ fct1.DefineConst("radius", Ri);
+ fct1.DefineConst("interfaceThickness", interfaceThickness * dx);
+
+ //MultiphaseVelocityFormInitDistributionsBlockVisitor initVisitor;
+ //initVisitor.setPhi(fct1);
+ //grid->accept(initVisitor);
+
+ InitDistributionsBlockVisitor initVisitor;
+ grid->accept(initVisitor);
+
+ // boundary conditions grid
+ {
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ SPtr<WriteBoundaryConditionsSimulationObserver> ppgeo(new WriteBoundaryConditionsSimulationObserver(grid, geoSch, outputPath, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppgeo->update(0);
+ ppgeo.reset();
+ }
+
+ grid->accept(bcVisitor);
+
+ OneDistributionSetConnectorsBlockVisitor setConnsVisitor(comm);
+ //TwoDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ // ThreeDistributionsDoubleGhostLayerSetConnectorsBlockVisitor setConnsVisitor(comm);
+ grid->accept(setConnsVisitor);
+
+ int numOfThreads = 1;
+ omp_set_num_threads(numOfThreads);
+
+ SPtr<UbScheduler> nupsSch = std::make_shared<UbScheduler>(10, 10, 100);
+ SPtr<NUPSCounterSimulationObserver> nupsSimulationObserver = make_shared<NUPSCounterSimulationObserver>(grid, nupsSch, numOfThreads, comm);
+
+ //// write data for visualization of macroscopic quantities
+ SPtr<UbScheduler> visSch(new UbScheduler(vtkSteps));
+ // SPtr<UbScheduler> visSch(new UbScheduler(1, 8700, 8800));
+ // visSch->addSchedule(1, 8700, 8800);
+ //SPtr<WriteSharpInterfaceQuantitiesSimulationObserver> writeMQSimulationObserver(new WriteSharpInterfaceQuantitiesSimulationObserver(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+ //writeMQSimulationObserver->update(0);
+
+ SPtr<WriteMacroscopicQuantitiesSimulationObserver> writeMQSimulationObserver(new WriteMacroscopicQuantitiesSimulationObserver(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(), SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+ writeMQSimulationObserver->update(0);
+
+ int endTime = 10000000;
+ SPtr<Simulation> simulation(new Simulation(grid, lScheduler, endTime));
+ simulation->addSimulationObserver(nupsSimulationObserver);
+ simulation->addSimulationObserver(lcSimulationObserver);
+ simulation->addSimulationObserver(writeMQSimulationObserver);
+
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
+ simulation->run();
+ 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;
+ }
+ return 0;
+}
diff --git a/src/cpu/LiggghtsCoupling/3rdParty/LiggghtsCouplingWrapper.cpp b/src/cpu/LiggghtsCoupling/3rdParty/LiggghtsCouplingWrapper.cpp
index 9be87887a26d654d03dc8a32ed9e456ec352fef2..c745e8a67b5a628321e27ea1132b55fba4bebd95 100644
--- a/src/cpu/LiggghtsCoupling/3rdParty/LiggghtsCouplingWrapper.cpp
+++ b/src/cpu/LiggghtsCoupling/3rdParty/LiggghtsCouplingWrapper.cpp
@@ -70,7 +70,7 @@ void LiggghtsCouplingWrapper::setVariable(char const *name, std::string &value)
void LiggghtsCouplingWrapper::run(int nSteps)
{
std::stringstream cmd;
- cmd << "run " << nSteps;
+ cmd << "run " << nSteps << " pre no post no";
execCommand(cmd);
}
void LiggghtsCouplingWrapper::runUpto(int nSteps)
diff --git a/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.cpp b/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.cpp
index 3879f11c945cbedd2d253cd0544b25a26904bebf..efe93f4c951684e41e6dd1c4d2b81c57c711727f 100644
--- a/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.cpp
+++ b/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.cpp
@@ -134,9 +134,9 @@ void IBcumulantK17LBMKernel::calculate(int step)
}
/////////////////////////////////////
- localDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
- nonLocalDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
- restDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
+ localDistributionsF = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
+ nonLocalDistributionsF = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
+ restDistributionsF = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
SPtr<BCArray3D> bcArray = this->getBCSet()->getBCArray();
@@ -191,35 +191,35 @@ void IBcumulantK17LBMKernel::calculate(int step)
// a b c
//-1 0 1
- LBMReal mfcbb = (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3);
- LBMReal mfbcb = (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3);
- LBMReal mfbbc = (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3);
- LBMReal mfccb = (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3);
- LBMReal mfacb = (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3);
- LBMReal mfcbc = (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3);
- LBMReal mfabc = (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3);
- LBMReal mfbcc = (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3);
- LBMReal mfbac = (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3);
- LBMReal mfccc = (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3);
- LBMReal mfacc = (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3);
- LBMReal mfcac = (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3);
- LBMReal mfaac = (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3);
-
- LBMReal mfabb = (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3);
- LBMReal mfbab = (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3);
- LBMReal mfbba = (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p);
- LBMReal mfaab = (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3);
- LBMReal mfcab = (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3);
- LBMReal mfaba = (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p);
- LBMReal mfcba = (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p);
- LBMReal mfbaa = (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p);
- LBMReal mfbca = (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p);
- LBMReal mfaaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p);
- LBMReal mfcaa = (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p);
- LBMReal mfaca = (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p);
- LBMReal mfcca = (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p);
-
- LBMReal mfbbb = (*this->restDistributions)(x1, x2, x3);
+ LBMReal mfcbb = (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3);
+ LBMReal mfbcb = (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3);
+ LBMReal mfbbc = (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3);
+ LBMReal mfccb = (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3);
+ LBMReal mfacb = (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3);
+ LBMReal mfcbc = (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3);
+ LBMReal mfabc = (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3);
+ LBMReal mfbcc = (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3);
+ LBMReal mfbac = (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3);
+ LBMReal mfccc = (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3);
+ LBMReal mfacc = (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3);
+ LBMReal mfcac = (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3);
+ LBMReal mfaac = (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3);
+
+ LBMReal mfabb = (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3);
+ LBMReal mfbab = (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3);
+ LBMReal mfbba = (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p);
+ LBMReal mfaab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3);
+ LBMReal mfcab = (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3);
+ LBMReal mfaba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p);
+ LBMReal mfcba = (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p);
+ LBMReal mfbaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p);
+ LBMReal mfbca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p);
+ LBMReal mfaaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p);
+ LBMReal mfcaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p);
+ LBMReal mfaca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p);
+ LBMReal mfcca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p);
+
+ LBMReal mfbbb = (*this->restDistributionsF)(x1, x2, x3);
LBMReal f[D3Q27System::ENDF + 1];
LBMReal fEq[D3Q27System::ENDF + 1];
@@ -753,35 +753,35 @@ void IBcumulantK17LBMKernel::calculate(int step)
//! href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017),
//! DOI:10.3390/computation5020019 ]</b></a>
//!
- (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = mfabb;
- (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = mfbab;
- (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = mfbba;
- (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = mfaab;
- (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab;
- (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = mfaba;
- (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba;
- (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa;
- (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca;
- (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa;
- (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa;
- (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca;
- (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca;
-
- (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb;
- (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb;
- (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc;
- (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb;
- (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb;
- (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac;
- (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc;
- (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac;
- (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
-
- (*this->restDistributions)(x1, x2, x3) = mfbbb;
+ (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) = mfabb;
+ (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) = mfbab;
+ (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) = mfbba;
+ (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) = mfaab;
+ (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) = mfcab;
+ (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) = mfaba;
+ (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) = mfcba;
+ (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) = mfbaa;
+ (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) = mfbca;
+ (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) = mfaaa;
+ (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) = mfcaa;
+ (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) = mfaca;
+ (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) = mfcca;
+
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) = mfcbb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) = mfbcb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) = mfbbc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) = mfccb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) = mfacb;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) = mfcbc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) = mfabc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) = mfbcc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) = mfbac;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) = mfccc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) = mfacc;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac;
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
+
+ (*this->restDistributionsF)(x1, x2, x3) = mfbbb;
//////////////////////////////////////////////////////////////////////////
f[vf::lbm::dir::DIR_000] = mfbbb;
@@ -813,91 +813,93 @@ void IBcumulantK17LBMKernel::calculate(int step)
f[vf::lbm::dir::DIR_MPM] = mfaca;
f[vf::lbm::dir::DIR_PPM] = mfcca;
}
- if ((*particleData)(x1, x2, x3)->solidFraction < SOLFRAC_MIN)
- continue;
-
- LBMReal vx1, vx2, vx3, drho;
- D3Q27System::calcCompMacroscopicValues(f, drho, vx1, vx2, vx3);
- D3Q27System::calcCompFeq(fEq, drho, vx1, vx2, vx3);
-
- std::array<double, 3> uPart;
- uPart[0] = (*particleData)(x1, x2, x3)->uPart[0] * (1. + drho);
- uPart[1] = (*particleData)(x1, x2, x3)->uPart[1] * (1. + drho);
- uPart[2] = (*particleData)(x1, x2, x3)->uPart[2] * (1. + drho);
-
- D3Q27System::calcCompFeq(fEqSolid, drho, uPart[0], uPart[1], uPart[2]);
-
- if ((*particleData)(x1, x2, x3)->solidFraction > SOLFRAC_MAX) {
- double const bb0 = fEq[vf::lbm::dir::DIR_000] - fEqSolid[vf::lbm::dir::DIR_000];
- f[vf::lbm::dir::DIR_000] = fPre[vf::lbm::dir::DIR_000] + bb0;
- for (int iPop = D3Q27System::FSTARTDIR; iPop <= D3Q27System::FENDDIR; iPop++) {
- const int iOpp = D3Q27System::INVDIR[iPop];
- double const bb = ((fPre[iOpp] - fEq[iOpp]) - (fPre[iPop] - fEqSolid[iPop]));
- double const bbOpp = ((fPre[iPop] - fEq[iPop]) - (fPre[iOpp] - fEqSolid[iOpp]));
-
-
- f[iPop] = fPre[iPop] + bb;
- f[iOpp] = fPre[iOpp] + bbOpp;
-
- (*particleData)(x1, x2, x3)->hydrodynamicForce[0] -= D3Q27System::DX1[iPop] * (bb - bbOpp);
- (*particleData)(x1, x2, x3)->hydrodynamicForce[1] -= D3Q27System::DX2[iPop] * (bb - bbOpp);
- (*particleData)(x1, x2, x3)->hydrodynamicForce[2] -= D3Q27System::DX3[iPop] * (bb - bbOpp);
- }
- } else { /* particleData.solidFraction < SOLFRAC_MAX */
-//#ifdef LBDEM_USE_WEIGHING
- double const ooo = 1. / omega - 0.5;
- double const B = (*particleData)(x1, x2, x3)->solidFraction * ooo / ((1. - (*particleData)(x1, x2, x3)->solidFraction) + ooo);
-//#else
-// T const B = particleData.solidFraction;
-//#endif
- double const oneMinB = 1. - B;
-
- double const bb0 = fEq[vf::lbm::dir::DIR_000] - fEqSolid[vf::lbm::dir::DIR_000];
- f[vf::lbm::dir::DIR_000] = fPre[vf::lbm::dir::DIR_000] + oneMinB * (f[vf::lbm::dir::DIR_000] - fPre[vf::lbm::dir::DIR_000]) + B * bb0;
-
- for (int iPop = D3Q27System::FSTARTDIR; iPop <= D3Q27System::FENDDIR; iPop++) {
- int const iOpp = D3Q27System::INVDIR[iPop];
- double const bb = B * ((fPre[iOpp] - fEq[iOpp]) - (fPre[iPop] - fEqSolid[iPop]));
- double const bbOpp = B * ((fPre[iPop] - fEq[iPop]) - (fPre[iOpp] - fEqSolid[iOpp]));
-
- f[iPop] = fPre[iPop] + oneMinB * (f[iPop] - fPre[iPop]) + bb;
- f[iOpp] = fPre[iOpp] + oneMinB * (f[iOpp] - fPre[iOpp]) + bbOpp;
-
- (*particleData)(x1, x2, x3)->hydrodynamicForce[0] -= D3Q27System::DX1[iPop] * (bb - bbOpp);
- (*particleData)(x1, x2, x3)->hydrodynamicForce[1] -= D3Q27System::DX2[iPop] * (bb - bbOpp);
- (*particleData)(x1, x2, x3)->hydrodynamicForce[2] -= D3Q27System::DX3[iPop] * (bb - bbOpp);
- }
- } /* if solidFraction > SOLFRAC_MAX */
-
- (*this->restDistributions)(x1, x2, x3) = f[vf::lbm::dir::DIR_000];
-
- (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = f[vf::lbm::dir::DIR_M00];
- (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = f[vf::lbm::dir::DIR_0M0];
- (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = f[vf::lbm::dir::DIR_00M];
- (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = f[vf::lbm::dir::DIR_MM0];
- (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = f[vf::lbm::dir::DIR_PM0];
- (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = f[vf::lbm::dir::DIR_M0M];
- (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = f[vf::lbm::dir::DIR_P0M];
- (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = f[vf::lbm::dir::DIR_0MM];
- (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = f[vf::lbm::dir::DIR_0PM];
- (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = f[vf::lbm::dir::DIR_MMM];
- (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = f[vf::lbm::dir::DIR_PMM];
- (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = f[vf::lbm::dir::DIR_MPM];
- (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = f[vf::lbm::dir::DIR_PPM];
-
- (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = f[vf::lbm::dir::DIR_P00];
- (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = f[vf::lbm::dir::DIR_0P0];
- (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = f[vf::lbm::dir::DIR_00P];
- (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = f[vf::lbm::dir::DIR_PP0];
- (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = f[vf::lbm::dir::DIR_MP0];
- (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = f[vf::lbm::dir::DIR_P0P];
- (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = f[vf::lbm::dir::DIR_M0P];
- (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = f[vf::lbm::dir::DIR_0PP];
- (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = f[vf::lbm::dir::DIR_0MP];
- (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = f[vf::lbm::dir::DIR_PPP];
- (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = f[vf::lbm::dir::DIR_MPP];
- (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = f[vf::lbm::dir::DIR_PMP];
- (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = f[vf::lbm::dir::DIR_MMP];
+ if ((*particleData)(x1, x2, x3)->solidFraction >= SOLFRAC_MIN) {
+ // if ((*particleData)(x1, x2, x3)->solidFraction < SOLFRAC_MIN)
+ // continue;
+ //
+ LBMReal vx1, vx2, vx3, drho;
+ D3Q27System::calcCompMacroscopicValues(f, drho, vx1, vx2, vx3);
+ D3Q27System::calcCompFeq(fEq, drho, vx1, vx2, vx3);
+
+ std::array<double, 3> uPart;
+ uPart[0] = (*particleData)(x1, x2, x3)->uPart[0] * (1. + drho);
+ uPart[1] = (*particleData)(x1, x2, x3)->uPart[1] * (1. + drho);
+ uPart[2] = (*particleData)(x1, x2, x3)->uPart[2] * (1. + drho);
+
+ D3Q27System::calcCompFeq(fEqSolid, drho, uPart[0], uPart[1], uPart[2]);
+
+ if ((*particleData)(x1, x2, x3)->solidFraction > SOLFRAC_MAX) {
+ double const bb0 = fEq[vf::lbm::dir::DIR_000] - fEqSolid[vf::lbm::dir::DIR_000];
+ f[vf::lbm::dir::DIR_000] = fPre[vf::lbm::dir::DIR_000] + bb0;
+ for (int iPop = D3Q27System::FSTARTDIR; iPop <= D3Q27System::FENDDIR; iPop++) {
+ const int iOpp = D3Q27System::INVDIR[iPop];
+ double const bb = ((fPre[iOpp] - fEq[iOpp]) - (fPre[iPop] - fEqSolid[iPop]));
+ double const bbOpp = ((fPre[iPop] - fEq[iPop]) - (fPre[iOpp] - fEqSolid[iOpp]));
+
+
+ f[iPop] = fPre[iPop] + bb;
+ f[iOpp] = fPre[iOpp] + bbOpp;
+
+ (*particleData)(x1, x2, x3)->hydrodynamicForce[0] -= D3Q27System::DX1[iPop] * (bb - bbOpp);
+ (*particleData)(x1, x2, x3)->hydrodynamicForce[1] -= D3Q27System::DX2[iPop] * (bb - bbOpp);
+ (*particleData)(x1, x2, x3)->hydrodynamicForce[2] -= D3Q27System::DX3[iPop] * (bb - bbOpp);
+ }
+ } else { /* particleData.solidFraction < SOLFRAC_MAX */
+ //#ifdef LBDEM_USE_WEIGHING
+ double const ooo = 1. / omega - 0.5;
+ double const B = (*particleData)(x1, x2, x3)->solidFraction * ooo / ((1. - (*particleData)(x1, x2, x3)->solidFraction) + ooo);
+ //#else
+ // T const B = particleData.solidFraction;
+ //#endif
+ double const oneMinB = 1. - B;
+
+ double const bb0 = fEq[vf::lbm::dir::DIR_000] - fEqSolid[vf::lbm::dir::DIR_000];
+ f[vf::lbm::dir::DIR_000] = fPre[vf::lbm::dir::DIR_000] + oneMinB * (f[vf::lbm::dir::DIR_000] - fPre[vf::lbm::dir::DIR_000]) + B * bb0;
+
+ for (int iPop = D3Q27System::FSTARTDIR; iPop <= D3Q27System::FENDDIR; iPop++) {
+ int const iOpp = D3Q27System::INVDIR[iPop];
+ double const bb = B * ((fPre[iOpp] - fEq[iOpp]) - (fPre[iPop] - fEqSolid[iPop]));
+ double const bbOpp = B * ((fPre[iPop] - fEq[iPop]) - (fPre[iOpp] - fEqSolid[iOpp]));
+
+ f[iPop] = fPre[iPop] + oneMinB * (f[iPop] - fPre[iPop]) + bb;
+ f[iOpp] = fPre[iOpp] + oneMinB * (f[iOpp] - fPre[iOpp]) + bbOpp;
+
+ (*particleData)(x1, x2, x3)->hydrodynamicForce[0] -= D3Q27System::DX1[iPop] * (bb - bbOpp);
+ (*particleData)(x1, x2, x3)->hydrodynamicForce[1] -= D3Q27System::DX2[iPop] * (bb - bbOpp);
+ (*particleData)(x1, x2, x3)->hydrodynamicForce[2] -= D3Q27System::DX3[iPop] * (bb - bbOpp);
+ }
+ } /* if solidFraction > SOLFRAC_MAX */
+
+ (*this->restDistributionsF)(x1, x2, x3) = f[vf::lbm::dir::DIR_000];
+
+ (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) = f[vf::lbm::dir::DIR_M00];
+ (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) = f[vf::lbm::dir::DIR_0M0];
+ (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) = f[vf::lbm::dir::DIR_00M];
+ (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) = f[vf::lbm::dir::DIR_MM0];
+ (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) = f[vf::lbm::dir::DIR_PM0];
+ (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) = f[vf::lbm::dir::DIR_M0M];
+ (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) = f[vf::lbm::dir::DIR_P0M];
+ (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) = f[vf::lbm::dir::DIR_0MM];
+ (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) = f[vf::lbm::dir::DIR_0PM];
+ (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) = f[vf::lbm::dir::DIR_MMM];
+ (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) = f[vf::lbm::dir::DIR_PMM];
+ (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) = f[vf::lbm::dir::DIR_MPM];
+ (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) = f[vf::lbm::dir::DIR_PPM];
+
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) = f[vf::lbm::dir::DIR_P00];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) = f[vf::lbm::dir::DIR_0P0];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) = f[vf::lbm::dir::DIR_00P];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) = f[vf::lbm::dir::DIR_PP0];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) = f[vf::lbm::dir::DIR_MP0];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) = f[vf::lbm::dir::DIR_P0P];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) = f[vf::lbm::dir::DIR_M0P];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) = f[vf::lbm::dir::DIR_0PP];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) = f[vf::lbm::dir::DIR_0MP];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) = f[vf::lbm::dir::DIR_PPP];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) = f[vf::lbm::dir::DIR_MPP];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) = f[vf::lbm::dir::DIR_PMP];
+ (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = f[vf::lbm::dir::DIR_MMP];
+ }
}
}
}
diff --git a/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.h b/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.h
index 985ba75df5f0d726868e4ff854c384ad62eeb630..d675e72aabd3122a01e092e99fa22eb491270805 100644
--- a/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.h
+++ b/src/cpu/LiggghtsCoupling/LBM/IBcumulantK17LBMKernel.h
@@ -34,7 +34,7 @@
#ifndef IBcumulantK17LBMKernel_h__
#define IBcumulantK17LBMKernel_h__
-#include "LBMKernel.h"
+#include "LiggghtsCouplingLBMKernel.h"
#include "BCSet.h"
#include "D3Q27System.h"
#include "UbTiming.h"
@@ -50,7 +50,7 @@
//! <a href="http://dx.doi.org/10.1016/j.jcp.2017.05.040"><b>[ Geier et al., (2017), 10.1016/j.jcp.2017.05.040]</b></a>,
//! <a href="http://dx.doi.org/10.1016/j.jcp.2017.07.004"><b>[ Geier et al., (2017), 10.1016/j.jcp.2017.07.004]</b></a>
//!
-class IBcumulantK17LBMKernel : public LBMKernel
+class IBcumulantK17LBMKernel : public LiggghtsCouplingLBMKernel
{
public:
IBcumulantK17LBMKernel();
@@ -58,11 +58,6 @@ public:
void calculate(int step) override;
SPtr<LBMKernel> clone() override;
double getCalculationTime() override { return .0; }
- CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr getParticleData() { return particleData; };
- void setParticleData(CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData)
- {
- this->particleData = particleData;
- };
protected:
inline void forwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K);
@@ -72,9 +67,9 @@ protected:
virtual void initDataSet();
- CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr localDistributions;
- CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributions;
- CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr restDistributions;
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr localDistributionsF;
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributionsF;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr restDistributionsF;
mu::value_type muX1, muX2, muX3;
mu::value_type muDeltaT;
@@ -83,7 +78,6 @@ protected:
LBMReal forcingX2;
LBMReal forcingX3;
- CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData;
};
////////////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.cpp b/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.cpp
index e5934c275ba41db03c4c55be69a304e62bedb684..c502c4b48150e558aa785fa1bdc8c16bbecc1986 100644
--- a/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.cpp
+++ b/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.cpp
@@ -184,11 +184,11 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
localDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
nonLocalDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
- zeroDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
+ restDistributionsF = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
localDistributionsH1 = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getLocalDistributions();
nonLocalDistributionsH1 = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getNonLocalDistributions();
- zeroDistributionsH1 = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getZeroDistributions();
+ restDistributionsH1 = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getHdistributions())->getZeroDistributions();
CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr pressure = dataSet->getPressureField();
@@ -239,7 +239,7 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
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);
+ real mfbbb = (*this->restDistributionsH1)(x1, x2, x3);
// omegaDRho = 2.0;// 1.5;
// real phiOld = (*phaseField)(x1, x2, x3);
@@ -680,12 +680,12 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
real mfcaa = (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p);
real mfaca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p);
real mfcca = (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p);
- real mfbbb = (*this->zeroDistributionsF)(x1, x2, x3);
+ real mfbbb = (*this->restDistributionsF)(x1, x2, x3);
- LBMReal f[D3Q27System::ENDF + 1];
- LBMReal fEq[D3Q27System::ENDF + 1];
- LBMReal fEqSolid[D3Q27System::ENDF + 1];
- LBMReal fPre[D3Q27System::ENDF + 1];
+ real f[D3Q27System::ENDF + 1];
+ real fEq[D3Q27System::ENDF + 1];
+ real fEqSolid[D3Q27System::ENDF + 1];
+ real fPre[D3Q27System::ENDF + 1];
f[vf::lbm::dir::DIR_000] = mfbbb;
@@ -1169,7 +1169,7 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
(*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac; //* rho * c1o3;
(*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac; //* rho * c1o3;
- (*this->zeroDistributionsF)(x1, x2, x3) = mfbbb; // *rho* c1o3;
+ (*this->restDistributionsF)(x1, x2, x3) = mfbbb; // *rho* c1o3;
f[vf::lbm::dir::DIR_000] = mfbbb;
@@ -1202,7 +1202,7 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
f[vf::lbm::dir::DIR_PPM] = mfcca;
}
if ((*particleData)(x1, x2, x3)->solidFraction >= SOLFRAC_MIN) {
- LBMReal vx1, vx2, vx3, drho;
+ real vx1, vx2, vx3, drho;
D3Q27System::calcIncompMacroscopicValues(f, drho, vx1, vx2, vx3);
D3Q27System::calcIncompFeq(fEq, drho, vx1, vx2, vx3);
@@ -1254,7 +1254,7 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
}
} /* if solidFraction > SOLFRAC_MAX */
- (*this->zeroDistributionsF)(x1, x2, x3) = f[vf::lbm::dir::DIR_000];
+ (*this->restDistributionsF)(x1, x2, x3) = f[vf::lbm::dir::DIR_000];
(*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) = f[vf::lbm::dir::DIR_M00];
(*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) = f[vf::lbm::dir::DIR_0M0];
@@ -1318,7 +1318,7 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
mfcaa = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BSE, x1, x2p, x3p);
mfaca = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
mfcca = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
- mfbbb = (*this->zeroDistributionsH1)(x1, x2, x3);
+ mfbbb = (*this->restDistributionsH1)(x1, x2, x3);
////////////////////////////////////////////////////////////////////////////////////
//! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq. (J1)-(J3) \ref
@@ -1529,7 +1529,7 @@ void IBsharpInterfaceLBMKernel::calculate(int step)
(*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p) = mfcac;
(*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p) = mfaac;
- (*this->zeroDistributionsH1)(x1, x2, x3) = mfbbb;
+ (*this->restDistributionsH1)(x1, x2, x3) = mfbbb;
}
}
@@ -1715,7 +1715,7 @@ void IBsharpInterfaceLBMKernel::computePhasefield()
h[DIR_MPM] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNW, x1p, x2, x3p);
h[DIR_PPM] = (*this->nonLocalDistributionsH1)(D3Q27System::ET_BNE, x1, x2, x3p);
- h[DIR_000] = (*this->zeroDistributionsH1)(x1, x2, x3);
+ h[DIR_000] = (*this->restDistributionsH1)(x1, x2, x3);
}
}
}
diff --git a/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.h b/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.h
index c402529083b478be51c913145c9dd76b2566342b..a24131488a642a0a44712dd73fe4188cfeafaab2 100644
--- a/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.h
+++ b/src/cpu/LiggghtsCoupling/LBM/IBsharpInterfaceLBMKernel.h
@@ -36,7 +36,7 @@
#include "BCSet.h"
#include "D3Q27System.h"
-#include "LBMKernel.h"
+#include "LiggghtsCouplingLBMKernel.h"
#include "basics/container/CbArray3D.h"
#include "basics/container/CbArray4D.h"
#include "basics/utilities/UbTiming.h"
@@ -45,7 +45,7 @@
//! \brief Multiphase Cascaded Cumulant LBM kernel.
//! \details CFD solver that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model
//! \author M. Geier, K. Kutscher, Hesameddin Safari
-class IBsharpInterfaceLBMKernel : public LBMKernel
+class IBsharpInterfaceLBMKernel : public LiggghtsCouplingLBMKernel
{
public:
IBsharpInterfaceLBMKernel();
@@ -61,15 +61,6 @@ public:
return .0;
}
- CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr getParticleData()
- {
- return particleData;
- };
- void setParticleData(CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData)
- {
- this->particleData = particleData;
- };
-
protected:
virtual void initDataSet();
void swapDistributions() override;
@@ -85,11 +76,11 @@ protected:
CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr localDistributionsF;
CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributionsF;
- CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr zeroDistributionsF;
+ CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr restDistributionsF;
CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr localDistributionsH1;
CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributionsH1;
- CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr zeroDistributionsH1;
+ CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr restDistributionsH1;
CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr pressureOld;
CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr p1Old;
@@ -137,7 +128,6 @@ protected:
real forcingX2;
real forcingX3;
- CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData;
};
#endif
diff --git a/src/cpu/LiggghtsCoupling/LBM/LiggghtsCouplingLBMKernel.cpp b/src/cpu/LiggghtsCoupling/LBM/LiggghtsCouplingLBMKernel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ff7a3bf1b19d5efd626f638d47dc05727fc2bc89
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/LBM/LiggghtsCouplingLBMKernel.cpp
@@ -0,0 +1,92 @@
+#include "LiggghtsCouplingLBMKernel.h"
+#include "D3Q27System.h"
+
+//void LiggghtsCouplingLBMKernel::collisionOperator(int x1, int x2, int x3, real collFactorM, real fPre[])
+//{
+// //if ((*particleData)(x1, x2, x3)->solidFraction >= SOLFRAC_MIN) {
+// LBMReal f[D3Q27System::ENDF + 1];
+// LBMReal fEq[D3Q27System::ENDF + 1];
+// LBMReal fEqSolid[D3Q27System::ENDF + 1];
+// LBMReal vx1, vx2, vx3, drho;
+// D3Q27System::calcIncompMacroscopicValues(f, drho, vx1, vx2, vx3);
+// D3Q27System::calcIncompFeq(fEq, drho, vx1, vx2, vx3);
+//
+// std::array<double, 3> uPart;
+// uPart[0] = (*particleData)(x1, x2, x3)->uPart[0];
+// uPart[1] = (*particleData)(x1, x2, x3)->uPart[1];
+// uPart[2] = (*particleData)(x1, x2, x3)->uPart[2];
+//
+// D3Q27System::calcIncompFeq(fEqSolid, drho, uPart[0], uPart[1], uPart[2]);
+// real rhoPhaseField = (phi[DIR_000] > c1o2) ? c1o1 : c1o1 / densityRatio;
+// if ((*particleData)(x1, x2, x3)->solidFraction > SOLFRAC_MAX) {
+// double const bb0 = fEq[vf::lbm::dir::DIR_000] - fEqSolid[vf::lbm::dir::DIR_000];
+// f[vf::lbm::dir::DIR_000] = fPre[vf::lbm::dir::DIR_000] + bb0;
+// for (int iPop = D3Q27System::FSTARTDIR; iPop <= D3Q27System::FENDDIR; iPop++) {
+// const int iOpp = D3Q27System::INVDIR[iPop];
+// double const bb = ((fPre[iOpp] - fEq[iOpp]) - (fPre[iPop] - fEqSolid[iPop]));
+// double const bbOpp = ((fPre[iPop] - fEq[iPop]) - (fPre[iOpp] - fEqSolid[iOpp]));
+//
+// f[iPop] = fPre[iPop] + bb;
+// f[iOpp] = fPre[iOpp] + bbOpp;
+//
+// (*particleData)(x1, x2, x3)->hydrodynamicForce[0] -= D3Q27System::DX1[iPop] * (bb - bbOpp) * rhoPhaseField;
+// (*particleData)(x1, x2, x3)->hydrodynamicForce[1] -= D3Q27System::DX2[iPop] * (bb - bbOpp) * rhoPhaseField;
+// (*particleData)(x1, x2, x3)->hydrodynamicForce[2] -= D3Q27System::DX3[iPop] * (bb - bbOpp) * rhoPhaseField;
+// }
+// } else { /* particleData.solidFraction < SOLFRAC_MAX */
+// // #ifdef LBDEM_USE_WEIGHING
+// double const ooo = 1. / collFactorM - 0.5;
+// double const B = (*particleData)(x1, x2, x3)->solidFraction * ooo / ((1. - (*particleData)(x1, x2, x3)->solidFraction) + ooo);
+// // #else
+// // T const B = particleData.solidFraction;
+// // #endif
+// double const oneMinB = 1. - B;
+//
+// double const bb0 = fEq[vf::lbm::dir::DIR_000] - fEqSolid[vf::lbm::dir::DIR_000];
+// f[vf::lbm::dir::DIR_000] = fPre[vf::lbm::dir::DIR_000] + oneMinB * (f[vf::lbm::dir::DIR_000] - fPre[vf::lbm::dir::DIR_000]) + B * bb0;
+//
+// for (int iPop = D3Q27System::FSTARTDIR; iPop <= D3Q27System::FENDDIR; iPop++) {
+// int const iOpp = D3Q27System::INVDIR[iPop];
+// double const bb = B * ((fPre[iOpp] - fEq[iOpp]) - (fPre[iPop] - fEqSolid[iPop]));
+// double const bbOpp = B * ((fPre[iPop] - fEq[iPop]) - (fPre[iOpp] - fEqSolid[iOpp]));
+//
+// f[iPop] = fPre[iPop] + oneMinB * (f[iPop] - fPre[iPop]) + bb;
+// f[iOpp] = fPre[iOpp] + oneMinB * (f[iOpp] - fPre[iOpp]) + bbOpp;
+//
+// (*particleData)(x1, x2, x3)->hydrodynamicForce[0] -= D3Q27System::DX1[iPop] * (bb - bbOpp) * rhoPhaseField;
+// (*particleData)(x1, x2, x3)->hydrodynamicForce[1] -= D3Q27System::DX2[iPop] * (bb - bbOpp) * rhoPhaseField;
+// (*particleData)(x1, x2, x3)->hydrodynamicForce[2] -= D3Q27System::DX3[iPop] * (bb - bbOpp) * rhoPhaseField;
+// }
+// } /* if solidFraction > SOLFRAC_MAX */
+//
+// // (*this->restDistributionsF)(x1, x2, x3) = f[vf::lbm::dir::DIR_000];
+//
+// // (*this->localDistributionsF)(D3Q27System::ET_E, x1, x2, x3) = f[vf::lbm::dir::DIR_M00];
+// // (*this->localDistributionsF)(D3Q27System::ET_N, x1, x2, x3) = f[vf::lbm::dir::DIR_0M0];
+// // (*this->localDistributionsF)(D3Q27System::ET_T, x1, x2, x3) = f[vf::lbm::dir::DIR_00M];
+// // (*this->localDistributionsF)(D3Q27System::ET_NE, x1, x2, x3) = f[vf::lbm::dir::DIR_MM0];
+// // (*this->localDistributionsF)(D3Q27System::ET_NW, x1p, x2, x3) = f[vf::lbm::dir::DIR_PM0];
+// // (*this->localDistributionsF)(D3Q27System::ET_TE, x1, x2, x3) = f[vf::lbm::dir::DIR_M0M];
+// // (*this->localDistributionsF)(D3Q27System::ET_TW, x1p, x2, x3) = f[vf::lbm::dir::DIR_P0M];
+// // (*this->localDistributionsF)(D3Q27System::ET_TN, x1, x2, x3) = f[vf::lbm::dir::DIR_0MM];
+// // (*this->localDistributionsF)(D3Q27System::ET_TS, x1, x2p, x3) = f[vf::lbm::dir::DIR_0PM];
+// // (*this->localDistributionsF)(D3Q27System::ET_TNE, x1, x2, x3) = f[vf::lbm::dir::DIR_MMM];
+// // (*this->localDistributionsF)(D3Q27System::ET_TNW, x1p, x2, x3) = f[vf::lbm::dir::DIR_PMM];
+// // (*this->localDistributionsF)(D3Q27System::ET_TSE, x1, x2p, x3) = f[vf::lbm::dir::DIR_MPM];
+// // (*this->localDistributionsF)(D3Q27System::ET_TSW, x1p, x2p, x3) = f[vf::lbm::dir::DIR_PPM];
+//
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_W, x1p, x2, x3) = f[vf::lbm::dir::DIR_P00];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_S, x1, x2p, x3) = f[vf::lbm::dir::DIR_0P0];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_B, x1, x2, x3p) = f[vf::lbm::dir::DIR_00P];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_SW, x1p, x2p, x3) = f[vf::lbm::dir::DIR_PP0];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_SE, x1, x2p, x3) = f[vf::lbm::dir::DIR_MP0];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BW, x1p, x2, x3p) = f[vf::lbm::dir::DIR_P0P];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BE, x1, x2, x3p) = f[vf::lbm::dir::DIR_M0P];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BS, x1, x2p, x3p) = f[vf::lbm::dir::DIR_0PP];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BN, x1, x2, x3p) = f[vf::lbm::dir::DIR_0MP];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BSW, x1p, x2p, x3p) = f[vf::lbm::dir::DIR_PPP];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BSE, x1, x2p, x3p) = f[vf::lbm::dir::DIR_MPP];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BNW, x1p, x2, x3p) = f[vf::lbm::dir::DIR_PMP];
+// // (*this->nonLocalDistributionsF)(D3Q27System::ET_BNE, x1, x2, x3p) = f[vf::lbm::dir::DIR_MMP];
+// //}
+//}
\ No newline at end of file
diff --git a/src/cpu/LiggghtsCoupling/LBM/LiggghtsCouplingLBMKernel.h b/src/cpu/LiggghtsCoupling/LBM/LiggghtsCouplingLBMKernel.h
new file mode 100644
index 0000000000000000000000000000000000000000..86ac895486d357e9facc30c210b325af527dfbf2
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/LBM/LiggghtsCouplingLBMKernel.h
@@ -0,0 +1,66 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// This file is part of VirtualFluids. VirtualFluids is free software: you can
+// redistribute it and/or modify it under the terms of the GNU General Public
+// License as published by the Free Software Foundation, either version 3 of
+// the License, or (at your option) any later version.
+//
+// VirtualFluids is distributed in the hope that it will be useful, but WITHOUT
+// ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
+// FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
+// for more details.
+//
+// You should have received a copy of the GNU General Public License along
+// with VirtualFluids (see COPYING.txt). If not, see <http://www.gnu.org/licenses/>.
+//
+//! \file LiggghtsCouplingLBMKernel.h
+//! \ingroup LBM
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#ifndef LiggghtsCouplingLBMKernel_h
+#define LiggghtsCouplingLBMKernel_h
+
+#include "LBMKernel.h"
+#include "IBdynamicsParticleData.h"
+#include "basics/container/CbArray3D.h"
+#include "basics/container/CbArray4D.h"
+
+class LiggghtsCouplingLBMKernel : public LBMKernel
+{
+public:
+ virtual ~LiggghtsCouplingLBMKernel() = default;
+
+ CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr getParticleData()
+ {
+ return particleData;
+ };
+ void setParticleData(CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData)
+ {
+ this->particleData = particleData;
+ };
+
+
+ protected:
+ //void collisionOperator();
+ CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData;
+
+ //CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr localDistributionsF;
+ //CbArray4D<real, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributionsF;
+ //CbArray3D<real, IndexerX3X2X1>::CbArray3DPtr restDistributionsF;
+};
+
+#endif
\ No newline at end of file
diff --git a/src/cpu/LiggghtsCoupling/SimulationObserver/LiggghtsCouplingSimulationObserver.cpp b/src/cpu/LiggghtsCoupling/SimulationObserver/LiggghtsCouplingSimulationObserver.cpp
index 93be6707ac56cfd93b57b288b0100ee92f214beb..651421bd11312fcfa99cd002c5ce2a4a4a44a63f 100644
--- a/src/cpu/LiggghtsCoupling/SimulationObserver/LiggghtsCouplingSimulationObserver.cpp
+++ b/src/cpu/LiggghtsCoupling/SimulationObserver/LiggghtsCouplingSimulationObserver.cpp
@@ -105,32 +105,39 @@ void LiggghtsCouplingSimulationObserver::setSingleSphere3D(double *x, double *v,
SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData =
- //dynamicPointerCast<IBcumulantK17LBMKernel>(kernel)->getParticleData();
- dynamicPointerCast<IBsharpInterfaceLBMKernel>(kernel)->getParticleData();
+ dynamicPointerCast<LiggghtsCouplingLBMKernel>(kernel)->getParticleData();
if (!particleData)
continue;
- //int minX1 = 1;
- //int minX2 = 1;
- //int minX3 = 1;
+ int minX1b = 1;
+ int minX2b = 1;
+ int minX3b = 1;
- //int maxX1 = (int)(distributions->getNX1()) - 1;
- //int maxX2 = (int)(distributions->getNX2()) - 1;
- //int maxX3 = (int)(distributions->getNX3()) - 1;
+ int maxX1b = (int)(distributions->getNX1()) - 2;
+ int maxX2b = (int)(distributions->getNX2()) - 2;
+ int maxX3b = (int)(distributions->getNX3()) - 2;
real deltax = grid->getDeltaX(block);
UbTupleInt3 nodesMin = grid->getNodeIndexes(block, x[0] - r - deltax, x[1] - r - deltax, x[2] - r - deltax);
UbTupleInt3 nodesMax = grid->getNodeIndexes(block, x[0] + r + deltax, x[1] + r + deltax, x[2] + r + deltax);
- int minX1 = val<1>(nodesMin);
- int minX2 = val<2>(nodesMin);
- int minX3 = val<3>(nodesMin);
+ int minX1 = (val<1>(nodesMin) < minX1b) ? minX1b : val<1>(nodesMin);
+ int minX2 = (val<2>(nodesMin) < minX2b) ? minX2b : val<2>(nodesMin);
+ int minX3 = (val<3>(nodesMin) < minX3b) ? minX3b : val<3>(nodesMin);
- int maxX1 = val<1>(nodesMax);
- int maxX2 = val<2>(nodesMax);
- int maxX3 = val<3>(nodesMin);
+ int maxX1 = (val<1>(nodesMax) > maxX1b) ? maxX1b : val<1>(nodesMax);
+ int maxX2 = (val<2>(nodesMax) > maxX2b) ? maxX2b : val<2>(nodesMax);
+ int maxX3 = (val<3>(nodesMax) > maxX3b) ? maxX3b : val<3>(nodesMax);
+
+ //int minX1 = minX1b;
+ //int minX2 = minX2b;
+ //int minX3 = minX3b;
+
+ //int maxX1 = maxX1b;
+ //int maxX2 = maxX2b;
+ //int maxX3 = maxX3b;
for (int ix3 = minX3; ix3 <= maxX3; ix3++) {
@@ -350,8 +357,7 @@ void LiggghtsCouplingSimulationObserver::SumForceTorque3D(ParticleData::Particle
SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData =
- //dynamicPointerCast<IBcumulantK17LBMKernel>(kernel)->getParticleData();
- dynamicPointerCast<IBsharpInterfaceLBMKernel>(kernel)->getParticleData();
+ dynamicPointerCast<LiggghtsCouplingLBMKernel>(kernel)->getParticleData();
if (!particleData)
continue;
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NoSlipBCStrategy.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NoSlipBCStrategy.cpp
index cabd41b30dd31ac8751c60f00f46c899cbfa2334..bb98e499421328abc36f8911a14bcb8c90ebabd1 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NoSlipBCStrategy.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NoSlipBCStrategy.cpp
@@ -54,23 +54,22 @@ void NoSlipBCStrategy::addDistributions(SPtr<DistributionArray3D> distributions)
//////////////////////////////////////////////////////////////////////////
void NoSlipBCStrategy::applyBC()
{
- real f[D3Q27System::ENDF + 1];
- real feq[D3Q27System::ENDF + 1];
+ using namespace vf::basics::constant;
+ using namespace D3Q27System;
+ real f[ENDF + 1];
+ real feq[ENDF + 1];
distributions->getDistributionInv(f, x1, x2, x3);
real rho, vx1, vx2, vx3;
calcMacrosFct(f, rho, vx1, vx2, vx3);
calcFeqFct(feq, rho, vx1, vx2, vx3);
- for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
+ for (int fdir = FSTARTDIR; fdir <= FENDDIR; fdir++) {
if (bcPtr->hasNoSlipBoundaryFlag(fdir)) {
// quadratic bounce back
- const int invDir = D3Q27System::INVDIR[fdir];
- real q = bcPtr->getQ(invDir);
- real fReturn = ((vf::basics::constant::c1o1 - q) / (vf::basics::constant::c1o1 + q)) * ((f[invDir] - feq[invDir]) / (vf::basics::constant::c1o1 - collFactor) + feq[invDir]) +
- ((q / (vf::basics::constant::c1o1 + q)) * (f[invDir] + f[fdir]));
- distributions->setDistributionForDirection(fReturn, x1 + D3Q27System::DX1[invDir],
- x2 + D3Q27System::DX2[invDir], x3 + D3Q27System::DX3[invDir],
- fdir);
+ const int invDir = INVDIR[fdir];
+ real q = bcPtr->getQ(invDir);
+ real fReturn = ((c1o1 - q) / (c1o1 + q)) * ((f[invDir] - feq[invDir]) / (c1o1 - collFactor) + feq[invDir]) + ((q / (c1o1 + q)) * (f[invDir] + f[fdir]));
+ distributions->setDistributionForDirection(fReturn, x1 + DX1[invDir], x2 + DX2[invDir], x3 + DX3[invDir], fdir);
}
}
}