diff --git a/apps/cpu/CouetteFlow/cflow.cpp b/apps/cpu/CouetteFlow/cflow.cpp
index f5a3abca7007e5907cd567460474105ba620d556..72fd03b2ffbdeb91eb1af9faaf22ec10064ab69d 100644
--- a/apps/cpu/CouetteFlow/cflow.cpp
+++ b/apps/cpu/CouetteFlow/cflow.cpp
@@ -34,8 +34,8 @@ void bflow(string configname)
//double tau0 = config.getValue<double>("tau0");
double velocity = config.getValue<double>("velocity");
double n = config.getValue<double>("n");
- double Re = config.getValue<double>("Re");
- double Bn = config.getValue<double>("Bn");
+// double Re = config.getValue<double>("Re");
+// double Bn = config.getValue<double>("Bn");
SPtr<Communicator> comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
@@ -81,8 +81,8 @@ void bflow(string configname)
double blockLength = 3.0 * deltax;
- double h = (g_maxX2) / 2.0;
- double dex = g_maxX1;
+// double h = (g_maxX2) / 2.0;
+// double dex = g_maxX1;
//LBMReal tau0 = 1.2e-7;
//LBMReal k = nuLB;
diff --git a/apps/cpu/FlowAroundCylinder/cylinder.cpp b/apps/cpu/FlowAroundCylinder/cylinder.cpp
index 8e0073b779fac7f678cc8688e66bb52c9a4584c4..044ef765fd5af2a25c3719d9cb4dac403c9f0ee0 100644
--- a/apps/cpu/FlowAroundCylinder/cylinder.cpp
+++ b/apps/cpu/FlowAroundCylinder/cylinder.cpp
@@ -73,7 +73,7 @@ void run(string configname)
SPtr<LBMUnitConverter> conv = SPtr<LBMUnitConverter>(new LBMUnitConverter());
- const int baseLevel = 0;
+// const int baseLevel = 0;
SPtr<Grid3D> grid(new Grid3D(comm));
diff --git a/apps/cpu/HerschelBulkleyModel/hbflow.cpp b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
index d83f0f6d0494e2da4e1771273d14b447c13850f9..10381afe2bd95ca7df043e20ab92d4d40553244d 100644
--- a/apps/cpu/HerschelBulkleyModel/hbflow.cpp
+++ b/apps/cpu/HerschelBulkleyModel/hbflow.cpp
@@ -34,8 +34,8 @@ void bflow(string configname)
double tau0 = config.getValue<double>("tau0");
double velocity = config.getValue<double>("velocity");
double n = config.getValue<double>("n");
- double Re = config.getValue<double>("Re");
- double Bn = config.getValue<double>("Bn");
+// double Re = config.getValue<double>("Re");
+// double Bn = config.getValue<double>("Bn");
double scaleFactor = config.getValue<double>("scaleFactor");
SPtr<Communicator> comm = MPICommunicator::getInstance();
@@ -84,8 +84,8 @@ void bflow(string configname)
double blockLength = 3.0 * deltax;
- double h = (g_maxX2) / 2.0;
- double dex = g_maxX1;
+// double h = (g_maxX2) / 2.0;
+// double dex = g_maxX1;
//LBMReal tau0 = 1.2e-7;
//LBMReal k = nuLB;
diff --git a/apps/cpu/rheometer/rheometer.cpp b/apps/cpu/rheometer/rheometer.cpp
index 96c50e5c3c2640ac76b8582ab24432316c130397..fba3d72356ec2e37bace1f7ad46c80ec1f45301c 100644
--- a/apps/cpu/rheometer/rheometer.cpp
+++ b/apps/cpu/rheometer/rheometer.cpp
@@ -105,7 +105,7 @@ void bflow(string configname)
//double g_maxX2 = boundingBox[1]/2.0;
//double g_maxX3 = boundingBox[2]/2.0;
- double blockLength = 3.0 * deltax;
+// double blockLength = 3.0 * deltax;
// double d = 2.0 * radius;
// double U = uLB;
diff --git a/apps/cpu/sphere/sphere.cpp b/apps/cpu/sphere/sphere.cpp
index 388d36718312e74958ef54e6eadd071ef11507bc..f7f8481d11a14fb458dfde70d72ddd0bcf098276 100644
--- a/apps/cpu/sphere/sphere.cpp
+++ b/apps/cpu/sphere/sphere.cpp
@@ -12,8 +12,6 @@ void run(string configname)
SPtr<Communicator> comm = MPICommunicator::getInstance();
int myid = comm->getProcessID();
- int mybundle = comm->getBundleID();
- int root = comm->getRoot();
//ConfigurationFile config;
//config.load(configname);
@@ -41,7 +39,7 @@ void run(string configname)
LBMReal nuLB = (uLB*2.0*radius)/Re;
double dp_LB = 1e-6;
- double rhoLBinflow = dp_LB*3.0;
+// double rhoLBinflow = dp_LB*3.0;
SPtr<BCAdapter> noSlipBCAdapter(new NoSlipBCAdapter());
noSlipBCAdapter->setBcAlgorithm(SPtr<BCAlgorithm>(new NoSlipBCAlgorithm()));
@@ -98,9 +96,6 @@ void run(string configname)
if (true)
{
-
- const int baseLevel = 0;
-
//bounding box
double d_minX1 = 0.0;
double d_minX2 = 0.0;
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp
index 006f45619d46051cea296584e48440acdddd7291..3c96343077f1030f7e6477bb352635c870ece4ce 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowAndThixotropyBCAlgorithm.cpp
@@ -64,15 +64,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
switch (direction)
{
case E:
- f[E] = ftemp[E] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[E];
- f[NE] = ftemp[NE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[NE];
- f[SE] = ftemp[SE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[SE];
- f[TE] = ftemp[TE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TE];
- f[BE] = ftemp[BE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BE];
- f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TNE];
- f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TSE];
- f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BNE];
- f[BSE] = ftemp[BSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BSE];
+ f[E] = ftemp[E] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[E];
+ f[NE] = ftemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[NE];
+ f[SE] = ftemp[SE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[SE];
+ f[TE] = ftemp[TE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TE];
+ f[BE] = ftemp[BE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BE];
+ f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TNE];
+ f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TSE];
+ f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BNE];
+ f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BSE];
distributions->setDistributionInvForDirection(f[E], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -85,15 +85,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BSE], x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW], TNW);
break;
case W:
- f[W] = ftemp[W] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[W];
- f[NW] = ftemp[NW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[NW];
- f[SW] = ftemp[SW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[SW];
- f[TW] = ftemp[TW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TW];
- f[BW] = ftemp[BW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BW];
- f[TNW] = ftemp[TNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TNW];
- f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TSW];
- f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BNW];
- f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BSW];
+ f[W] = ftemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[W];
+ f[NW] = ftemp[NW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[NW];
+ f[SW] = ftemp[SW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[SW];
+ f[TW] = ftemp[TW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TW];
+ f[BW] = ftemp[BW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BW];
+ f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TNW];
+ f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TSW];
+ f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BNW];
+ f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BSW];
distributions->setDistributionInvForDirection(f[W], x1 + DX1[E], x2 + DX2[E], x3 + DX3[E], E);
distributions->setDistributionInvForDirection(f[NW], x1 + DX1[SE], x2 + DX2[SE], x3 + DX3[SE], SE);
@@ -106,15 +106,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
break;
case N:
- f[N] = ftemp[N] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[N];
- f[NE] = ftemp[NE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NE];
- f[NW] = ftemp[NW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NW];
- f[TN] = ftemp[TN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TN];
- f[BN] = ftemp[BN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BN];
- f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNE];
- f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNW];
- f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNE];
- f[BNW] = ftemp[BNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNW];
+ f[N] = ftemp[N] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[N];
+ f[NE] = ftemp[NE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NE];
+ f[NW] = ftemp[NW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NW];
+ f[TN] = ftemp[TN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TN];
+ f[BN] = ftemp[BN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BN];
+ f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNE];
+ f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNW];
+ f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNE];
+ f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNW];
distributions->setDistributionInvForDirection(f[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -127,15 +127,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BNW], x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE], TSE);
break;
case S:
- f[S] = ftemp[S] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[S];
- f[SE] = ftemp[SE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SE];
- f[SW] = ftemp[SW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SW];
- f[TS] = ftemp[TS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TS];
- f[BS] = ftemp[BS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BS];
- f[TSE] = ftemp[TSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSE];
- f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSW];
- f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSE];
- f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSW];
+ f[S] = ftemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[S];
+ f[SE] = ftemp[SE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SE];
+ f[SW] = ftemp[SW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SW];
+ f[TS] = ftemp[TS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TS];
+ f[BS] = ftemp[BS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BS];
+ f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSE];
+ f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSW];
+ f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSE];
+ f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSW];
distributions->setDistributionInvForDirection(f[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
distributions->setDistributionInvForDirection(f[SE], x1 + DX1[NW], x2 + DX2[NW], x3 + DX3[NW], NW);
@@ -148,15 +148,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
break;
case T:
- f[T] = ftemp[T] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[T];
- f[TE] = ftemp[TE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TE];
- f[TW] = ftemp[TW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TW];
- f[TN] = ftemp[TN] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TN];
- f[TS] = ftemp[TS] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TS];
- f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNE];
- f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNW];
- f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSE];
- f[TSW] = ftemp[TSW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSW];
+ f[T] = ftemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[T];
+ f[TE] = ftemp[TE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TE];
+ f[TW] = ftemp[TW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TW];
+ f[TN] = ftemp[TN] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TN];
+ f[TS] = ftemp[TS] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TS];
+ f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNE];
+ f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNW];
+ f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSE];
+ f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSW];
distributions->setDistributionInvForDirection(f[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
distributions->setDistributionInvForDirection(f[TE], x1 + DX1[BW], x2 + DX2[BW], x3 + DX3[BW], BW);
@@ -169,15 +169,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[TSW], x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE], BNE);
break;
case B:
- f[B] = ftemp[B] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[B];
- f[BE] = ftemp[BE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BE];
- f[BW] = ftemp[BW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BW];
- f[BN] = ftemp[BN] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BN];
- f[BS] = ftemp[BS] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BS];
- f[BNE] = ftemp[BNE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNE];
- f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNW];
- f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSE];
- f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSW];
+ f[B] = ftemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[B];
+ f[BE] = ftemp[BE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BE];
+ f[BW] = ftemp[BW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BW];
+ f[BN] = ftemp[BN] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BN];
+ f[BS] = ftemp[BS] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BS];
+ f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNE];
+ f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNW];
+ f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSE];
+ f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSW];
distributions->setDistributionInvForDirection(f[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
distributions->setDistributionInvForDirection(f[BE], x1 + DX1[TW], x2 + DX2[TW], x3 + DX3[TW], TW);
@@ -208,15 +208,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
switch (direction)
{
case E:
- h[E] = htemp[E] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[E];
- h[NE] = htemp[NE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[NE];
- h[SE] = htemp[SE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[SE];
- h[TE] = htemp[TE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[TE];
- h[BE] = htemp[BE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[BE];
- h[TNE] = htemp[TNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[TNE];
- h[TSE] = htemp[TSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[TSE];
- h[BNE] = htemp[BNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[BNE];
- h[BSE] = htemp[BSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * h[BSE];
+ h[E] = htemp[E] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[E];
+ h[NE] = htemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[NE];
+ h[SE] = htemp[SE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[SE];
+ h[TE] = htemp[TE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[TE];
+ h[BE] = htemp[BE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[BE];
+ h[TNE] = htemp[TNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[TNE];
+ h[TSE] = htemp[TSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[TSE];
+ h[BNE] = htemp[BNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[BNE];
+ h[BSE] = htemp[BSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * h[BSE];
distributionsH->setDistributionInvForDirection(h[E], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
distributionsH->setDistributionInvForDirection(h[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -229,15 +229,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributionsH->setDistributionInvForDirection(h[BSE], x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW], TNW);
break;
case W:
- h[W] = htemp[W] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[W];
- h[NW] = htemp[NW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[NW];
- h[SW] = htemp[SW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[SW];
- h[TW] = htemp[TW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[TW];
- h[BW] = htemp[BW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[BW];
- h[TNW] = htemp[TNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[TNW];
- h[TSW] = htemp[TSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[TSW];
- h[BNW] = htemp[BNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[BNW];
- h[BSW] = htemp[BSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * h[BSW];
+ h[W] = htemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[W];
+ h[NW] = htemp[NW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[NW];
+ h[SW] = htemp[SW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[SW];
+ h[TW] = htemp[TW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[TW];
+ h[BW] = htemp[BW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[BW];
+ h[TNW] = htemp[TNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[TNW];
+ h[TSW] = htemp[TSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[TSW];
+ h[BNW] = htemp[BNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[BNW];
+ h[BSW] = htemp[BSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * h[BSW];
distributionsH->setDistributionInvForDirection(h[W], x1 + DX1[E], x2 + DX2[E], x3 + DX3[E], E);
distributionsH->setDistributionInvForDirection(h[NW], x1 + DX1[SE], x2 + DX2[SE], x3 + DX3[SE], SE);
@@ -250,15 +250,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributionsH->setDistributionInvForDirection(h[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
break;
case N:
- h[N] = htemp[N] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[N];
- h[NE] = htemp[NE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[NE];
- h[NW] = htemp[NW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[NW];
- h[TN] = htemp[TN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[TN];
- h[BN] = htemp[BN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[BN];
- h[TNE] = htemp[TNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[TNE];
- h[TNW] = htemp[TNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[TNW];
- h[BNE] = htemp[BNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[BNE];
- h[BNW] = htemp[BNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * h[BNW];
+ h[N] = htemp[N] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[N];
+ h[NE] = htemp[NE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[NE];
+ h[NW] = htemp[NW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[NW];
+ h[TN] = htemp[TN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[TN];
+ h[BN] = htemp[BN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[BN];
+ h[TNE] = htemp[TNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[TNE];
+ h[TNW] = htemp[TNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[TNW];
+ h[BNE] = htemp[BNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[BNE];
+ h[BNW] = htemp[BNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * h[BNW];
distributionsH->setDistributionInvForDirection(h[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
distributionsH->setDistributionInvForDirection(h[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -271,15 +271,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributionsH->setDistributionInvForDirection(h[BNW], x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE], TSE);
break;
case S:
- h[S] = htemp[S] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[S];
- h[SE] = htemp[SE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[SE];
- h[SW] = htemp[SW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[SW];
- h[TS] = htemp[TS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[TS];
- h[BS] = htemp[BS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[BS];
- h[TSE] = htemp[TSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[TSE];
- h[TSW] = htemp[TSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[TSW];
- h[BSE] = htemp[BSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[BSE];
- h[BSW] = htemp[BSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * h[BSW];
+ h[S] = htemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[S];
+ h[SE] = htemp[SE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[SE];
+ h[SW] = htemp[SW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[SW];
+ h[TS] = htemp[TS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[TS];
+ h[BS] = htemp[BS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[BS];
+ h[TSE] = htemp[TSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[TSE];
+ h[TSW] = htemp[TSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[TSW];
+ h[BSE] = htemp[BSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[BSE];
+ h[BSW] = htemp[BSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * h[BSW];
distributionsH->setDistributionInvForDirection(h[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
distributionsH->setDistributionInvForDirection(h[SE], x1 + DX1[NW], x2 + DX2[NW], x3 + DX3[NW], NW);
@@ -292,15 +292,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributionsH->setDistributionInvForDirection(h[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
break;
case T:
- h[T] = htemp[T] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[T];
- h[TE] = htemp[TE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TE];
- h[TW] = htemp[TW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TW];
- h[TN] = htemp[TN] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TN];
- h[TS] = htemp[TS] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TS];
- h[TNE] = htemp[TNE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TNE];
- h[TNW] = htemp[TNW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TNW];
- h[TSE] = htemp[TSE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TSE];
- h[TSW] = htemp[TSW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * h[TSW];
+ h[T] = htemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[T];
+ h[TE] = htemp[TE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TE];
+ h[TW] = htemp[TW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TW];
+ h[TN] = htemp[TN] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TN];
+ h[TS] = htemp[TS] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TS];
+ h[TNE] = htemp[TNE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TNE];
+ h[TNW] = htemp[TNW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TNW];
+ h[TSE] = htemp[TSE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TSE];
+ h[TSW] = htemp[TSW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * h[TSW];
distributionsH->setDistributionInvForDirection(h[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
distributionsH->setDistributionInvForDirection(h[TE], x1 + DX1[BW], x2 + DX2[BW], x3 + DX3[BW], BW);
@@ -313,15 +313,15 @@ void NonReflectingOutflowAndThixotropyBCAlgorithm::applyBC()
distributionsH->setDistributionInvForDirection(h[TSW], x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE], BNE);
break;
case B:
- h[B] = htemp[B] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[B];
- h[BE] = htemp[BE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BE];
- h[BW] = htemp[BW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BW];
- h[BN] = htemp[BN] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BN];
- h[BS] = htemp[BS] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BS];
- h[BNE] = htemp[BNE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BNE];
- h[BNW] = htemp[BNW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BNW];
- h[BSE] = htemp[BSE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BSE];
- h[BSW] = htemp[BSW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * h[BSW];
+ h[B] = htemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[B];
+ h[BE] = htemp[BE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BE];
+ h[BW] = htemp[BW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BW];
+ h[BN] = htemp[BN] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BN];
+ h[BS] = htemp[BS] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BS];
+ h[BNE] = htemp[BNE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BNE];
+ h[BNW] = htemp[BNW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BNW];
+ h[BSE] = htemp[BSE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BSE];
+ h[BSW] = htemp[BSW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * h[BSW];
distributionsH->setDistributionInvForDirection(h[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
distributionsH->setDistributionInvForDirection(h[BE], x1 + DX1[TW], x2 + DX2[TW], x3 + DX3[TW], TW);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
index 0b5884d12229cfccab6647bd0ec52add2e238bbc..0a4280c75832c4eb6c71a62682f346aed4fba64f 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowBCAlgorithm.cpp
@@ -66,15 +66,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
switch (direction) {
case E:
- f[E] = ftemp[E] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[E];
- f[NE] = ftemp[NE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[NE];
- f[SE] = ftemp[SE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[SE];
- f[TE] = ftemp[TE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TE];
- f[BE] = ftemp[BE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BE];
- f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TNE];
- f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[TSE];
- f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BNE];
- f[BSE] = ftemp[BSE] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[BSE];
+ f[E] = ftemp[E] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[E];
+ f[NE] = ftemp[NE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[NE];
+ f[SE] = ftemp[SE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[SE];
+ f[TE] = ftemp[TE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TE];
+ f[BE] = ftemp[BE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BE];
+ f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TNE];
+ f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[TSE];
+ f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BNE];
+ f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 + vx1) + (1.0 - UbMath::one_over_sqrt3 - vx1) * f[BSE];
distributions->setDistributionInvForDirection(f[E], x1 + DX1[W], x2 + DX2[W], x3 + DX3[W], W);
distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -87,15 +87,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BSE], x1 + DX1[TNW], x2 + DX2[TNW], x3 + DX3[TNW], TNW);
break;
case W:
- f[W] = ftemp[W] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[W];
- f[NW] = ftemp[NW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[NW];
- f[SW] = ftemp[SW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[SW];
- f[TW] = ftemp[TW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TW];
- f[BW] = ftemp[BW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BW];
- f[TNW] = ftemp[TNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TNW];
- f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[TSW];
- f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BNW];
- f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[BSW];
+ f[W] = ftemp[W] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[W];
+ f[NW] = ftemp[NW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[NW];
+ f[SW] = ftemp[SW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[SW];
+ f[TW] = ftemp[TW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TW];
+ f[BW] = ftemp[BW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BW];
+ f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TNW];
+ f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[TSW];
+ f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BNW];
+ f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx1) + (1.0 - UbMath::one_over_sqrt3 + vx1) * f[BSW];
distributions->setDistributionInvForDirection(f[W], x1 + DX1[E], x2 + DX2[E], x3 + DX3[E], E);
distributions->setDistributionInvForDirection(f[NW], x1 + DX1[SE], x2 + DX2[SE], x3 + DX3[SE], SE);
@@ -108,15 +108,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
break;
case N:
- f[N] = ftemp[N] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[N];
- f[NE] = ftemp[NE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NE];
- f[NW] = ftemp[NW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[NW];
- f[TN] = ftemp[TN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TN];
- f[BN] = ftemp[BN] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BN];
- f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNE];
- f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[TNW];
- f[BNE] = ftemp[BNE] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNE];
- f[BNW] = ftemp[BNW] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[BNW];
+ f[N] = ftemp[N] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[N];
+ f[NE] = ftemp[NE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NE];
+ f[NW] = ftemp[NW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[NW];
+ f[TN] = ftemp[TN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TN];
+ f[BN] = ftemp[BN] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BN];
+ f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNE];
+ f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[TNW];
+ f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNE];
+ f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 + vx2) + (1.0 - UbMath::one_over_sqrt3 - vx2) * f[BNW];
distributions->setDistributionInvForDirection(f[N], x1 + DX1[S], x2 + DX2[S], x3 + DX3[S], S);
distributions->setDistributionInvForDirection(f[NE], x1 + DX1[SW], x2 + DX2[SW], x3 + DX3[SW], SW);
@@ -129,15 +129,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BNW], x1 + DX1[TSE], x2 + DX2[TSE], x3 + DX3[TSE], TSE);
break;
case S:
- f[S] = ftemp[S] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[S];
- f[SE] = ftemp[SE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SE];
- f[SW] = ftemp[SW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[SW];
- f[TS] = ftemp[TS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TS];
- f[BS] = ftemp[BS] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BS];
- f[TSE] = ftemp[TSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSE];
- f[TSW] = ftemp[TSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[TSW];
- f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSE];
- f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[BSW];
+ f[S] = ftemp[S] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[S];
+ f[SE] = ftemp[SE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SE];
+ f[SW] = ftemp[SW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[SW];
+ f[TS] = ftemp[TS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TS];
+ f[BS] = ftemp[BS] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BS];
+ f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSE];
+ f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[TSW];
+ f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSE];
+ f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx2) + (1.0 - UbMath::one_over_sqrt3 + vx2) * f[BSW];
distributions->setDistributionInvForDirection(f[S], x1 + DX1[N], x2 + DX2[N], x3 + DX3[N], N);
distributions->setDistributionInvForDirection(f[SE], x1 + DX1[NW], x2 + DX2[NW], x3 + DX3[NW], NW);
@@ -150,15 +150,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[BSW], x1 + DX1[TNE], x2 + DX2[TNE], x3 + DX3[TNE], TNE);
break;
case T:
- f[T] = ftemp[T] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[T];
- f[TE] = ftemp[TE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TE];
- f[TW] = ftemp[TW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TW];
- f[TN] = ftemp[TN] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TN];
- f[TS] = ftemp[TS] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TS];
- f[TNE] = ftemp[TNE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNE];
- f[TNW] = ftemp[TNW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TNW];
- f[TSE] = ftemp[TSE] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSE];
- f[TSW] = ftemp[TSW] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[TSW];
+ f[T] = ftemp[T] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[T];
+ f[TE] = ftemp[TE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TE];
+ f[TW] = ftemp[TW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TW];
+ f[TN] = ftemp[TN] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TN];
+ f[TS] = ftemp[TS] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TS];
+ f[TNE] = ftemp[TNE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNE];
+ f[TNW] = ftemp[TNW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TNW];
+ f[TSE] = ftemp[TSE] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSE];
+ f[TSW] = ftemp[TSW] * (UbMath::one_over_sqrt3 + vx3) + (1.0 - UbMath::one_over_sqrt3 - vx3) * f[TSW];
distributions->setDistributionInvForDirection(f[T], x1 + DX1[B], x2 + DX2[B], x3 + DX3[B], B);
distributions->setDistributionInvForDirection(f[TE], x1 + DX1[BW], x2 + DX2[BW], x3 + DX3[BW], BW);
@@ -171,15 +171,15 @@ void NonReflectingOutflowBCAlgorithm::applyBC()
distributions->setDistributionInvForDirection(f[TSW], x1 + DX1[BNE], x2 + DX2[BNE], x3 + DX3[BNE], BNE);
break;
case B:
- f[B] = ftemp[B] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[B];
- f[BE] = ftemp[BE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BE];
- f[BW] = ftemp[BW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BW];
- f[BN] = ftemp[BN] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BN];
- f[BS] = ftemp[BS] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BS];
- f[BNE] = ftemp[BNE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNE];
- f[BNW] = ftemp[BNW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BNW];
- f[BSE] = ftemp[BSE] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSE];
- f[BSW] = ftemp[BSW] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[BSW];
+ f[B] = ftemp[B] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[B];
+ f[BE] = ftemp[BE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BE];
+ f[BW] = ftemp[BW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BW];
+ f[BN] = ftemp[BN] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BN];
+ f[BS] = ftemp[BS] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BS];
+ f[BNE] = ftemp[BNE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNE];
+ f[BNW] = ftemp[BNW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BNW];
+ f[BSE] = ftemp[BSE] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSE];
+ f[BSW] = ftemp[BSW] * (UbMath::one_over_sqrt3 - vx3) + (1.0 - UbMath::one_over_sqrt3 + vx3) * f[BSW];
distributions->setDistributionInvForDirection(f[B], x1 + DX1[T], x2 + DX2[T], x3 + DX3[T], T);
distributions->setDistributionInvForDirection(f[BE], x1 + DX1[TW], x2 + DX2[TW], x3 + DX3[TW], TW);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
index 2614c703b4da29f666392edbe121b86d3ad78393..4ea802eb4f798526b822fc69c3b9111207973cc8 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.cpp
@@ -26,10 +26,10 @@ void RheologicalVelocityBCAlgorithm::applyBC()
calcMacrosFct(f, drho, vx1, vx2, vx3);
calcFeqFct(feq, drho, vx1, vx2, vx3);
- LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
- LBMReal collFactorF = getThyxotropyCollFactor(collFactor, shearRate, drho);
+ LBMReal shearRate = D3Q27System::getShearRate(f, collFactor);
+ LBMReal collFactorF = getThyxotropyCollFactor(collFactor, shearRate, drho);
- rho = 1.0+drho*compressibleFactor;
+ rho = 1.0+drho*compressibleFactor;
for (int fdir = D3Q27System::FSTARTDIR; fdir<=D3Q27System::FENDDIR; fdir++)
{
@@ -38,7 +38,7 @@ void RheologicalVelocityBCAlgorithm::applyBC()
const int invDir = D3Q27System::INVDIR[fdir];
LBMReal q = bcPtr->getQ(invDir);// m+m q=0 stabiler
LBMReal velocity = bcPtr->getBoundaryVelocity(invDir);
- LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactor)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity*rho)/(1.0+q));
+ LBMReal fReturn = ((1.0-q)/(1.0+q))*((f[invDir]-feq[invDir])/(1.0-collFactorF)+feq[invDir])+((q*(f[invDir]+f[fdir])-velocity*rho)/(1.0+q));
distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
}
}
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h
index 60237e3db6be7c1a77291146fe0f70db69018f0f..179bf57c3450f0cdaf77a9985fd31c7f1f28f9a6 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/RheologicalVelocityBCAlgorithm.h
@@ -12,7 +12,7 @@ public:
RheologicalVelocityBCAlgorithm();
~RheologicalVelocityBCAlgorithm();
virtual SPtr<BCAlgorithm> clone() override { UB_THROW(UbException("LBMReal clone() - belongs in the derived class")); }
- void addDistributions(SPtr<DistributionArray3D> distributions);
+ void addDistributions(SPtr<DistributionArray3D> distributions) override;
void applyBC() override;
protected:
virtual LBMReal getThyxotropyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho) const = 0; // { UB_THROW(UbException("LBMReal getThyxotropyCollFactor() - belongs in the derived class")); }
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.h
index 3da68623200d51d39e7c69b6ffb21554c9979eab..bfe4a766bc4fe1f05cf55cb679cfb3900f574c39 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleSlipBCAlgorithm.h
@@ -44,8 +44,8 @@ class SimpleSlipBCAlgorithm : public BCAlgorithm
public:
SimpleSlipBCAlgorithm();
virtual ~SimpleSlipBCAlgorithm();
- SPtr<BCAlgorithm> clone();
- void addDistributions(SPtr<DistributionArray3D> distributions);
+ SPtr<BCAlgorithm> clone() override;
+ void addDistributions(SPtr<DistributionArray3D> distributions) override;
void applyBC() override;
};
#endif // SimpleSlipBCAlgorithm_h__
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.cpp
index 27e0615e1879b6e52569ff354a152e132395aae8..9e136e0153a52ebb8b8c598701d4746fc21d4ab3 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.cpp
@@ -72,7 +72,6 @@ void SimpleVelocityBCAlgorithm::applyBC()
if (bcPtr->hasVelocityBoundaryFlag(fdir))
{
const int invDir = D3Q27System::INVDIR[fdir];
- LBMReal q = bcPtr->getQ(invDir);
LBMReal velocity = bcPtr->getBoundaryVelocity(invDir);
LBMReal fReturn = f[invDir] - velocity;
distributions->setDistributionForDirection(fReturn, x1+D3Q27System::DX1[invDir], x2+D3Q27System::DX2[invDir], x3+D3Q27System::DX3[invDir], fdir);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.h
index 714276760e379b93a433255a3a45981250f4bc91..b751d220c0cd15b597969fd9980184c16c5de36c 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/SimpleVelocityBCAlgorithm.h
@@ -45,8 +45,8 @@ class SimpleVelocityBCAlgorithm : public BCAlgorithm
public:
SimpleVelocityBCAlgorithm();
~SimpleVelocityBCAlgorithm();
- SPtr<BCAlgorithm> clone();
- void addDistributions(SPtr<DistributionArray3D> distributions);
+ SPtr<BCAlgorithm> clone() override;
+ void addDistributions(SPtr<DistributionArray3D> distributions) override;
void applyBC() override;
};
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp
index 0ee6d281d78998066a4f02e63a9aea0c5845369a..9fae576b70615b56e54bd673c9c26b3def83ad26 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.cpp
@@ -2,13 +2,6 @@
#include "DistributionArray3D.h"
#include "BoundaryConditions.h"
-ThixotropyNoSlipBCAlgorithm::ThixotropyNoSlipBCAlgorithm()
-{
-}
-//////////////////////////////////////////////////////////////////////////
-ThixotropyNoSlipBCAlgorithm::~ThixotropyNoSlipBCAlgorithm()
-{
-}
//////////////////////////////////////////////////////////////////////////
void ThixotropyNoSlipBCAlgorithm::addDistributions(SPtr<DistributionArray3D> distributions)
{
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.h b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.h
index 284cbbf2f6105c61d97784a0d215662a4fcb460e..51448a836607080fd34d97fae00546d67ae0ff6e 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.h
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/ThixotropyNoSlipBCAlgorithm.h
@@ -9,10 +9,10 @@ class DistributionArray3D;
class ThixotropyNoSlipBCAlgorithm : public BCAlgorithm
{
public:
- ThixotropyNoSlipBCAlgorithm();
- ~ThixotropyNoSlipBCAlgorithm();
+ ThixotropyNoSlipBCAlgorithm() = default;
+ ~ThixotropyNoSlipBCAlgorithm() = default;
virtual SPtr<BCAlgorithm> clone() override { UB_THROW(UbException("LBMReal clone() - belongs in the derived class")); }
- void addDistributions(SPtr<DistributionArray3D> distributions);
+ void addDistributions(SPtr<DistributionArray3D> distributions) override;
void applyBC() override;
protected:
virtual LBMReal getThyxotropyCollFactor(LBMReal omegaInf, LBMReal shearRate, LBMReal drho) const = 0; // { UB_THROW(UbException("LBMReal getThyxotropyCollFactor() - belongs in the derived class")); }
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityAndThixotropyBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityAndThixotropyBCAlgorithm.cpp
index d08d112d468229bc6bd948ab41472dd594045322..8f2436b387c1a38e35985738f6b17d45cdc431ef 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityAndThixotropyBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityAndThixotropyBCAlgorithm.cpp
@@ -81,7 +81,6 @@ void VelocityWithDensityAndThixotropyBCAlgorithm::applyBC()
if (bcArray->isSolid(nX1,nX2,nX3))
{
const int invDir = D3Q27System::INVDIR[fdir];
- LBMReal q =1.0;// bcPtr->getQ(invDir);// m+m q=0 stabiler
LBMReal velocity = bcPtr->getBoundaryVelocity(fdir);
LBMReal fReturn = (f[fdir] + f[invDir] - velocity*rho) / 2.0 - drho*D3Q27System::WEIGTH[invDir];
@@ -91,7 +90,6 @@ void VelocityWithDensityAndThixotropyBCAlgorithm::applyBC()
if (bcPtr->hasVelocityBoundaryFlag(fdir))
{
- const int invDir = D3Q27System::INVDIR[fdir];
LBMReal htemp = D3Q27System::getCompFeqForDirection(fdir, lambda, vx1, vx2, vx3);
htemp = D3Q27System::getCompFeqForDirection(fdir, lambdaBC, vx1, vx2, vx3) + h[fdir] - htemp;
distributionsH->setDistributionForDirection(htemp, nx1, nx2, nx3, fdir);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCAlgorithm.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCAlgorithm.cpp
index 96c4654394ca156a126bc0272ba92d087566f563..45615466c3ba40881dcf03dff27bbcb3f11adc80 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCAlgorithm.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCAlgorithm.cpp
@@ -25,7 +25,7 @@ void VelocityWithDensityBCAlgorithm::applyBC()
{
//velocity bc for non reflecting pressure bc
LBMReal f[D3Q27System::ENDF+1];
- LBMReal feq[D3Q27System::ENDF+1];
+ //LBMReal feq[D3Q27System::ENDF+1];
distributions->getDistributionInv(f, x1, x2, x3);
LBMReal rho, vx1, vx2, vx3, drho;
calcMacrosFct(f, drho, vx1, vx2, vx3);
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
index af7533ef5bffee3871f669adb2682412faf3b7bc..0976157630976f4cf429e29837e6c8aa6ffc5a10 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/CalculateTorqueCoProcessor.cpp
@@ -86,7 +86,6 @@ void CalculateTorqueCoProcessor::calculateForces()
forceX1global = 0.0;
forceX2global = 0.0;
forceX3global = 0.0;
- int counter = 0;
for(SPtr<D3Q27Interactor> interactor : interactors)
{
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
index 43687af9c025aa66b9db0afc5222716054dd488f..04191a9927f34d8616ae5bcb1d9670432e5ddb0e 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteMacroscopicQuantitiesCoProcessor.cpp
@@ -146,7 +146,6 @@ void WriteMacroscopicQuantitiesCoProcessor::clearData()
void WriteMacroscopicQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
{
double level = (double)block->getLevel();
- double blockID = (double)block->getGlobalID();
// Diese Daten werden geschrieben:
datanames.resize(0);
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
index 72e7114e6d13504727ea8d86d7641627df6056f6..27ec00a767b9ca19981f2606b657e0d551e19023 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.cpp
@@ -110,13 +110,12 @@ void WriteThixotropyQuantitiesCoProcessor::clearData()
//////////////////////////////////////////////////////////////////////////
void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
{
- UbTupleDouble3 org = grid->getBlockWorldCoordinates(block);
- UbTupleDouble3 blockLengths = grid->getBlockLengths(block);
+ UbTupleDouble3 org = grid->getBlockWorldCoordinates(block);;
UbTupleDouble3 nodeOffset = grid->getNodeOffset(block);
double dx = grid->getDeltaX(block);
- double level = (double)block->getLevel();
- double blockID = (double)block->getGlobalID();
+ //double level = (double)block->getLevel();
+ //double blockID = (double)block->getGlobalID();
//Diese Daten werden geschrieben:
datanames.resize(0);
@@ -217,7 +216,7 @@ void WriteThixotropyQuantitiesCoProcessor::addDataMQ(SPtr<Block3D> block)
//LBMReal omega = Thixotropy::getHerschelBulkleyCollFactor(collFactor, shearRate, rho);
//LBMReal omega = Thixotropy::getPowellEyringCollFactor(collFactor, shearRate, rho);
LBMReal omega = Thixotropy::getBinghamCollFactor(collFactor, shearRate, rho);
- LBMReal viscosity = (omega == 0) ? 0 : c1o3 * (c1/omega-c1o2);
+ LBMReal viscosity = (omega == 0) ? 0 : UbMath::c1o3 * (UbMath::c1/omega-UbMath::c1o2);
data[index++].push_back(viscosity);
diff --git a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.h b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.h
index 58430fd9ece771d9e9dacf496860787f7b1364f5..269bb4a3b9be8072f33793e6ce105481ebd3b243 100644
--- a/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.h
+++ b/src/cpu/VirtualFluidsCore/CoProcessors/WriteThixotropyQuantitiesCoProcessor.h
@@ -13,7 +13,7 @@ class WriteThixotropyQuantitiesCoProcessor : public CoProcessor
public:
WriteThixotropyQuantitiesCoProcessor();
WriteThixotropyQuantitiesCoProcessor(SPtr<Grid3D> grid, SPtr<UbScheduler> s, const std::string& path, WbWriter* const writer, SPtr<LBMUnitConverter> conv, SPtr<Communicator> comm);
- ~WriteThixotropyQuantitiesCoProcessor() {}
+ ~WriteThixotropyQuantitiesCoProcessor() = default;
void process(double step) override;
@@ -31,12 +31,12 @@ private:
std::string path;
WbWriter* writer;
SPtr<LBMUnitConverter> conv;
- bool bcInformation;
+// bool bcInformation;
std::vector<std::vector<SPtr<Block3D> > > blockVector;
int minInitLevel;
int maxInitLevel;
int gridRank;
SPtr<Communicator> comm;
- double ConcentrationSum;
+// double ConcentrationSum;
};
#endif
diff --git a/src/cpu/VirtualFluidsCore/Connectors/ThixotropyFullDirectConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/ThixotropyFullDirectConnector.cpp
index 3fb819ab7b90b10964d07758ff8fbcb8a8603f38..39561bd819e86e8beece92e52766396c2b031d91 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/ThixotropyFullDirectConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/ThixotropyFullDirectConnector.cpp
@@ -66,7 +66,7 @@ void ThixotropyFullDirectConnector::sendVectors()
nonLocalDistributionsToh = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(this->hTo)->getNonLocalDistributions();
zeroDistributionsToh = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(this->hTo)->getZeroDistributions();
- bool con = /*(from.lock()->getGlobalID()==11) &&*/ (to.lock()->getGlobalID() == 1);
+ //bool con = /*(from.lock()->getGlobalID()==11) &&*/ (to.lock()->getGlobalID() == 1);
//EAST
if (sendDir == D3Q27System::E)
diff --git a/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.cpp
index 553c49961fe2ef66df9ff3ef5511ff495f051766..12c4b6e891214a683e3e4ac23b52e14e72e093a2 100644
--- a/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.cpp
@@ -49,10 +49,6 @@ void BasicLBMKernel::calculate(int step)
{
if (!bcArray->isSolid(x1, x2, x3) && !bcArray->isUndefined(x1, x2, x3))
{
- int x1p = x1 + 1;
- int x2p = x2 + 1;
- int x3p = x3 + 1;
-
nodeCollision(step, x1, x2, x3);
}
}
diff --git a/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.h
index 4a6b4632f3648e2cc7e92e2e8314c9ea4fe2c616..2e8a6c5b330fd3ad3672249099096cc67a557c03 100644
--- a/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/BasicLBMKernel.h
@@ -1,4 +1,4 @@
-#ifndef BasicLBMKernell_h__
+#ifndef BasicLBMKernel_h__
#define BasicLBMKernel_h__
#include "LBMKernel.h"
diff --git a/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp
index 57e26c7e0f4d33ad2cf9983daab85dca37d8ed41..1ecfc5a4ce6e4106750fad71b9d63ac7e5dd0fc9 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.cpp
@@ -8,6 +8,8 @@
#define PROOF_CORRECTNESS
+using namespace UbMath;
+
//////////////////////////////////////////////////////////////////////////
CumulantLBMKernel::CumulantLBMKernel()
{
@@ -16,11 +18,6 @@ CumulantLBMKernel::CumulantLBMKernel()
this->OxyyMxzz = 1.0;
this->bulkOmegaToOmega = false;
this->OxxPyyPzz = 1.0;
-}
-//////////////////////////////////////////////////////////////////////////
-CumulantLBMKernel::~CumulantLBMKernel(void)
-{
-
}
//////////////////////////////////////////////////////////////////////////
void CumulantLBMKernel::initDataSet()
@@ -60,7 +57,7 @@ SPtr<LBMKernel> CumulantLBMKernel::clone()
}
else
{
- dynamicPointerCast<CumulantLBMKernel>(kernel)->OxxPyyPzz = one;
+ dynamicPointerCast<CumulantLBMKernel>(kernel)->OxxPyyPzz = UbMath::one;
}
return kernel;
}
@@ -1067,10 +1064,6 @@ void CumulantLBMKernel::initData()
muForcingX1.DefineVar("nu", &muNu);
muForcingX2.DefineVar("nu", &muNu);
muForcingX3.DefineVar("nu", &muNu);
-
- LBMReal forcingX1 = 0;
- LBMReal forcingX2 = 0;
- LBMReal forcingX3 = 0;
}
localDistributions = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getLocalDistributions();
nonLocalDistributions = dynamicPointerCast<D3Q27EsoTwist3DSplittedVector>(dataSet->getFdistributions())->getNonLocalDistributions();
@@ -1139,7 +1132,7 @@ void CumulantLBMKernel::nodeCollision(int step, int x1, int x2, int x3)
(((mfbac + mfbca) + (mfbaa + mfbcc)) + ((mfabc + mfcba) + (mfaba + mfcbc)) + ((mfacb + mfcab) + (mfaab + mfccb))) +
((mfabb + mfcbb) + (mfbab + mfbcb)) + (mfbba + mfbbc)) + mfbbb;
- LBMReal rho = one + drho;
+ LBMReal rho = UbMath::one + drho;
////////////////////////////////////////////////////////////////////////////////////
LBMReal vvx = ((((mfccc - mfaaa) + (mfcac - mfaca)) + ((mfcaa - mfacc) + (mfcca - mfaac))) +
(((mfcba - mfabc) + (mfcbc - mfaba)) + ((mfcab - mfacb) + (mfccb - mfaab))) +
@@ -1429,7 +1422,7 @@ void CumulantLBMKernel::nodeCollision(int step, int x1, int x2, int x3)
//////////////////////////////
LBMReal OxyyPxzz = one;//three * (two - omega) / (three - omega);//
//LBMReal OxyyMxzz = one;//six * (two - omega) / (six - omega);//
- LBMReal Oxyz = one;//twelve * (two - omega) / (twelve + omega);//
+ //LBMReal Oxyz = one;//twelve * (two - omega) / (twelve + omega);//
//////////////////////////////
//LBMReal OxyyPxzz = two-omega;//
//LBMReal OxyyMxzz = two-omega;//
diff --git a/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.h
index b75b13031fbfc0c112f2230607cb56928bcba24a..995ce63d877d833e7907d6335e609cfa7478aebd 100644
--- a/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/CumulantLBMKernel.h
@@ -18,10 +18,10 @@ public:
enum Parameter { NORMAL, MAGIC };
public:
CumulantLBMKernel();
- virtual ~CumulantLBMKernel(void);
+ virtual ~CumulantLBMKernel() = default;
//virtual void calculate(int step);
- SPtr<LBMKernel> clone();
- double getCalculationTime();
+ SPtr<LBMKernel> clone() override;
+ double getCalculationTime() override;
void setBulkOmegaToOmega(bool value);
void setRelaxationParameter(Parameter p);
protected:
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
index 6c2f84f06d701dc5cfe137e6a2cedf103a438efc..874a24fbb0e18dcd963048ca3dbe0548131acf56 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.cpp
@@ -1,4 +1,5 @@
#include "D3Q27System.h"
+
namespace D3Q27System
{
using namespace UbMath;
@@ -28,85 +29,85 @@ const int INVDIR[] = { INV_E, INV_W, INV_N, INV_S, INV_T, INV_B, INV
INV_NW, INV_TE, INV_BW, INV_BE, INV_TW, INV_TN, INV_BS, INV_BN, INV_TS,
INV_TNE, INV_TNW, INV_TSE, INV_TSW, INV_BNE, INV_BNW, INV_BSE, INV_BSW };
-// The x,y,z component for each normalized direction
-const double cNorm[3][ENDDIR] = { { double(DX1[0]),
- double(DX1[1]),
- double(DX1[2]),
- double(DX1[3]),
- double(DX1[4]),
- double(DX1[5]),
- double(DX1[6]) / std::sqrt(double(2)),
- double(DX1[7]) / std::sqrt(double(2)),
- double(DX1[8]) / std::sqrt(double(2)),
- double(DX1[9]) / std::sqrt(double(2)),
- double(DX1[10]) / std::sqrt(double(2)),
- double(DX1[11]) / std::sqrt(double(2)),
- double(DX1[12]) / std::sqrt(double(2)),
- double(DX1[13]) / std::sqrt(double(2)),
- double(DX1[14]),
- double(DX1[15]),
- double(DX1[16]),
- double(DX1[17]),
- double(DX1[18]) / std::sqrt(double(3)),
- double(DX1[19]) / std::sqrt(double(3)),
- double(DX1[20]) / std::sqrt(double(3)),
- double(DX1[21]) / std::sqrt(double(3)),
- double(DX1[22]) / std::sqrt(double(3)),
- double(DX1[23]) / std::sqrt(double(3)),
- double(DX1[24]) / std::sqrt(double(3)),
- double(DX1[25]) / std::sqrt(double(3)) },
- { double(DX2[0]),
- double(DX2[1]),
- double(DX2[2]),
- double(DX2[3]),
- double(DX2[4]),
- double(DX2[5]),
- double(DX2[6]) / std::sqrt(double(2)),
- double(DX2[7]) / std::sqrt(double(2)),
- double(DX2[8]) / std::sqrt(double(2)),
- double(DX2[9]) / std::sqrt(double(2)),
- double(DX2[10]),
- double(DX2[11]),
- double(DX2[12]),
- double(DX2[13]),
- double(DX2[14]) / std::sqrt(double(2)),
- double(DX2[15]) / std::sqrt(double(2)),
- double(DX2[16]) / std::sqrt(double(2)),
- double(DX2[17]) / std::sqrt(double(2)),
- double(DX2[18]) / std::sqrt(double(3)),
- double(DX2[19]) / std::sqrt(double(3)),
- double(DX2[20]) / std::sqrt(double(3)),
- double(DX2[21]) / std::sqrt(double(3)),
- double(DX2[22]) / std::sqrt(double(3)),
- double(DX2[23]) / std::sqrt(double(3)),
- double(DX2[24]) / std::sqrt(double(3)),
- double(DX2[25]) / std::sqrt(double(3)) },
- { double(DX3[0]),
- double(DX3[1]),
- double(DX3[2]),
- double(DX3[3]),
- double(DX3[4]),
- double(DX3[5]),
- double(DX3[6]),
- double(DX3[7]),
- double(DX3[8]),
- double(DX3[9]),
- double(DX3[10]) / std::sqrt(double(2)),
- double(DX3[11]) / std::sqrt(double(2)),
- double(DX3[12]) / std::sqrt(double(2)),
- double(DX3[13]) / std::sqrt(double(2)),
- double(DX3[14]) / std::sqrt(double(2)),
- double(DX3[15]) / std::sqrt(double(2)),
- double(DX3[16]) / std::sqrt(double(2)),
- double(DX3[17]) / std::sqrt(double(2)),
- double(DX3[18]) / std::sqrt(double(3)),
- double(DX3[19]) / std::sqrt(double(3)),
- double(DX3[20]) / std::sqrt(double(3)),
- double(DX3[21]) / std::sqrt(double(3)),
- double(DX3[22]) / std::sqrt(double(3)),
- double(DX3[23]) / std::sqrt(double(3)),
- double(DX3[24]) / std::sqrt(double(3)),
- double(DX3[25]) / std::sqrt(double(3)) } };
+//// The x,y,z component for each normalized direction
+//const double cNorm[3][ENDDIR] = { { double(DX1[0]),
+// double(DX1[1]),
+// double(DX1[2]),
+// double(DX1[3]),
+// double(DX1[4]),
+// double(DX1[5]),
+// double(DX1[6]) / std::sqrt(double(2)),
+// double(DX1[7]) / std::sqrt(double(2)),
+// double(DX1[8]) / std::sqrt(double(2)),
+// double(DX1[9]) / std::sqrt(double(2)),
+// double(DX1[10]) / std::sqrt(double(2)),
+// double(DX1[11]) / std::sqrt(double(2)),
+// double(DX1[12]) / std::sqrt(double(2)),
+// double(DX1[13]) / std::sqrt(double(2)),
+// double(DX1[14]),
+// double(DX1[15]),
+// double(DX1[16]),
+// double(DX1[17]),
+// double(DX1[18]) / std::sqrt(double(3)),
+// double(DX1[19]) / std::sqrt(double(3)),
+// double(DX1[20]) / std::sqrt(double(3)),
+// double(DX1[21]) / std::sqrt(double(3)),
+// double(DX1[22]) / std::sqrt(double(3)),
+// double(DX1[23]) / std::sqrt(double(3)),
+// double(DX1[24]) / std::sqrt(double(3)),
+// double(DX1[25]) / std::sqrt(double(3)) },
+// { double(DX2[0]),
+// double(DX2[1]),
+// double(DX2[2]),
+// double(DX2[3]),
+// double(DX2[4]),
+// double(DX2[5]),
+// double(DX2[6]) / std::sqrt(double(2)),
+// double(DX2[7]) / std::sqrt(double(2)),
+// double(DX2[8]) / std::sqrt(double(2)),
+// double(DX2[9]) / std::sqrt(double(2)),
+// double(DX2[10]),
+// double(DX2[11]),
+// double(DX2[12]),
+// double(DX2[13]),
+// double(DX2[14]) / std::sqrt(double(2)),
+// double(DX2[15]) / std::sqrt(double(2)),
+// double(DX2[16]) / std::sqrt(double(2)),
+// double(DX2[17]) / std::sqrt(double(2)),
+// double(DX2[18]) / std::sqrt(double(3)),
+// double(DX2[19]) / std::sqrt(double(3)),
+// double(DX2[20]) / std::sqrt(double(3)),
+// double(DX2[21]) / std::sqrt(double(3)),
+// double(DX2[22]) / std::sqrt(double(3)),
+// double(DX2[23]) / std::sqrt(double(3)),
+// double(DX2[24]) / std::sqrt(double(3)),
+// double(DX2[25]) / std::sqrt(double(3)) },
+// { double(DX3[0]),
+// double(DX3[1]),
+// double(DX3[2]),
+// double(DX3[3]),
+// double(DX3[4]),
+// double(DX3[5]),
+// double(DX3[6]),
+// double(DX3[7]),
+// double(DX3[8]),
+// double(DX3[9]),
+// double(DX3[10]) / std::sqrt(double(2)),
+// double(DX3[11]) / std::sqrt(double(2)),
+// double(DX3[12]) / std::sqrt(double(2)),
+// double(DX3[13]) / std::sqrt(double(2)),
+// double(DX3[14]) / std::sqrt(double(2)),
+// double(DX3[15]) / std::sqrt(double(2)),
+// double(DX3[16]) / std::sqrt(double(2)),
+// double(DX3[17]) / std::sqrt(double(2)),
+// double(DX3[18]) / std::sqrt(double(3)),
+// double(DX3[19]) / std::sqrt(double(3)),
+// double(DX3[20]) / std::sqrt(double(3)),
+// double(DX3[21]) / std::sqrt(double(3)),
+// double(DX3[22]) / std::sqrt(double(3)),
+// double(DX3[23]) / std::sqrt(double(3)),
+// double(DX3[24]) / std::sqrt(double(3)),
+// double(DX3[25]) / std::sqrt(double(3)) } };
} // namespace D3Q27System
diff --git a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
index a15a7a82ba59f851df5a345989b607b8779e9ac2..7cd494823ca9d5335a3029609f5c3a83fd60b523 100644
--- a/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
+++ b/src/cpu/VirtualFluidsCore/LBM/D3Q27System.h
@@ -67,11 +67,6 @@ class for global system-functions
usage: ...
*/
-
-#ifndef SWIG
- using namespace UbMath;
-#endif
-
namespace D3Q27System
{
//enum COLLISIONMODEL { UNDEFINED, INCOMPLBGKMODEL, COMPLBGKMODEL, COMPLBGKWTMODEL, INCOMPLBGKLESMODEL, INCOMPLBGKNONNEWTONIANMODEL
@@ -173,8 +168,6 @@ usage: ...
extern const int DX3[ENDDIR + 1];
extern const double WEIGTH[ENDDIR + 1];
- extern const double cNorm[3][ENDDIR];
-
//static const int ZERO /*f0 */ = 0;
//static const int E /*f1 */ = 1;
//static const int W /*f2 */ = 2;
@@ -287,30 +280,6 @@ usage: ...
static const int ET_TSW = 12;
static const int ET_BNE = 12;
- static const int M_RHO = 0;
- static const int M_EN = 1;
- static const int M_EPS = 2;
- static const int M_JX1 = 3;
- static const int M_QX1 = 4;
- static const int M_JX2 = 5;
- static const int M_QX2 = 6;
- static const int M_JX3 = 7;
- static const int M_QX3 = 8;
- static const int M_3PX1X1 = 9;
- static const int M_3PIX1X1 = 10;
- static const int M_PWW = 11;
- static const int M_PIWW = 12;
- static const int M_PX1X2 = 13;
- static const int M_PX2X3 = 14;
- static const int M_PX1X3 = 15;
- static const int M_MX1 = 16;
- static const int M_MX2 = 17;
- static const int M_MX3 = 18;
-
- static const int STARTM = 0;
- static const int ENDM = 18; //D3Q27
-
-
//////////////////////////////////////////////////////////////////////////
//MACROSCOPIC VALUES
@@ -326,7 +295,7 @@ usage: ...
//ACHTUNG: gilt nicht fuer alle modelle -> praedikat verwenden anstelle static! toDo
static LBMReal getPressure(const LBMReal* const& f/*[27]*/)
{
- return REAL_CAST(c1o3) * getDensity(f);
+ return REAL_CAST(UbMath::c1o3) * getDensity(f);
}
/*=====================================================================*/
static LBMReal getIncompVelocityX1(const LBMReal* const& f/*[27]*/)
@@ -437,7 +406,7 @@ usage: ...
D3Q27System::calcIncompVelocityX1(f, vx1);
D3Q27System::calcIncompVelocityX2(f, vx2);
D3Q27System::calcIncompVelocityX3(f, vx3);
- LBMReal rho = drho + one;
+ LBMReal rho = drho + UbMath::one;
vx1 /= rho;
vx2 /= rho;
vx3 /= rho;
@@ -449,68 +418,68 @@ usage: ...
//switch(direction)
//{
- // case ZERO : return REAL_CAST( c8o27*rho*(1.0-cu_sq));
- // case E : return REAL_CAST( c2o27*rho*(1.0+3.0*( vx1 )+c9o2*( vx1 )*( vx1 )-cu_sq));
- // case W : return REAL_CAST( c2o27*rho*(1.0+3.0*(-vx1 )+c9o2*(-vx1 )*(-vx1 )-cu_sq));
- // case N : return REAL_CAST( c2o27*rho*(1.0+3.0*( vx2)+c9o2*( vx2)*( vx2)-cu_sq));
- // case S : return REAL_CAST( c2o27*rho*(1.0+3.0*( -vx2)+c9o2*( -vx2)*( -vx2)-cu_sq));
- // case T : return REAL_CAST( c2o27*rho*(1.0+3.0*( vx3 )+c9o2*( vx3)*( vx3)-cu_sq));
- // case B : return REAL_CAST( c2o27*rho*(1.0+3.0*( -vx3)+c9o2*( -vx3)*( -vx3)-cu_sq));
- // case NE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1+vx2)+c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq));
- // case SW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1-vx2)+c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq));
- // case SE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1-vx2)+c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq));
- // case NW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1+vx2)+c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq));
- // case TE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1+vx3)+c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq));
- // case BW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1-vx3)+c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq));
- // case BE : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx1-vx3)+c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq));
- // case TW : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx1+vx3)+c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq));
- // case TN : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx2+vx3)+c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq));
- // case BS : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx2-vx3)+c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq));
- // case BN : return REAL_CAST( c1o54*rho*(1.0+3.0*( vx2-vx3)+c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq));
- // case TS : return REAL_CAST( c1o54*rho*(1.0+3.0*(-vx2+vx3)+c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq));
- // case TNE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq));
- // case BSW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq));
- // case BNE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq));
- // case TSW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq));
- // case TSE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq));
- // case BNW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq));
- // case BSE : return REAL_CAST(c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq));
- // case TNW : return REAL_CAST(c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq));
+ // case ZERO : return REAL_CAST( UbMath::c8o27*rho*(1.0-cu_sq));
+ // case E : return REAL_CAST( UbMath::c2o27*rho*(1.0+3.0*( vx1 )+UbMath::c9o2*( vx1 )*( vx1 )-cu_sq));
+ // case W : return REAL_CAST( UbMath::c2o27*rho*(1.0+3.0*(-vx1 )+UbMath::c9o2*(-vx1 )*(-vx1 )-cu_sq));
+ // case N : return REAL_CAST( UbMath::c2o27*rho*(1.0+3.0*( vx2)+UbMath::c9o2*( vx2)*( vx2)-cu_sq));
+ // case S : return REAL_CAST( UbMath::c2o27*rho*(1.0+3.0*( -vx2)+UbMath::c9o2*( -vx2)*( -vx2)-cu_sq));
+ // case T : return REAL_CAST( UbMath::c2o27*rho*(1.0+3.0*( vx3 )+UbMath::c9o2*( vx3)*( vx3)-cu_sq));
+ // case B : return REAL_CAST( UbMath::c2o27*rho*(1.0+3.0*( -vx3)+UbMath::c9o2*( -vx3)*( -vx3)-cu_sq));
+ // case NE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1+vx2)+UbMath::c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq));
+ // case SW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx2)+UbMath::c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq));
+ // case SE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1-vx2)+UbMath::c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq));
+ // case NW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx2)+UbMath::c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq));
+ // case TE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1+vx3)+UbMath::c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq));
+ // case BW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx3)+UbMath::c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq));
+ // case BE : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx1-vx3)+UbMath::c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq));
+ // case TW : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx3)+UbMath::c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq));
+ // case TN : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx2+vx3)+UbMath::c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq));
+ // case BS : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx2-vx3)+UbMath::c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq));
+ // case BN : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*( vx2-vx3)+UbMath::c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq));
+ // case TS : return REAL_CAST( UbMath::c1o54*rho*(1.0+3.0*(-vx2+vx3)+UbMath::c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq));
+ // case TNE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+UbMath::c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq));
+ // case BSW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+UbMath::c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq));
+ // case BNE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+UbMath::c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq));
+ // case TSW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+UbMath::c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq));
+ // case TSE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+UbMath::c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq));
+ // case BNW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+UbMath::c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq));
+ // case BSE : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+UbMath::c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq));
+ // case TNW : return REAL_CAST(UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+UbMath::c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq));
// default: throw UbException(UB_EXARGS,"unknown dir");
//}
////-----
- LBMReal rho = drho + one;
+ LBMReal rho = drho + UbMath::one;
switch (direction)
{
- case ZERO: return REAL_CAST(c8o27 * (drho + rho * (-cu_sq)));
- case E: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq)));
- case W: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq)));
- case N: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq)));
- case S: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq)));
- case T: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (vx3)+c9o2 * (vx3) * (vx3)-cu_sq)));
- case B: return REAL_CAST(c2o27 * (drho + rho * (3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq)));
- case NE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq)));
- case SW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq)));
- case SE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq)));
- case NW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq)));
- case TE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq)));
- case BW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq)));
- case BE: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq)));
- case TW: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq)));
- case TN: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq)));
- case BS: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq)));
- case BN: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq)));
- case TS: return REAL_CAST(c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq)));
- case TNE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq)));
- case BSW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq)));
- case BNE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq)));
- case TSW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq)));
- case TSE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq)));
- case BNW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq)));
- case BSE: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq)));
- case TNW: return REAL_CAST(c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq)));
+ case ZERO: return REAL_CAST(UbMath::c8o27 * (drho + rho * (-cu_sq)));
+ case E: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq)));
+ case W: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq)));
+ case N: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq)));
+ case S: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq)));
+ case T: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq)));
+ case B: return REAL_CAST(UbMath::c2o27 * (drho + rho * (3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq)));
+ case NE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq)));
+ case SW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq)));
+ case SE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq)));
+ case NW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq)));
+ case TE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq)));
+ case BW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq)));
+ case BE: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq)));
+ case TW: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq)));
+ case TN: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq)));
+ case BS: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq)));
+ case BN: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq)));
+ case TS: return REAL_CAST(UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq)));
+ case TNE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq)));
+ case BSW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq)));
+ case BNE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq)));
+ case TSW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq)));
+ case TSE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq)));
+ case BNW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq)));
+ case BSE: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq)));
+ case TNW: return REAL_CAST(UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq)));
default: throw UbException(UB_EXARGS, "unknown dir");
}
@@ -520,66 +489,66 @@ usage: ...
{
//LBMReal cu_sq=1.5*(vx1*vx1+vx2*vx2+vx3*vx3);
- //feq[ZERO] = c8o27*rho*(1.0-cu_sq);
- //feq[E] = c2o27*rho*(1.0+3.0*( vx1 )+c9o2*( vx1 )*( vx1 )-cu_sq);
- //feq[W] = c2o27*rho*(1.0+3.0*(-vx1 )+c9o2*(-vx1 )*(-vx1 )-cu_sq);
- //feq[N] = c2o27*rho*(1.0+3.0*( vx2)+c9o2*( vx2)*( vx2)-cu_sq);
- //feq[S] = c2o27*rho*(1.0+3.0*( -vx2)+c9o2*( -vx2)*( -vx2)-cu_sq);
- //feq[T] = c2o27*rho*(1.0+3.0*( vx3 )+c9o2*( vx3)*( vx3)-cu_sq);
- //feq[B] = c2o27*rho*(1.0+3.0*( -vx3)+c9o2*( -vx3)*( -vx3)-cu_sq);
- //feq[NE] = c1o54*rho*(1.0+3.0*( vx1+vx2)+c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq);
- //feq[SW] = c1o54*rho*(1.0+3.0*(-vx1-vx2)+c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq);
- //feq[SE] = c1o54*rho*(1.0+3.0*( vx1-vx2)+c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq);
- //feq[NW] = c1o54*rho*(1.0+3.0*(-vx1+vx2)+c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq);
- //feq[TE] = c1o54*rho*(1.0+3.0*( vx1+vx3)+c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq);
- //feq[BW] = c1o54*rho*(1.0+3.0*(-vx1-vx3)+c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq);
- //feq[BE] = c1o54*rho*(1.0+3.0*( vx1-vx3)+c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq);
- //feq[TW] = c1o54*rho*(1.0+3.0*(-vx1+vx3)+c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq);
- //feq[TN] = c1o54*rho*(1.0+3.0*( vx2+vx3)+c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq);
- //feq[BS] = c1o54*rho*(1.0+3.0*(-vx2-vx3)+c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq);
- //feq[BN] = c1o54*rho*(1.0+3.0*( vx2-vx3)+c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq);
- //feq[TS] = c1o54*rho*(1.0+3.0*(-vx2+vx3)+c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq);
- //feq[TNE] = c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
- //feq[BSW] = c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
- //feq[BNE] = c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
- //feq[TSW] = c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
- //feq[TSE] = c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
- //feq[BNW] = c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
- //feq[BSE] = c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
- //feq[TNW] = c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
+ //feq[ZERO] = UbMath::c8o27*rho*(1.0-cu_sq);
+ //feq[E] = UbMath::c2o27*rho*(1.0+3.0*( vx1 )+UbMath::c9o2*( vx1 )*( vx1 )-cu_sq);
+ //feq[W] = UbMath::c2o27*rho*(1.0+3.0*(-vx1 )+UbMath::c9o2*(-vx1 )*(-vx1 )-cu_sq);
+ //feq[N] = UbMath::c2o27*rho*(1.0+3.0*( vx2)+UbMath::c9o2*( vx2)*( vx2)-cu_sq);
+ //feq[S] = UbMath::c2o27*rho*(1.0+3.0*( -vx2)+UbMath::c9o2*( -vx2)*( -vx2)-cu_sq);
+ //feq[T] = UbMath::c2o27*rho*(1.0+3.0*( vx3 )+UbMath::c9o2*( vx3)*( vx3)-cu_sq);
+ //feq[B] = UbMath::c2o27*rho*(1.0+3.0*( -vx3)+UbMath::c9o2*( -vx3)*( -vx3)-cu_sq);
+ //feq[NE] = UbMath::c1o54*rho*(1.0+3.0*( vx1+vx2)+UbMath::c9o2*( vx1+vx2)*( vx1+vx2)-cu_sq);
+ //feq[SW] = UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx2)+UbMath::c9o2*(-vx1-vx2)*(-vx1-vx2)-cu_sq);
+ //feq[SE] = UbMath::c1o54*rho*(1.0+3.0*( vx1-vx2)+UbMath::c9o2*( vx1-vx2)*( vx1-vx2)-cu_sq);
+ //feq[NW] = UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx2)+UbMath::c9o2*(-vx1+vx2)*(-vx1+vx2)-cu_sq);
+ //feq[TE] = UbMath::c1o54*rho*(1.0+3.0*( vx1+vx3)+UbMath::c9o2*( vx1+vx3)*( vx1+vx3)-cu_sq);
+ //feq[BW] = UbMath::c1o54*rho*(1.0+3.0*(-vx1-vx3)+UbMath::c9o2*(-vx1-vx3)*(-vx1-vx3)-cu_sq);
+ //feq[BE] = UbMath::c1o54*rho*(1.0+3.0*( vx1-vx3)+UbMath::c9o2*( vx1-vx3)*( vx1-vx3)-cu_sq);
+ //feq[TW] = UbMath::c1o54*rho*(1.0+3.0*(-vx1+vx3)+UbMath::c9o2*(-vx1+vx3)*(-vx1+vx3)-cu_sq);
+ //feq[TN] = UbMath::c1o54*rho*(1.0+3.0*( vx2+vx3)+UbMath::c9o2*( vx2+vx3)*( vx2+vx3)-cu_sq);
+ //feq[BS] = UbMath::c1o54*rho*(1.0+3.0*(-vx2-vx3)+UbMath::c9o2*(-vx2-vx3)*(-vx2-vx3)-cu_sq);
+ //feq[BN] = UbMath::c1o54*rho*(1.0+3.0*( vx2-vx3)+UbMath::c9o2*( vx2-vx3)*( vx2-vx3)-cu_sq);
+ //feq[TS] = UbMath::c1o54*rho*(1.0+3.0*(-vx2+vx3)+UbMath::c9o2*(-vx2+vx3)*(-vx2+vx3)-cu_sq);
+ //feq[TNE] = UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2+vx3)+UbMath::c9o2*( vx1+vx2+vx3)*( vx1+vx2+vx3)-cu_sq);
+ //feq[BSW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2-vx3)+UbMath::c9o2*(-vx1-vx2-vx3)*(-vx1-vx2-vx3)-cu_sq);
+ //feq[BNE] = UbMath::c1o216*rho*(1.0+3.0*( vx1+vx2-vx3)+UbMath::c9o2*( vx1+vx2-vx3)*( vx1+vx2-vx3)-cu_sq);
+ //feq[TSW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1-vx2+vx3)+UbMath::c9o2*(-vx1-vx2+vx3)*(-vx1-vx2+vx3)-cu_sq);
+ //feq[TSE] = UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2+vx3)+UbMath::c9o2*( vx1-vx2+vx3)*( vx1-vx2+vx3)-cu_sq);
+ //feq[BNW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2-vx3)+UbMath::c9o2*(-vx1+vx2-vx3)*(-vx1+vx2-vx3)-cu_sq);
+ //feq[BSE] = UbMath::c1o216*rho*(1.0+3.0*( vx1-vx2-vx3)+UbMath::c9o2*( vx1-vx2-vx3)*( vx1-vx2-vx3)-cu_sq);
+ //feq[TNW] = UbMath::c1o216*rho*(1.0+3.0*(-vx1+vx2+vx3)+UbMath::c9o2*(-vx1+vx2+vx3)*(-vx1+vx2+vx3)-cu_sq);
//////////////////////////////////////////////////////////////////////////
LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
- LBMReal rho = drho + one;
-
- feq[ZERO] = c8o27 * (drho + rho * (-cu_sq));
- feq[E] = c2o27 * (drho + rho * (3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq));
- feq[W] = c2o27 * (drho + rho * (3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq));
- feq[N] = c2o27 * (drho + rho * (3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq));
- feq[S] = c2o27 * (drho + rho * (3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq));
- feq[T] = c2o27 * (drho + rho * (3.0 * (vx3)+c9o2 * (vx3) * (vx3)-cu_sq));
- feq[B] = c2o27 * (drho + rho * (3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq));
- feq[NE] = c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
- feq[SW] = c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
- feq[SE] = c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
- feq[NW] = c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
- feq[TE] = c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
- feq[BW] = c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
- feq[BE] = c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
- feq[TW] = c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
- feq[TN] = c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
- feq[BS] = c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
- feq[BN] = c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
- feq[TS] = c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
- feq[TNE] = c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
- feq[BSW] = c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
- feq[BNE] = c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
- feq[TSW] = c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
- feq[TSE] = c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
- feq[BNW] = c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
- feq[BSE] = c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
- feq[TNW] = c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
+ LBMReal rho = drho + UbMath::one;
+
+ feq[ZERO] = UbMath::c8o27 * (drho + rho * (-cu_sq));
+ feq[E] = UbMath::c2o27 * (drho + rho * (3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq));
+ feq[W] = UbMath::c2o27 * (drho + rho * (3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq));
+ feq[N] = UbMath::c2o27 * (drho + rho * (3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq));
+ feq[S] = UbMath::c2o27 * (drho + rho * (3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq));
+ feq[T] = UbMath::c2o27 * (drho + rho * (3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq));
+ feq[B] = UbMath::c2o27 * (drho + rho * (3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq));
+ feq[NE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
+ feq[SW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
+ feq[SE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
+ feq[NW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
+ feq[TE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
+ feq[BW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
+ feq[BE] = UbMath::c1o54 * (drho + rho * (3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
+ feq[TW] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
+ feq[TN] = UbMath::c1o54 * (drho + rho * (3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
+ feq[BS] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
+ feq[BN] = UbMath::c1o54 * (drho + rho * (3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
+ feq[TS] = UbMath::c1o54 * (drho + rho * (3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
+ feq[TNE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
+ feq[BSW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
+ feq[BNE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
+ feq[TSW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
+ feq[TSE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
+ feq[BNW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
+ feq[BSE] = UbMath::c1o216 * (drho + rho * (3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
+ feq[TNW] = UbMath::c1o216 * (drho + rho * (3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
}
//////////////////////////////////////////////////////////////////////////
static LBMReal getIncompFeqForDirection(const int& direction, const LBMReal& drho, const LBMReal& vx1, const LBMReal& vx2, const LBMReal& vx3)
@@ -588,33 +557,33 @@ usage: ...
switch (direction)
{
- case ZERO: return REAL_CAST(c8o27 * (drho - cu_sq));
- case E: return REAL_CAST(c2o27 * (drho + 3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq));
- case W: return REAL_CAST(c2o27 * (drho + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq));
- case N: return REAL_CAST(c2o27 * (drho + 3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq));
- case S: return REAL_CAST(c2o27 * (drho + 3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq));
- case T: return REAL_CAST(c2o27 * (drho + 3.0 * (vx3)+c9o2 * (vx3) * (vx3)-cu_sq));
- case B: return REAL_CAST(c2o27 * (drho + 3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq));
- case NE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
- case SW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
- case SE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
- case NW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
- case TE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
- case BW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
- case BE: return REAL_CAST(c1o54 * (drho + 3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
- case TW: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
- case TN: return REAL_CAST(c1o54 * (drho + 3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
- case BS: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
- case BN: return REAL_CAST(c1o54 * (drho + 3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
- case TS: return REAL_CAST(c1o54 * (drho + 3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
- case TNE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
- case BSW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
- case BNE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
- case TSW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
- case TSE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
- case BNW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
- case BSE: return REAL_CAST(c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
- case TNW: return REAL_CAST(c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
+ case ZERO: return REAL_CAST(UbMath::c8o27 * (drho - cu_sq));
+ case E: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq));
+ case W: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq));
+ case N: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq));
+ case S: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq));
+ case T: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq));
+ case B: return REAL_CAST(UbMath::c2o27 * (drho + 3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq));
+ case NE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq));
+ case SW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq));
+ case SE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq));
+ case NW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq));
+ case TE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq));
+ case BW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq));
+ case BE: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq));
+ case TW: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq));
+ case TN: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq));
+ case BS: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq));
+ case BN: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq));
+ case TS: return REAL_CAST(UbMath::c1o54 * (drho + 3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq));
+ case TNE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq));
+ case BSW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq));
+ case BNE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq));
+ case TSW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq));
+ case TSE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq));
+ case BNW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq));
+ case BSE: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq));
+ case TNW: return REAL_CAST(UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq));
default: throw UbException(UB_EXARGS, "unknown dir");
}
}
@@ -623,33 +592,33 @@ usage: ...
{
LBMReal cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
- feq[ZERO] = c8o27 * (drho - cu_sq);
- feq[E] = c2o27 * (drho + 3.0 * (vx1)+c9o2 * (vx1) * (vx1)-cu_sq);
- feq[W] = c2o27 * (drho + 3.0 * (-vx1) + c9o2 * (-vx1) * (-vx1) - cu_sq);
- feq[N] = c2o27 * (drho + 3.0 * (vx2)+c9o2 * (vx2) * (vx2)-cu_sq);
- feq[S] = c2o27 * (drho + 3.0 * (-vx2) + c9o2 * (-vx2) * (-vx2) - cu_sq);
- feq[T] = c2o27 * (drho + 3.0 * (vx3)+c9o2 * (vx3) * (vx3)-cu_sq);
- feq[B] = c2o27 * (drho + 3.0 * (-vx3) + c9o2 * (-vx3) * (-vx3) - cu_sq);
- feq[NE] = c1o54 * (drho + 3.0 * (vx1 + vx2) + c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq);
- feq[SW] = c1o54 * (drho + 3.0 * (-vx1 - vx2) + c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq);
- feq[SE] = c1o54 * (drho + 3.0 * (vx1 - vx2) + c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq);
- feq[NW] = c1o54 * (drho + 3.0 * (-vx1 + vx2) + c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq);
- feq[TE] = c1o54 * (drho + 3.0 * (vx1 + vx3) + c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq);
- feq[BW] = c1o54 * (drho + 3.0 * (-vx1 - vx3) + c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq);
- feq[BE] = c1o54 * (drho + 3.0 * (vx1 - vx3) + c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq);
- feq[TW] = c1o54 * (drho + 3.0 * (-vx1 + vx3) + c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq);
- feq[TN] = c1o54 * (drho + 3.0 * (vx2 + vx3) + c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq);
- feq[BS] = c1o54 * (drho + 3.0 * (-vx2 - vx3) + c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq);
- feq[BN] = c1o54 * (drho + 3.0 * (vx2 - vx3) + c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq);
- feq[TS] = c1o54 * (drho + 3.0 * (-vx2 + vx3) + c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq);
- feq[TNE] = c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq);
- feq[BSW] = c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq);
- feq[BNE] = c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq);
- feq[TSW] = c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq);
- feq[TSE] = c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq);
- feq[BNW] = c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq);
- feq[BSE] = c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq);
- feq[TNW] = c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
+ feq[ZERO] = UbMath::c8o27 * (drho - cu_sq);
+ feq[E] = UbMath::c2o27 * (drho + 3.0 * (vx1)+UbMath::c9o2 * (vx1) * (vx1)-cu_sq);
+ feq[W] = UbMath::c2o27 * (drho + 3.0 * (-vx1) + UbMath::c9o2 * (-vx1) * (-vx1) - cu_sq);
+ feq[N] = UbMath::c2o27 * (drho + 3.0 * (vx2)+UbMath::c9o2 * (vx2) * (vx2)-cu_sq);
+ feq[S] = UbMath::c2o27 * (drho + 3.0 * (-vx2) + UbMath::c9o2 * (-vx2) * (-vx2) - cu_sq);
+ feq[T] = UbMath::c2o27 * (drho + 3.0 * (vx3)+UbMath::c9o2 * (vx3) * (vx3)-cu_sq);
+ feq[B] = UbMath::c2o27 * (drho + 3.0 * (-vx3) + UbMath::c9o2 * (-vx3) * (-vx3) - cu_sq);
+ feq[NE] = UbMath::c1o54 * (drho + 3.0 * (vx1 + vx2) + UbMath::c9o2 * (vx1 + vx2) * (vx1 + vx2) - cu_sq);
+ feq[SW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx2) + UbMath::c9o2 * (-vx1 - vx2) * (-vx1 - vx2) - cu_sq);
+ feq[SE] = UbMath::c1o54 * (drho + 3.0 * (vx1 - vx2) + UbMath::c9o2 * (vx1 - vx2) * (vx1 - vx2) - cu_sq);
+ feq[NW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx2) + UbMath::c9o2 * (-vx1 + vx2) * (-vx1 + vx2) - cu_sq);
+ feq[TE] = UbMath::c1o54 * (drho + 3.0 * (vx1 + vx3) + UbMath::c9o2 * (vx1 + vx3) * (vx1 + vx3) - cu_sq);
+ feq[BW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 - vx3) + UbMath::c9o2 * (-vx1 - vx3) * (-vx1 - vx3) - cu_sq);
+ feq[BE] = UbMath::c1o54 * (drho + 3.0 * (vx1 - vx3) + UbMath::c9o2 * (vx1 - vx3) * (vx1 - vx3) - cu_sq);
+ feq[TW] = UbMath::c1o54 * (drho + 3.0 * (-vx1 + vx3) + UbMath::c9o2 * (-vx1 + vx3) * (-vx1 + vx3) - cu_sq);
+ feq[TN] = UbMath::c1o54 * (drho + 3.0 * (vx2 + vx3) + UbMath::c9o2 * (vx2 + vx3) * (vx2 + vx3) - cu_sq);
+ feq[BS] = UbMath::c1o54 * (drho + 3.0 * (-vx2 - vx3) + UbMath::c9o2 * (-vx2 - vx3) * (-vx2 - vx3) - cu_sq);
+ feq[BN] = UbMath::c1o54 * (drho + 3.0 * (vx2 - vx3) + UbMath::c9o2 * (vx2 - vx3) * (vx2 - vx3) - cu_sq);
+ feq[TS] = UbMath::c1o54 * (drho + 3.0 * (-vx2 + vx3) + UbMath::c9o2 * (-vx2 + vx3) * (-vx2 + vx3) - cu_sq);
+ feq[TNE] = UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 + vx3) + UbMath::c9o2 * (vx1 + vx2 + vx3) * (vx1 + vx2 + vx3) - cu_sq);
+ feq[BSW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 - vx3) + UbMath::c9o2 * (-vx1 - vx2 - vx3) * (-vx1 - vx2 - vx3) - cu_sq);
+ feq[BNE] = UbMath::c1o216 * (drho + 3.0 * (vx1 + vx2 - vx3) + UbMath::c9o2 * (vx1 + vx2 - vx3) * (vx1 + vx2 - vx3) - cu_sq);
+ feq[TSW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 - vx2 + vx3) + UbMath::c9o2 * (-vx1 - vx2 + vx3) * (-vx1 - vx2 + vx3) - cu_sq);
+ feq[TSE] = UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 + vx3) + UbMath::c9o2 * (vx1 - vx2 + vx3) * (vx1 - vx2 + vx3) - cu_sq);
+ feq[BNW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 - vx3) + UbMath::c9o2 * (-vx1 + vx2 - vx3) * (-vx1 + vx2 - vx3) - cu_sq);
+ feq[BSE] = UbMath::c1o216 * (drho + 3.0 * (vx1 - vx2 - vx3) + UbMath::c9o2 * (vx1 - vx2 - vx3) * (vx1 - vx2 - vx3) - cu_sq);
+ feq[TNW] = UbMath::c1o216 * (drho + 3.0 * (-vx1 + vx2 + vx3) + UbMath::c9o2 * (-vx1 + vx2 + vx3) * (-vx1 + vx2 + vx3) - cu_sq);
}
//////////////////////////////////////////////////////////////////////////
static inline float getBoundaryVelocityForDirection(const int& direction, const float& bcVelocityX1, const float& bcVelocityX2, const float& bcVelocityX3)
@@ -737,22 +706,19 @@ usage: ...
static inline void calcDistanceToNeighbors(std::vector<double>& distNeigh, const double& deltaX1)
{
//distNeigh.resize(FENDDIR+1, UbMath::sqrt2*deltaX1);
- double sqrt3 = UbMath::sqrt3;
- double sqrt2 = UbMath::sqrt2;
+
distNeigh[E] = distNeigh[W] = distNeigh[N] = deltaX1;
distNeigh[S] = distNeigh[T] = distNeigh[B] = deltaX1;
- distNeigh[NE] = distNeigh[NW] = distNeigh[SW] = distNeigh[SE] = sqrt2 * deltaX1;
- distNeigh[TE] = distNeigh[TN] = distNeigh[TW] = distNeigh[TS] = sqrt2 * deltaX1;
- distNeigh[BE] = distNeigh[BN] = distNeigh[BW] = distNeigh[BS] = sqrt2 * deltaX1;
- distNeigh[TNE] = distNeigh[TNW] = distNeigh[TSE] = distNeigh[TSW] = sqrt3 * deltaX1;
- distNeigh[BNE] = distNeigh[BNW] = distNeigh[BSE] = distNeigh[BSW] = sqrt3 * deltaX1;
+ distNeigh[NE] = distNeigh[NW] = distNeigh[SW] = distNeigh[SE] = UbMath::sqrt2 * deltaX1;
+ distNeigh[TE] = distNeigh[TN] = distNeigh[TW] = distNeigh[TS] = UbMath::sqrt2 * deltaX1;
+ distNeigh[BE] = distNeigh[BN] = distNeigh[BW] = distNeigh[BS] = UbMath::sqrt2 * deltaX1;
+ distNeigh[TNE] = distNeigh[TNW] = distNeigh[TSE] = distNeigh[TSW] = UbMath::sqrt3 * deltaX1;
+ distNeigh[BNE] = distNeigh[BNW] = distNeigh[BSE] = distNeigh[BSW] = UbMath::sqrt3 * deltaX1;
}
//////////////////////////////////////////////////////////////////////////
static inline void calcDistanceToNeighbors(std::vector<double>& distNeigh, const double& deltaX1, const double& deltaX2, const double& deltaX3)
{
//distNeigh.resize(FENDDIR+1, UbMath::sqrt2*deltaX1);
- double sqrt3 = UbMath::sqrt3;
- double sqrt2 = UbMath::sqrt2;
distNeigh[E] = distNeigh[W] = deltaX1;
distNeigh[N] = distNeigh[S] = deltaX2;
distNeigh[T] = distNeigh[B] = deltaX3;
@@ -800,7 +766,7 @@ usage: ...
{
LBMReal op = 1.0;
return ((f[E] + f[W] + f[N] + f[S] + f[T] + f[B] + 2. * (f[NE] + f[SW] + f[SE] + f[NW] + f[TE] + f[BW] + f[BE] + f[TW] + f[TN] + f[BS] + f[BN] + f[TS]) +
- 3. * (f[TNE] + f[TSW] + f[TSE] + f[TNW] + f[BNE] + f[BSW] + f[BSE] + f[BNW]) - (vx1 * vx1 + vx2 * vx2 + vx3 * vx3)) * (1 - 0.5 * op) + op * 0.5 * (rho)) * c1o3;
+ 3. * (f[TNE] + f[TSW] + f[TSE] + f[TNW] + f[BNE] + f[BSW] + f[BSE] + f[BNW]) - (vx1 * vx1 + vx2 * vx2 + vx3 * vx3)) * (1 - 0.5 * op) + op * 0.5 * (rho)) * UbMath::c1o3;
}
//////////////////////////////////////////////////////////////////////////
@@ -889,8 +855,6 @@ usage: ...
vy2 = vvy * vvy;
vz2 = vvz * vvz;
////////////////////////////////////////////////////////////////////////////////////
- LBMReal qudricLimit = 0.01;
- ////////////////////////////////////////////////////////////////////////////////////
//Hin
////////////////////////////////////////////////////////////////////////////////////
// mit 1/36, 1/9, 1/36, 1/9, 4/9, 1/9, 1/36, 1/9, 1/36 Konditionieren
@@ -900,7 +864,7 @@ usage: ...
m1 = mfaac - mfaaa;
m0 = m2 + mfaab;
mfaaa = m0;
- m0 += c1o36 * oMdrho;
+ m0 += UbMath::c1o36 * oMdrho;
mfaab = m1 - m0 * vvz;
mfaac = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -908,7 +872,7 @@ usage: ...
m1 = mfabc - mfaba;
m0 = m2 + mfabb;
mfaba = m0;
- m0 += c1o9 * oMdrho;
+ m0 += UbMath::c1o9 * oMdrho;
mfabb = m1 - m0 * vvz;
mfabc = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -916,7 +880,7 @@ usage: ...
m1 = mfacc - mfaca;
m0 = m2 + mfacb;
mfaca = m0;
- m0 += c1o36 * oMdrho;
+ m0 += UbMath::c1o36 * oMdrho;
mfacb = m1 - m0 * vvz;
mfacc = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -925,7 +889,7 @@ usage: ...
m1 = mfbac - mfbaa;
m0 = m2 + mfbab;
mfbaa = m0;
- m0 += c1o9 * oMdrho;
+ m0 += UbMath::c1o9 * oMdrho;
mfbab = m1 - m0 * vvz;
mfbac = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -933,7 +897,7 @@ usage: ...
m1 = mfbbc - mfbba;
m0 = m2 + mfbbb;
mfbba = m0;
- m0 += c4o9 * oMdrho;
+ m0 += UbMath::c4o9 * oMdrho;
mfbbb = m1 - m0 * vvz;
mfbbc = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -941,7 +905,7 @@ usage: ...
m1 = mfbcc - mfbca;
m0 = m2 + mfbcb;
mfbca = m0;
- m0 += c1o9 * oMdrho;
+ m0 += UbMath::c1o9 * oMdrho;
mfbcb = m1 - m0 * vvz;
mfbcc = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -950,7 +914,7 @@ usage: ...
m1 = mfcac - mfcaa;
m0 = m2 + mfcab;
mfcaa = m0;
- m0 += c1o36 * oMdrho;
+ m0 += UbMath::c1o36 * oMdrho;
mfcab = m1 - m0 * vvz;
mfcac = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -958,7 +922,7 @@ usage: ...
m1 = mfcbc - mfcba;
m0 = m2 + mfcbb;
mfcba = m0;
- m0 += c1o9 * oMdrho;
+ m0 += UbMath::c1o9 * oMdrho;
mfcbb = m1 - m0 * vvz;
mfcbc = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -966,7 +930,7 @@ usage: ...
m1 = mfccc - mfcca;
m0 = m2 + mfccb;
mfcca = m0;
- m0 += c1o36 * oMdrho;
+ m0 += UbMath::c1o36 * oMdrho;
mfccb = m1 - m0 * vvz;
mfccc = m2 - 2. * m1 * vvz + vz2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -978,7 +942,7 @@ usage: ...
m1 = mfaca - mfaaa;
m0 = m2 + mfaba;
mfaaa = m0;
- m0 += c1o6 * oMdrho;
+ m0 += UbMath::c1o6 * oMdrho;
mfaba = m1 - m0 * vvy;
mfaca = m2 - 2. * m1 * vvy + vy2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -993,7 +957,7 @@ usage: ...
m1 = mfacc - mfaac;
m0 = m2 + mfabc;
mfaac = m0;
- m0 += c1o18 * oMdrho;
+ m0 += UbMath::c1o18 * oMdrho;
mfabc = m1 - m0 * vvy;
mfacc = m2 - 2. * m1 * vvy + vy2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -1002,7 +966,7 @@ usage: ...
m1 = mfbca - mfbaa;
m0 = m2 + mfbba;
mfbaa = m0;
- m0 += c2o3 * oMdrho;
+ m0 += UbMath::c2o3 * oMdrho;
mfbba = m1 - m0 * vvy;
mfbca = m2 - 2. * m1 * vvy + vy2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -1017,7 +981,7 @@ usage: ...
m1 = mfbcc - mfbac;
m0 = m2 + mfbbc;
mfbac = m0;
- m0 += c2o9 * oMdrho;
+ m0 += UbMath::c2o9 * oMdrho;
mfbbc = m1 - m0 * vvy;
mfbcc = m2 - 2. * m1 * vvy + vy2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -1026,7 +990,7 @@ usage: ...
m1 = mfcca - mfcaa;
m0 = m2 + mfcba;
mfcaa = m0;
- m0 += c1o6 * oMdrho;
+ m0 += UbMath::c1o6 * oMdrho;
mfcba = m1 - m0 * vvy;
mfcca = m2 - 2. * m1 * vvy + vy2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -1041,7 +1005,7 @@ usage: ...
m1 = mfccc - mfcac;
m0 = m2 + mfcbc;
mfcac = m0;
- m0 += c1o18 * oMdrho;
+ m0 += UbMath::c1o18 * oMdrho;
mfcbc = m1 - m0 * vvy;
mfccc = m2 - 2. * m1 * vvy + vy2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -1068,7 +1032,7 @@ usage: ...
m1 = mfcca - mfaca;
m0 = m2 + mfbca;
mfaca = m0;
- m0 += c1o3 * oMdrho;
+ m0 += UbMath::c1o3 * oMdrho;
mfbca = m1 - m0 * vvx;
mfcca = m2 - 2. * m1 * vvx + vx2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -1099,7 +1063,7 @@ usage: ...
m1 = mfcac - mfaac;
m0 = m2 + mfbac;
mfaac = m0;
- m0 += c1o3 * oMdrho;
+ m0 += UbMath::c1o3 * oMdrho;
mfbac = m1 - m0 * vvx;
mfcac = m2 - 2. * m1 * vvx + vx2 * m0;
////////////////////////////////////////////////////////////////////////////////////
@@ -1114,66 +1078,63 @@ usage: ...
m1 = mfccc - mfacc;
m0 = m2 + mfbcc;
mfacc = m0;
- m0 += c1o9 * oMdrho;
+ m0 += UbMath::c1o9 * oMdrho;
mfbcc = m1 - m0 * vvx;
mfccc = m2 - 2. * m1 * vvx + vx2 * m0;
////////////////////////////////////////////////////////////////////////////////////
// Cumulants
////////////////////////////////////////////////////////////////////////////////////
LBMReal OxxPyyPzz = 1.; //omega2 or bulk viscosity
- LBMReal OxyyPxzz = 1.;//-s9;//2+s9;//
- //LBMReal OxyyMxzz = 1.;//2+s9;//
- LBMReal O4 = 1.;
- LBMReal O5 = 1.;
- LBMReal O6 = 1.;
+ // LBMReal OxyyPxzz = 1.;//-s9;//2+s9;//
+ // LBMReal OxyyMxzz = 1.;//2+s9;//
//Cum 4.
- //LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
- //LBMReal CUMbcb = mfbcb - ((mfaca + c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
- //LBMReal CUMbbc = mfbbc - ((mfaac + c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015
+ //LBMReal CUMcbb = mfcbb - ((mfcaa + UbMath::c1o3 * oMdrho) * mfabb + 2. * mfbba * mfbab); // till 18.05.2015
+ //LBMReal CUMbcb = mfbcb - ((mfaca + UbMath::c1o3 * oMdrho) * mfbab + 2. * mfbba * mfabb); // till 18.05.2015
+ //LBMReal CUMbbc = mfbbc - ((mfaac + UbMath::c1o3 * oMdrho) * mfbba + 2. * mfbab * mfabb); // till 18.05.2015
- LBMReal CUMcbb = mfcbb - ((mfcaa + c1o3) * mfabb + 2. * mfbba * mfbab);
- LBMReal CUMbcb = mfbcb - ((mfaca + c1o3) * mfbab + 2. * mfbba * mfabb);
- LBMReal CUMbbc = mfbbc - ((mfaac + c1o3) * mfbba + 2. * mfbab * mfabb);
+ // LBMReal CUMcbb = mfcbb - ((mfcaa + UbMath::c1o3) * mfabb + 2. * mfbba * mfbab);
+ // LBMReal CUMbcb = mfbcb - ((mfaca + UbMath::c1o3) * mfbab + 2. * mfbba * mfabb);
+ // LBMReal CUMbbc = mfbbc - ((mfaac + UbMath::c1o3) * mfbba + 2. * mfbab * mfabb);
- LBMReal CUMcca = mfcca - ((mfcaa * mfaca + 2. * mfbba * mfbba) + c1o3 * (mfcaa + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
- LBMReal CUMcac = mfcac - ((mfcaa * mfaac + 2. * mfbab * mfbab) + c1o3 * (mfcaa + mfaac) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
- LBMReal CUMacc = mfacc - ((mfaac * mfaca + 2. * mfabb * mfabb) + c1o3 * (mfaac + mfaca) * oMdrho + c1o9 * (oMdrho - 1) * oMdrho);
+ // LBMReal CUMcca = mfcca - ((mfcaa * mfaca + 2. * mfbba * mfbba) + UbMath::c1o3 * (mfcaa + mfaca) * oMdrho + UbMath::c1o9 * (oMdrho - 1) * oMdrho);
+ // LBMReal CUMcac = mfcac - ((mfcaa * mfaac + 2. * mfbab * mfbab) + UbMath::c1o3 * (mfcaa + mfaac) * oMdrho + UbMath::c1o9 * (oMdrho - 1) * oMdrho);
+ // LBMReal CUMacc = mfacc - ((mfaac * mfaca + 2. * mfabb * mfabb) + UbMath::c1o3 * (mfaac + mfaca) * oMdrho + UbMath::c1o9 * (oMdrho - 1) * oMdrho);
//Cum 5.
- LBMReal CUMbcc = mfbcc - (mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb + 2. * (mfbab * mfacb + mfbba * mfabc)) - c1o3 * (mfbca + mfbac) * oMdrho;
- LBMReal CUMcbc = mfcbc - (mfaac * mfcba + mfcaa * mfabc + 4. * mfbab * mfbbb + 2. * (mfabb * mfcab + mfbba * mfbac)) - c1o3 * (mfcba + mfabc) * oMdrho;
- LBMReal CUMccb = mfccb - (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) - c1o3 * (mfacb + mfcab) * oMdrho;
+ // LBMReal CUMbcc = mfbcc - (mfaac * mfbca + mfaca * mfbac + 4. * mfabb * mfbbb + 2. * (mfbab * mfacb + mfbba * mfabc)) - UbMath::c1o3 * (mfbca + mfbac) * oMdrho;
+ // LBMReal CUMcbc = mfcbc - (mfaac * mfcba + mfcaa * mfabc + 4. * mfbab * mfbbb + 2. * (mfabb * mfcab + mfbba * mfbac)) - UbMath::c1o3 * (mfcba + mfabc) * oMdrho;
+ // LBMReal CUMccb = mfccb - (mfcaa * mfacb + mfaca * mfcab + 4. * mfbba * mfbbb + 2. * (mfbab * mfbca + mfabb * mfcba)) - UbMath::c1o3 * (mfacb + mfcab) * oMdrho;
//Cum 6.
- LBMReal CUMccc = mfccc + ((-4. * mfbbb * mfbbb
- - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
- - 4. * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
- - 2. * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb))
- + (4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
- + 2. * (mfcaa * mfaca * mfaac)
- + 16. * mfbba * mfbab * mfabb)
- - c1o3 * (mfacc + mfcac + mfcca) * oMdrho - c1o9 * oMdrho * oMdrho
- - c1o9 * (mfcaa + mfaca + mfaac) * oMdrho * (1. - 2. * oMdrho) - c1o27 * oMdrho * oMdrho * (-2. * oMdrho)
- + (2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
- + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3 * oMdrho) + c1o27 * oMdrho;
+// LBMReal CUMccc = mfccc + ((-4. * mfbbb * mfbbb
+// - (mfcaa * mfacc + mfaca * mfcac + mfaac * mfcca)
+// - 4. * (mfabb * mfcbb + mfbab * mfbcb + mfbba * mfbbc)
+// - 2. * (mfbca * mfbac + mfcba * mfabc + mfcab * mfacb))
+// + (4. * (mfbab * mfbab * mfaca + mfabb * mfabb * mfcaa + mfbba * mfbba * mfaac)
+// + 2. * (mfcaa * mfaca * mfaac)
+// + 16. * mfbba * mfbab * mfabb)
+// - UbMath::c1o3 * (mfacc + mfcac + mfcca) * oMdrho - UbMath::c1o9 * oMdrho * oMdrho
+// - UbMath::c1o9 * (mfcaa + mfaca + mfaac) * oMdrho * (1. - 2. * oMdrho) - c1o27 * oMdrho * oMdrho * (-2. * oMdrho)
+// + (2. * (mfbab * mfbab + mfabb * mfabb + mfbba * mfbba)
+// + (mfaac * mfaca + mfaac * mfcaa + mfaca * mfcaa)) * c2o3 * oMdrho) + c1o27 * oMdrho;
LBMReal mxxPyyPzz = mfcaa + mfaca + mfaac;
LBMReal mxxMyy = mfcaa - mfaca;
LBMReal mxxMzz = mfcaa - mfaac;
- LBMReal dxux = -c1o2 * collFactorF * (mxxMyy + mxxMzz) + c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
- LBMReal dyuy = dxux + collFactorF * c3o2 * mxxMyy;
- LBMReal dzuz = dxux + collFactorF * c3o2 * mxxMzz;
+ LBMReal dxux = -UbMath::c1o2 * collFactorF * (mxxMyy + mxxMzz) + UbMath::c1o2 * OxxPyyPzz * (mfaaa - mxxPyyPzz);
+ LBMReal dyuy = dxux + collFactorF * UbMath::c3o2 * mxxMyy;
+ LBMReal dzuz = dxux + collFactorF * UbMath::c3o2 * mxxMzz;
- LBMReal Dxy = -three * collFactorF * mfbba;
- LBMReal Dxz = -three * collFactorF * mfbab;
- LBMReal Dyz = -three * collFactorF * mfabb;
+ LBMReal Dxy = -UbMath::three * collFactorF * mfbba;
+ LBMReal Dxz = -UbMath::three * collFactorF * mfbab;
+ LBMReal Dyz = -UbMath::three * collFactorF * mfabb;
//TODO: may be factor 2
- return sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+ return sqrt(UbMath::c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + UbMath::one);
}
}
#endif
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
index 42c3240077e4b6cc3e3a279608bca0325b304b8c..de780bc0574ce169432007dd9eb257de987249a9 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.cpp
@@ -10,6 +10,8 @@
#define PROOF_CORRECTNESS
+using namespace UbMath;
+
//////////////////////////////////////////////////////////////////////////
RheologyK17LBMKernel::RheologyK17LBMKernel()
{
@@ -83,10 +85,6 @@ void RheologyK17LBMKernel::calculate(int step)
muForcingX1.DefineVar("nu", &muNu);
muForcingX2.DefineVar("nu", &muNu);
muForcingX3.DefineVar("nu", &muNu);
-
- LBMReal forcingX1 = 0;
- LBMReal forcingX2 = 0;
- LBMReal forcingX3 = 0;
}
/////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h
index f2c211dffa5da0ca5e189759c75bdaad5a9a5006..f9838ba2e73e8d49a5d637b50a871faf751b4e7c 100644
--- a/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h
+++ b/src/cpu/VirtualFluidsCore/LBM/RheologyK17LBMKernel.h
@@ -19,8 +19,8 @@ public:
public:
RheologyK17LBMKernel();
virtual ~RheologyK17LBMKernel(void);
- virtual void calculate(int step);
- virtual SPtr<LBMKernel> clone();
+ virtual void calculate(int step) override;
+ virtual SPtr<LBMKernel> clone() override;
double getCalculationTime() override;
//! The value should not be equal to a shear viscosity
void setBulkViscosity(LBMReal value);
diff --git a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
index a5ad23dabe7f4d048e854ac9d42a1674d73458a5..d54e2ac5b598eb2b01e27beb9c9e9410592cf266 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
+++ b/src/cpu/VirtualFluidsCore/LBM/Thixotropy.h
@@ -61,7 +61,7 @@ inline LBMReal Thixotropy::getBinghamCollFactor(LBMReal omegaInf, LBMReal shearR
//LBMReal omega = cs2 * cs2 * shearRate * shearRate * omegaInf * rho * rho / (cs2 * cs2 * shearRate * shearRate * rho * rho + cs2 * shearRate * omegaInf * rho * tau0+omegaInf*omegaInf*tau0*tau0);
- LBMReal a = omegaInf * tau0 / (cs2 * shearRate * rho);
+ // LBMReal a = omegaInf * tau0 / (cs2 * shearRate * rho);
//10 iterations
//LBMReal omega = omegaInf / (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a * (1 + a))))))))));
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp
index adbca6c6c57eab9fd520fe7e7fe1f8d7c055053c..a13e4fc716725156adc28841da22c4b2516dc24f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyExpLBMKernel.cpp
@@ -8,6 +8,8 @@
#define PROOF_CORRECTNESS
+using namespace UbMath;
+
//////////////////////////////////////////////////////////////////////////
ThixotropyExpLBMKernel::ThixotropyExpLBMKernel()
{
@@ -84,10 +86,6 @@ void ThixotropyExpLBMKernel::calculate(int step)
muForcingX1.DefineVar("nu", &muNu);
muForcingX2.DefineVar("nu", &muNu);
muForcingX3.DefineVar("nu", &muNu);
-
- LBMReal forcingX1 = 0;
- LBMReal forcingX2 = 0;
- LBMReal forcingX3 = 0;
}
/////////////////////////////////////
@@ -233,8 +231,8 @@ void ThixotropyExpLBMKernel::calculate(int step)
(mfbbc - mfbba));
- LBMReal eta0 = (1/collFactor-c1o2)*c1o3;
- LBMReal eta = (1 + lambda)* eta0;
+// LBMReal eta0 = (1/collFactor-c1o2)*c1o3;
+// LBMReal eta = (1 + lambda)* eta0;
//collFactorF = one/(3*eta/(rho+one)+c1o2);
collFactorF = collFactor;
@@ -1761,9 +1759,9 @@ void ThixotropyExpLBMKernel::calculate(int step)
//proof correctness
//////////////////////////////////////////////////////////////////////////
//#ifdef PROOF_CORRECTNESS
- LBMReal drho_post = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca)
- + (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
- + (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
+// LBMReal 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;
//UBLOG(logINFO, "lambda ="<<drho_post);
// //LBMReal dif = fabs(rho - rho_post);
// dif = drho - drho_post;
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp
index 0f1edd355acde4f072c0bdfd36b25fc541046ced..a673f473e13715e1a90db22e502446066e695ca0 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.cpp
@@ -82,29 +82,29 @@ void ThixotropyInterpolationProcessor::setOffsets(LBMReal xoff, LBMReal yoff, LB
//////////////////////////////////////////////////////////////////////////
void ThixotropyInterpolationProcessor::interpolateCoarseToFine(D3Q27ICell& icellC, D3Q27ICell& icellF, LBMReal xoff, LBMReal yoff, LBMReal zoff)
{
- setOffsets(xoff, yoff, zoff);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, -0.25, -0.25, -1, -1, -1);
- calcInterpolatedNode(icellF.BSW, /*omegaF,*/ -0.25, -0.25, -0.25, calcPressBSW(), -1, -1, -1);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, 0.25, -0.25, 1, 1, -1);
- calcInterpolatedNode(icellF.BNE, /*omegaF,*/ 0.25, 0.25, -0.25, calcPressBNE(), 1, 1, -1);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, -0.25, 0.25, 0.25, -1, 1, 1);
- calcInterpolatedNode(icellF.TNW, /*omegaF,*/ -0.25, 0.25, 0.25, calcPressTNW(), -1, 1, 1);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, -0.25, 0.25, 1, -1, 1);
- calcInterpolatedNode(icellF.TSE, /*omegaF,*/ 0.25, -0.25, 0.25, calcPressTSE(), 1, -1, 1);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, -0.25, 0.25, -0.25, -1, 1, -1);
- calcInterpolatedNode(icellF.BNW, /*omegaF,*/ -0.25, 0.25, -0.25, calcPressBNW(), -1, 1, -1);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, -0.25, -0.25, 1, -1, -1);
- calcInterpolatedNode(icellF.BSE, /*omegaF,*/ 0.25, -0.25, -0.25, calcPressBSE(), 1, -1, -1);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, -0.25, -0.25, 0.25, -1, -1, 1);
- calcInterpolatedNode(icellF.TSW, /*omegaF,*/ -0.25, -0.25, 0.25, calcPressTSW(), -1, -1, 1);
- calcInterpolatedCoefficiets(icellC, omegaC, 0.5, 0.25, 0.25, 0.25, 1, 1, 1);
- calcInterpolatedNode(icellF.TNE, /*omegaF,*/ 0.25, 0.25, 0.25, calcPressTNE(), 1, 1, 1);
+ setOffsets(xoff, yoff, zoff);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, 0.25, -0.25, -0.25, -1, -1, -1);
+ calcInterpolatedNode(icellF.BSW, /*omegaF,*/ -0.25, -0.25, -0.25, calcPressBSW(), -1, -1, -1);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, 0.25, 0.25, -0.25, 1, 1, -1);
+ calcInterpolatedNode(icellF.BNE, /*omegaF,*/ 0.25, 0.25, -0.25, calcPressBNE(), 1, 1, -1);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, -0.25, 0.25, 0.25, -1, 1, 1);
+ calcInterpolatedNode(icellF.TNW, /*omegaF,*/ -0.25, 0.25, 0.25, calcPressTNW(), -1, 1, 1);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, 0.25, -0.25, 0.25, 1, -1, 1);
+ calcInterpolatedNode(icellF.TSE, /*omegaF,*/ 0.25, -0.25, 0.25, calcPressTSE(), 1, -1, 1);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, -0.25, 0.25, -0.25, -1, 1, -1);
+ calcInterpolatedNode(icellF.BNW, /*omegaF,*/ -0.25, 0.25, -0.25, calcPressBNW(), -1, 1, -1);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, 0.25, -0.25, -0.25, 1, -1, -1);
+ calcInterpolatedNode(icellF.BSE, /*omegaF,*/ 0.25, -0.25, -0.25, calcPressBSE(), 1, -1, -1);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, -0.25, -0.25, 0.25, -1, -1, 1);
+ calcInterpolatedNode(icellF.TSW, /*omegaF,*/ -0.25, -0.25, 0.25, calcPressTSW(), -1, -1, 1);
+ calcInterpolatedCoefficiets_intern(icellC, omegaC, 0.5, 0.25, 0.25, 0.25, 1, 1, 1);
+ calcInterpolatedNode(icellF.TNE, /*omegaF,*/ 0.25, 0.25, 0.25, calcPressTNE(), 1, 1, 1);
}
//////////////////////////////////////////////////////////////////////////
void ThixotropyInterpolationProcessor::interpolateFineToCoarse(D3Q27ICell& icellF, LBMReal* icellC, LBMReal xoff, LBMReal yoff, LBMReal zoff)
{
setOffsets(xoff, yoff, zoff);
- calcInterpolatedCoefficiets(icellF, omegaF, 2.0, 0, 0, 0, 0, 0, 0);
+ calcInterpolatedCoefficiets_intern(icellF, omegaF, 2.0, 0, 0, 0, 0, 0, 0);
calcInterpolatedNodeFC(icellC, omegaC);
}
//////////////////////////////////////////////////////////////////////////
@@ -128,7 +128,15 @@ void ThixotropyInterpolationProcessor::calcMoments(const LBMReal* const f, LBMRe
kxxMzz = -3./2.*omega*((((f[NW]+f[SE])-(f[BS]+f[TN]))+((f[SW]+f[NE])-(f[TS]+f[BN])))+((f[W]+f[E])-(f[B]+f[T]))-(vx1*vx1-vx3*vx3));
}
//////////////////////////////////////////////////////////////////////////
-void ThixotropyInterpolationProcessor::calcInterpolatedCoefficiets(const D3Q27ICell& icell, LBMReal omega, LBMReal eps_new, LBMReal x, LBMReal y, LBMReal z, LBMReal xs, LBMReal ys, LBMReal zs)
+void ThixotropyInterpolationProcessor::calcInterpolatedCoefficiets_intern(const D3Q27ICell& icell,
+ LBMReal omega,
+ LBMReal eps_new,
+ LBMReal x,
+ LBMReal y,
+ LBMReal z,
+ LBMReal xs,
+ LBMReal ys,
+ LBMReal zs)
{
LBMReal vx1_SWT,vx2_SWT,vx3_SWT;
LBMReal vx1_NWT,vx2_NWT,vx3_NWT;
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h
index 9f55b8f5fec27ffa116f32f8b96df3c44e1dd77e..daaa6449ea8a37d16aa9d359cf72ab58d6d872f0 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyInterpolationProcessor.h
@@ -70,8 +70,6 @@ private:
LBMReal kxyAverage, kyzAverage, kxzAverage, kxxMyyAverage, kxxMzzAverage;
- LBMReal a,b,c;
-
LBMReal rho;
LBMReal shearRate;
@@ -80,7 +78,7 @@ private:
void setOffsets(LBMReal xoff, LBMReal yoff, LBMReal zoff);
void calcMoments(const LBMReal* const f, LBMReal omegaInf, LBMReal& rho, LBMReal& vx1, LBMReal& vx2, LBMReal& vx3,
LBMReal& kxy, LBMReal& kyz, LBMReal& kxz, LBMReal& kxxMyy, LBMReal& kxxMzz);
- void calcInterpolatedCoefficiets(const D3Q27ICell& icell, LBMReal omega, LBMReal eps_new, LBMReal x, LBMReal y, LBMReal z, LBMReal xs, LBMReal ys, LBMReal zs);
+ void calcInterpolatedCoefficiets_intern(const D3Q27ICell& icell, LBMReal omega, LBMReal eps_new, LBMReal x, LBMReal y, LBMReal z, LBMReal xs, LBMReal ys, LBMReal zs);
void calcInterpolatedNode(LBMReal* f, /*LBMReal omega,*/ LBMReal x, LBMReal y, LBMReal z, LBMReal press, LBMReal xs, LBMReal ys, LBMReal zs);
LBMReal calcPressBSW();
LBMReal calcPressTSW();
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp
index 949e7777d644c5b9051611fbcc4b8277b37d630d..b369b45a6c7b10efb91716634443c88aa520a8cf 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyLBMKernel.cpp
@@ -8,6 +8,8 @@
#define PROOF_CORRECTNESS
+using namespace UbMath;
+
//////////////////////////////////////////////////////////////////////////
ThixotropyLBMKernel::ThixotropyLBMKernel()
{
@@ -84,10 +86,6 @@ void ThixotropyLBMKernel::calculate(int step)
muForcingX1.DefineVar("nu", &muNu);
muForcingX2.DefineVar("nu", &muNu);
muForcingX3.DefineVar("nu", &muNu);
-
- LBMReal forcingX1 = 0;
- LBMReal forcingX2 = 0;
- LBMReal forcingX3 = 0;
}
/////////////////////////////////////
@@ -1755,9 +1753,9 @@ void ThixotropyLBMKernel::calculate(int step)
//proof correctness
//////////////////////////////////////////////////////////////////////////
//#ifdef PROOF_CORRECTNESS
- LBMReal drho_post = (mfaaa + mfaac + mfaca + mfcaa + mfacc + mfcac + mfccc + mfcca)
- + (mfaab + mfacb + mfcab + mfccb) + (mfaba + mfabc + mfcba + mfcbc) + (mfbaa + mfbac + mfbca + mfbcc)
- + (mfabb + mfcbb) + (mfbab + mfbcb) + (mfbba + mfbbc) + mfbbb;
+// LBMReal 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;
//UBLOG(logINFO, "lambda ="<<drho_post);
// //LBMReal dif = fabs(rho - rho_post);
// dif = drho - drho_post;
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
index d0d0542b63ba68136c0c23ec9cbc86bb37354ea5..1685eb6e1393055f0e04f7f4289fd295e805c33f 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel.cpp
@@ -8,6 +8,8 @@
#define PROOF_CORRECTNESS
+using namespace UbMath;
+
ThixotropyModelLBMKernel::ThixotropyModelLBMKernel() : forcingX1(0), forcingX2(0), forcingX3(0)
{
compressible = false;
diff --git a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
index 8839818a70892d30b57de2938360b0199af3e234..44decadb7a4fb720ce0785373f9acd00106cbe10 100644
--- a/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/ThixotropyModelLBMKernel2.cpp
@@ -9,6 +9,9 @@
#define PROOF_CORRECTNESS
+using namespace UbMath;
+
+
ThixotropyModelLBMKernel2::ThixotropyModelLBMKernel2() : forcingX1(0), forcingX2(0), forcingX3(0)
{
compressible = false;
@@ -476,12 +479,12 @@ void ThixotropyModelLBMKernel2::calculate(int step)
LBMReal dyuy = dxux + collFactorF * c3o2 * mxxMyy;
LBMReal dzuz = dxux + collFactorF * c3o2 * mxxMzz;
- LBMReal Dxy = -three * collFactorF * mfbba;
- LBMReal Dxz = -three * collFactorF * mfbab;
- LBMReal Dyz = -three * collFactorF * mfabb;
+// LBMReal Dxy = -three * collFactorF * mfbba;
+// LBMReal Dxz = -three * collFactorF * mfbab;
+// LBMReal Dyz = -three * collFactorF * mfabb;
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
//non Newtonian fluid collision factor
- LBMReal shearRate = sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
+// LBMReal shearRate = sqrt(c2 * (dxux * dxux + dyuy * dyuy + dzuz * dzuz) + Dxy * Dxy + Dxz * Dxz + Dyz * Dyz) / (rho + one);
LBMReal shearFactor = sqrt(c1o2 * ((mfcaa - mfaaa * c1o3) * (mfcaa - mfaaa * c1o3) + (mfaca - mfaaa * c1o3) * (mfaca - mfaaa * c1o3) + (mfaac - mfaaa * c1o3) * (mfaac - mfaaa * c1o3)) + mfbba * mfbba + mfbab * mfbab + mfabb * mfabb) + UbMath::Epsilon<LBMReal>::val();
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
index c17d8f659546e3a9b2ceabbdd9397bdbb10a5075..8704355286f145efb926e78e8af61e03a310c26c 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
@@ -163,10 +163,6 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
LBMReal f[D3Q27System::ENDF+1];
- size_t nx1 = distributions->getNX1();
- size_t nx2 = distributions->getNX2();
- size_t nx3 = distributions->getNX3();
-
for(int ix3=0; ix3<bcArray->getNX3(); ix3++)
for(int ix2=0; ix2<bcArray->getNX2(); ix2++)
for(int ix1=0; ix1<bcArray->getNX1(); ix1++)
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
index e63ad0e74d3a0c8f1ab2adb55e73c054399131e9..17668ba25330dadeb05cc8769f84565b2f5cbe2e 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
@@ -221,7 +221,7 @@ void InitThixotropyBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D> bl
if(!block) UB_THROW( UbException(UB_EXARGS,"block is not exist") );
- double dx = grid->getDeltaX(block);
+// double dx = grid->getDeltaX(block);
//define vars for functions
mu::value_type x1,x2,x3;
@@ -235,8 +235,6 @@ void InitThixotropyBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D> bl
//Funktionszeiger
typedef void (*CalcFeqsFct)(LBMReal* const& /*feq[27]*/,const LBMReal& /*(d)rho*/,const LBMReal& /*vx1*/,const LBMReal& /*vx2*/,const LBMReal& /*vx3*/);
CalcFeqsFct calcFeqsFct = NULL;
-
- LBMReal vx1,vx2,vx3,rho;
int gridRank = grid->getRank();
int blockRank = block->getRank();
@@ -255,14 +253,8 @@ void InitThixotropyBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D> bl
SPtr<BCArray3D> bcArray = kernel->getBCProcessor()->getBCArray();
SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getHdistributions();
- LBMReal o = kernel->getCollisionFactor();
-
LBMReal h[D3Q27System::ENDF+1];
- size_t nx1 = distributions->getNX1();
- size_t nx2 = distributions->getNX2();
- size_t nx3 = distributions->getNX3();
-
for(int ix3=0; ix3<bcArray->getNX3(); ix3++)
for(int ix2=0; ix2<bcArray->getNX2(); ix2++)
for(int ix1=0; ix1<bcArray->getNX1(); ix1++)
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp
index 9ea340b1618bcc2ad31a21a958424d2b256cf357..c5a2bdebc92ea478942a2a26c3c1f8ebfd8c48bf 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetConnectorsBlockVisitor.cpp
@@ -143,11 +143,6 @@ void SetConnectorsBlockVisitor::setRemoteConnectors(SPtr<Block3D> sblock, SPtr<B
void SetConnectorsBlockVisitor::setInterpolationConnectors(SPtr<Grid3D> grid, SPtr<Block3D> block)
{
UBLOG(logDEBUG5, "D3Q27SetConnectorsBlockVisitor::setInterpolationConnectors() - start");
- int blockRank = block->getRank();
- if (block->getGlobalID()==394)
- {
- int test=0;
- }
//search for all blocks with different ranks
if (block->hasInterpolationFlagCF() && block->isActive())
@@ -490,8 +485,6 @@ void SetConnectorsBlockVisitor::createTransmitters(SPtr<Block3D> cBlock, SPtr<Bl
{
UBLOG(logDEBUG5, "D3Q27SetConnectorsBlockVisitor::createTransmitters(...) - start");
CreateTransmittersHelper helper;
- bool MPIpool = true;
- bool orthogonal = false;
int fBlockRank = fBlock->getRank();
int cBlockRank = cBlock->getRank();
if(fBlockRank == cBlockRank && fBlockRank == gridRank)