From 3a2b0d33041860de577f542370f8e1e775e96170 Mon Sep 17 00:00:00 2001
From: Kutscher <kutscher@irmb.tu-bs.de>
Date: Fri, 14 Oct 2022 13:57:13 +0200
Subject: [PATCH] adds apps: Falling Sphere, Nozzle
---
apps/cpu/Applications.cmake | 2 +
apps/cpu/FallingSphere/CMakeLists.txt | 3 +
apps/cpu/FallingSphere/FallingSphere.cpp | 175 ++++++++++++
apps/cpu/FallingSphere/in.lbdem | 75 +++++
apps/cpu/FallingSphere/in2.lbdem | 25 ++
apps/cpu/LiggghtsApp/LiggghtsApp.cpp | 35 ++-
apps/cpu/LiggghtsApp/in.lbdem | 8 +-
apps/cpu/LiggghtsApp/in2.lbdem | 8 +-
apps/cpu/Nozzle/CMakeLists.txt | 3 +
apps/cpu/Nozzle/in.nozzle | 125 ++++++++
apps/cpu/Nozzle/nozzle.cpp | 268 ++++++++++++++++++
apps/cpu/ViskomatXL/viskomat.cfg | 3 +-
apps/cpu/ViskomatXL/viskomat.cpp | 2 +-
.../LiggghtsCouplingCoProcessor.cpp | 8 -
.../VirtualFluidsCore/Utilities/MemoryUtil.h | 45 +++
.../Visitors/SetKernelBlockVisitor.cpp | 4 +-
.../Visitors/SetKernelBlockVisitor.h | 2 +-
17 files changed, 757 insertions(+), 34 deletions(-)
create mode 100644 apps/cpu/FallingSphere/CMakeLists.txt
create mode 100644 apps/cpu/FallingSphere/FallingSphere.cpp
create mode 100644 apps/cpu/FallingSphere/in.lbdem
create mode 100644 apps/cpu/FallingSphere/in2.lbdem
create mode 100644 apps/cpu/Nozzle/CMakeLists.txt
create mode 100644 apps/cpu/Nozzle/in.nozzle
create mode 100644 apps/cpu/Nozzle/nozzle.cpp
diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 37fc4de29..143f20a9e 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -10,6 +10,8 @@ add_subdirectory(${APPS_ROOT_CPU}/FlowAroundCylinder)
add_subdirectory(${APPS_ROOT_CPU}/LaminarTubeFlow)
add_subdirectory(${APPS_ROOT_CPU}/MultiphaseDropletTest)
add_subdirectory(${APPS_ROOT_CPU}/LiggghtsApp)
+add_subdirectory(${APPS_ROOT_CPU}/FallingSphere)
+add_subdirectory(${APPS_ROOT_CPU}/Nozzle)
#add_subdirectory(tests)
#add_subdirectory(Applications/gridRf)
diff --git a/apps/cpu/FallingSphere/CMakeLists.txt b/apps/cpu/FallingSphere/CMakeLists.txt
new file mode 100644
index 000000000..94eab3ae0
--- /dev/null
+++ b/apps/cpu/FallingSphere/CMakeLists.txt
@@ -0,0 +1,3 @@
+PROJECT(FallingSphere)
+
+vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} LiggghtsCoupling FILES FallingSphere.cpp )
diff --git a/apps/cpu/FallingSphere/FallingSphere.cpp b/apps/cpu/FallingSphere/FallingSphere.cpp
new file mode 100644
index 000000000..88b5373c9
--- /dev/null
+++ b/apps/cpu/FallingSphere/FallingSphere.cpp
@@ -0,0 +1,175 @@
+#include <iostream>
+#include <string>
+#include <memory>
+
+#include "VirtualFluids.h"
+
+#include "LiggghtsCouplingCoProcessor.h"
+#include "LiggghtsCouplingWrapper.h"
+#include "IBcumulantK17LBMKernel.h"
+
+using namespace std;
+
+
+int main(int argc, char *argv[])
+{
+ SPtr<Communicator> comm = MPICommunicator::getInstance();
+ int myid = comm->getProcessID();
+
+
+ // bounding box
+ double g_minX1 = 0;
+ double g_minX2 = 0;
+ double g_minX3 = 0;
+
+ double g_maxX1 = 1;
+ double g_maxX2 = 1;
+ double g_maxX3 = 2;
+
+ //int blockNX[3] = { 16, 16, 16 };
+ //double dx = 1./32.;
+
+ double d_part = 0.1;
+ double r_p = d_part / 2.0;
+
+ int blockNX[3] = { 10, 10, 10 };
+ double dx = 0.05;
+
+
+ double nuLB = 1e-2;
+
+ 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);
+
+ SPtr<Grid3D> grid = make_shared<Grid3D>(comm);
+ grid->setPeriodicX1(true);
+ grid->setPeriodicX2(true);
+ grid->setPeriodicX3(false);
+ grid->setDeltaX(dx);
+ grid->setBlockNX(blockNX[0], blockNX[1], blockNX[2]);
+
+ string outputPath = "d:/temp/FallingSpheres";
+
+ UbSystem::makeDirectory(outputPath);
+ UbSystem::makeDirectory(outputPath + "/liggghts");
+
+ 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)
+ GbSystem3D::writeGeoObject(gridCube.get(), outputPath + "/geo/gridCube", WbWriterVtkXmlBinary::getInstance());
+
+ GenBlocksGridVisitor genBlocks(gridCube);
+ grid->accept(genBlocks);
+
+ SPtr<CoProcessor> ppblocks =
+ make_shared <WriteBlocksCoProcessor>(grid, SPtr<UbScheduler>(new UbScheduler(1)), outputPath,
+ WbWriterVtkXmlBinary::getInstance(), comm);
+ ppblocks->process(0);
+ ppblocks.reset();
+
+ 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);
+
+ SPtr<UbScheduler> lScheduler = make_shared<UbScheduler>(1);
+ string inFile1 = "d:/Projects/VirtualFluids_LIGGGHTS_coupling/apps/cpu/FallingSphere/in.lbdem";
+ string inFile2 = "d:/Projects/VirtualFluids_LIGGGHTS_coupling/apps/cpu/FallingSphere/in2.lbdem";
+ MPI_Comm mpi_comm = *(MPI_Comm*)(comm->getNativeCommunicator());
+ LiggghtsCouplingWrapper wrapper(argv, mpi_comm);
+
+
+
+ // 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>(r_p, 0.1, 1000, r_p / dx, 0.01);
+ std::cout << units->toString() << std::endl;
+
+ double v_frac = 0.1;
+ double dt_phys = units->getFactorTimeLbToW();
+ int demSubsteps = 10;
+ double dt_dem = dt_phys / (double)demSubsteps;
+ int vtkSteps = 100;
+ string demOutDir = outputPath;
+
+ wrapper.execCommand("echo none");
+
+ //wrapper.setVariable("d_part", d_part);
+ wrapper.setVariable("r_part", d_part/2.);
+ wrapper.setVariable("v_frac", v_frac);
+
+ //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);
+
+ 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, WbWriterVtkXmlBinary::getInstance(),
+ SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+
+ int endTime = 3000; //20;
+ SPtr<Calculator> calculator(new BasicCalculator(grid, lScheduler, endTime));
+ calculator->addCoProcessor(lcCoProcessor);
+ calculator->addCoProcessor(writeMQCoProcessor);
+
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
+ calculator->calculate();
+ if (myid == 0) UBLOG(logINFO, "Simulation-end");
+
+
+ return 0;
+}
diff --git a/apps/cpu/FallingSphere/in.lbdem b/apps/cpu/FallingSphere/in.lbdem
new file mode 100644
index 000000000..b47a85c99
--- /dev/null
+++ b/apps/cpu/FallingSphere/in.lbdem
@@ -0,0 +1,75 @@
+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. 10. units box
+create_box 1 box
+
+variable skin equal 0.01
+neighbor ${skin} bin
+neigh_modify delay 0 binsize 0.01 one 1000
+
+fix grav all gravity 0.981 vector 0 0 -1
+
+
+fix m1 all property/global youngsModulus peratomtype 1e8
+fix m2 all property/global poissonsRatio peratomtype 0.4
+fix m3 all property/global coefficientRestitution peratomtypepair 1 0.95
+fix m4 all property/global coefficientFriction peratomtypepair 1 0.45
+fix m5 all property/global coefficientRollingFriction peratomtypepair 1 0.020
+
+# lb coupling fix
+fix lbcoupling all couple/lb/onetoone
+
+
+pair_style gran model hertz tangential history rolling_friction cdt
+pair_coeff * *
+
+fix 1 all nve/sphere
+
+fix xwalls1 all wall/gran model hertz tangential history primitive type 1 xplane 0.
+fix xwalls2 all wall/gran model hertz tangential history primitive type 1 xplane 1.
+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.
+
+create_atoms 1 single 0.5 0.5 9.75
+#create_atoms 1 single 0.38 0.05 0.05
+
+set group all diameter ${d_part} density 2400
+
+atom_modify sort 0 0.0
+
+#fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 1000 radius constant 0.015
+#fix pts2 all particletemplate/sphere 1 atom_type 1 density constant 1000 radius constant 0.01
+#fix pts1 all particletemplate/sphere 1 atom_type 1 density constant 1100 radius constant ${r_part}
+
+# fix pdd1 all particledistribution/discrete 6778 1 pts1 1.0
+# #fix pdd2 all particledistribution/discrete 6778 2 pts2 0.2 pts3 0.8
+
+# # region insreg block 0.1 0.9 0.1 0.9 1.3 1.9 units box
+
+
+# #fix ins all insert/pack seed 1001 distributiontemplate pdd1 insert_every once &
+# # overlapcheck yes particles_in_region 350 region insreg ntry_mc 10000
+# fix ins all insert/pack seed 1001 distributiontemplate pdd1 insert_every once &
+ # overlapcheck yes volumefraction_region ${v_frac} region insreg ntry_mc 10000
+# #fix ins all insert/pack seed 1001 distributiontemplate pdd2 insert_every once &
+# # overlapcheck yes volumefraction_region 0.05 region insreg ntry_mc 10000
+# #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/FallingSphere/in2.lbdem b/apps/cpu/FallingSphere/in2.lbdem
new file mode 100644
index 000000000..f11767f12
--- /dev/null
+++ b/apps/cpu/FallingSphere/in2.lbdem
@@ -0,0 +1,25 @@
+
+echo none
+
+timestep ${t_step}
+
+# thermo settings
+fix ts all check/timestep/gran 10000 0.1 0.1
+compute 1 all erotate/sphere
+thermo_style custom step atoms ke c_1 f_ts[1] f_ts[2] cpu
+thermo 10000
+thermo_modify lost ignore norm no flush yes
+compute_modify thermo_temp dynamic yes
+
+# particle dump
+variable dmp_fname string ${dmp_dir}d_*.liggghts
+
+# dump dmp all custom ${dmp_stp} ${dmp_fname} &
+# id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
+
+# 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 ${dmp_stp} ${dmp_dir}/liggghts/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/apps/cpu/LiggghtsApp/LiggghtsApp.cpp b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
index 978ae376f..6eb27d6dc 100644
--- a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
+++ b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
@@ -30,14 +30,25 @@ int main(int argc, char *argv[])
double g_maxX1 = 1;
double g_maxX2 = 1;
- double g_maxX3 = 2;
+ double g_maxX3 = 10;
int blockNX[3] = { 16, 16, 16 };
double dx = 1./32.;
- double Re = 300;
- double nuLB = 5e-5;
+
+ 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::AIR_20C, r_p / dx);
+ // SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(r_p, 0.1, 1000, r_p / dx, 0.01);
+ SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(r_p, 0.1, 1000, r_p / dx, 0.01);
+ //SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(r_p, LBMUnitConverter::OIL, r_p / dx);
+ std::cout << units->toString() << std::endl;
+
+ //double Re = 300;
+ double nuLB = 1e-2; // 5e-5;
SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
SPtr<BCProcessor> bcProc = make_shared<BCProcessor>();
@@ -60,7 +71,7 @@ int main(int argc, char *argv[])
grid->setDeltaX(dx);
grid->setBlockNX(blockNX[0], blockNX[1], blockNX[2]);
- string outputPath = "d:/temp/lll4";
+ string outputPath = "d:/temp/lll8";
SPtr<Grid3DVisitor> metisVisitor(new MetisPartitioningGridVisitor(comm, MetisPartitioningGridVisitor::LevelBased,
D3Q27System::BSW, MetisPartitioner::RECURSIVE));
@@ -110,13 +121,9 @@ int main(int argc, char *argv[])
LiggghtsCouplingWrapper wrapper(argv, mpi_comm);
- 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::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;
@@ -124,10 +131,10 @@ int main(int argc, char *argv[])
double dt_phys = units->getFactorTimeLbToW();
int demSubsteps = 10;
double dt_dem = dt_phys / (double)demSubsteps;
- int vtkSteps = 1;
+ int vtkSteps = 100;
string demOutDir = outputPath; // "d:/temp/lll2/";
- wrapper.execCommand("echo log");
+ //wrapper.execCommand("echo none");
wrapper.setVariable("r_part", d_part / 2);
wrapper.setVariable("v_frac", v_frac);
@@ -167,7 +174,7 @@ int main(int argc, char *argv[])
new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(),
SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
- int endTime = 1000; //20;
+ int endTime = 3000; //20;
SPtr<Calculator> calculator(new BasicCalculator(grid, lScheduler, endTime));
calculator->addCoProcessor(lcCoProcessor);
calculator->addCoProcessor(writeMQCoProcessor);
diff --git a/apps/cpu/LiggghtsApp/in.lbdem b/apps/cpu/LiggghtsApp/in.lbdem
index ae2baa373..c356a2f6b 100644
--- a/apps/cpu/LiggghtsApp/in.lbdem
+++ b/apps/cpu/LiggghtsApp/in.lbdem
@@ -1,4 +1,4 @@
-echo none
+#verbose no
units si
atom_style granular
@@ -12,7 +12,7 @@ boundary f f f
newton off
processors * * 1
-region box block 0. 1. 0. 1. 0. 2. units box
+region box block 0. 1. 0. 1. 0. 10. units box
create_box 1 box
variable skin equal 0.01
@@ -44,7 +44,7 @@ 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 1.75
+create_atoms 1 single 0.5 0.5 9.75
#create_atoms 1 single 0.38 0.05 0.05
set group all diameter 0.25 density 2400
@@ -64,7 +64,7 @@ atom_modify sort 0 0.0
# #fix ins all insert/pack seed 1001 distributiontemplate pdd1 insert_every once &
# # overlapcheck yes particles_in_region 350 region insreg ntry_mc 10000
# fix ins all insert/pack seed 1001 distributiontemplate pdd1 insert_every once &
- # overlapcheck yes volumefraction_region ${v_frac} region insreg ntry_mc 10000
+# overlapcheck yes volumefraction_region ${v_frac} region insreg ntry_mc 10000
# #fix ins all insert/pack seed 1001 distributiontemplate pdd2 insert_every once &
# # overlapcheck yes volumefraction_region 0.05 region insreg ntry_mc 10000
# #fix ins all insert/pack seed 1001 distributiontemplate pdd1 insert_every once &
diff --git a/apps/cpu/LiggghtsApp/in2.lbdem b/apps/cpu/LiggghtsApp/in2.lbdem
index 8d6a0748b..229fe747f 100644
--- a/apps/cpu/LiggghtsApp/in2.lbdem
+++ b/apps/cpu/LiggghtsApp/in2.lbdem
@@ -14,11 +14,13 @@ compute_modify thermo_temp dynamic yes
# particle dump
variable dmp_fname string ${dmp_dir}d_*.liggghts
-# dump dmp all custom ${dmp_stp} ${dmp_fname} &
- # id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
+#dump dmp all custom ${dmp_stp} ${dmp_fname} &
+# id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
-# dump dmp all custom ${dmp_stp} ${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 ${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
diff --git a/apps/cpu/Nozzle/CMakeLists.txt b/apps/cpu/Nozzle/CMakeLists.txt
new file mode 100644
index 000000000..b653be79c
--- /dev/null
+++ b/apps/cpu/Nozzle/CMakeLists.txt
@@ -0,0 +1,3 @@
+PROJECT(Nozzle)
+
+vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} LiggghtsCoupling FILES nozzle.cpp )
diff --git a/apps/cpu/Nozzle/in.nozzle b/apps/cpu/Nozzle/in.nozzle
new file mode 100644
index 000000000..4479b5e3c
--- /dev/null
+++ b/apps/cpu/Nozzle/in.nozzle
@@ -0,0 +1,125 @@
+# shotcrete nozzle simulation
+
+atom_style granular
+atom_modify map array
+boundary f f f
+newton off
+
+communicate single vel yes
+
+units si
+
+#region domain block -1.35 -1.25 0.34 0.44 -0.25 0.25 units box
+region domain block -1.35 -1.25 0.34 0.44 -0.4 0.25 units box
+create_box 1 domain
+
+neighbor 0.002 bin
+neigh_modify delay 0
+
+
+#Material properties required for new pair styles
+
+fix m1 all property/global youngsModulus peratomtype 5.e6
+fix m2 all property/global poissonsRatio peratomtype 0.45
+fix m3 all property/global coefficientRestitution peratomtypepair 1 0.3
+fix m4 all property/global coefficientFriction peratomtypepair 1 0.5
+fix m5 all property/global k_finnie peratomtypepair 1 1.0
+
+# lb coupling fix
+fix lbcoupling all couple/lb/onetoone
+
+#New pair style
+pair_style gran model hertz tangential history #Hertzian without cohesion
+pair_coeff * *
+
+timestep ${t_step}
+
+fix gravi all gravity 9.81 vector 0.0 0.0 -1.0
+
+#the chute
+#variable meshes_dir string d:/Projects/TRR277/Project/WP4/Liggghts/
+variable meshes_dir string d:/Projects/TRR277/Project/WP4/Liggghts/A04/
+
+fix cad1 all mesh/surface file ${meshes_dir}Duese_Acc_Einlass.stl type 1 scale 0.001
+fix cad2 all mesh/surface file ${meshes_dir}Duese_Acc_Verteiler.stl type 1 scale 0.001
+fix cad3 all mesh/surface file ${meshes_dir}Duese_Air_Einlass.stl type 1 scale 0.001
+fix cad4 all mesh/surface file ${meshes_dir}Duese_Air_Verteiler.stl type 1 scale 0.001
+fix cad5 all mesh/surface file ${meshes_dir}Duese_Volcan_Duese.stl type 1 scale 0.001
+fix cad6 all mesh/surface file ${meshes_dir}Duese_Zwischenstueck.stl type 1 element_exclusion_list read list.file scale 0.001 curvature_tolerant yes
+
+fix inface all mesh/surface file ${meshes_dir}InsertDisk2.stl type 1 scale 0.001
+fix wallTop all mesh/surface file ${meshes_dir}InsertDisk3.stl type 1 scale 0.001
+
+#fix granwalls all wall/gran model hertz tangential history mesh n_meshes 7 meshes cad1 cad2 cad3 cad4 cad5 cad6 wallTop
+fix granwalls all wall/gran model hertz tangential history mesh n_meshes 8 meshes cad1 cad2 cad3 cad4 cad5 cad6 wallTop inface
+
+#distributions for insertion
+
+fix pts1 all particletemplate/sphere 15485863 atom_type 1 density constant 2500 radius constant 0.001
+fix pts2 all particletemplate/sphere 15485867 atom_type 1 density constant 2500 radius constant 0.002
+fix pdd1 all particledistribution/discrete 32452843 2 pts1 0.3 pts2 0.7
+
+#region and insertion
+group nve_group region domain
+#region bc cylinder z 0.0 0.0 0.015 0.201 0.23 units box
+#region bc cylinder z 0.0 0.0 10 213 220 units box
+
+region bc cylinder z -1.3013105 0.388582 0.01275005 0.18055 0.20105 units box
+
+#particle insertion
+# fix ins nve_group insert/stream seed 32452867 distributiontemplate pdd1 &
+ # nparticles 6000 massrate 0.1 insert_every 1000 overlapcheck yes all_in no vel constant 0.0 0.0 -1.0 &
+ # insertion_face inface
+
+
+# fix ins nve_group insert/stream seed 32452867 distributiontemplate pdd1 &
+ # nparticles 6000 massrate 0.1 insert_every 1000 overlapcheck yes all_in no vel constant 0.0 0.0 -1.0 &
+ # insertion_face inface
+
+# fix ins nve_group insert/stream seed 32452867 distributiontemplate pdd1 &
+ # nparticles 6000 massrate 0.1 insert_every ones overlapcheck yes all_in no vel constant 0.0 0.0 -1.0 &
+ # insertion_face inface
+
+fix ins nve_group insert/pack seed 32452867 distributiontemplate pdd1 insert_every 1000 &
+ overlapcheck yes volumefraction_region 0.1 region bc ntry_mc 1001
+
+
+# fix ins all insert/stream seed 32452867 distributiontemplate pdd1 &
+ # nparticles INF massrate 0.1 overlapcheck yes all_in yes vel constant 0.0 0.0 -1.0 &
+ # insertion_face inface extrude_length 0.25
+
+
+
+#apply nve integration to all particles that are inserted as single particles
+fix integr nve_group nve/sphere
+
+#output settings, include total thermal energy
+compute 1 all erotate/sphere
+thermo_style custom step atoms ke c_1 vol
+thermo 1000
+thermo_modify lost ignore norm no
+
+variable dmp_time_cad equal 100000000
+
+dump dumpcad1 all mesh/stl ${dmp_time_cad} ${dmp_dir}/cad1_*.stl cad1
+dump dumpcad2 all mesh/stl ${dmp_time_cad} ${dmp_dir}/cad2_*.stl cad2
+dump dumpcad3 all mesh/stl ${dmp_time_cad} ${dmp_dir}/cad3_*.stl cad3
+dump dumpcad4 all mesh/stl ${dmp_time_cad} ${dmp_dir}/cad4_*.stl cad4
+dump dumpcad5 all mesh/stl ${dmp_time_cad} ${dmp_dir}/cad5_*.stl cad5
+dump dumpcad6 all mesh/stl ${dmp_time_cad} ${dmp_dir}/cad6_*.stl cad6
+dump dumpinface all mesh/stl ${dmp_time_cad} ${dmp_dir}/inface_*.stl inface
+dump dumpwallTop all mesh/stl ${dmp_time_cad} ${dmp_dir}/wallTop_*.stl wallTop
+
+#insert the first particles so that dump is not empty
+run 1
+dump dmp all custom/vtk ${dmp_stp} ${dmp_dir}/particles_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
+
+
+#run 1
+#run 1
+
+#insert particles
+#run 100000 upto
+#unfix ins
+
+
diff --git a/apps/cpu/Nozzle/nozzle.cpp b/apps/cpu/Nozzle/nozzle.cpp
new file mode 100644
index 000000000..ea3616be7
--- /dev/null
+++ b/apps/cpu/Nozzle/nozzle.cpp
@@ -0,0 +1,268 @@
+#include <iostream>
+#include <string>
+#include <memory>
+
+#include "VirtualFluids.h"
+
+#include "LiggghtsCouplingCoProcessor.h"
+#include "LiggghtsCouplingWrapper.h"
+#include "IBcumulantK17LBMKernel.h"
+
+using namespace std;
+
+
+int main(int argc, char *argv[])
+{
+ SPtr<Communicator> comm = 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 dx = 1e-3;
+
+ double nuLB = 1e-3;
+ double uLB = -0.01;
+ double rhoLB = 0.0;
+
+ SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
+ //SPtr<LBMKernel> kernel = make_shared<CumulantK17LBMKernel>();
+ SPtr<BCProcessor> bcProc = make_shared<BCProcessor>();
+ kernel->setBCProcessor(bcProc);
+
+ SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
+ noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
+
+ mu::Parser fct;
+ fct.SetExpr("U");
+ fct.DefineConst("U", uLB);
+ SPtr<BCAdapter> inflowBCAdapter(new VelocityBCAdapter(false, false, true, fct, 0, BCFunction::INFCONST));
+ inflowBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new VelocityBCAlgorithm()));
+
+ SPtr<BCAdapter> outflowBCAdapter(new DensityBCAdapter(rhoLB));
+ outflowBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NonEqDensityBCAlgorithm()));
+ //////////////////////////////////////////////////////////////////////////////////
+ // BC visitor
+ BoundaryConditionsBlockVisitor bcVisitor;
+ bcVisitor.addBC(noSlipBCAdapter);
+ bcVisitor.addBC(inflowBCAdapter);
+ bcVisitor.addBC(outflowBCAdapter);
+
+ SPtr<Grid3D> grid = make_shared<Grid3D>(comm);
+ grid->setPeriodicX1(true);
+ grid->setPeriodicX2(true);
+ grid->setPeriodicX3(false);
+ grid->setDeltaX(dx);
+ grid->setBlockNX(blockNX[0], blockNX[1], blockNX[2]);
+
+ string geoPath = "d:/Projects/TRR277/Project/WP4/NozzleGeo";
+
+ string outputPath = "d:/temp/NozzleFlow";
+ UbSystem::makeDirectory(outputPath);
+ UbSystem::makeDirectory(outputPath + "/liggghts");
+
+ 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)
+ 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, noSlipBCAdapter, 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, noSlipBCAdapter, 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, noSlipBCAdapter, 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, noSlipBCAdapter, 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, noSlipBCAdapter, 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, noSlipBCAdapter, 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, noSlipBCAdapter, Interactor3D::SOLID, Interactor3D::EDGES);
+ ///////////////////////////////////////////////////////////
+ //box
+ SPtr<D3Q27Interactor> intrBox = SPtr<D3Q27Interactor>(new D3Q27Interactor(gridCube, grid, noSlipBCAdapter, 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));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoInflow.get(), outputPath + "/geo/geoInflow", WbWriterVtkXmlASCII::getInstance());
+ SPtr<D3Q27Interactor> intrInflow = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoInflow, grid, inflowBCAdapter, 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));
+ if (myid == 0) GbSystem3D::writeGeoObject(geoOutflow.get(), outputPath + "/geo/geoOutflow", WbWriterVtkXmlASCII::getInstance());
+ SPtr<D3Q27Interactor> intrOutflow = SPtr<D3Q27Interactor>(new D3Q27Interactor(geoOutflow, grid, inflowBCAdapter, Interactor3D::SOLID));
+ ///////////////////////////////////////////////////////////
+
+ InteractorsHelper intHelper(grid, metisVisitor, true);
+ intHelper.addInteractor(intrBox);
+ intHelper.addInteractor(intrInflow);
+ intHelper.addInteractor(intrOutflow);
+ intHelper.addInteractor(intrNozzleAirDistributor);
+ intHelper.addInteractor(intrNozzleAirInlet);
+ intHelper.addInteractor(intrNozzleSpacer);
+ intHelper.addInteractor(intrNozzleAccDistributor);
+ intHelper.addInteractor(intrNozzleAccInlet);
+ intHelper.addInteractor(intrNozzleVolcanNozzle1);
+ intHelper.addInteractor(intrNozzleVolcanNozzle2);
+
+
+ intHelper.selectBlocks();
+
+ SPtr<CoProcessor> ppblocks = make_shared<WriteBlocksCoProcessor>(
+ grid, SPtr<UbScheduler>(new UbScheduler(1)), outputPath, WbWriterVtkXmlBinary::getInstance(), comm);
+ ppblocks->process(0);
+ ppblocks.reset();
+
+ if (myid == 0) UBLOG(logINFO, Utilities::toString(grid, comm->getNumberOfProcesses()));
+
+
+ SetKernelBlockVisitor kernelVisitor(kernel, nuLB, comm->getNumberOfProcesses());
+ grid->accept(kernelVisitor);
+
+ intHelper.setBC();
+
+ InitDistributionsBlockVisitor initVisitor;
+ grid->accept(initVisitor);
+
+
+ string inFile1 = "d:/Projects/VirtualFluids_LIGGGHTS_coupling/apps/cpu/Nozzle/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 d_part = 1e-3;
+
+ // 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, 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 = 100;
+ 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);
+
+ SPtr<UbScheduler> lScheduler = make_shared<UbScheduler>(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);
+
+ int numOfThreads = 18;
+ omp_set_num_threads(numOfThreads);
+
+ SPtr<UbScheduler> nupsSch = std::make_shared<UbScheduler>(10, 10, 100);
+ SPtr<NUPSCounterCoProcessor> nupsCoProcessor = make_shared<NUPSCounterCoProcessor>(grid, nupsSch, numOfThreads, comm);
+
+ //// write data for visualization of macroscopic quantities
+ SPtr<UbScheduler> visSch(new UbScheduler(vtkSteps));
+ SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(
+ new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlBinary::getInstance(),
+ SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+ writeMQCoProcessor->process(0);
+
+ int endTime = 1000000; //20;
+ SPtr<Calculator> calculator(new BasicCalculator(grid, lScheduler, endTime));
+ calculator->addCoProcessor(nupsCoProcessor);
+ calculator->addCoProcessor(lcCoProcessor);
+ calculator->addCoProcessor(writeMQCoProcessor);
+
+ if (myid == 0) UBLOG(logINFO, "Simulation-start");
+ calculator->calculate();
+ if (myid == 0) UBLOG(logINFO, "Simulation-end");
+
+
+ return 0;
+}
diff --git a/apps/cpu/ViskomatXL/viskomat.cfg b/apps/cpu/ViskomatXL/viskomat.cfg
index 4227ba9f8..575c244ee 100644
--- a/apps/cpu/ViskomatXL/viskomat.cfg
+++ b/apps/cpu/ViskomatXL/viskomat.cfg
@@ -2,7 +2,7 @@ outputPath = d:/temp/viskomatCylinderRestartTest3_Migration
geoPath = d:/Projects/TRR277/Project/WP1/Rheometer/Aileen
geoFile = fishbone.stl
-numOfThreads = 4
+numOfThreads = 18
availMem = 8e9
logToFile = false
@@ -36,6 +36,7 @@ refineLevel = 0
#nuLB = 1.5e-4
OmegaLB = 1e-4
tau0 = 20e-7
+N = 30
resolution = 32
scaleFactor = 1
diff --git a/apps/cpu/ViskomatXL/viskomat.cpp b/apps/cpu/ViskomatXL/viskomat.cpp
index 3cbc797f8..fa75c3a97 100644
--- a/apps/cpu/ViskomatXL/viskomat.cpp
+++ b/apps/cpu/ViskomatXL/viskomat.cpp
@@ -366,7 +366,7 @@ void bflow(string configname)
wallXminInt->initInteractor();
}
- omp_set_num_threads(numOfThreads);
+ //omp_set_num_threads(numOfThreads);
//set connectors
//InterpolationProcessorPtr iProcessor(new ThixotropyInterpolationProcessor());
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
index da8c93296..e6588e766 100644
--- a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
@@ -16,15 +16,7 @@ LiggghtsCouplingCoProcessor::LiggghtsCouplingCoProcessor(SPtr<Grid3D> grid, SPtr
SPtr<Communicator> comm, LiggghtsCouplingWrapper &wrapper, int demSteps, SPtr<LBMUnitConverter> units)
: CoProcessor(grid, s), comm(comm), wrapper(wrapper), demSteps(demSteps), units(units)
{
- //gridRank = comm->getProcessID();
- //minInitLevel = this->grid->getCoarsestInitializedLevel();
- //maxInitLevel = this->grid->getFinestInitializedLevel();
- //blockVector.resize(maxInitLevel + 1);
-
- //for (int level = minInitLevel; level <= maxInitLevel; level++) {
- // grid->getBlocks(level, gridRank, true, blockVector[level]);
- //}
}
LiggghtsCouplingCoProcessor::~LiggghtsCouplingCoProcessor()
diff --git a/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h b/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h
index 670a597cb..3466d9730 100644
--- a/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h
+++ b/src/cpu/VirtualFluidsCore/Utilities/MemoryUtil.h
@@ -60,6 +60,14 @@
#if defined(__CYGWIN__)
#define MEMORYUTIL_CYGWIN
#endif
+
+#include <iostream>
+#include <sstream>
+#include <string>
+#include <vector>
+#include "Grid3D.h"
+#include "Communicator.h"
+
//////////////////////////////////////////////////////////////////////////
// MemoryUtil
//////////////////////////////////////////////////////////////////////////
@@ -159,6 +167,43 @@ static long long getPhysMemUsedByMe()
}
//////////////////////////////////////////////////////////////////////////
+static std::string toString(SPtr<Grid3D> grid, int numberOfProcesses)
+{
+ unsigned long long numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
+ int ghostLayer = grid->getGhostLayerWidth()*2+1;
+ UbTupleInt3 blockNx = grid->getBlockNX();
+
+ unsigned long long numberOfNodesPerBlock = (unsigned long long)(val<1>(blockNx)) *
+ (unsigned long long)(val<2>(blockNx)) *
+ (unsigned long long)(val<3>(blockNx));
+ unsigned long long numberOfNodes = numberOfBlocks * numberOfNodesPerBlock;
+ unsigned long long numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (val<1>(blockNx) + ghostLayer) *
+ (val<2>(blockNx) + ghostLayer) *
+ (val<3>(blockNx) + ghostLayer);
+ double needMemAll = double(numberOfNodesPerBlockWithGhostLayer*(27*sizeof(double)+sizeof(int)+sizeof(float)*4));
+ double needMem = needMemAll / double(numberOfProcesses);
+
+ std::ostringstream out;
+ out << "Grid information:" << std::endl;
+ out << "###################################################" << std::endl;
+ out << "# Number of blocks = " << numberOfBlocks << std::endl;
+ out << "# Number of nodes = " << numberOfNodes << std::endl;
+ int minInitLevel = grid->getCoarsestInitializedLevel();
+ int maxInitLevel = grid->getFinestInitializedLevel();
+ for (int level = minInitLevel; level<=maxInitLevel; level++)
+ {
+ int nobl = grid->getNumberOfBlocks(level);
+ out << "# Number of blocks for level " << level << " = " << nobl << std::endl;
+ out << "# Number of nodes for level " << level << " = " << nobl * numberOfNodesPerBlock << std::endl;
+ }
+ out << "# Necessary memory = " << needMemAll << " bytes" << std::endl;
+ out << "# Necessary memory per process = " << needMem << " bytes" << std::endl;
+ out << "# Available memory per process = " << (double)getTotalPhysMem() << " bytes" << std::endl;
+ out << "###################################################" << std::endl;
+
+ return out.str();
+}
+
} // namespace Utilities
#endif
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
index 5c813d289..b55c405d6 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
@@ -53,7 +53,7 @@ SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue
}
}
-SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, int &numberOfProcesses,
+SetKernelBlockVisitor::SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, int numberOfProcesses,
SetKernelBlockVisitor::Action action)
: Block3DVisitor(0, Grid3DSystem::MAXLEVEL), kernel(std::move(kernel)), nue(nue), action(action), dataSetFlag(true),
numberOfProcesses(numberOfProcesses)
@@ -127,7 +127,7 @@ double SetKernelBlockVisitor::getRequiredPhysicalMemory(const SPtr<Grid3D> &grid
unsigned long long numberOfNodesPerBlockWithGhostLayer;
auto numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
auto blockNx = grid->getBlockNX();
- int ghostLayer = 3;
+ int ghostLayer = grid->getGhostLayerWidth() * 2 + 1;
numberOfNodesPerBlockWithGhostLayer = numberOfBlocks * (val<1>(blockNx) + ghostLayer) *
(val<2>(blockNx) + ghostLayer) * (val<3>(blockNx) + ghostLayer);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
index 51cbc256c..1e0621f22 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.h
@@ -52,7 +52,7 @@ public:
SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, double availMem, double needMem,
SetKernelBlockVisitor::Action action = SetKernelBlockVisitor::NewKernel);
- SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, int &numberOfProcesses,
+ SetKernelBlockVisitor(SPtr<LBMKernel> kernel, LBMReal nue, int numberOfProcesses,
SetKernelBlockVisitor::Action action = SetKernelBlockVisitor::NewKernel);
~SetKernelBlockVisitor() override = default;
--
GitLab