diff --git a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
index 634ca28dd4a21a39ae56d57d31e225ba2aba6399..fd892b6b9288b496fc11007ea034a145484d1c5e 100644
--- a/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
+++ b/apps/cpu/LiggghtsApp/LiggghtsApp.cpp
@@ -11,9 +11,8 @@
//#include "fix_lb_coupling_onetoone.h"
#include "LiggghtsCouplingCoProcessor.h"
-
#include "LiggghtsCouplingWrapper.h"
-
+#include "IBcumulantK17LBMKernel.h"
using namespace std;
@@ -31,9 +30,9 @@ int main(int argc, char *argv[])
double g_maxX1 = 1;
double g_maxX2 = 1;
- double g_maxX3 = 2;
+ double g_maxX3 = 1;
- int blockNX[3] = { 10, 10, 20 };
+ int blockNX[3] = { 10, 10, 10 };
double dx = 0.1;
@@ -61,7 +60,7 @@ int main(int argc, char *argv[])
ppblocks->process(0);
ppblocks.reset();
- SPtr<LBMKernel> kernel = make_shared<IncompressibleCumulantLBMKernel>();
+ SPtr<LBMKernel> kernel = make_shared<IBcumulantK17LBMKernel>();
SPtr<BCProcessor> bcProc = make_shared<BCProcessor>();
kernel->setBCProcessor(bcProc);
@@ -78,11 +77,17 @@ int main(int argc, char *argv[])
MPI_Comm mpi_comm = *(MPI_Comm*)(comm->getNativeCommunicator());
LiggghtsCouplingWrapper wrapper(argv, mpi_comm);
- double d_part = 0.1;
+
+ double d_part = 0.3;
+ double r_p = d_part / 2.0;
+
+ // SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(r_p, 1.480, 2060, r_p/dx);
+ SPtr<LBMUnitConverter> units = std::make_shared<LBMUnitConverter>(r_p, LBMUnitConverter::WATER, r_p / dx);
+
double v_frac = 0.1;
- double dt_phys = 1; // units.getPhysTime(1);
- int demSubsteps = 1;
- double dt_dem = 1e-1; //dt_phys / (double)demSubsteps;
+ double dt_phys = units->getFactorTimeLbToW();
+ int demSubsteps = 10;
+ double dt_dem = dt_phys / (double)demSubsteps;
int vtkSteps = 1;
string demOutDir = "d:/temp/lll2/";
@@ -99,8 +104,9 @@ int main(int argc, char *argv[])
wrapper.execFile((char *)inFile2.c_str());
//wrapper.runUpto(demSubsteps - 1);
+
SPtr<LiggghtsCouplingCoProcessor> lcCoProcessor =
- make_shared<LiggghtsCouplingCoProcessor>(grid, lScheduler, comm, wrapper, demSubsteps);
+ make_shared<LiggghtsCouplingCoProcessor>(grid, lScheduler, comm, wrapper, demSubsteps, units);
// write data for visualization of macroscopic quantities
SPtr<UbScheduler> visSch(new UbScheduler(vtkSteps));
@@ -108,7 +114,7 @@ int main(int argc, char *argv[])
new WriteMacroscopicQuantitiesCoProcessor(grid, visSch, outputPath, WbWriterVtkXmlASCII::getInstance(),
SPtr<LBMUnitConverter>(new LBMUnitConverter()), comm));
- int endTime = 20;
+ int endTime = 1000; //20;
SPtr<Calculator> calculator(new BasicCalculator(grid, lScheduler, endTime));
calculator->addCoProcessor(lcCoProcessor);
calculator->addCoProcessor(writeMQCoProcessor);
diff --git a/apps/cpu/LiggghtsApp/in.lbdem b/apps/cpu/LiggghtsApp/in.lbdem
index f15eab6f821c7b20e88d77077a926ec7317a276e..366f7afd134f76492d3d9756bf6743c5bde0c212 100644
--- a/apps/cpu/LiggghtsApp/in.lbdem
+++ b/apps/cpu/LiggghtsApp/in.lbdem
@@ -9,7 +9,7 @@ boundary f f f
newton off
processors * * 1
-region box block 0. 1. 0. 1. 0. 2. units box
+region box block 0. 1. 0. 1. 0. 1. units box
create_box 1 box
variable skin equal 0.01
@@ -41,7 +41,7 @@ fix ywalls2 all wall/gran model hertz tangential history primitive type 1 yplane
fix zwalls1 all wall/gran model hertz tangential history primitive type 1 zplane 0.
fix zwalls2 all wall/gran model hertz tangential history primitive type 1 zplane 2.
-create_atoms 1 single 0.5 0.5 1.75
+create_atoms 1 single 0.5 0.5 0.75
#create_atoms 1 single 0.38 0.05 0.05
set group all diameter 0.3 density 2400
diff --git a/apps/cpu/LiggghtsApp/in2.lbdem b/apps/cpu/LiggghtsApp/in2.lbdem
index 3e2fa90d2a25e3d4e01e73f32cb2f6f314c9f1bf..7098890cc28ea5f22ef442a388b1f79b4708107c 100644
--- a/apps/cpu/LiggghtsApp/in2.lbdem
+++ b/apps/cpu/LiggghtsApp/in2.lbdem
@@ -18,4 +18,4 @@ variable dmp_fname string ${dmp_dir}d_*.liggghts
# dump dmp all custom ${dmp_stp} ${dmp_dir}d_*.liggghts &
# id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
-dump dmp all custom/vtk 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
+dump dmp all custom/vtk 10 ${dmp_dir}post/atom_*.vtk id type type x y z ix iy iz vx vy vz fx fy fz omegax omegay omegaz radius
diff --git a/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.cpp
similarity index 55%
rename from src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp
rename to src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.cpp
index 346a61009b789b2e6bbc043d1e425da49eb1829f..ac0bc35844cb32544d901cbf5a93503d49de4070 100644
--- a/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp
+++ b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.cpp
@@ -55,14 +55,14 @@ void IBcumulantK17LBMKernel::initDataSet()
dataSet->setFdistributions(d);
particleData =
- std::make_shared<CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>>(nx[0], nx[1], nx[2]);
+ std::make_shared<CbArray3D<SPtr<IBdynamicsParticleData>, IndexerX3X2X1>>(nx[0] + 2, nx[1] + 2, nx[2] + 2);
- int minX1 = ghostLayerWidth;
- int minX2 = ghostLayerWidth;
- int minX3 = ghostLayerWidth;
- int maxX1 = nx[0];
- int maxX2 = nx[1];
- int maxX3 = nx[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;
@@ -150,37 +150,40 @@ void IBcumulantK17LBMKernel::calculate(int step)
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))
- {
+ 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>
+ //! - 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
+ // E N T
+ // c c c
//////////
- //W S B
- //a a a
+ // W S B
+ // a a a
- //Rest is b
+ // Rest is b
- //mfxyz
- //a - negative
- //b - null
- //c - positive
+ // mfxyz
+ // a - negative
+ // b - null
+ // c - positive
// a b c
//-1 0 1
@@ -215,31 +218,107 @@ void IBcumulantK17LBMKernel::calculate(int step)
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>
+ //! - 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;
+ (((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) +
+ ((mfacb + mfcab) + (mfaab + mfccb))) +
+ ((mfabb + mfcbb) + (mfbab + mfbcb)) + (mfbba + mfbbc)) +
+ mfbbb;
- LBMReal rho = c1 + drho;
+ 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;
+ (mfcbb - mfabb)) /
+ rho;
LBMReal vvy = ((((mfccc - mfaaa) + (mfaca - mfcac)) + ((mfacc - mfcaa) + (mfcca - mfaac))) +
(((mfbca - mfbac) + (mfbcc - mfbaa)) + ((mfacb - mfcab) + (mfccb - mfaab))) +
- (mfbcb - mfbab)) / rho;
+ (mfbcb - mfbab)) /
+ rho;
LBMReal vvz = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfacc - mfcaa) + (mfaac - mfcca))) +
(((mfbac - mfbca) + (mfbcc - mfbaa)) + ((mfabc - mfcba) + (mfcbc - mfaba))) +
- (mfbbc - mfbba)) / rho;
+ (mfbbc - mfbba)) /
+ rho;
////////////////////////////////////////////////////////////////////////////////////
- //forcing
+ // forcing
///////////////////////////////////////////////////////////////////////////////////////////
- if (withForcing)
- {
+ 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);
@@ -250,7 +329,8 @@ void IBcumulantK17LBMKernel::calculate(int step)
////////////////////////////////////////////////////////////////////////////////////
//! - 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>
+ //! <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
@@ -262,18 +342,20 @@ void IBcumulantK17LBMKernel::calculate(int step)
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>
+ //! - 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>
+ //! - 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
@@ -312,86 +394,116 @@ void IBcumulantK17LBMKernel::calculate(int step)
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>
+ //! - 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$.
+ //! - 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.
+ // 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);
+ // 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.
+ // 4.
LBMReal O4 = c1;
////////////////////////////////////////////////////////////
- //5.
+ // 5.
LBMReal O5 = c1;
////////////////////////////////////////////////////////////
- //6.
+ // 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).
+ //! - 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);
+ 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>
+ //! <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.
+ // 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));
+ 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;
+ // 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));
+ // 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.
+ // 2.
LBMReal mxxPyyPzz = mfcaa + mfaca + mfaac;
- LBMReal mxxMyy = mfcaa - mfaca;
- LBMReal mxxMzz = mfcaa - mfaac;
+ LBMReal mxxMyy = mfcaa - mfaca;
+ LBMReal mxxMzz = mfcaa - mfaac;
////////////////////////////////////////////////////////////
- //3.
+ // 3.
LBMReal mxxyPyzz = mfcba + mfabc;
LBMReal mxxyMyzz = mfcba - mfabc;
@@ -402,44 +514,48 @@ void IBcumulantK17LBMKernel::calculate(int step)
LBMReal mxyyMxzz = mfbca - mfbac;
////////////////////////////////////////////////////////////////////////////////////
- //incl. correction
+ // 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.
+ //! <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 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>
+ //! <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);
+ 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);
+ // mxxPyyPzz += OxxPyyPzz*(mfaaa - mxxPyyPzz);
+ // mxxMyy += -(-omega) * (-mxxMyy);
+ // mxxMzz += -(-omega) * (-mxxMzz);
/////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
mfabb += omega * (-mfabb);
mfbab += omega * (-mfbab);
mfbba += omega * (-mfbba);
////////////////////////////////////////////////////////////////////////////////////
- //relax
+ // 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>
+ //! <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);
@@ -457,13 +573,13 @@ void IBcumulantK17LBMKernel::calculate(int step)
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);
+ // 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
@@ -481,10 +597,11 @@ void IBcumulantK17LBMKernel::calculate(int step)
//////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
- //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>
+ // 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);
@@ -494,55 +611,75 @@ void IBcumulantK17LBMKernel::calculate(int step)
CUMcbb = -O4 * (c1 / omega - c1o2) * Dyz * c1o3 * B + (c1 - O4) * (CUMcbb);
//////////////////////////////////////////////////////////////////////////
- //5.
+ // 5.
CUMbcc += O5 * (-CUMbcc);
CUMcbc += O5 * (-CUMcbc);
CUMccb += O5 * (-CUMccb);
//////////////////////////////////////////////////////////////////////////
- //6.
+ // 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>
+ //! <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.
+ // 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;
+ 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;
+ // 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));
-
+ // 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>
+ //! <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;
@@ -550,10 +687,11 @@ void IBcumulantK17LBMKernel::calculate(int step)
////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////
- //! - 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>
+ //! - 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
@@ -593,12 +731,13 @@ void IBcumulantK17LBMKernel::calculate(int step)
////////////////////////////////////////////////////////////////////////////////////
//////////////////////////////////////////////////////////////////////////
- //proof correctness
+ // 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;
+#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)
@@ -606,51 +745,165 @@ void IBcumulantK17LBMKernel::calculate(int step)
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)));
+ 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>
+ //! - 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_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_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->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;
-
-
+ 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/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.h
similarity index 99%
rename from src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h
rename to src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.h
index 1d3fb0f568cf5c71191c3ee3981e9a454d6b5d3a..503d1970979afa83cf6d1d65ccd43356a5123720 100644
--- a/src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h
+++ b/src/cpu/LiggghtsCoupling/IBcumulantK17LBMKernel.h
@@ -71,7 +71,6 @@ protected:
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;
diff --git a/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h b/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
index 8bf11ca7a052d7f403c4b147975f1086a0402e1f..28f11bf2eb84a509d59e6c098c0fbb69874f9517 100644
--- a/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
+++ b/src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h
@@ -41,10 +41,21 @@ 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. };
+ 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
index ae3335374161bff8ff8c189a54d634cb53e64bce..d4926e0b11bd96b4d097ccf79060b437de6ce8da 100644
--- a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
@@ -9,11 +9,12 @@
#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)
- : CoProcessor(grid, s), comm(comm), wrapper(wrapper), demSteps(demSteps)
+ 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();
@@ -32,9 +33,11 @@ LiggghtsCouplingCoProcessor::~LiggghtsCouplingCoProcessor()
void LiggghtsCouplingCoProcessor::process(double actualTimeStep)
{
- setSpheresOnLattice();
- wrapper.run(demSteps);
std::cout << "step: " << actualTimeStep << "\n";
+
+ wrapper.run(demSteps);
+
+ setSpheresOnLattice();
}
void LiggghtsCouplingCoProcessor::setSpheresOnLattice()
@@ -64,12 +67,12 @@ void LiggghtsCouplingCoProcessor::setSpheresOnLattice()
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];
+ 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];
+ 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;
@@ -80,8 +83,6 @@ void LiggghtsCouplingCoProcessor::setSpheresOnLattice()
}
}
-void LiggghtsCouplingCoProcessor::getForcesFromLattice() {}
-
void LiggghtsCouplingCoProcessor::setSingleSphere3D(double *x, double *v, double *omega, /* double *com,*/ double r,
int id /*, bool initVelFlag*/)
{
@@ -98,7 +99,6 @@ void LiggghtsCouplingCoProcessor::setSingleSphere3D(double *x, double *v, double
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 =
@@ -111,21 +111,27 @@ void LiggghtsCouplingCoProcessor::setSingleSphere3D(double *x, double *v, double
int minX2 = 1;
int minX3 = 1;
- int maxX1 = (int)(distributions->getNX1()) - 1;
- int maxX2 = (int)(distributions->getNX2()) - 1;
- int maxX3 = (int)(distributions->getNX3()) - 1;
+ int maxX1 = (int)(distributions->getNX1()) - 2;
+ int maxX2 = (int)(distributions->getNX2()) - 2;
+ int maxX3 = (int)(distributions->getNX3()) - 2;
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];
+ //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);
+ 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;
@@ -137,7 +143,7 @@ void LiggghtsCouplingCoProcessor::setSingleSphere3D(double *x, double *v, double
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);
+ 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
@@ -145,12 +151,17 @@ void LiggghtsCouplingCoProcessor::setSingleSphere3D(double *x, double *v, double
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);
+ 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;
+ //}
}
}
}
@@ -206,10 +217,12 @@ double LiggghtsCouplingCoProcessor::calcSolidFraction(double const dx_, double c
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)
+ 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.from_cArray(v);
+ 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;
@@ -227,4 +240,172 @@ void LiggghtsCouplingCoProcessor::setToZero(IBdynamicsParticleData &p)
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()) - 2;
+ int maxX2 = (int)(distributions->getNX2()) - 2;
+ int maxX3 = (int)(distributions->getNX3()) - 2;
+
+ 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
index 509e83b4a5ed3150d40c539084d9c854870f8cbe..dcaa6e16c0a46a69d64998285119e86077500b6b 100644
--- a/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
+++ b/src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
@@ -51,12 +51,18 @@ 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);
+ LiggghtsCouplingWrapper &wrapper, int demSteps, SPtr<LBMUnitConverter> units);
virtual ~LiggghtsCouplingCoProcessor();
void process(double actualTimeStep) override;
@@ -64,24 +70,34 @@ public:
protected:
void setSpheresOnLattice();
- void getForcesFromLattice();
- void setSingleSphere3D(double *x, double *v, double *omega, /* double *com,*/ double r,
- int id /*, bool initVelFlag*/);
+
+ 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 setValues(IBdynamicsParticleData &p, double const sf, double const dx, double const dy, double const dz,
- double *omega, int id);
+ void addForce(int const partId, int const coord, double const value, double *force);
- void setToZero(IBdynamicsParticleData &p);
+ 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;
+ //std::vector<std::vector<SPtr<Block3D>>> blockVector;
+ //int minInitLevel;
+ //int maxInitLevel;
int gridRank;
+
+ double *force, *torque;
};
#endif
diff --git a/src/cpu/VirtualFluidsCore/CMakeLists.txt b/src/cpu/VirtualFluidsCore/CMakeLists.txt
index ac6302dba2f89114906356d58221f9dcae3d4bad..260c9d50eb67b471d6c6396c6b4dd45e8d685f35 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()
@@ -25,12 +24,13 @@ if(BUILD_USE_OPENMP)
list(APPEND VF_LIBRARIES OpenMP::OpenMP_CXX)
endif()
-vf_add_library(BUILDTYPE static PUBLIC_LINK basics muparser MPI::MPI_CXX ${VF_LIBRARIES} liggghts)
+IF(${USE_LIGGGHTS})
+ list(APPEND VF_LIBRARIES optimized ${LIGGGHTS_RELEASE_LIBRARY} debug ${LIGGGHTS_DEBUG_LIBRARY})
+ENDIF()
-vf_get_library_name(library_name)
+vf_add_library(BUILDTYPE static PUBLIC_LINK basics muparser MPI::MPI_CXX ${VF_LIBRARIES})
-target_link_directories(${library_name} PUBLIC d:/Tools/LIGGGHTS/build/Debug)
-target_include_directories(${library_name} PUBLIC d:/Tools/LIGGGHTS/src)
+vf_get_library_name(library_name)
target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/BoundaryConditions)
target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/Connectors)
@@ -42,8 +42,6 @@ target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/Gr
target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/Visitors)
target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/CoProcessors)
target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/Utilities)
-target_include_directories(${library_name} PUBLIC ${CMAKE_CURRENT_SOURCE_DIR}/LiggghtsCoupling)
-target_include_directories(${library_name} PUBLIC src/LiggghtsCoupling/Rconstructor)
IF(${USE_METIS} AND METIS_INCLUDEDIR)
@@ -51,6 +49,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 8d5cec10521830ba2aec9f2c06ef2a796da6b954..48d3a5c7f69c230622e42be2722d47751259cfef 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