From 75cd008bc6691a2e8344118070a97ea4f06c3765 Mon Sep 17 00:00:00 2001
From: Konstantin Kutscher <kutscher@irmb.tu-bs.de>
Date: Mon, 18 Jul 2022 23:21:18 +0200
Subject: [PATCH] add Immersed Boundary: it doesn't work yet
---
apps/cpu/LiggghtsApp/LiggghtsApp.cpp | 101 ++-
apps/cpu/LiggghtsApp/in.lbdem | 72 ++
apps/cpu/LiggghtsApp/in2.lbdem | 21 +
.../LiggghtsCoupling/IBdynamicsParticleData.h | 50 ++
.../LiggghtsCouplingCoProcessor.cpp | 219 +++++-
.../LiggghtsCouplingCoProcessor.h | 36 +-
.../LBM/IBcumulantK17LBMKernel.cpp | 660 ++++++++++++++++++
.../LBM/IBcumulantK17LBMKernel.h | 156 +++++
8 files changed, 1311 insertions(+), 4 deletions(-)
create mode 100644 apps/cpu/LiggghtsApp/in.lbdem
create mode 100644 apps/cpu/LiggghtsApp/in2.lbdem
create mode 100644 src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
create mode 100644 src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp
create mode 100644 src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h
diff --git a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
index 61f835128..634ca28dd 100644
--- a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
+++ b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
@@ -21,10 +21,109 @@ using namespace std;
int main(int argc, char *argv[])
{
SPtr<Communicator> comm = MPICommunicator::getInstance();
- //MPI_Init(&argc, &argv);
+ 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] = { 10, 10, 20 };
+
+ double dx = 0.1;
+
+ double nuLB = 0.005;
+
+ SPtr<Grid3D> grid = make_shared<Grid3D>(comm);
+ grid->setPeriodicX1(true);
+ grid->setPeriodicX2(true);
+ grid->setPeriodicX3(true);
+ grid->setDeltaX(dx);
+ grid->setBlockNX(blockNX[0], blockNX[1], blockNX[2]);
+
+ string outputPath = "d:/temp/lll2";
+
+ 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();
+
+ SPtr<LBMKernel> kernel = make_shared<IncompressibleCumulantLBMKernel>();
+ SPtr<BCProcessor> bcProc = make_shared<BCProcessor>();
+ kernel->setBCProcessor(bcProc);
+
+ SetKernelBlockVisitor kernelVisitor(kernel, nuLB, 1e9, 1e9);
+ grid->accept(kernelVisitor);
+
+ 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.1;
+ double v_frac = 0.1;
+ double dt_phys = 1; // units.getPhysTime(1);
+ int demSubsteps = 1;
+ double dt_dem = 1e-1; //dt_phys / (double)demSubsteps;
+ int vtkSteps = 1;
+ string demOutDir = "d:/temp/lll2/";
+
+ 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);
+
+ // write data for visualization of macroscopic quantities
+ SPtr<UbScheduler> visSch(new UbScheduler(vtkSteps));
+ SPtr<WriteMacroscopicQuantitiesCoProcessor> writeMQCoProcessor(
+ new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlASCII::getInstance(),
+ SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
+
+ int endTime = 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"};
diff --git a/apps/cpu/LiggghtsApp/in.lbdem b/apps/cpu/LiggghtsApp/in.lbdem
new file mode 100644
index 000000000..f15eab6f8
--- /dev/null
+++ b/apps/cpu/LiggghtsApp/in.lbdem
@@ -0,0 +1,72 @@
+
+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.3 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
+
+
+
+run 1
diff --git a/apps/cpu/LiggghtsApp/in2.lbdem b/apps/cpu/LiggghtsApp/in2.lbdem
new file mode 100644
index 000000000..3e2fa90d2
--- /dev/null
+++ b/apps/cpu/LiggghtsApp/in2.lbdem
@@ -0,0 +1,21 @@
+
+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 1 ${dmp_dir}post/atom_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
diff --git a/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h b/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
new file mode 100644
index 000000000..8bf11ca7a
--- /dev/null
+++ b/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
@@ -0,0 +1,50 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// 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:
+ int partId{0};
+ double solidFraction{0};
+ std::array<double, 3> uPart{ 0., 0., 0. };
+ std::array<double, 3> hydrodynamicForce{ 0., 0., 0. };
+};
+
+#endif
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
index 6dc8539c6..ae3335374 100644
--- a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
@@ -1,7 +1,29 @@
#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"
-LiggghtsCouplingCoProcessor::LiggghtsCouplingCoProcessor()
+LiggghtsCouplingCoProcessor::LiggghtsCouplingCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s,
+ SPtr<Communicator> comm, LiggghtsCouplingWrapper &wrapper,
+ int demSteps)
+ : CoProcessor(grid, s), comm(comm), wrapper(wrapper), demSteps(demSteps)
{
+ //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()
@@ -9,5 +31,200 @@ LiggghtsCouplingCoProcessor::~LiggghtsCouplingCoProcessor()
}
void LiggghtsCouplingCoProcessor::process(double actualTimeStep)
+{
+ setSpheresOnLattice();
+ wrapper.run(demSteps);
+ std::cout << "step: " << actualTimeStep << "\n";
+}
+
+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];
+ v[i] = wrapper.lmp->atom->v[iS][i];
+ omega[i] = wrapper.lmp->atom->omega[iS][i];
+ }
+
+ r = wrapper.lmp->atom->radius[iS];
+
+ 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::getForcesFromLattice() {}
+
+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<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
+ 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++) {
+
+ Vector3D nX = grid->getNodeCoordinates(block, ix1, ix2, ix3);
+
+ double const dx = nX[0] - x[0];
+ double const dy = nX[1] - x[1];
+ double const dz = nX[2] - x[2];
+
+ double const sf = calcSolidFraction(dx, dy, dz, r);
+
+ 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(), sf, dx, dy, dz, omega, id);
+ 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(), sf, dx, dy, dz, omega, id);
+ break; // else do nothing
+ }
+ // if desired, initialize interior of sphere with sphere velocity
+ // if (initVelFlag && sf > SOLFRAC_MAX)
+ // cell.defineVelocity(particleData->uPart);
+ }
+ }
+ }
+ }
+ }
+
+}
+
+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, double const sf, double const dx,
+ double const dy, double const dz, double* omega, int id)
+{
+ //p.uPart.from_cArray(v);
+ 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;
}
\ No newline at end of file
diff --git a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
index d186aa750..509e83b4a 100644
--- a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
@@ -40,16 +40,48 @@
#include "input.h"
#include "atom.h"
#include "modify.h"
-#include "fix_lb_coupling_onetoone.h"
+
+#include <memory>
+#include <vector>
+
+
+class CoProcessor;
+class Communicator;
+class LiggghtsCouplingWrapper;
+class Grid3D;
+class Block3D;
+struct IBdynamicsParticleData;
class LiggghtsCouplingCoProcessor : public CoProcessor
{
public:
- LiggghtsCouplingCoProcessor();
+ LiggghtsCouplingCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, SPtr<Communicator> comm,
+ LiggghtsCouplingWrapper &wrapper, int demSteps);
virtual ~LiggghtsCouplingCoProcessor();
void process(double actualTimeStep) override;
+
+protected:
+ void setSpheresOnLattice();
+ void getForcesFromLattice();
+ void setSingleSphere3D(double *x, double *v, double *omega, /* double *com,*/ double r,
+ int id /*, bool initVelFlag*/);
+ double calcSolidFraction(double const dx_, double const dy_, double const dz_, double const r_);
+
+ void setValues(IBdynamicsParticleData &p, double const sf, double const dx, double const dy, double const dz,
+ double *omega, int id);
+
+ void setToZero(IBdynamicsParticleData &p);
+
+private:
+ SPtr<Communicator> comm;
+ LiggghtsCouplingWrapper &wrapper;
+ int demSteps;
+ std::vector<std::vector<SPtr<Block3D>>> blockVector;
+ int minInitLevel;
+ int maxInitLevel;
+ int gridRank;
};
#endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp
new file mode 100644
index 000000000..346a61009
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp
@@ -0,0 +1,660 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// 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], nx[1], nx[2]);
+
+ int minX1 = ghostLayerWidth;
+ int minX2 = ghostLayerWidth;
+ int minX3 = ghostLayerWidth;
+ int maxX1 = nx[0];
+ int maxX2 = nx[1];
+ int maxX3 = nx[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);
+
+ ////////////////////////////////////////////////////////////////////////////////////
+ //! - 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;
+ //////////////////////////////////////////////////////////////////////////
+
+ if ((*particleData)(ix1, ix2, ix3)->solidFraction < SOLFRAC_MIN)
+ return;
+
+
+ }
+ }
+ }
+ }
+}
+//////////////////////////////////////////////////////////////////////////
+
diff --git a/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h
new file mode 100644
index 000000000..1d3fb0f56
--- /dev/null
+++ b/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h
@@ -0,0 +1,156 @@
+//=======================================================================================
+// ____ ____ __ ______ __________ __ __ __ __
+// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
+// \ \ | | | | | |_) | | | | | | | / \ | |
+// \ \ | | | | | _ / | | | | | | / /\ \ | |
+// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
+// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
+// \ \ | | ________________________________________________________________
+// \ \ | | | ______________________________________________________________|
+// \ \| | | | __ __ __ __ ______ _______
+// \ | | |_____ | | | | | | | | | _ \ / _____)
+// \ | | _____| | | | | | | | | | | \ \ \_______
+// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
+// \ _____| |__| |________| \_______/ |__| |______/ (_______/
+//
+// 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();
+ LBMReal f[D3Q27System::ENDF + 1];
+
+ 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
--
GitLab