diff --git a/CMake/cmake_config_files/BOMBADIL.config.cmake b/CMake/cmake_config_files/BOMBADIL.config.cmake
index c0c5cf2f08b2cc925be248441c99a336160fd1bd..9c4bd4ecffab1e63161343ecc493eb9d9bc951a4 100644
--- a/CMake/cmake_config_files/BOMBADIL.config.cmake
+++ b/CMake/cmake_config_files/BOMBADIL.config.cmake
@@ -26,6 +26,15 @@ SET(BOOST_LIBRARYDIR ${BOOST_ROOT}"/stageMSVC64/lib")
#################################################################################
set(VTK_DIR "d:/Tools/VTK/build/VTK-8.0.0")
+#################################################################################
+# LIGGGHTS
+#################################################################################
+set(LIGGGHTS_SOURCE_DIR "d:/Tools/LIGGGHTS/src")
+set(LIGGGHTS_DEBUG_LIBRARY "d:/Tools/LIGGGHTS/build/Debug/liggghts.lib")
+set(LIGGGHTS_RELEASE_LIBRARY "d:/Tools/LIGGGHTS/build/Release/liggghts.lib")
+
+
+
#################################################################################
# METIS
#################################################################################
diff --git a/apps/cpu/Applications.cmake b/apps/cpu/Applications.cmake
index 436db3c5d3cc2f9c3bc867ce633bb6d9585ef0a8..b63689a9dc0d178b7d82c32a0c67a1b4b7946328 100644
--- a/apps/cpu/Applications.cmake
+++ b/apps/cpu/Applications.cmake
@@ -11,6 +11,9 @@ add_subdirectory(${APPS_ROOT_CPU}/LaminarTubeFlow)
add_subdirectory(${APPS_ROOT_CPU}/MultiphaseDropletTest)
add_subdirectory(${APPS_ROOT_CPU}/RisingBubble2D)
add_subdirectory(${APPS_ROOT_CPU}/JetBreakup)
+add_subdirectory(${APPS_ROOT_CPU}/LiggghtsApp)
+add_subdirectory(${APPS_ROOT_CPU}/FallingSphere)
+add_subdirectory(${APPS_ROOT_CPU}/Nozzle)
#add_subdirectory(tests)
#add_subdirectory(Applications/gridRf)
#add_subdirectory(Applications/greenvortex)
diff --git a/apps/cpu/LiggghtsApp/CMakeLists.txt b/apps/cpu/LiggghtsApp/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..f3a2d925f1d082c8f2e9e52e31d8179fe82c9235
--- /dev/null
+++ b/apps/cpu/LiggghtsApp/CMakeLists.txt
@@ -0,0 +1,3 @@
+PROJECT(LiggghtsApp)
+
+vf_add_library(BUILDTYPE binary PRIVATE_LINK VirtualFluidsCore basics ${MPI_CXX_LIBRARIES} LiggghtsCoupling FILES LiggghtsApp.cpp )
diff --git a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..978ae376f1c842b0d108730fb9d2ffe43fa65537
--- /dev/null
+++ b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
@@ -0,0 +1,259 @@
+#include <iostream>
+#include <string>
+#include <memory>
+
+#include "VirtualFluids.h"
+
+//#include "lammps.h"
+//#include "input.h"
+//#include "atom.h"
+//#include "modify.h"
+//#include "fix_lb_coupling_onetoone.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 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);
+
+
+
+
+ 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/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)
+ 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/LiggghtsApp/in.lbdem";
+ //string inFile1 = "d:/Tools/LIGGGHTS/examples/LIGGGHTS/Tutorials_public/chute_wear/in.chute_wear2";
+ 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 = 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;
+
+ 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 = outputPath; // "d:/temp/lll2/";
+
+ wrapper.execCommand("echo log");
+
+ 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 *)inFile2.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 = 1000; //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");
+
+ //MPI_Init(&argc, &argv);
+ //MPI_Comm mpi_comm = *(MPI_Comm*)(comm->getNativeCommunicator());
+ //LiggghtsCouplingWrapper wrapper(argv, mpi_comm);
+
+ //wrapper.execFile("in2.lbdem");
+ //wrapper.runUpto(demSubsteps - 1);
+
+ //LAMMPS_NS::LAMMPS *lmp;
+ // // custom argument vector for LAMMPS library
+ // const char *lmpargv[] {"liblammps", "-log", "none"};
+ // int lmpargc = sizeof(lmpargv)/sizeof(const char *);
+
+ // // explicitly initialize MPI
+ // MPI_Init(&argc, &argv);
+
+ // // create LAMMPS instance
+ // lmp = new LAMMPS_NS::LAMMPS(lmpargc, (char **)lmpargv, MPI_COMM_WORLD);
+ // lmp->input->file("in.lbdem");
+ // //lmp->input->one("run 1");
+ //
+ // //# Try extracting a global value
+ // // print("")
+ // // print("Attempting to get the number of atoms in simulation")
+ // // numAtoms = lmp.extract_global("natoms", 0)
+ // // print("natoms =", numAtoms)
+
+ // // # Try extracting atom's positions
+ // // print("")
+ // // print("Attempting to get the atom's positions")
+ // // pos = lmp.extract_atom("x",3)
+ // // for k in range(0,numAtoms):
+ // // print("Pos[%i] = [%f, %f, %f]" % (k, pos[k][0], pos[k][1], pos[k][2]))
+
+ // LAMMPS_NS::FixLbCouplingOnetoone
+ // *couplingFix
+ // = dynamic_cast<LAMMPS_NS::FixLbCouplingOnetoone*>
+ // (lmp->modify->find_fix_style("couple/lb/onetoone",0));
+
+ // cout << "test1\n";
+ //
+ // //double **t_liggghts = couplingFix->get_torque_ptr();
+ // cout << "test2\n";
+
+ // lmp->input->one("run 9 upto");
+
+ // for (int step = 0; step < 10; step++)
+ // {
+ //
+
+ // int numAtoms = lmp->atom->natoms;
+
+ // //double** pos = (double**)lmp->atom->extract("x");
+ // double** pos = lmp->atom->x;
+ //
+ // //double* forceX = lmp->atom->fx;
+
+ // for (int i = 0; i < numAtoms; i++)
+ // {
+ // double **f_liggghts = couplingFix->get_force_ptr();
+ // double** force = lmp->atom->f;
+ // cout << "Pos[" << i << "] = [" << pos[i][0] << ", " << pos[i][1] << ", " << pos[i][2] << "]\n";
+ // cout << "Force1[" << i << "] = [" << f_liggghts[i][0] << ", " << f_liggghts[i][1] << ", " << f_liggghts[i][2] << "]\n";
+ // f_liggghts[i][0] += 0;
+ // f_liggghts[i][1] += 0;
+ // f_liggghts[i][2] += 500;
+ // cout << "Force2[" << i << "] = [" << force[i][0] << ", " << force[i][1] << ", " << force[i][2] << "]\n";
+ // }
+
+ // couplingFix->comm_force_torque();
+
+ // lmp->input->one("run 10000");
+ //
+ // }
+
+ // // delete LAMMPS instance
+ // delete lmp;
+
+ // // stop MPI environment
+ //MPI_Finalize();
+ return 0;
+}
diff --git a/apps/cpu/LiggghtsApp/in.lbdem b/apps/cpu/LiggghtsApp/in.lbdem
new file mode 100644
index 0000000000000000000000000000000000000000..ae2baa373dc75edec41634522e53ade31f459b4d
--- /dev/null
+++ b/apps/cpu/LiggghtsApp/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. 2. 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 1.75
+#create_atoms 1 single 0.38 0.05 0.05
+
+set group all diameter 0.25 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/LiggghtsApp/in2.lbdem b/apps/cpu/LiggghtsApp/in2.lbdem
new file mode 100644
index 0000000000000000000000000000000000000000..8d6a0748befeea97d6faa7efa87ffd60e5c56cfd
--- /dev/null
+++ b/apps/cpu/LiggghtsApp/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}/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/cpu.cmake b/cpu.cmake
index bad575aa3854fb2be777e90bba9ea71fb2cbba5a..e5b77c8a43c40c4987f3d2d00d34833c7cb27c61 100644
--- a/cpu.cmake
+++ b/cpu.cmake
@@ -29,6 +29,8 @@ SET(USE_CATALYST OFF CACHE BOOL "include Paraview Catalyst support")
SET(USE_HLRN_LUSTRE OFF CACHE BOOL "include HLRN Lustre support")
SET(USE_DEM_COUPLING OFF CACHE BOOL "PE plugin")
+SET(USE_LIGGGHTS ON CACHE BOOL "include LIGGGHTS library support")
+
#MPI
IF((NOT ${CMAKE_CXX_COMPILER} MATCHES mpicxx) AND (NOT ${CMAKE_CXX_COMPILER} MATCHES mpiicpc))# OR NOT ${CMAKE_CXX_COMPILER} MATCHES cc OR NOT ${CMAKE_CXX_COMPILER} MATCHES mpiCC)
FIND_PACKAGE(MPI REQUIRED)
@@ -85,5 +87,9 @@ if(BUILD_VF_PYTHON_BINDINGS)
add_subdirectory(src/cpu/pythonbindings)
endif()
+if(USE_LIGGGHTS)
+ add_subdirectory(src/cpu/LiggghtsCoupling)
+endif()
+
set (APPS_ROOT_CPU "${VF_ROOT_DIR}/apps/cpu/")
include(${APPS_ROOT_CPU}/Applications.cmake)
\ No newline at end of file
diff --git a/src/cpu/LiggghtsCoupling/CMakeLists.txt b/src/cpu/LiggghtsCoupling/CMakeLists.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ed9d1f0e2bd8d0302f1d37aa3bde8a120ceb2312
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/CMakeLists.txt
@@ -0,0 +1,10 @@
+
+set(LIGGGHTS_LIBRARIES optimized ${LIGGGHTS_RELEASE_LIBRARY} debug ${LIGGGHTS_DEBUG_LIBRARY})
+
+vf_add_library(BUILDTYPE static PUBLIC_LINK basics muparser MPI::MPI_CXX VirtualFluidsCore ${LIGGGHTS_LIBRARIES})
+
+vf_get_library_name(library_name)
+
+#target_link_directories(${library_name} PUBLIC ${LIGGGHTS_BINARY_DIR})
+target_include_directories(${library_name} PUBLIC ${LIGGGHTS_SOURCE_DIR})
+
diff --git a/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.cpp b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..ac0bc35844cb32544d901cbf5a93503d49de4070
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.cpp
@@ -0,0 +1,913 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// 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 IBcumulantK17LBMKernel.cpp
+//! \ingroup LBM
+//! \author Konstantin Kutscher, Martin Geier
+//=======================================================================================
+#include "IBcumulantK17LBMKernel.h"
+#include "D3Q27System.h"
+#include "D3Q27EsoTwist3DSplittedVector.h"
+#include <cmath>
+#include "DataSet3D.h"
+#include "LBMKernel.h"
+#include "Block3D.h"
+#include "BCArray3D.h"
+
+#define PROOF_CORRECTNESS
+
+using namespace UbMath;
+
+//////////////////////////////////////////////////////////////////////////
+IBcumulantK17LBMKernel::IBcumulantK17LBMKernel()
+{
+ this->compressible = true;
+}
+//////////////////////////////////////////////////////////////////////////
+void IBcumulantK17LBMKernel::initDataSet()
+{
+ SPtr<DistributionArray3D> d(new D3Q27EsoTwist3DSplittedVector(nx[0] + 2, nx[1] + 2, nx[2] + 2, -999.9));
+ dataSet->setFdistributions(d);
+
+ particleData =
+ std::make_shared<CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>>(nx[0] + 2, nx[1] + 2, nx[2] + 2);
+
+ int minX1 = 0;
+ int minX2 = 0;
+ int minX3 = 0;
+ int maxX1 = nx[0]+2;
+ int maxX2 = nx[1]+2;
+ int maxX3 = nx[2]+2;
+
+ LBMReal omega = collFactor;
+
+ for (int x3 = minX3; x3 < maxX3; x3++) {
+ for (int x2 = minX2; x2 < maxX2; x2++) {
+ for (int x1 = minX1; x1 < maxX1; x1++) {
+ (*particleData)(x1, x2, x3) = std::make_shared<IBdynamicsParticleData>();
+ }
+ }
+ }
+
+}
+//////////////////////////////////////////////////////////////////////////
+SPtr<LBMKernel> IBcumulantK17LBMKernel::clone()
+{
+ SPtr<LBMKernel> kernel(new IBcumulantK17LBMKernel());
+ kernel->setNX(nx);
+ std::dynamic_pointer_cast<IBcumulantK17LBMKernel>(kernel)->initDataSet();
+ kernel->setCollisionFactor(this->collFactor);
+ kernel->setBCProcessor(bcProcessor->clone(kernel));
+ kernel->setWithForcing(withForcing);
+ kernel->setForcingX1(muForcingX1);
+ kernel->setForcingX2(muForcingX2);
+ kernel->setForcingX3(muForcingX3);
+ kernel->setIndex(ix1, ix2, ix3);
+ kernel->setDeltaT(deltaT);
+ kernel->setBlock(block.lock());
+
+ return kernel;
+}
+//////////////////////////////////////////////////////////////////////////
+void IBcumulantK17LBMKernel::calculate(int step)
+{
+ //////////////////////////////////////////////////////////////////////////
+ //! Cumulant K17 Kernel is based on
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //! and
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.07.004"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.07.004 ]</b></a>
+ //!
+ //! The cumulant kernel is executed in the following steps
+ //!
+ ////////////////////////////////////////////////////////////////////////////////
+ //! - Get node index coordinates from thredIdx, blockIdx, blockDim and gridDim.
+ //!
+
+ using namespace std;
+
+ //initializing of forcing stuff
+ if (withForcing)
+ {
+ muForcingX1.DefineVar("x1", &muX1); muForcingX1.DefineVar("x2", &muX2); muForcingX1.DefineVar("x3", &muX3);
+ muForcingX2.DefineVar("x1", &muX1); muForcingX2.DefineVar("x2", &muX2); muForcingX2.DefineVar("x3", &muX3);
+ muForcingX3.DefineVar("x1", &muX1); muForcingX3.DefineVar("x2", &muX2); muForcingX3.DefineVar("x3", &muX3);
+
+ muDeltaT = deltaT;
+
+ muForcingX1.DefineVar("dt", &muDeltaT);
+ muForcingX2.DefineVar("dt", &muDeltaT);
+ muForcingX3.DefineVar("dt", &muDeltaT);
+
+ muNu = (1.0 / 3.0) * (1.0 / collFactor - 1.0 / 2.0);
+
+ muForcingX1.DefineVar("nu", &muNu);
+ muForcingX2.DefineVar("nu", &muNu);
+ muForcingX3.DefineVar("nu", &muNu);
+ }
+ /////////////////////////////////////
+
+ localDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
+ nonLocalDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
+ restDistributions = dynamic_pointer_cast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getZeroDistributions();
+
+ SPtr<BCArray3D> bcArray = this->getBCProcessor()->getBCArray();
+
+ const int bcArrayMaxX1 = (int)bcArray->getNX1();
+ const int bcArrayMaxX2 = (int)bcArray->getNX2();
+ const int bcArrayMaxX3 = (int)bcArray->getNX3();
+
+ int minX1 = ghostLayerWidth;
+ int minX2 = ghostLayerWidth;
+ int minX3 = ghostLayerWidth;
+ int maxX1 = bcArrayMaxX1 - ghostLayerWidth;
+ int maxX2 = bcArrayMaxX2 - ghostLayerWidth;
+ int maxX3 = bcArrayMaxX3 - ghostLayerWidth;
+
+ LBMReal omega = collFactor;
+
+
+
+ for (int x3 = minX3; x3 < maxX3; x3++)
+ {
+ for (int x2 = minX2; x2 < maxX2; x2++)
+ {
+ for (int x1 = minX1; x1 < maxX1; x1++)
+ {
+ if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3)) {
+ int x1p = x1 + 1;
+ int x2p = x2 + 1;
+ int x3p = x3 + 1;
+ //////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+ //! - Read distributions: style of reading and writing the distributions from/to stored arrays
+ //! dependent on timestep is based on the esoteric twist algorithm <a
+ //! href="https://doi.org/10.3390/computation5020019"><b>[ M. Geier et al. (2017),
+ //! DOI:10.3390/computation5020019 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////
+ //////////////////////////////////////////////////////////////////////////
+
+ // E N T
+ // c c c
+ //////////
+ // W S B
+ // a a a
+
+ // Rest is b
+
+ // mfxyz
+ // a - negative
+ // b - null
+ // c - positive
+
+ // 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 f[D3Q27System::ENDF + 1];
+ LBMReal fEq[D3Q27System::ENDF + 1];
+ LBMReal fEqSolid[D3Q27System::ENDF + 1];
+ LBMReal fPre[D3Q27System::ENDF + 1];
+
+ f[D3Q27System::REST] = mfbbb;
+
+ f[D3Q27System::E] = mfcbb;
+ f[D3Q27System::N] = mfbcb;
+ f[D3Q27System::T] = mfbbc;
+ f[D3Q27System::NE] = mfccb;
+ f[D3Q27System::NW] = mfacb;
+ f[D3Q27System::TE] = mfcbc;
+ f[D3Q27System::TW] = mfabc;
+ f[D3Q27System::TN] = mfbcc;
+ f[D3Q27System::TS] = mfbac;
+ f[D3Q27System::TNE] = mfccc;
+ f[D3Q27System::TNW] = mfacc;
+ f[D3Q27System::TSE] = mfcac;
+ f[D3Q27System::TSW] = mfaac;
+
+ f[D3Q27System::W] = mfabb;
+ f[D3Q27System::S] = mfbab;
+ f[D3Q27System::B] = mfbba;
+ f[D3Q27System::SW] = mfaab;
+ f[D3Q27System::SE] = mfcab;
+ f[D3Q27System::BW] = mfaba;
+ f[D3Q27System::BE] = mfcba;
+ f[D3Q27System::BS] = mfbaa;
+ f[D3Q27System::BN] = mfbca;
+ f[D3Q27System::BSW] = mfaaa;
+ f[D3Q27System::BSE] = mfcaa;
+ f[D3Q27System::BNW] = mfaca;
+ f[D3Q27System::BNE] = mfcca;
+
+ if ((*particleData)(x1, x2, x3)->solidFraction > SOLFRAC_MIN) {
+ fPre[D3Q27System::REST] = mfbbb;
+
+ fPre[D3Q27System::E] = mfcbb;
+ fPre[D3Q27System::N] = mfbcb;
+ fPre[D3Q27System::T] = mfbbc;
+ fPre[D3Q27System::NE] = mfccb;
+ fPre[D3Q27System::NW] = mfacb;
+ fPre[D3Q27System::TE] = mfcbc;
+ fPre[D3Q27System::TW] = mfabc;
+ fPre[D3Q27System::TN] = mfbcc;
+ fPre[D3Q27System::TS] = mfbac;
+ fPre[D3Q27System::TNE] = mfccc;
+ fPre[D3Q27System::TNW] = mfacc;
+ fPre[D3Q27System::TSE] = mfcac;
+ fPre[D3Q27System::TSW] = mfaac;
+
+ fPre[D3Q27System::W] = mfabb;
+ fPre[D3Q27System::S] = mfbab;
+ fPre[D3Q27System::B] = mfbba;
+ fPre[D3Q27System::SW] = mfaab;
+ fPre[D3Q27System::SE] = mfcab;
+ fPre[D3Q27System::BW] = mfaba;
+ fPre[D3Q27System::BE] = mfcba;
+ fPre[D3Q27System::BS] = mfbaa;
+ fPre[D3Q27System::BN] = mfbca;
+ fPre[D3Q27System::BSW] = mfaaa;
+ fPre[D3Q27System::BSE] = mfcaa;
+ fPre[D3Q27System::BNW] = mfaca;
+ fPre[D3Q27System::BNE] = mfcca;
+ }
+
+ (*particleData)(x1, x2, x3)->hydrodynamicForce.fill(0.0);
+
+ if ((*particleData)(x1, x2, x3)->solidFraction <= SOLFRAC_MAX) {
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Calculate density and velocity using pyramid summation for low round-off errors as in Eq.
+ //! (J1)-(J3) <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
+ //! DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
+ //!
+ LBMReal drho = ((((mfccc + mfaaa) + (mfaca + mfcac)) + ((mfacc + mfcaa) + (mfaac + mfcca))) +
+ (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) +
+ ((mfacb + mfcab) + (mfaab + mfccb))) +
+ ((mfabb + mfcbb) + (mfbab + mfbcb)) + (mfbba + mfbbc)) +
+ mfbbb;
+
+ LBMReal rho = c1 + drho;
+ LBMReal OOrho = c1 / rho;
+ ////////////////////////////////////////////////////////////////////////////////////
+ LBMReal vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
+ (((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
+ (mfcbb - mfabb)) /
+ rho;
+ LBMReal vvy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
+ (((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
+ (mfbcb - mfbab)) /
+ rho;
+ LBMReal vvz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
+ (((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
+ (mfbbc - mfbba)) /
+ rho;
+ ////////////////////////////////////////////////////////////////////////////////////
+ // forcing
+ ///////////////////////////////////////////////////////////////////////////////////////////
+ if (withForcing) {
+ muX1 = static_cast<double>(x1 - 1 + ix1 * maxX1);
+ muX2 = static_cast<double>(x2 - 1 + ix2 * maxX2);
+ muX3 = static_cast<double>(x3 - 1 + ix3 * maxX3);
+
+ forcingX1 = muForcingX1.Eval();
+ forcingX2 = muForcingX2.Eval();
+ forcingX3 = muForcingX3.Eval();
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Add half of the acceleration (body force) to the velocity as in Eq. (42)
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
+ //! DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
+ //!
+ vvx += forcingX1 * deltaT * c1o2; // X
+ vvy += forcingX2 * deltaT * c1o2; // Y
+ vvz += forcingX3 * deltaT * c1o2; // Z
+ }
+ ////////////////////////////////////////////////////////////////////////////////////
+ // calculate the square of velocities for this lattice node
+ LBMReal vx2 = vvx * vvx;
+ LBMReal vy2 = vvy * vvy;
+ LBMReal vz2 = vvz * vvz;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Set relaxation limiters for third order cumulants to default value \f$ \lambda=0.001 \f$
+ //! according to section 6 in <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et
+ //! al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ LBMReal wadjust;
+ LBMReal qudricLimitP = c1o100;
+ LBMReal qudricLimitM = c1o100;
+ LBMReal qudricLimitD = c1o100;
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from well conditioned distributions to central moments as defined in
+ //! Appendix J in <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
+ //! DOI:10.1016/j.camwa.2015.05.001 ]</b></a> see also Eq. (6)-(14) in <a
+ //! href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ forwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36, c1o36);
+ forwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36, c1o36);
+ forwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
+ forwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36, c1o36);
+ forwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9, c1o9);
+ forwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36, c1o36);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ forwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6, c1o6);
+ forwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
+ forwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18, c1o18);
+ forwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
+ forwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
+ forwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
+ forwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6, c1o6);
+ forwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
+ forwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18, c1o18);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ forwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1, c1);
+ forwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
+ forwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3, c1o3);
+ forwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
+ forwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
+ forwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
+ forwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3, c1o3);
+ forwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
+ forwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9, c1o9);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Setting relaxation rates for non-hydrodynamic cumulants (default values). Variable names and
+ //! equations according to <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al.
+ //! (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //! => [NAME IN PAPER]=[NAME IN CODE]=[DEFAULT VALUE].
+ //! - Trace of second order cumulants \f$ C_{200}+C_{020}+C_{002} \f$ used to adjust bulk
+ //! viscosity:\f$\omega_2=OxxPyyPzz=1.0 \f$.
+ //! - Third order cumulants \f$ C_{120}+C_{102} \f$, \f$ C_{210}+C_{012} \f$, \f$ C_{201}+C_{021}
+ //! \f$: \f$\omega_3=OxyyPxzz\f$ set according to Eq. (111) with simplifications assuming
+ //! \f$\omega_2=1.0\f$.
+ //! - Third order cumulants \f$ C_{120}-C_{102} \f$, \f$ C_{210}-C_{012} \f$, \f$ C_{201}-C_{021}
+ //! \f$: \f$\omega_4 = OxyyMxzz\f$ set according to Eq. (112) with simplifications assuming
+ //! \f$\omega_2 = 1.0\f$.
+ //! - Third order cumulants \f$ C_{111} \f$: \f$\omega_5 = Oxyz\f$ set according to Eq. (113) with
+ //! simplifications assuming \f$\omega_2 = 1.0\f$ (modify for different bulk viscosity).
+ //! - Fourth order cumulants \f$ C_{220} \f$, \f$ C_{202} \f$, \f$ C_{022} \f$, \f$ C_{211} \f$,
+ //! \f$ C_{121} \f$, \f$ C_{112} \f$: for simplification all set to the same default value \f$
+ //! \omega_6=\omega_7=\omega_8=O4=1.0 \f$.
+ //! - Fifth order cumulants \f$ C_{221}\f$, \f$C_{212}\f$, \f$C_{122}\f$: \f$\omega_9=O5=1.0\f$.
+ //! - Sixth order cumulant \f$ C_{222}\f$: \f$\omega_{10}=O6=1.0\f$.
+ //!
+ ////////////////////////////////////////////////////////////
+ // 2.
+ LBMReal OxxPyyPzz = c1;
+ ////////////////////////////////////////////////////////////
+ // 3.
+ LBMReal OxyyPxzz =
+ c8 * (-c2 + omega) * (c1 + c2 * omega) / (-c8 - c14 * omega + c7 * omega * omega);
+ LBMReal OxyyMxzz =
+ c8 * (-c2 + omega) * (-c7 + c4 * omega) / (c56 - c50 * omega + c9 * omega * omega);
+ LBMReal Oxyz = c24 * (-c2 + omega) * (-c2 - c7 * omega + c3 * omega * omega) /
+ (c48 + c152 * omega - c130 * omega * omega + c29 * omega * omega * omega);
+ ////////////////////////////////////////////////////////////
+ // 4.
+ LBMReal O4 = c1;
+ ////////////////////////////////////////////////////////////
+ // 5.
+ LBMReal O5 = c1;
+ ////////////////////////////////////////////////////////////
+ // 6.
+ LBMReal O6 = c1;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - A and B: parameters for fourth order convergence of the diffusion term according to Eq. (114)
+ //! and (115) <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a> with simplifications assuming \f$\omega_2 = 1.0\f$
+ //! (modify for different bulk viscosity).
+ //!
+ LBMReal A = (c4 + c2 * omega - c3 * omega * omega) / (c2 - c7 * omega + c5 * omega * omega);
+ LBMReal B = (c4 + c28 * omega - c14 * omega * omega) / (c6 - c21 * omega + c15 * omega * omega);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute cumulants from central moments according to Eq. (20)-(23) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////
+ // 4.
+ LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + c2 * mfbba * mfbab) * OOrho;
+ LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + c2 * mfbba * mfabb) * OOrho;
+ LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + c2 * mfbab * mfabb) * OOrho;
+
+ LBMReal CUMcca = mfcca - (((mfcaa * mfaca + c2 * mfbba * mfbba) + c1o3 * (mfcaa + mfaca)) * OOrho -
+ c1o9 * (drho * OOrho));
+ LBMReal CUMcac = mfcac - (((mfcaa * mfaac + c2 * mfbab * mfbab) + c1o3 * (mfcaa + mfaac)) * OOrho -
+ c1o9 * (drho * OOrho));
+ LBMReal CUMacc = mfacc - (((mfaac * mfaca + c2 * mfabb * mfabb) + c1o3 * (mfaac + mfaca)) * OOrho -
+ c1o9 * (drho * OOrho));
+ ////////////////////////////////////////////////////////////
+ // 5.
+ LBMReal CUMbcc =
+ mfbcc -
+ ((mfaac * mfbca + mfaca * mfbac + c4 * mfabb * mfbbb + c2 * (mfbab * mfacb + mfbba * mfabc)) +
+ c1o3 * (mfbca + mfbac)) *
+ OOrho;
+ LBMReal CUMcbc =
+ mfcbc -
+ ((mfaac * mfcba + mfcaa * mfabc + c4 * mfbab * mfbbb + c2 * (mfabb * mfcab + mfbba * mfbac)) +
+ c1o3 * (mfcba + mfabc)) *
+ OOrho;
+ LBMReal CUMccb =
+ mfccb -
+ ((mfcaa * mfacb + mfaca * mfcab + c4 * mfbba * mfbbb + c2 * (mfbab * mfbca + mfabb * mfcba)) +
+ c1o3 * (mfacb + mfcab)) *
+ OOrho;
+ ////////////////////////////////////////////////////////////
+ // 6.
+ LBMReal CUMccc =
+ mfccc + ((-c4 * mfbbb * mfbbb - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca) -
+ c4 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc) -
+ c2 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) *
+ OOrho +
+ (c4 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac) +
+ c2 * (mfcaa * mfaca * mfaac) + c16 * mfbba * mfbab * mfabb) *
+ OOrho * OOrho -
+ c1o3 * (mfacc + mfcac + mfcca) * OOrho - c1o9 * (mfcaa + mfaca + mfaac) * OOrho +
+ (c2 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) +
+ (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 * (mfaac + mfaca + mfcaa)) *
+ OOrho * OOrho * c2o3 +
+ c1o27 * ((drho * drho - drho) * OOrho * OOrho));
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute linear combinations of second and third order cumulants
+ //!
+ ////////////////////////////////////////////////////////////
+ // 2.
+ LBMReal mxxPyyPzz = mfcaa + mfaca + mfaac;
+ LBMReal mxxMyy = mfcaa - mfaca;
+ LBMReal mxxMzz = mfcaa - mfaac;
+ ////////////////////////////////////////////////////////////
+ // 3.
+ LBMReal mxxyPyzz = mfcba + mfabc;
+ LBMReal mxxyMyzz = mfcba - mfabc;
+
+ LBMReal mxxzPyyz = mfcab + mfacb;
+ LBMReal mxxzMyyz = mfcab - mfacb;
+
+ LBMReal mxyyPxzz = mfbca + mfbac;
+ LBMReal mxyyMxzz = mfbca - mfbac;
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // incl. correction
+ ////////////////////////////////////////////////////////////
+ //! - Compute velocity gradients from second order cumulants according to Eq. (27)-(32)
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a> Further explanations of the correction in viscosity in
+ //! Appendix H of <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
+ //! DOI:10.1016/j.camwa.2015.05.001 ]</b></a> Note that the division by rho is omitted here as we
+ //! need rho times the gradients later.
+ //!
+ LBMReal Dxy = -c3 * omega * mfbba;
+ LBMReal Dxz = -c3 * omega * mfbab;
+ LBMReal Dyz = -c3 * omega * mfabb;
+ LBMReal dxux = c1o2 * (-omega) * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
+ LBMReal dyuy = dxux + omega * c3o2 * mxxMyy;
+ LBMReal dzuz = dxux + omega * c3o2 * mxxMzz;
+ ////////////////////////////////////////////////////////////
+ //! - Relaxation of second order cumulants with correction terms according to Eq. (33)-(35) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ mxxPyyPzz += OxxPyyPzz * (mfaaa - mxxPyyPzz) -
+ c3 * (c1 - c1o2 * OxxPyyPzz) * (vx2 * dxux + vy2 * dyuy + vz2 * dzuz);
+ mxxMyy += omega * (-mxxMyy) - c3 * (c1 + c1o2 * (-omega)) * (vx2 * dxux - vy2 * dyuy);
+ mxxMzz += omega * (-mxxMzz) - c3 * (c1 + c1o2 * (-omega)) * (vx2 * dxux - vz2 * dzuz);
+
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ ////no correction
+ // mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz);
+ // mxxMyy += -(-omega) * (-mxxMyy);
+ // mxxMzz += -(-omega) * (-mxxMzz);
+ /////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
+ mfabb += omega * (-mfabb);
+ mfbab += omega * (-mfbab);
+ mfbba += omega * (-mfbba);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // relax
+ //////////////////////////////////////////////////////////////////////////
+ // incl. limiter
+ //! - Relaxation of third order cumulants including limiter according to Eq. (116)-(123)
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ wadjust = Oxyz + (c1 - Oxyz) * abs(mfbbb) / (abs(mfbbb) + qudricLimitD);
+ mfbbb += wadjust * (-mfbbb);
+ wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxxyPyzz) / (abs(mxxyPyzz) + qudricLimitP);
+ mxxyPyzz += wadjust * (-mxxyPyzz);
+ wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxxyMyzz) / (abs(mxxyMyzz) + qudricLimitM);
+ mxxyMyzz += wadjust * (-mxxyMyzz);
+ wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxxzPyyz) / (abs(mxxzPyyz) + qudricLimitP);
+ mxxzPyyz += wadjust * (-mxxzPyyz);
+ wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxxzMyyz) / (abs(mxxzMyyz) + qudricLimitM);
+ mxxzMyyz += wadjust * (-mxxzMyyz);
+ wadjust = OxyyPxzz + (c1 - OxyyPxzz) * abs(mxyyPxzz) / (abs(mxyyPxzz) + qudricLimitP);
+ mxyyPxzz += wadjust * (-mxyyPxzz);
+ wadjust = OxyyMxzz + (c1 - OxyyMxzz) * abs(mxyyMxzz) / (abs(mxyyMxzz) + qudricLimitM);
+ mxyyMxzz += wadjust * (-mxyyMxzz);
+ //////////////////////////////////////////////////////////////////////////
+ // no limiter
+ // mfbbb += OxyyMxzz * (-mfbbb);
+ // mxxyPyzz += OxyyPxzz * (-mxxyPyzz);
+ // mxxyMyzz += OxyyMxzz * (-mxxyMyzz);
+ // mxxzPyyz += OxyyPxzz * (-mxxzPyyz);
+ // mxxzMyyz += OxyyMxzz * (-mxxzMyyz);
+ // mxyyPxzz += OxyyPxzz * (-mxyyPxzz);
+ // mxyyMxzz += OxyyMxzz * (-mxyyMxzz);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute inverse linear combinations of second and third order cumulants
+ //!
+ mfcaa = c1o3 * (mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaca = c1o3 * (-c2 * mxxMyy + mxxMzz + mxxPyyPzz);
+ mfaac = c1o3 * (mxxMyy - c2 * mxxMzz + mxxPyyPzz);
+
+ mfcba = (mxxyMyzz + mxxyPyzz) * c1o2;
+ mfabc = (-mxxyMyzz + mxxyPyzz) * c1o2;
+ mfcab = (mxxzMyyz + mxxzPyyz) * c1o2;
+ mfacb = (-mxxzMyyz + mxxzPyyz) * c1o2;
+ mfbca = (mxyyMxzz + mxyyPxzz) * c1o2;
+ mfbac = (-mxyyMxzz + mxyyPxzz) * c1o2;
+ //////////////////////////////////////////////////////////////////////////
+
+ //////////////////////////////////////////////////////////////////////////
+ // 4.
+ // no limiter
+ //! - Relax fourth order cumulants to modified equilibrium for fourth order convergence of diffusion
+ //! according to Eq. (43)-(48) <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et
+ //! al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ CUMacc = -O4 * (c1 / omega - c1o2) * (dyuy + dzuz) * c2o3 * A + (c1 - O4) * (CUMacc);
+ CUMcac = -O4 * (c1 / omega - c1o2) * (dxux + dzuz) * c2o3 * A + (c1 - O4) * (CUMcac);
+ CUMcca = -O4 * (c1 / omega - c1o2) * (dyuy + dxux) * c2o3 * A + (c1 - O4) * (CUMcca);
+ CUMbbc = -O4 * (c1 / omega - c1o2) * Dxy * c1o3 * B + (c1 - O4) * (CUMbbc);
+ CUMbcb = -O4 * (c1 / omega - c1o2) * Dxz * c1o3 * B + (c1 - O4) * (CUMbcb);
+ CUMcbb = -O4 * (c1 / omega - c1o2) * Dyz * c1o3 * B + (c1 - O4) * (CUMcbb);
+
+ //////////////////////////////////////////////////////////////////////////
+ // 5.
+ CUMbcc += O5 * (-CUMbcc);
+ CUMcbc += O5 * (-CUMcbc);
+ CUMccb += O5 * (-CUMccb);
+
+ //////////////////////////////////////////////////////////////////////////
+ // 6.
+ CUMccc += O6 * (-CUMccc);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Compute central moments from post collision cumulants according to Eq. (53)-(56) in
+ //! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+
+ //////////////////////////////////////////////////////////////////////////
+ // 4.
+ mfcbb = CUMcbb + c1o3 * ((c3 * mfcaa + c1) * mfabb + c6 * mfbba * mfbab) * OOrho;
+ mfbcb = CUMbcb + c1o3 * ((c3 * mfaca + c1) * mfbab + c6 * mfbba * mfabb) * OOrho;
+ mfbbc = CUMbbc + c1o3 * ((c3 * mfaac + c1) * mfbba + c6 * mfbab * mfabb) * OOrho;
+
+ mfcca = CUMcca + (((mfcaa * mfaca + c2 * mfbba * mfbba) * c9 + c3 * (mfcaa + mfaca)) * OOrho -
+ (drho * OOrho)) *
+ c1o9;
+ mfcac = CUMcac + (((mfcaa * mfaac + c2 * mfbab * mfbab) * c9 + c3 * (mfcaa + mfaac)) * OOrho -
+ (drho * OOrho)) *
+ c1o9;
+ mfacc = CUMacc + (((mfaac * mfaca + c2 * mfabb * mfabb) * c9 + c3 * (mfaac + mfaca)) * OOrho -
+ (drho * OOrho)) *
+ c1o9;
+
+ //////////////////////////////////////////////////////////////////////////
+ // 5.
+ mfbcc = CUMbcc + c1o3 *
+ (c3 * (mfaac * mfbca + mfaca * mfbac + c4 * mfabb * mfbbb +
+ c2 * (mfbab * mfacb + mfbba * mfabc)) +
+ (mfbca + mfbac)) *
+ OOrho;
+ mfcbc = CUMcbc + c1o3 *
+ (c3 * (mfaac * mfcba + mfcaa * mfabc + c4 * mfbab * mfbbb +
+ c2 * (mfabb * mfcab + mfbba * mfbac)) +
+ (mfcba + mfabc)) *
+ OOrho;
+ mfccb = CUMccb + c1o3 *
+ (c3 * (mfcaa * mfacb + mfaca * mfcab + c4 * mfbba * mfbbb +
+ c2 * (mfbab * mfbca + mfabb * mfcba)) +
+ (mfacb + mfcab)) *
+ OOrho;
+
+ //////////////////////////////////////////////////////////////////////////
+ // 6.
+ mfccc =
+ CUMccc - ((-c4 * mfbbb * mfbbb - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca) -
+ c4 * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc) -
+ c2 * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb)) *
+ OOrho +
+ (c4 * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac) +
+ c2 * (mfcaa * mfaca * mfaac) + c16 * mfbba * mfbab * mfabb) *
+ OOrho * OOrho -
+ c1o3 * (mfacc + mfcac + mfcca) * OOrho - c1o9 * (mfcaa + mfaca + mfaac) * OOrho +
+ (c2 * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba) +
+ (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa) + c1o3 * (mfaac + mfaca + mfcaa)) *
+ OOrho * OOrho * c2o3 +
+ c1o27 * ((drho * drho - drho) * OOrho * OOrho));
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Add acceleration (body force) to first order cumulants according to Eq. (85)-(87) in
+ //! <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
+ //! DOI:10.1016/j.camwa.2015.05.001 ]</b></a>
+ //!
+ mfbaa = -mfbaa;
+ mfaba = -mfaba;
+ mfaab = -mfaab;
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Chimera transform from central moments to well conditioned distributions as defined in
+ //! Appendix J in <a href="https://doi.org/10.1016/j.camwa.2015.05.001"><b>[ M. Geier et al. (2015),
+ //! DOI:10.1016/j.camwa.2015.05.001 ]</b></a> see also Eq. (88)-(96) in <a
+ //! href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017),
+ //! DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+ //!
+ ////////////////////////////////////////////////////////////////////////////////////
+ // X - Dir
+ backwardInverseChimeraWithK(mfaaa, mfbaa, mfcaa, vvx, vx2, c1, c1);
+ backwardChimera(mfaba, mfbba, mfcba, vvx, vx2);
+ backwardInverseChimeraWithK(mfaca, mfbca, mfcca, vvx, vx2, c3, c1o3);
+ backwardChimera(mfaab, mfbab, mfcab, vvx, vx2);
+ backwardChimera(mfabb, mfbbb, mfcbb, vvx, vx2);
+ backwardChimera(mfacb, mfbcb, mfccb, vvx, vx2);
+ backwardInverseChimeraWithK(mfaac, mfbac, mfcac, vvx, vx2, c3, c1o3);
+ backwardChimera(mfabc, mfbbc, mfcbc, vvx, vx2);
+ backwardInverseChimeraWithK(mfacc, mfbcc, mfccc, vvx, vx2, c9, c1o9);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Y - Dir
+ backwardInverseChimeraWithK(mfaaa, mfaba, mfaca, vvy, vy2, c6, c1o6);
+ backwardChimera(mfaab, mfabb, mfacb, vvy, vy2);
+ backwardInverseChimeraWithK(mfaac, mfabc, mfacc, vvy, vy2, c18, c1o18);
+ backwardInverseChimeraWithK(mfbaa, mfbba, mfbca, vvy, vy2, c3o2, c2o3);
+ backwardChimera(mfbab, mfbbb, mfbcb, vvy, vy2);
+ backwardInverseChimeraWithK(mfbac, mfbbc, mfbcc, vvy, vy2, c9o2, c2o9);
+ backwardInverseChimeraWithK(mfcaa, mfcba, mfcca, vvy, vy2, c6, c1o6);
+ backwardChimera(mfcab, mfcbb, mfccb, vvy, vy2);
+ backwardInverseChimeraWithK(mfcac, mfcbc, mfccc, vvy, vy2, c18, c1o18);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ // Z - Dir
+ backwardInverseChimeraWithK(mfaaa, mfaab, mfaac, vvz, vz2, c36, c1o36);
+ backwardInverseChimeraWithK(mfaba, mfabb, mfabc, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfaca, mfacb, mfacc, vvz, vz2, c36, c1o36);
+ backwardInverseChimeraWithK(mfbaa, mfbab, mfbac, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfbba, mfbbb, mfbbc, vvz, vz2, c9o4, c4o9);
+ backwardInverseChimeraWithK(mfbca, mfbcb, mfbcc, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfcaa, mfcab, mfcac, vvz, vz2, c36, c1o36);
+ backwardInverseChimeraWithK(mfcba, mfcbb, mfcbc, vvz, vz2, c9, c1o9);
+ backwardInverseChimeraWithK(mfcca, mfccb, mfccc, vvz, vz2, c36, c1o36);
+ ////////////////////////////////////////////////////////////////////////////////////
+
+ //////////////////////////////////////////////////////////////////////////
+ // proof correctness
+ //////////////////////////////////////////////////////////////////////////
+#ifdef PROOF_CORRECTNESS
+ LBMReal drho_post = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca) +
+ (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) +
+ (mfbaa + mfbac + mfbca + mfbcc) + (mfabb + mfcbb) + (mfbab + mfbcb) +
+ (mfbba + mfbbc) + mfbbb;
+ LBMReal dif = drho - drho_post;
+#ifdef SINGLEPRECISION
+ if (dif > 10.0E-7 || dif < -10.0E-7)
+#else
+ if (dif > 10.0E-15 || dif < -10.0E-15)
+#endif
+ {
+ UB_THROW(UbException(
+ UB_EXARGS,
+ "rho=" + UbSystem::toString(drho) + ", rho_post=" + UbSystem::toString(drho_post) +
+ " dif=" + UbSystem::toString(dif) + " rho is not correct for node " +
+ UbSystem::toString(x1) + "," + UbSystem::toString(x2) + "," + UbSystem::toString(x3) +
+ " in " + block.lock()->toString() + " step = " + UbSystem::toString(step)));
+ }
+#endif
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - Write distributions: style of reading and writing the distributions from/to stored arrays
+ //! dependent on timestep is based on the esoteric twist algorithm <a
+ //! 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;
+ //////////////////////////////////////////////////////////////////////////
+ f[D3Q27System::REST] = mfbbb;
+
+ f[D3Q27System::E] = mfcbb;
+ f[D3Q27System::N] = mfbcb;
+ f[D3Q27System::T] = mfbbc;
+ f[D3Q27System::NE] = mfccb;
+ f[D3Q27System::NW] = mfacb;
+ f[D3Q27System::TE] = mfcbc;
+ f[D3Q27System::TW] = mfabc;
+ f[D3Q27System::TN] = mfbcc;
+ f[D3Q27System::TS] = mfbac;
+ f[D3Q27System::TNE] = mfccc;
+ f[D3Q27System::TNW] = mfacc;
+ f[D3Q27System::TSE] = mfcac;
+ f[D3Q27System::TSW] = mfaac;
+
+ f[D3Q27System::W] = mfabb;
+ f[D3Q27System::S] = mfbab;
+ f[D3Q27System::B] = mfbba;
+ f[D3Q27System::SW] = mfaab;
+ f[D3Q27System::SE] = mfcab;
+ f[D3Q27System::BW] = mfaba;
+ f[D3Q27System::BE] = mfcba;
+ f[D3Q27System::BS] = mfbaa;
+ f[D3Q27System::BN] = mfbca;
+ f[D3Q27System::BSW] = mfaaa;
+ f[D3Q27System::BSE] = mfcaa;
+ f[D3Q27System::BNW] = mfaca;
+ f[D3Q27System::BNE] = 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[D3Q27System::REST] - fEqSolid[D3Q27System::REST];
+ f[D3Q27System::REST] = fPre[D3Q27System::REST] + 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[D3Q27System::REST] - fEqSolid[D3Q27System::REST];
+ f[D3Q27System::REST] = fPre[D3Q27System::REST] + oneMinB * (f[D3Q27System::REST] - fPre[D3Q27System::REST]) + 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[D3Q27System::REST];
+
+ (*this->localDistributions)(D3Q27System::ET_E, x1, x2, x3) = f[D3Q27System::W] ;
+ (*this->localDistributions)(D3Q27System::ET_N, x1, x2, x3) = f[D3Q27System::S] ;
+ (*this->localDistributions)(D3Q27System::ET_T, x1, x2, x3) = f[D3Q27System::B] ;
+ (*this->localDistributions)(D3Q27System::ET_NE, x1, x2, x3) = f[D3Q27System::SW] ;
+ (*this->localDistributions)(D3Q27System::ET_NW, x1p, x2, x3) = f[D3Q27System::SE] ;
+ (*this->localDistributions)(D3Q27System::ET_TE, x1, x2, x3) = f[D3Q27System::BW] ;
+ (*this->localDistributions)(D3Q27System::ET_TW, x1p, x2, x3) = f[D3Q27System::BE] ;
+ (*this->localDistributions)(D3Q27System::ET_TN, x1, x2, x3) = f[D3Q27System::BS] ;
+ (*this->localDistributions)(D3Q27System::ET_TS, x1, x2p, x3) = f[D3Q27System::BN] ;
+ (*this->localDistributions)(D3Q27System::ET_TNE, x1, x2, x3) = f[D3Q27System::BSW] ;
+ (*this->localDistributions)(D3Q27System::ET_TNW, x1p, x2, x3) = f[D3Q27System::BSE] ;
+ (*this->localDistributions)(D3Q27System::ET_TSE, x1, x2p, x3) = f[D3Q27System::BNW] ;
+ (*this->localDistributions)(D3Q27System::ET_TSW, x1p, x2p, x3) = f[D3Q27System::BNE] ;
+
+ (*this->nonLocalDistributions)(D3Q27System::ET_W, x1p, x2, x3) = f[D3Q27System::E] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_S, x1, x2p, x3) = f[D3Q27System::N] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_B, x1, x2, x3p) = f[D3Q27System::T] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_SW, x1p, x2p, x3) = f[D3Q27System::NE] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_SE, x1, x2p, x3) = f[D3Q27System::NW] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BW, x1p, x2, x3p) = f[D3Q27System::TE] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BE, x1, x2, x3p) = f[D3Q27System::TW] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BS, x1, x2p, x3p) = f[D3Q27System::TN] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BN, x1, x2, x3p) = f[D3Q27System::TS] ;
+ (*this->nonLocalDistributions)(D3Q27System::ET_BSW, x1p, x2p, x3p) = f[D3Q27System::TNE];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BSE, x1, x2p, x3p) = f[D3Q27System::TNW];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BNW, x1p, x2, x3p) = f[D3Q27System::TSE];
+ (*this->nonLocalDistributions)(D3Q27System::ET_BNE, x1, x2, x3p) = f[D3Q27System::TSW];
+ }
+ }
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+
diff --git a/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.h b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.h
new file mode 100644
index 0000000000000000000000000000000000000000..503d1970979afa83cf6d1d65ccd43356a5123720
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.h
@@ -0,0 +1,155 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// 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 IBcumulantK17LBMKernel.h
+//! \ingroup LBM
+//! \author Konstantin Kutscher, Martin Geier
+//=======================================================================================
+
+#ifndef IBcumulantK17LBMKernel_h__
+#define IBcumulantK17LBMKernel_h__
+
+#include "LBMKernel.h"
+#include "BCProcessor.h"
+#include "D3Q27System.h"
+#include "UbTiming.h"
+#include "CbArray4D.h"
+#include "CbArray3D.h"
+#include "IBdynamicsParticleData.h"
+
+
+//! \brief Compressible cumulant LBM kernel.
+//! \details LBM implementation that use Cascaded Cumulant Lattice Boltzmann method for D3Q27 model
+//!
+//! The model is publisched in
+//! <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
+{
+public:
+ IBcumulantK17LBMKernel();
+ ~IBcumulantK17LBMKernel() = default;
+ 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);
+ inline void backwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K);
+ inline void forwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
+ inline void backwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2);
+
+ virtual void initDataSet();
+
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr localDistributions;
+ CbArray4D<LBMReal, IndexerX4X3X2X1>::CbArray4DPtr nonLocalDistributions;
+ CbArray3D<LBMReal, IndexerX3X2X1>::CbArray3DPtr restDistributions;
+
+ mu::value_type muX1, muX2, muX3;
+ mu::value_type muDeltaT;
+ mu::value_type muNu;
+ LBMReal forcingX1;
+ LBMReal forcingX2;
+ LBMReal forcingX3;
+
+ CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData;
+};
+
+////////////////////////////////////////////////////////////////////////////////
+//! \brief forward chimera transformation \ref forwardInverseChimeraWithK
+//! Transformation from distributions to central moments according to Eq. (6)-(14) in
+//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+//! Modified for lower round-off errors.
+////////////////////////////////////////////////////////////////////////////////
+inline void IBcumulantK17LBMKernel::forwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K)
+{
+ using namespace UbMath;
+ LBMReal m2 = mfa + mfc;
+ LBMReal m1 = mfc - mfa;
+ LBMReal m0 = m2 + mfb;
+ mfa = m0;
+ m0 *= Kinverse;
+ m0 += c1;
+ mfb = (m1 * Kinverse - m0 * vv) * K;
+ mfc = ((m2 - c2 * m1 * vv) * Kinverse + v2 * m0) * K;
+}
+////////////////////////////////////////////////////////////////////////////////
+//! \brief backward chimera transformation \ref backwardInverseChimeraWithK
+//! Transformation from central moments to distributions according to Eq. (57)-(65) in
+//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+//! ] Modified for lower round-off errors.
+////////////////////////////////////////////////////////////////////////////////
+inline void IBcumulantK17LBMKernel::backwardInverseChimeraWithK(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2, LBMReal Kinverse, LBMReal K)
+{
+ using namespace UbMath;
+ LBMReal m0 = (((mfc - mfb) * c1o2 + mfb * vv) * Kinverse + (mfa * Kinverse + c1) * (v2 - vv) * c1o2) * K;
+ LBMReal m1 = (((mfa - mfc) - c2 * mfb * vv) * Kinverse + (mfa * Kinverse + c1) * (-v2)) * K;
+ mfc = (((mfc + mfb) * c1o2 + mfb * vv) * Kinverse + (mfa * Kinverse + c1) * (v2 + vv) * c1o2) * K;
+ mfa = m0;
+ mfb = m1;
+}
+////////////////////////////////////////////////////////////////////////////////
+//! \brief forward chimera transformation \ref forwardChimera
+//! Transformation from distributions to central moments according to Eq. (6)-(14) in
+//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+//! for \f$ K_{abc}=0 \f$. This is to avoid unnessary floating point operations.
+//! Modified for lower round-off errors.
+////////////////////////////////////////////////////////////////////////////////
+inline void IBcumulantK17LBMKernel::forwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2)
+{
+ using namespace UbMath;
+ LBMReal m1 = (mfa + mfc) + mfb;
+ LBMReal m2 = mfc - mfa;
+ mfc = (mfc + mfa) + (v2 * m1 - c2 * vv * m2);
+ mfb = m2 - vv * m1;
+ mfa = m1;
+}
+////////////////////////////////////////////////////////////////////////////////
+//! \brief backward chimera transformation \ref backwardChimera
+//! Transformation from central moments to distributions according to Eq. (57)-(65) in
+//! <a href="https://doi.org/10.1016/j.jcp.2017.05.040"><b>[ M. Geier et al. (2017), DOI:10.1016/j.jcp.2017.05.040 ]</b></a>
+//! for \f$ K_{abc}=0 \f$. This is to avoid unnessary floating point operations.
+//! Modified for lower round-off errors.
+////////////////////////////////////////////////////////////////////////////////
+inline void IBcumulantK17LBMKernel::backwardChimera(LBMReal& mfa, LBMReal& mfb, LBMReal& mfc, LBMReal vv, LBMReal v2)
+{
+ using namespace UbMath;
+ LBMReal ma = (mfc + mfa * (v2 - vv)) * c1o2 + mfb * (vv - c1o2);
+ LBMReal mb = ((mfa - mfc) - mfa * v2) - c2 * mfb * vv;
+ mfc = (mfc + mfa * (v2 + vv)) * c1o2 + mfb * (vv + c1o2);
+ mfb = mb;
+ mfa = ma;
+}
+
+#endif // IBcumulantK17LBMKernel_h__
\ No newline at end of file
diff --git a/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h b/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
new file mode 100644
index 0000000000000000000000000000000000000000..28f11bf2eb84a509d59e6c098c0fbb69874f9517
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
@@ -0,0 +1,61 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// 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 DataSet3D.h
+//! \ingroup Data
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#ifndef IBdynamicsParticleData_h
+#define IBdynamicsParticleData_h
+
+#include<array>
+
+constexpr auto SOLFRAC_MIN = 0.001;
+constexpr auto SOLFRAC_MAX = 0.999;
+
+struct IBdynamicsParticleData {
+public:
+ IBdynamicsParticleData()
+ : partId(0), solidFraction(0.)
+ {
+ uPart[0] = 0.;
+ uPart[1] = 0.;
+ uPart[2] = 0.;
+
+ hydrodynamicForce[0] = 0.;
+ hydrodynamicForce[1] = 0.;
+ hydrodynamicForce[2] = 0.;
+ };
+ int partId;
+ double solidFraction;
+ std::array<double, 3> uPart;
+ std::array<double, 3> hydrodynamicForce;
+};
+
+#endif
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..da8c93296a333050c036f7fe084bbe23aa39ca23
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
@@ -0,0 +1,413 @@
+#include "LiggghtsCouplingCoProcessor.h"
+#include "GbSphere3D.h"
+#include "MPICommunicator.h"
+#include "CoProcessor.h"
+#include "LiggghtsCouplingWrapper.h"
+#include "Grid3D.h"
+#include "Block3D.h"
+#include "LBMKernel.h"
+#include "DistributionArray3D.h"
+#include "DataSet3D.h"
+#include "IBcumulantK17LBMKernel.h"
+#include "LBMUnitConverter.h"
+#include "fix_lb_coupling_onetoone.h"
+
+LiggghtsCouplingCoProcessor::LiggghtsCouplingCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
+ 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()
+{
+}
+
+void LiggghtsCouplingCoProcessor::process(double actualTimeStep)
+{
+ std::cout << "step: " << actualTimeStep << "\n";
+
+ getForcesFromLattice();
+
+ wrapper.run(demSteps);
+
+ setSpheresOnLattice();
+}
+
+void LiggghtsCouplingCoProcessor::setSpheresOnLattice()
+{
+ std::vector<int> excludeType;
+
+ int nPart = wrapper.lmp->atom->nlocal + wrapper.lmp->atom->nghost;
+
+ for (int iS = 0; iS < nPart; iS++)
+ {
+ int type = (int)wrapper.lmp->atom->type[iS];
+ bool excludeFlag(false);
+ for (int iT = 0; iT < excludeType.size(); iT++) {
+ //std::cout << iS << " " << type << " " << excludeType[iT] << std::endl;
+ if (type == excludeType[iT]) {
+ excludeFlag = true;
+ break;
+ }
+ }
+
+ if (excludeFlag)
+ continue;
+
+ double x[3], v[3], omega[3];
+ double r;
+ int id = wrapper.lmp->atom->tag[iS];
+
+ for (int i = 0; i < 3; i++)
+ {
+ x[i] = wrapper.lmp->atom->x[iS][i]; // * units->getFactorLentghWToLb(); // - 0.5; ????
+ v[i] = wrapper.lmp->atom->v[iS][i] * units->getFactorVelocityWToLb();
+ omega[i] = wrapper.lmp->atom->omega[iS][i] / units->getFactorTimeWToLb();
+ }
+
+ 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;
+
+ setSingleSphere3D(x, v, omega, r, id);
+ }
+}
+
+void LiggghtsCouplingCoProcessor::setSingleSphere3D(double *x, double *v, double *omega, /* double *com,*/ double r,
+ int id /*, bool initVelFlag*/)
+{
+ int level = 0;
+ //UbTupleInt3 bi = grid->getBlockIndexes(x[0], x[1], x[2], level);
+ //SPtr<Block3D> block = grid->getBlock(val<1>(bi), val<2>(bi), val<3>(bi), level);
+
+ std::vector<SPtr<Block3D>> blocks;
+ grid->getBlocksByCuboid(level, x[0] - r, x[1] - r, x[2] - r, x[0] + r, x[1] + r, x[2] + r, blocks);
+
+
+
+
+ for (SPtr<Block3D> block : blocks) {
+ if (block) {
+ SPtr<ILBMKernel> kernel = block->getKernel();
+ SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
+
+ CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData =
+ dynamicPointerCast<IBcumulantK17LBMKernel>(kernel)->getParticleData();
+
+ if (!particleData)
+ continue;
+
+ int minX1 = 1;
+ int minX2 = 1;
+ int minX3 = 1;
+
+ 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++) {
+ for (int ix1 = minX1; ix1 < maxX1; ix1++) {
+
+ //UbTupleInt3 blockNX = grid->getBlockNX();
+
+ //double const dx = val<1>(blockNX) * block->getX1() + ix1 - x[0];
+ //double const dy = val<2>(blockNX) * block->getX2() + ix2 - x[1];
+ //double const dz = val<3>(blockNX) * block->getX3() + ix3 - x[2];
+
+ Vector3D worldCoordinates = grid->getNodeCoordinates(block, ix1, ix2, ix3);
+
+ double const dx = (worldCoordinates[0] - x[0]) * units->getFactorLentghWToLb();
+ double const dy = (worldCoordinates[1] - x[1]) * units->getFactorLentghWToLb();
+ double const dz = (worldCoordinates[2] - x[2]) * units->getFactorLentghWToLb();
+
+ double const sf = calcSolidFraction(dx, dy, dz, r * units->getFactorLentghWToLb());
+
+ double const sf_old = (*particleData)(ix1,ix2,ix3)->solidFraction;
+ int const id_old = (int)(*particleData)(ix1,ix2,ix3)->partId;
+
+ int const decFlag = (sf > SOLFRAC_MIN) + 2 * (sf_old > SOLFRAC_MIN);
+
+ switch (decFlag) {
+ case 0: // sf == 0 && sf_old == 0
+ setToZero(*(*particleData)(ix1, ix2, ix3).get());
+ break; // do nothing
+ case 1: // sf > 0 && sf_old == 0
+ setValues(*(*particleData)(ix1, ix2, ix3).get(), id, sf, v, dx, dy, dz, omega);
+ break;
+ case 2: // sf == 0 && sf_old > 0
+ if (id_old == id) // then particle has left this cell
+ setToZero(*(*particleData)(ix1, ix2, ix3).get());
+ break; // else do nothing
+ case 3: // sf > 0 && sf_old > 0
+ if (sf > sf_old || id_old == id)
+ setValues(*(*particleData)(ix1, ix2, ix3).get(), id, sf, v, dx, dy, dz, omega);
+ break; // else do nothing
+ }
+ // if desired, initialize interior of sphere with sphere velocity
+ // if (initVelFlag && sf > SOLFRAC_MAX)
+ // cell.defineVelocity(particleData->uPart);
+
+ //if (sf > 0) {
+ // std::cout << "sf = " << sf << std::endl;
+ // std::cout << "ix1 = " << ix1 << ", ix2 = " << ix2 << ", ix3 = " << ix3 << std::endl;
+ //}
+ }
+ }
+ }
+ }
+ }
+
+}
+
+double LiggghtsCouplingCoProcessor::calcSolidFraction(double const dx_, double const dy_, double const dz_,
+ double const r_)
+{
+ static int const slicesPerDim = 5;
+ static double const sliceWidth = 1. / ((double)slicesPerDim);
+ static double const fraction = 1. / ((double)(slicesPerDim * slicesPerDim * slicesPerDim));
+
+ // should be sqrt(3.)/2.
+ // add a little to avoid roundoff errors
+ static const double sqrt3half = (double)sqrt(3.1) / 2.;
+
+ double const dist = dx_ * dx_ + dy_ * dy_ + dz_ * dz_;
+
+ double const r_p = r_ + sqrt3half;
+ if (dist > r_p * r_p)
+ return 0;
+
+ double const r_m = r_ - sqrt3half;
+ if (dist < r_m * r_m)
+ return 1;
+
+ double const r_sq = r_ * r_;
+ double dx_sq[slicesPerDim], dy_sq[slicesPerDim], dz_sq[slicesPerDim];
+
+ // pre-calculate d[xyz]_sq for efficiency
+ for (int i = 0; i < slicesPerDim; i++) {
+ double const delta = -0.5 + ((double)i + 0.5) * sliceWidth;
+ double const dx = dx_ + delta;
+ dx_sq[i] = dx * dx;
+ double const dy = dy_ + delta;
+ dy_sq[i] = dy * dy;
+ double const dz = dz_ + delta;
+ dz_sq[i] = dz * dz;
+ }
+
+ unsigned int n(0);
+ for (int i = 0; i < slicesPerDim; i++) {
+ for (int j = 0; j < slicesPerDim; j++) {
+ for (int k = 0; k < slicesPerDim; k++) {
+ n += (dx_sq[i] + dy_sq[j] + dz_sq[k] < r_sq);
+ }
+ }
+ }
+
+ return fraction * ((double)n);
+}
+
+ void LiggghtsCouplingCoProcessor::setValues(IBdynamicsParticleData &p, int const id, double const sf, double const *v, double const dx, double const dy, double const dz, double const *omega)
+{
+ p.uPart[0] = v[0];
+ p.uPart[1] = v[1];
+ p.uPart[2] = v[2];
+
+ if (omega != 0) {
+ p.uPart[0] += omega[1] * dz - omega[2] * dy;
+ p.uPart[1] += -omega[0] * dz + omega[2] * dx;
+ p.uPart[2] += omega[0] * dy - omega[1] * dx;
+ }
+ p.solidFraction = sf;
+ p.partId = id;
+}
+
+
+void LiggghtsCouplingCoProcessor::setToZero(IBdynamicsParticleData &p)
+{
+ p.uPart[0] = 0;
+ p.uPart[1] = 0;
+ p.uPart[2] = 0;
+ p.solidFraction = 0;
+ p.partId = 0;
+}
+
+void LiggghtsCouplingCoProcessor::getForcesFromLattice()
+{
+ static std::vector<double> force, torque;
+ static typename ParticleData::ParticleDataArrayVector x_lb;
+
+ int const nPart = wrapper.lmp->atom->nlocal + wrapper.lmp->atom->nghost;
+ int const n_force = nPart * 3;
+
+ if (nPart == 0)
+ return; // no particles - no work
+
+ if (nPart > x_lb.size()) {
+ for (int iPart = 0; iPart < x_lb.size(); iPart++) {
+ x_lb[iPart][0] = wrapper.lmp->atom->x[iPart][0];
+ x_lb[iPart][1] = wrapper.lmp->atom->x[iPart][1];
+ x_lb[iPart][2] = wrapper.lmp->atom->x[iPart][2];
+ }
+ for (int iPart = x_lb.size(); iPart < nPart; iPart++) {
+ std::array<double, 3> ar = {wrapper.lmp->atom->x[iPart][0],
+ wrapper.lmp->atom->x[iPart][1],
+ wrapper.lmp->atom->x[iPart][2]};
+ x_lb.push_back(ar);
+ }
+
+
+ } else {
+ for (int iPart = 0; iPart < nPart; iPart++) {
+ x_lb[iPart][0] = wrapper.lmp->atom->x[iPart][0];
+ x_lb[iPart][1] = wrapper.lmp->atom->x[iPart][1];
+ x_lb[iPart][2] = wrapper.lmp->atom->x[iPart][2];
+ }
+ }
+
+ if (n_force > force.size()) {
+ for (int i = 0; i < force.size(); i++) {
+ force[i] = 0;
+ torque[i] = 0;
+ }
+ for (int i = force.size(); i < n_force; i++) {
+ force.push_back(0.);
+ torque.push_back(0.);
+ }
+ } else {
+ for (int i = 0; i < n_force; i++) {
+ force[i] = 0;
+ torque[i] = 0;
+ }
+ }
+
+ SumForceTorque3D(x_lb, &force.front(), &torque.front());
+
+ LAMMPS_NS::FixLbCouplingOnetoone *couplingFix =
+ dynamic_cast<LAMMPS_NS::FixLbCouplingOnetoone *>(wrapper.lmp->modify->find_fix_style("couple/lb/onetoone", 0));
+
+ double **f_liggghts = couplingFix->get_force_ptr();
+ double **t_liggghts = couplingFix->get_torque_ptr();
+
+ for (int iPart = 0; iPart < nPart; iPart++)
+ for (int j = 0; j < 3; j++) {
+ f_liggghts[iPart][j] = 0;
+ t_liggghts[iPart][j] = 0;
+ }
+
+ for (int iPart = 0; iPart < nPart; iPart++) {
+ int tag = wrapper.lmp->atom->tag[iPart];
+ int liggghts_ind = wrapper.lmp->atom->map(tag);
+
+ for (int j = 0; j < 3; j++) {
+ f_liggghts[liggghts_ind][j] += force[3 * iPart + j] * units->getFactorForceLbToW();
+ t_liggghts[liggghts_ind][j] += torque[3 * iPart + j] * units->getFactorTorqueLbToW();
+ }
+ }
+ couplingFix->comm_force_torque();
+}
+
+void LiggghtsCouplingCoProcessor::SumForceTorque3D(ParticleData::ParticleDataArrayVector &x, double *force, double *torque)
+{
+ int nx = grid->getNX1(), ny = grid->getNX2(), nz = grid->getNX3();
+
+ std::vector < SPtr < Block3D > > blocks;
+ int level = 0;
+ grid->getBlocks(level, gridRank, true, blocks);
+
+
+ for (SPtr<Block3D> block : blocks) {
+ if (block) {
+ SPtr<ILBMKernel> kernel = block->getKernel();
+ SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
+
+ CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>::CbArray3DPtr particleData =
+ dynamicPointerCast<IBcumulantK17LBMKernel>(kernel)->getParticleData();
+
+ if (!particleData)
+ continue;
+
+ int minX1 = 1;
+ int minX2 = 1;
+ int minX3 = 1;
+
+ 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++) {
+ for (int ix1 = minX1; ix1 < maxX1; ix1++) {
+
+ // LIGGGHTS indices start at 1
+ int const id = (*particleData)(ix1, ix2, ix3)->partId;
+ if (id < 1)
+ continue; // no particle here
+
+ int const ind = wrapper.lmp->atom->map(id);
+
+ Vector3D worldCoordinates = grid->getNodeCoordinates(block, ix1, ix2, ix3);
+
+ double dx = (worldCoordinates[0] - x[ind][0]) * units->getFactorLentghWToLb();
+ double dy = (worldCoordinates[1] - x[ind][1]) * units->getFactorLentghWToLb();
+ double dz = (worldCoordinates[2] - x[ind][2]) * units->getFactorLentghWToLb();
+
+ // minimum image convention, needed if
+ // (1) PBC are used and
+ // (2) both ends of PBC lie on the same processor
+ if (dx > nx / 2)
+ dx -= nx;
+ else if (dx < -nx / 2)
+ dx += nx;
+ if (dy > ny / 2)
+ dy -= ny;
+ else if (dy < -ny / 2)
+ dy += ny;
+ if (dz > nz / 2)
+ dz -= nz;
+ else if (dz < -nz / 2)
+ dz += nz;
+
+ double const forceX = (*particleData)(ix1, ix2, ix3)->hydrodynamicForce[0];
+ double const forceY = (*particleData)(ix1, ix2, ix3)->hydrodynamicForce[1];
+ double const forceZ = (*particleData)(ix1, ix2, ix3)->hydrodynamicForce[2];
+
+ double const torqueX = dy * forceZ - dz * forceY;
+ double const torqueY = -dx * forceZ + dz * forceX;
+ double const torqueZ = dx * forceY - dy * forceX;
+
+ addForce(ind, 0, forceX, force);
+ addForce(ind, 1, forceY, force);
+ addForce(ind, 2, forceZ, force);
+
+ addTorque(ind, 0, torqueX, torque);
+ addTorque(ind, 1, torqueY, torque);
+ addTorque(ind, 2, torqueZ, torque);
+ }
+ }
+ }
+ }
+ }
+ }
+
+void LiggghtsCouplingCoProcessor::addForce(int const partId, int const coord, double const value, double *force)
+{
+ force[3 * partId + coord] += value;
+}
+
+void LiggghtsCouplingCoProcessor::addTorque(int const partId, int const coord, double const value, double *torque)
+{
+ torque[3 * partId + coord] += value;
+}
\ No newline at end of file
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
new file mode 100644
index 0000000000000000000000000000000000000000..dcaa6e16c0a46a69d64998285119e86077500b6b
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
@@ -0,0 +1,104 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// 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 LiggghtsCouplingCoProcessor.h
+//! \ingroup LiggghtsCoupling
+//! \author Konstantin Kutscher
+//=======================================================================================
+
+#ifndef LiggghtsCouplingCoProcessor_h
+#define LiggghtsCouplingCoProcessor_h
+
+#include "CoProcessor.h"
+
+#include "lammps.h"
+#include "input.h"
+#include "atom.h"
+#include "modify.h"
+
+#include <memory>
+#include <vector>
+
+
+class CoProcessor;
+class Communicator;
+class LiggghtsCouplingWrapper;
+class Grid3D;
+class Block3D;
+struct IBdynamicsParticleData;
+class LBMUnitConverter;
+
+struct ParticleData {
+ typedef typename std::vector<std::array<double, 3>> ParticleDataArrayVector;
+ typedef typename std::vector<double> ParticleDataScalarVector;
+};
+
+class LiggghtsCouplingCoProcessor : public CoProcessor
+{
+public:
+ LiggghtsCouplingCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<Communicator> comm,
+ LiggghtsCouplingWrapper &wrapper, int demSteps, SPtr<LBMUnitConverter> units);
+ virtual ~LiggghtsCouplingCoProcessor();
+
+ void process(double actualTimeStep) override;
+
+
+protected:
+ void setSpheresOnLattice();
+
+ void setSingleSphere3D(double *x, double *v, double *omega, double r, int id /*, bool initVelFlag*/);
+
+ double calcSolidFraction(double const dx_, double const dy_, double const dz_, double const r_);
+
+ void setValues(IBdynamicsParticleData &p, int const id, double const sf, double const *v, double const dx, double const dy, double const dz, double const *omega);
+
+ void setToZero(IBdynamicsParticleData &p);
+
+ void getForcesFromLattice();
+
+ void SumForceTorque3D(ParticleData::ParticleDataArrayVector &x, double *force, double *torque);
+
+ void addForce(int const partId, int const coord, double const value, double *force);
+
+ void addTorque(int const partId, int const coord, double const value, double *torque);
+
+private:
+ SPtr<Communicator> comm;
+ LiggghtsCouplingWrapper &wrapper;
+ SPtr<LBMUnitConverter> units;
+ int demSteps;
+ //std::vector<std::vector<SPtr<Block3D>>> blockVector;
+ //int minInitLevel;
+ //int maxInitLevel;
+ int gridRank;
+
+ double *force, *torque;
+};
+
+#endif
+
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingWrapper.cpp b/src/cpu/LiggghtsCoupling/LiggghtsCouplingWrapper.cpp
new file mode 100644
index 0000000000000000000000000000000000000000..9be87887a26d654d03dc8a32ed9e456ec352fef2
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingWrapper.cpp
@@ -0,0 +1,85 @@
+/*
+ * This file is part of the LBDEMcoupling software.
+ *
+ * LBDEMcoupling 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, version 3.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Copyright 2014 Johannes Kepler University Linz
+ *
+ * Author: Philippe Seil (philippe.seil@jku.at)
+ */
+
+#include "LiggghtsCouplingWrapper.h"
+#include "mpi.h"
+#include <iostream>
+#include <sstream>
+
+LiggghtsCouplingWrapper::LiggghtsCouplingWrapper(char **argv, MPI_Comm communicator)
+ : lmp(0)
+{
+ // todo: get LAMMPS to recognize command line options
+ int argc_lmp = 1;
+ char **argv_lmp = 0;
+ argv_lmp = new char*[1];
+ argv_lmp[0] = argv[0];
+
+ lmp = new LAMMPS_NS::LAMMPS(argc_lmp,argv_lmp,communicator);
+
+ // delete[] argv_lmp[0];
+ delete[] argv_lmp;
+}
+void LiggghtsCouplingWrapper::execFile(char* const fname)
+{
+ lmp->input->file(fname);
+}
+void LiggghtsCouplingWrapper::execCommand(std::stringstream const &cmd)
+{
+ lmp->input->one(cmd.str().c_str());
+}
+void LiggghtsCouplingWrapper::execCommand(char* const cmd)
+{
+ lmp->input->one(cmd);
+}
+int LiggghtsCouplingWrapper::getNumParticles()
+{
+ return lammps_get_natoms(lmp);
+}
+void LiggghtsCouplingWrapper::setVariable(char const *name, double value)
+{
+ std::stringstream cmd;
+ cmd << "variable " << name << " equal " << value;
+ std::cout << cmd.str() << std::endl;
+ execCommand(cmd);
+}
+void LiggghtsCouplingWrapper::setVariable(char const *name, std::string &value)
+{
+ std::stringstream cmd;
+ cmd << "variable " << name << " string " << value;
+ std::cout << cmd.str() << std::endl;
+ execCommand(cmd);
+}
+void LiggghtsCouplingWrapper::run(int nSteps)
+{
+ std::stringstream cmd;
+ cmd << "run " << nSteps;
+ execCommand(cmd);
+}
+void LiggghtsCouplingWrapper::runUpto(int nSteps)
+{
+ std::stringstream cmd;
+ cmd << "run " << nSteps << " upto";
+ execCommand(cmd);
+}
+
+
+
+
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingWrapper.h b/src/cpu/LiggghtsCoupling/LiggghtsCouplingWrapper.h
new file mode 100644
index 0000000000000000000000000000000000000000..a745a7e967ee7852a11cd37452012951e421c347
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingWrapper.h
@@ -0,0 +1,47 @@
+/*
+ * This file is part of the LBDEMcoupling software.
+ *
+ * LBDEMcoupling 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, version 3.
+ *
+ * This program 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 this program. If not, see <http://www.gnu.org/licenses/>.
+ *
+ * Copyright 2014 Johannes Kepler University Linz
+ *
+ * Author: Philippe Seil (philippe.seil@jku.at)
+ */
+
+#ifndef LIGGGHTSCOUPLINGWRAPPER_H
+#define LIGGGHTSCOUPLINGWRAPPER_H
+
+// necessary LAMMPS/LIGGGHTS includes
+
+#include "lammps.h"
+#include "input.h"
+#include "library.h"
+#include "library_cfd_coupling.h"
+
+class LiggghtsCouplingWrapper {
+public:
+ LiggghtsCouplingWrapper(char **argv, MPI_Comm communicator);
+ void execFile(char* const fname);
+ void execCommand(std::stringstream const &cmd);
+ void execCommand(char* const cmd);
+ void run(int nSteps);
+ void runUpto(int nSteps);
+ int getNumParticles();
+ void setVariable(char const *name, double value);
+ void setVariable(char const *name, std::string &value);
+
+ //private:
+ LAMMPS_NS::LAMMPS *lmp;
+};
+
+#endif /* LIGGGHTSCOUPLINGWRAPPER_H */
diff --git a/src/cpu/VirtualFluidsCore/CMakeLists.txt b/src/cpu/VirtualFluidsCore/CMakeLists.txt
index 15cdceffd99515f84d60c4b6169e2da7e74ecfc3..5abe76be9ed986d773785a544eb1e29b9954944d 100644
--- a/src/cpu/VirtualFluidsCore/CMakeLists.txt
+++ b/src/cpu/VirtualFluidsCore/CMakeLists.txt
@@ -16,7 +16,6 @@ IF(${USE_CATALYST})
list(APPEND VF_LIBRARIES optimized vtkParallelMPI debug vtkParallelMPI )
ENDIF()
-
IF(${USE_DEM_COUPLING})
INCLUDE(${CMAKE_CURRENT_SOURCE_DIR}/../DemCoupling/DemCoupling.cmake)
ENDIF()
@@ -27,6 +26,11 @@ endif()
vf_add_library(BUILDTYPE static PUBLIC_LINK basics muparser ${VF_LIBRARIES} PRIVATE_LINK lbm mpi logger)
+IF(${USE_LIGGGHTS})
+ list(APPEND VF_LIBRARIES optimized ${LIGGGHTS_RELEASE_LIBRARY} debug ${LIGGGHTS_DEBUG_LIBRARY})
+ENDIF()
+
+vf_add_library(BUILDTYPE static PUBLIC_LINK basics muparser MPI::MPI_CXX ${VF_LIBRARIES})
vf_get_library_name(library_name)
@@ -47,6 +51,12 @@ IF(${USE_METIS} AND METIS_INCLUDEDIR)
ENDIF()
target_include_directories(${library_name} PRIVATE ${ZOLTAN_INCLUDEDIR})
+
IF(${USE_VTK})
target_include_directories(${library_name} PRIVATE ${VTK_INCLUDE_DIRS})
ENDIF()
+
+IF(${USE_LIGGGHTS})
+ target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/../LiggghtsCoupling)
+ target_include_directories(${library_name} PUBLIC ${LIGGGHTS_SOURCE_DIR})
+ENDIF()
\ No newline at end of file
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/InSituVTKCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/InSituVTKCoProcessor.cpp
index 73c0a2325953994c337934347e872223ba18452a..ed258864a4a87b473ca276064abf60ad5910828d 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/InSituVTKCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/InSituVTKCoProcessor.cpp
@@ -226,11 +226,11 @@ void InSituVTKCoProcessor::addData(SPtr<Block3D> block)
UbSystem::toString(ix2) + "," + UbSystem::toString(ix3)));
// vx3=999.0;
- arrays[0]->InsertNextValue(rho * conv->getFactorDensityLbToW2());
- arrays[1]->InsertNextValue(vx1 * conv->getFactorVelocityLbToW2());
- arrays[2]->InsertNextValue(vx2 * conv->getFactorVelocityLbToW2());
- arrays[3]->InsertNextValue(vx3 * conv->getFactorVelocityLbToW2());
- arrays[4]->InsertNextValue(press * conv->getFactorPressureLbToW2());
+ arrays[0]->InsertNextValue(rho * conv->getFactorDensityLbToW());
+ arrays[1]->InsertNextValue(vx1 * conv->getFactorVelocityLbToW());
+ arrays[2]->InsertNextValue(vx2 * conv->getFactorVelocityLbToW());
+ arrays[3]->InsertNextValue(vx3 * conv->getFactorVelocityLbToW());
+ arrays[4]->InsertNextValue(press * conv->getFactorPressureLbToW());
}
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h b/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h
index 40570cc3847f71a1942791afa7e95145daafb53b..7c4ef7372c3126378c8fbe258b5914cafae5500b 100644
--- a/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h
+++ b/src/cpu/VirtualFluidsCore/LBM/LBMUnitConverter.h
@@ -97,14 +97,6 @@ public:
this->init(refLengthWorld, csWorld, rhoWorld, csWorld, refLengthLb, rhoLb, csLb);
}
- LBMUnitConverter(int /*dummy*/, double uReal, double uLB, double nuReal, double nuLB)
- {
- factorVelocityLbToW = uReal / uLB;
- factorViscosityLbToW = nuReal / nuLB;
- factorDensityLbToW = factorViscosityLbToW * factorVelocityLbToW * factorVelocityLbToW;
- factorPressureLbToW = factorDensityLbToW;
- }
-
virtual ~LBMUnitConverter() = default;
double getRefRhoLb() { return refRhoLb; }
@@ -124,10 +116,7 @@ public:
double getFactorDensityLbToW() { return this->factorMassLbToW / std::pow(factorLengthLbToW, 3.0); }
double getFactorDensityWToLb() { return 1.0 / this->getFactorDensityLbToW(); }
- double getFactorPressureLbToW()
- {
- return this->factorMassLbToW / (std::pow(factorTimeLbToW, 2.0) * factorLengthLbToW);
- }
+ double getFactorPressureLbToW(){ return this->factorMassLbToW / (factorLengthLbToW * factorTimeLbToW * factorTimeLbToW); }
double getFactorPressureWToLb() { return 1.0 / this->getFactorPressureLbToW(); }
double getFactorMassLbToW() { return this->factorMassLbToW; }
@@ -136,14 +125,14 @@ public:
double getFactorForceLbToW() { return factorMassLbToW * factorLengthLbToW / (factorTimeLbToW * factorTimeLbToW); }
double getFactorForceWToLb() { return 1.0 / this->getFactorForceLbToW(); }
+ double getFactorTorqueLbToW() { return factorMassLbToW * factorLengthLbToW * factorLengthLbToW / (factorTimeLbToW * factorTimeLbToW);}
+ double getFactorTorqueWToLb() { return 1.0 / this->getFactorTorqueWToLb(); }
+
double getFactorAccLbToW() { return factorLengthLbToW / (factorTimeLbToW * factorTimeLbToW); }
double getFactorAccWToLb() { return 1.0 / this->getFactorAccLbToW(); }
double getFactorTimeLbToW(double deltaX) const { return factorTimeWithoutDx * deltaX; }
- //////////////////////////////////////////////////////////////////////////
- double getFactorVelocityLbToW2() { return factorVelocityLbToW; }
- double getFactorDensityLbToW2() { return factorDensityLbToW; }
- double getFactorPressureLbToW2() { return factorPressureLbToW; }
+
/*==========================================================*/
friend inline std::ostream &operator<<(std::ostream &os, LBMUnitConverter c)
@@ -212,11 +201,6 @@ protected:
double factorMassLbToW{ 1.0 };
double refRhoLb{ 1.0 };
double factorTimeWithoutDx{ 0.0 };
-
- double factorVelocityLbToW{ 1.0 };
- double factorViscosityLbToW{ 1.0 };
- double factorDensityLbToW{ 1.0 };
- double factorPressureLbToW{ 1.0 };
};
#endif // LBMUNITCONVERTER_H