Skip to content
Snippets Groups Projects
Commit 75cd008b authored by Konstantin Kutscher's avatar Konstantin Kutscher
Browse files

add Immersed Boundary: it doesn't work yet

parent 4291d958
No related branches found
No related tags found
2 merge requests!171Newest Update,!83Fix MPICommunicator
Expand all Hide whitespace changes
Inline Side-by-side
Showing
with 1311 additions and 4 deletions
apps/cpu/LiggghtsApp/LiggghtsApp.cpp
+ 100
− 1
View file @ 75cd008b
......@@ -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"};
......
apps/cpu/LiggghtsApp/in.lbdem 0 → 100644
+ 72
− 0
View file @ 75cd008b
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
apps/cpu/LiggghtsApp/in2.lbdem 0 → 100644
+ 21
− 0
View file @ 75cd008b
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
src/cpu/LiggghtsCoupling/IBdynamicsParticleData.h 0 → 100644
+ 50
− 0
View file @ 75cd008b
//=======================================================================================
// ____ ____ __ ______ __________ __ __ __ __
// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
// \ \ | | | | | |_) | | | | | | | / \ | |
// \ \ | | | | | _ / | | | | | | / /\ \ | |
// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
// \ \ | | ________________________________________________________________
// \ \ | | | ______________________________________________________________|
// \ \| | | | __ __ __ __ ______ _______
// \ | | |_____ | | | | | | | | | _ \ / _____)
// \ | | _____| | | | | | | | | | | \ \ \_______
// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
// \ _____| |__| |________| \_______/ |__| |______/ (_______/
//
// 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
src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.cpp
+ 218
− 1
View file @ 75cd008b
#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
src/cpu/LiggghtsCoupling/LiggghtsCouplingCoProcessor.h
+ 34
− 2
View file @ 75cd008b
......@@ -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
......
src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.cpp 0 → 100644
+ 660
− 0
View file @ 75cd008b
This diff is collapsed.
src/cpu/VirtualFluidsCore/LBM/IBcumulantK17LBMKernel.h 0 → 100644
+ 156
− 0
View file @ 75cd008b
//=======================================================================================
// ____ ____ __ ______ __________ __ __ __ __
// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
// \ \ | | | | | |_) | | | | | | | / \ | |
// \ \ | | | | | _ / | | | | | | / /\ \ | |
// \ \ | | | | | | \ \ | | | \__/ | / ____ \ | |____
// \ \ | | |__| |__| \__\ |__| \________/ /__/ \__\ |_______|
// \ \ | | ________________________________________________________________
// \ \ | | | ______________________________________________________________|
// \ \| | | | __ __ __ __ ______ _______
// \ | | |_____ | | | | | | | | | _ \ / _____)
// \ | | _____| | | | | | | | | | | \ \ \_______
// \ | | | | |_____ | \_/ | | | | |_/ / _____ |
// \ _____| |__| |________| \_______/ |__| |______/ (_______/
//
// 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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment