diff --git a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
index fd892b6b9288b496fc11007ea034a145484d1c5e..978ae376f1c842b0d108730fb9d2ffe43fa65537 100644
--- a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
+++ b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
@@ -30,22 +30,40 @@ int main(int argc, char *argv[])
double g_maxX1 = 1;
double g_maxX2 = 1;
- double g_maxX3 = 1;
+ double g_maxX3 = 2;
+
+ int blockNX[3] = { 16, 16, 16 };
+
+ double dx = 1./32.;
+
+ double Re = 300;
+ double nuLB = 5e-5;
+
+ SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
+ SPtr<BCProcessor> bcProc = make_shared<BCProcessor>();
+ kernel->setBCProcessor(bcProc);
+
+ SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
+ //////////////////////////////////////////////////////////////////////////////////
+ // BC visitor
+ BoundaryConditionsBlockVisitor bcVisitor;
+ bcVisitor.addBC(noSlipBCAdapter);
- int blockNX[3] = { 10, 10, 10 };
- double dx = 0.1;
- double nuLB = 0.005;
SPtr<Grid3D> grid = make_shared<Grid3D>(comm);
grid->setPeriodicX1(true);
grid->setPeriodicX2(true);
- grid->setPeriodicX3(true);
+ grid->setPeriodicX3(false);
grid->setDeltaX(dx);
grid->setBlockNX(blockNX[0], blockNX[1], blockNX[2]);
- string outputPath = "d:/temp/lll2";
+ string outputPath = "d:/temp/lll4";
+
+ SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased,
+ D3Q27System::BSW, MetisPartitioner::RECURSIVE));
SPtr<GbObject3D> gridCube = make_shared <GbCuboid3D>(g_minX1, g_minX2, g_minX3, g_maxX1, g_maxX2, g_maxX3);
if (myid == 0)
@@ -60,13 +78,27 @@ int main(int argc, char *argv[])
ppblocks->process(0);
ppblocks.reset();
- SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
- SPtr<BCProcessor> bcProc = make_shared<BCProcessor>();
- kernel->setBCProcessor(bcProc);
+ double dx2 = 2.0 * dx;
+ GbCuboid3DPtr wallZmin(
+ new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_minX3 - dx2, g_maxX1 + dx2, g_maxX2 + dx2, g_minX3));
+ GbSystem3D::writeGeoObject(wallZmin.get(), outputPath + "/geo/wallZmin", WbWriterVtkXmlASCII::getInstance());
+ GbCuboid3DPtr wallZmax(
+ new GbCuboid3D(g_minX1 - dx2, g_minX2 - dx2, g_maxX3, g_maxX1 + dx2, g_maxX2 + dx2, g_maxX3 + dx2));
+ GbSystem3D::writeGeoObject(wallZmax.get(), outputPath + "/geo/wallZmax", WbWriterVtkXmlASCII::getInstance());
+
+ SPtr<D3Q27Interactor> wallZminInt(new D3Q27Interactor(wallZmin, grid, noSlipBCAdapter, Interactor3D::SOLID));
+ SPtr<D3Q27Interactor> wallZmaxInt(new D3Q27Interactor(wallZmax, grid, noSlipBCAdapter, Interactor3D::SOLID));
+
+ InteractorsHelper intHelper(grid, metisVisitor, true);
+ intHelper.addInteractor(wallZminInt);
+ intHelper.addInteractor(wallZmaxInt);
+ intHelper.selectBlocks();
SetKernelBlockVisitor kernelVisitor(kernel, nuLB, 1e9, 1e9);
grid->accept(kernelVisitor);
+ intHelper.setBC();
+
InitDistributionsBlockVisitor initVisitor;
grid->accept(initVisitor);
@@ -78,18 +110,24 @@ int main(int argc, char *argv[])
LiggghtsCouplingWrapper wrapper(argv, mpi_comm);
- double d_part = 0.3;
+ double d_part = 0.25;
double r_p = d_part / 2.0;
// 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::WATER, 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>(r_p, 0.1, 1000, r_p / dx, 0.01);
+ std::cout << units->toString() << std::endl;
+
+ //return 0;
double v_frac = 0.1;
double dt_phys = units->getFactorTimeLbToW();
int demSubsteps = 10;
double dt_dem = dt_phys / (double)demSubsteps;
int vtkSteps = 1;
- string demOutDir = "d:/temp/lll2/";
+ string demOutDir = outputPath; // "d:/temp/lll2/";
+
+ wrapper.execCommand("echo log");
wrapper.setVariable("r_part", d_part / 2);
wrapper.setVariable("v_frac", v_frac);
@@ -102,16 +140,31 @@ int main(int argc, char *argv[])
wrapper.setVariable("dmp_dir", demOutDir);
wrapper.execFile((char *)inFile2.c_str());
- //wrapper.runUpto(demSubsteps - 1);
+ wrapper.runUpto(demSubsteps - 1);
-
+
SPtr<LiggghtsCouplingCoProcessor> lcCoProcessor =
make_shared<LiggghtsCouplingCoProcessor>(grid, lScheduler, comm, wrapper, demSubsteps, units);
+ // boundary conditions grid
+ {
+ SPtr<UbScheduler> geoSch(new UbScheduler(1));
+ SPtr<WriteBoundaryConditionsCoProcessor> ppgeo(new WriteBoundaryConditionsCoProcessor(
+ grid, geoSch, outputPath, WbWriterVtkXmlBinary::getInstance(), comm));
+ ppgeo->process(0);
+ ppgeo.reset();
+ }
+
+ grid->accept(bcVisitor);
+
+ OneDistributionSetConnectorsBlockVisitor setConnsVisitor(comm);
+ grid->accept(setConnsVisitor);
+
+
// write data for visualization of macroscopic quantities
SPtr<UbScheduler> visSch(new UbScheduler(vtkSteps));
SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(
- new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlASCII::getInstance(),
+ new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(),
SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
int endTime = 1000; //20;
diff --git a/apps/cpu/LiggghtsApp/in.lbdem b/apps/cpu/LiggghtsApp/in.lbdem
index 366f7afd134f76492d3d9756bf6743c5bde0c212..ae2baa373dc75edec41634522e53ade31f459b4d 100644
--- a/apps/cpu/LiggghtsApp/in.lbdem
+++ b/apps/cpu/LiggghtsApp/in.lbdem
@@ -1,15 +1,18 @@
+echo none
units si
atom_style granular
atom_modify map array
+
+
communicate single vel yes
boundary f f f
newton off
processors * * 1
-region box block 0. 1. 0. 1. 0. 1. units box
+region box block 0. 1. 0. 1. 0. 2. units box
create_box 1 box
variable skin equal 0.01
@@ -41,10 +44,10 @@ fix ywalls2 all wall/gran model hertz tangential history primitive type 1 yplane
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.
-create_atoms 1 single 0.5 0.5 0.75
+create_atoms 1 single 0.5 0.5 1.75
#create_atoms 1 single 0.38 0.05 0.05
-set group all diameter 0.3 density 2400
+set group all diameter 0.25 density 2400
atom_modify sort 0 0.0
@@ -67,6 +70,6 @@ atom_modify sort 0 0.0
# #fix ins all insert/pack seed 1001 distributiontemplate pdd1 insert_every once &
# # overlapcheck yes particles_in_region 1 region insreg ntry_mc 10000
-
+echo none
run 1
diff --git a/apps/cpu/LiggghtsApp/in2.lbdem b/apps/cpu/LiggghtsApp/in2.lbdem
index 7098890cc28ea5f22ef442a388b1f79b4708107c..8d6a0748befeea97d6faa7efa87ffd60e5c56cfd 100644
--- a/apps/cpu/LiggghtsApp/in2.lbdem
+++ b/apps/cpu/LiggghtsApp/in2.lbdem
@@ -1,4 +1,6 @@
+echo none
+
timestep ${t_step}
# thermo settings
@@ -18,4 +20,6 @@ variable dmp_fname string ${dmp_dir}d_*.liggghts
# dump dmp all custom ${dmp_stp} ${dmp_dir}d_*.liggghts &
# id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
-dump dmp all custom/vtk 10 ${dmp_dir}post/atom_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
+dump dmp all custom/vtk ${dmp_stp} ${dmp_dir}/post/atom_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
+
+echo none
\ No newline at end of file
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
index d4926e0b11bd96b4d097ccf79060b437de6ce8da..da8c93296a333050c036f7fe084bbe23aa39ca23 100644
--- a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
@@ -35,6 +35,8 @@ void LiggghtsCouplingCoProcessor::process(double actualTimeStep)
{
std::cout << "step: " << actualTimeStep << "\n";
+ getForcesFromLattice();
+
wrapper.run(demSteps);
setSpheresOnLattice();
@@ -74,10 +76,10 @@ void LiggghtsCouplingCoProcessor::setSpheresOnLattice()
r = wrapper.lmp->atom->radius[iS]; // * units->getFactorLentghWToLb();
- std::cout << "x[0] = " << x[0] << ", x[1] = " << x[1] << ", x[2] = " << x[2] << std::endl;
- std::cout << "v[0] = " << v[0] << ", v[1] = " << v[1] << ", v[2] = " << v[2] << std::endl;
- std::cout << "omega[0] = " << omega[0] << ", omega[1] = " << omega[1] << ", omega[2] = " << omega[2] << std::endl;
- std::cout << "r = " << r << std::endl;
+ //std::cout << "x[0] = " << x[0] << ", x[1] = " << x[1] << ", x[2] = " << x[2] << std::endl;
+ //std::cout << "v[0] = " << v[0] << ", v[1] = " << v[1] << ", v[2] = " << v[2] << std::endl;
+ //std::cout << "omega[0] = " << omega[0] << ", omega[1] = " << omega[1] << ", omega[2] = " << omega[2] << std::endl;
+ //std::cout << "r = " << r << std::endl;
setSingleSphere3D(x, v, omega, r, id);
}
@@ -111,9 +113,9 @@ void LiggghtsCouplingCoProcessor::setSingleSphere3D(double *x, double *v, double
int minX2 = 1;
int minX3 = 1;
- int maxX1 = (int)(distributions->getNX1()) - 2;
- int maxX2 = (int)(distributions->getNX2()) - 2;
- int maxX3 = (int)(distributions->getNX3()) - 2;
+ int maxX1 = (int)(distributions->getNX1()) - 1;
+ int maxX2 = (int)(distributions->getNX2()) - 1;
+ int maxX3 = (int)(distributions->getNX3()) - 1;
for (int ix3 = minX3; ix3 < maxX3; ix3++) {
for (int ix2 = minX2; ix2 < maxX2; ix2++) {
@@ -341,9 +343,9 @@ void LiggghtsCouplingCoProcessor::SumForceTorque3D(ParticleData::ParticleDataArr
int minX2 = 1;
int minX3 = 1;
- int maxX1 = (int)(distributions->getNX1()) - 2;
- int maxX2 = (int)(distributions->getNX2()) - 2;
- int maxX3 = (int)(distributions->getNX3()) - 2;
+ int maxX1 = (int)(distributions->getNX1()) - 1;
+ int maxX2 = (int)(distributions->getNX2()) - 1;
+ int maxX3 = (int)(distributions->getNX3()) - 1;
for (int ix3 = minX3; ix3 < maxX3; ix3++) {
for (int ix2 = minX2; ix2 < maxX2; ix2++) {