diff --git a/src/cpu/VirtualFluidsCore/Connectors/CoarseToFineVectorConnector.h b/src/cpu/VirtualFluidsCore/Connectors/CoarseToFineVectorConnector.h
index d7f9f72c4f86f5d849fa4f2c0cd703333007091a..6dce407d05aa983767cfd9697d395533d90177a2 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/CoarseToFineVectorConnector.h
+++ b/src/cpu/VirtualFluidsCore/Connectors/CoarseToFineVectorConnector.h
@@ -47,6 +47,7 @@
#include "basics/transmitter/TbTransmitter.h"
#include "basics/transmitter/TbTransmitterLocal.h"
#include <PointerDefinitions.h>
+#include "basics/constants/NumericConstants.h"
#include "BCSet.h"
#include "FineToCoarseVectorConnector.h"
@@ -829,6 +830,8 @@ void CoarseToFineVectorConnector<VectorTransmitter>::fillSendVectorExt(SPtr<Dist
const int &lMaxX2, const int &lMaxX3,
vector_type &data, int &index)
{
+ using namespace vf::basics::constant;
+
if (data.size() == 0)
return;
int ix1, ix2, ix3;
@@ -845,9 +848,9 @@ void CoarseToFineVectorConnector<VectorTransmitter>::fillSendVectorExt(SPtr<Dist
if (howManySolids == 0 || howManySolids == 8) {
iprocessor->readICell(fFrom, icellC, ix1, ix2, ix3);
- xoff = 0.0;
- yoff = 0.0;
- zoff = 0.0;
+ xoff = c0o1;
+ yoff = c0o1;
+ zoff = c0o1;
} else {
if (!iprocessor->findNeighborICell(bcArray, fFrom, icellC, bMaxX1, bMaxX2, bMaxX3, ix1, ix2, ix3,
xoff, yoff, zoff)) {
@@ -1991,9 +1994,9 @@ void CoarseToFineVectorConnector<VectorTransmitter>::findCFnodes(SPtr<Distributi
if (howManySolids == 0 || howManySolids == 8) {
iprocessor->readICell(fFrom, icellC, ix1, ix2, ix3);
- xoff = 0.0;
- yoff = 0.0;
- zoff = 0.0;
+ xoff = c0o1;
+ yoff = c0o1;
+ zoff = c0o1;
} else {
if (!iprocessor->findNeighborICell(bcArray, fFrom, icellC, bMaxX1, bMaxX2, bMaxX3, ix1, ix2, ix3,
xoff, yoff, zoff)) {
diff --git a/src/cpu/VirtualFluidsCore/Connectors/FineToCoarseVectorConnector.h b/src/cpu/VirtualFluidsCore/Connectors/FineToCoarseVectorConnector.h
index ffab5008241e7f8f4457c7385c83872a296ee748..7a18e86791fed253c9a52a26709b6a1a9861ac3b 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/FineToCoarseVectorConnector.h
+++ b/src/cpu/VirtualFluidsCore/Connectors/FineToCoarseVectorConnector.h
@@ -45,6 +45,7 @@
#include "MathUtil.hpp"
#include "basics/transmitter/TbTransmitter.h"
#include <PointerDefinitions.h>
+#include "basics/constants/NumericConstants.h"
#include "BCSet.h"
#include "DataSet3D.h"
@@ -793,6 +794,8 @@ void FineToCoarseVectorConnector<VectorTransmitter>::fillSendVector(SPtr<Distrib
const int &lMaxX1, const int &lMaxX2,
const int &lMaxX3, vector_type &data, int &index)
{
+ using namespace vf::basics::constant;
+
int ix1, ix2, ix3;
real xoff, yoff, zoff;
SPtr<BCArray3D> bcArray = block.lock()->getKernel()->getBCSet()->getBCArray();
@@ -807,9 +810,9 @@ void FineToCoarseVectorConnector<VectorTransmitter>::fillSendVector(SPtr<Distrib
if (howManySolids == 0 || howManySolids == 8) {
iprocessor->readICell(fFrom, icellF, ix1, ix2, ix3);
- xoff = 0.0;
- yoff = 0.0;
- zoff = 0.0;
+ xoff = c0o1;
+ yoff = c0o1;
+ zoff = c0o1;
} else {
if (!iprocessor->findNeighborICell(bcArray, fFrom, icellF, bMaxX1, bMaxX2, bMaxX3, ix1, ix2, ix3,
xoff, yoff, zoff)) {
diff --git a/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp
index 1bdb92f6b0d51c3bfb8daf6e149be7c1be0fecf0..be66dfb5b0785198e826f3fe2017aa479b92b53a 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/OneDistributionFullVectorConnector.cpp
@@ -13,6 +13,7 @@ OneDistributionFullVectorConnector::OneDistributionFullVectorConnector(SPtr<Bloc
void OneDistributionFullVectorConnector::init()
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
FullVectorConnector::init();
@@ -25,36 +26,36 @@ void OneDistributionFullVectorConnector::init()
break;
case DIR_P00:
case DIR_M00:
- sender->getData().resize(maxX2 * maxX3 * anz, 0.0);
+ sender->getData().resize(maxX2 * maxX3 * anz, c0o1);
break;
case DIR_0P0:
case DIR_0M0:
- sender->getData().resize(maxX1 * maxX3 * anz, 0.0);
+ sender->getData().resize(maxX1 * maxX3 * anz, c0o1);
break;
case DIR_00P:
case DIR_00M:
- sender->getData().resize(maxX1 * maxX2 * anz, 0.0);
+ sender->getData().resize(maxX1 * maxX2 * anz, c0o1);
break;
case DIR_PP0:
case DIR_MM0:
case DIR_PM0:
case DIR_MP0:
- sender->getData().resize(maxX3 * anz, 0.0);
+ sender->getData().resize(maxX3 * anz, c0o1);
break;
case DIR_P0P:
case DIR_M0M:
case DIR_P0M:
case DIR_M0P:
- sender->getData().resize(maxX2 * anz, 0.0);
+ sender->getData().resize(maxX2 * anz, c0o1);
break;
case DIR_0PP:
case DIR_0MM:
case DIR_0PM:
case DIR_0MP:
- sender->getData().resize(maxX1 * anz, 0.0);
+ sender->getData().resize(maxX1 * anz, c0o1);
break;
case DIR_PPP:
@@ -65,7 +66,7 @@ void OneDistributionFullVectorConnector::init()
case DIR_MPM:
case DIR_PMM:
case DIR_MPP:
- sender->getData().resize(anz, 0.0);
+ sender->getData().resize(anz, c0o1);
break;
default:
diff --git a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp
index 8f6b88898a9da1cfca9aee49ae4cb084ee54217a..f8ff7d5c28bb04e41f51a5223dc9a64bcd0cf317 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsDoubleGhostLayerFullVectorConnector.cpp
@@ -51,6 +51,7 @@ TwoDistributionsDoubleGhostLayerFullVectorConnector::TwoDistributionsDoubleGhost
void TwoDistributionsDoubleGhostLayerFullVectorConnector::init()
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
FullVectorConnector::init();
@@ -63,26 +64,26 @@ void TwoDistributionsDoubleGhostLayerFullVectorConnector::init()
{
case DIR_000: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
case DIR_P00:
- case DIR_M00: sender->getData().resize(maxX2*maxX3*anz*2, 0.0); break;
+ case DIR_M00: sender->getData().resize(maxX2*maxX3*anz*2, c0o1); break;
case DIR_0P0:
- case DIR_0M0: sender->getData().resize(maxX1*maxX3*anz*2, 0.0); break;
+ case DIR_0M0: sender->getData().resize(maxX1*maxX3*anz*2, c0o1); break;
case DIR_00P:
- case DIR_00M: sender->getData().resize(maxX1*maxX2*anz*2, 0.0); break;
+ case DIR_00M: sender->getData().resize(maxX1*maxX2*anz*2, c0o1); break;
case DIR_PP0:
case DIR_MM0:
case DIR_PM0:
- case DIR_MP0: sender->getData().resize(maxX3*anz*4, 0.0); break;
+ case DIR_MP0: sender->getData().resize(maxX3*anz*4, c0o1); break;
case DIR_P0P:
case DIR_M0M:
case DIR_P0M:
- case DIR_M0P: sender->getData().resize(maxX2*anz*4, 0.0); break;
+ case DIR_M0P: sender->getData().resize(maxX2*anz*4, c0o1); break;
case DIR_0PP:
case DIR_0MM:
case DIR_0PM:
- case DIR_0MP: sender->getData().resize(maxX1*anz*4, 0.0); break;
+ case DIR_0MP: sender->getData().resize(maxX1*anz*4, c0o1); break;
case DIR_PPP:
case DIR_MMM:
@@ -91,7 +92,7 @@ void TwoDistributionsDoubleGhostLayerFullVectorConnector::init()
case DIR_PMP:
case DIR_MPM:
case DIR_PMM:
- case DIR_MPP: sender->getData().resize(anz*8, 0.0); break;
+ case DIR_MPP: sender->getData().resize(anz*8, c0o1); break;
default: UB_THROW(UbException(UB_EXARGS, "unknown sendDir"));
}
diff --git a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp
index 7987c2f6c8af52fbf897ff6bbcee47add3fc0056..36e34321d0bcaecec63b70d2b059542dc5514b95 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/TwoDistributionsFullVectorConnector.cpp
@@ -51,7 +51,8 @@ TwoDistributionsFullVectorConnector::TwoDistributionsFullVectorConnector(SPtr<Bl
void TwoDistributionsFullVectorConnector::init()
{
using namespace vf::lbm::dir;
-
+ using namespace vf::basics::constant;
+
FullVectorConnector::init();
fDis = dynamicPointerCast<EsoTwist3D>(block.lock()->getKernel()->getDataSet()->getFdistributions());
@@ -62,26 +63,26 @@ void TwoDistributionsFullVectorConnector::init()
{
case DIR_000: UB_THROW(UbException(UB_EXARGS, "ZERO not allowed")); break;
case DIR_P00:
- case DIR_M00: sender->getData().resize(maxX2*maxX3*anz, 0.0); break;
+ case DIR_M00: sender->getData().resize(maxX2*maxX3*anz, c0o1); break;
case DIR_0P0:
- case DIR_0M0: sender->getData().resize(maxX1*maxX3*anz, 0.0); break;
+ case DIR_0M0: sender->getData().resize(maxX1*maxX3*anz, c0o1); break;
case DIR_00P:
- case DIR_00M: sender->getData().resize(maxX1*maxX2*anz, 0.0); break;
+ case DIR_00M: sender->getData().resize(maxX1*maxX2*anz, c0o1); break;
case DIR_PP0:
case DIR_MM0:
case DIR_PM0:
- case DIR_MP0: sender->getData().resize(maxX3*anz, 0.0); break;
+ case DIR_MP0: sender->getData().resize(maxX3*anz, c0o1); break;
case DIR_P0P:
case DIR_M0M:
case DIR_P0M:
- case DIR_M0P: sender->getData().resize(maxX2*anz, 0.0); break;
+ case DIR_M0P: sender->getData().resize(maxX2*anz, c0o1); break;
case DIR_0PP:
case DIR_0MM:
case DIR_0PM:
- case DIR_0MP: sender->getData().resize(maxX1*anz, 0.0); break;
+ case DIR_0MP: sender->getData().resize(maxX1*anz, c0o1); break;
case DIR_PPP:
case DIR_MMM:
@@ -90,7 +91,7 @@ void TwoDistributionsFullVectorConnector::init()
case DIR_PMP:
case DIR_MPM:
case DIR_PMM:
- case DIR_MPP: sender->getData().resize(anz, 0.0); break;
+ case DIR_MPP: sender->getData().resize(anz, c0o1); break;
default: UB_THROW(UbException(UB_EXARGS, "unknown sendDir"));
}
diff --git a/src/cpu/VirtualFluidsCore/Simulation/Grid3D.cpp b/src/cpu/VirtualFluidsCore/Simulation/Grid3D.cpp
index a214b4bd0137b2bf319925b519f1dcb77fabded4..ed467597f473119fe1569fe2d198507ac0744909 100644
--- a/src/cpu/VirtualFluidsCore/Simulation/Grid3D.cpp
+++ b/src/cpu/VirtualFluidsCore/Simulation/Grid3D.cpp
@@ -64,10 +64,12 @@ Grid3D::Grid3D(std::shared_ptr<vf::mpi::Communicator> comm, int blockNx1, int bl
blockNx1(blockNx1), blockNx2(blockNx2), blockNx3(blockNx2), nx1(gridNx1), nx2(gridNx2), nx3(gridNx3)
{
+ using namespace vf::basics::constant;
+
levelSet.resize(D3Q27System::MAXLEVEL + 1);
bundle = comm->getBundleID();
rank = comm->getProcessID();
- trafo = std::make_shared<CoordinateTransformation3D>(0.0, 0.0, 0.0, (real)blockNx1, (real)blockNx2,
+ trafo = std::make_shared<CoordinateTransformation3D>(c0o1, c0o1, c0o1, (real)blockNx1, (real)blockNx2,
(real)blockNx3);
UbTupleInt3 minInd(0, 0, 0);
UbTupleInt3 maxInd(gridNx1, gridNx2, gridNx3);
@@ -481,7 +483,8 @@ UbTupleDouble3 Grid3D::getBlockLengths(const SPtr<Block3D> block) const
trafo->getX3CoordinateScaling() * delta);
}
//////////////////////////////////////////////////////////////////////////
-UbTupleDouble6 Grid3D::getBlockOversize() const { return makeUbTuple(0.0, 0.0, 0.0, 0.0, 0.0, 0.0); }
+using namespace vf::basics::constant;
+UbTupleDouble6 Grid3D::getBlockOversize() const { return makeUbTuple(c0o1, c0o1, c0o1, c0o1, c0o1, c0o1); }
//////////////////////////////////////////////////////////////////////////
void Grid3D::setCoordinateTransformator(SPtr<CoordinateTransformation3D> trafo) { this->trafo = trafo; }
//////////////////////////////////////////////////////////////////////////
@@ -2233,7 +2236,7 @@ int Grid3D::getGhostLayerWidth() const
//////////////////////////////////////////////////////////////////////////
void Grid3D::setGhostLayerWidth(int ghostLayerWidth)
{
- this->offset = static_cast<real>(ghostLayerWidth) - 0.5;
+ this->offset = static_cast<real>(ghostLayerWidth) - c1o2;
}
//////////////////////////////////////////////////////////////////////////
void Grid3D::setTimeStep(real step) { timeStep = step; }
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/AdjustForcingSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/AdjustForcingSimulationObserver.cpp
index 2254b9b02ea383e18c654a7569f0e5b2e973c839..3d6feaf15c6555cec3141c7bcc48d9a94db9c91c 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/AdjustForcingSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/AdjustForcingSimulationObserver.cpp
@@ -18,25 +18,27 @@ AdjustForcingSimulationObserver::AdjustForcingSimulationObserver(SPtr<Grid3D> gr
: SimulationObserver(grid, s), path(path), integrateValues(integrateValues), comm(comm), vx1Targed(vTarged)
{
+ using namespace vf::basics::constant;
+
// cnodes = integrateValues->getCNodes();
root = comm->isRoot();
Ta = scheduler->getMaxStep();
- Kpcrit = 3.0 / Ta; // 0.3;
- Tcrit = 3.0 * Ta; // 30.0;
- Tn = 0.5 * Tcrit;
- Tv = 0.12 * Tcrit;
+ Kpcrit = c3o1 / Ta; // 0.3;
+ Tcrit = c3o1 * Ta; // 30.0;
+ Tn = c1o2 * Tcrit;
+ Tv = c12o1 / c100o1 * Tcrit;
- Kp = 0.6 * Kpcrit;
+ Kp = c6o1 / c10o1 * Kpcrit;
Ki = Kp / Tn;
Kd = Kp * Tv;
- y = 0;
- e = 0;
- esum = 0;
- eold = 0;
- forcing = 0;
+ y = c0o1;
+ e = c0o1;
+ esum = c0o1;
+ eold = c0o1;
+ forcing = c0o1;
if (root) {
std::string fname = path + "/forcing/forcing.csv";
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/AverageValuesSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/AverageValuesSimulationObserver.cpp
index 1adf3ad9944a49c8065756988e95ab837e9f6d15..6213f5713c281b18dd516d304cd47f0d3195f9f3 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/AverageValuesSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/AverageValuesSimulationObserver.cpp
@@ -78,6 +78,8 @@ void AverageValuesSimulationObserver::update(real step)
void AverageValuesSimulationObserver::resetDataRMS(real step)
{
+ using namespace vf::basics::constant;
+
resetStepRMS = (int)step;
for (int level = minInitLevel; level <= maxInitLevel; level++) {
@@ -103,13 +105,13 @@ void AverageValuesSimulationObserver::resetDataRMS(real step)
//////////////////////////////////////////////////////////////////////////
// compute average values
//////////////////////////////////////////////////////////////////////////
- (*av)(AvVxx, ix1, ix2, ix3) = 0.0;
- (*av)(AvVyy, ix1, ix2, ix3) = 0.0;
- (*av)(AvVzz, ix1, ix2, ix3) = 0.0;
- (*av)(AvVxy, ix1, ix2, ix3) = 0.0;
- (*av)(AvVxz, ix1, ix2, ix3) = 0.0;
- (*av)(AvVyz, ix1, ix2, ix3) = 0.0;
- (*av)(AvPrms, ix1, ix2, ix3) = 0.0;
+ (*av)(AvVxx, ix1, ix2, ix3) = c0o1;
+ (*av)(AvVyy, ix1, ix2, ix3) = c0o1;
+ (*av)(AvVzz, ix1, ix2, ix3) = c0o1;
+ (*av)(AvVxy, ix1, ix2, ix3) = c0o1;
+ (*av)(AvVxz, ix1, ix2, ix3) = c0o1;
+ (*av)(AvVyz, ix1, ix2, ix3) = c0o1;
+ (*av)(AvPrms, ix1, ix2, ix3) = c0o1;
//////////////////////////////////////////////////////////////////////////
}
}
@@ -122,6 +124,8 @@ void AverageValuesSimulationObserver::resetDataRMS(real step)
//////////////////////////////////////////////////////////////////////////
void AverageValuesSimulationObserver::resetDataMeans(real step)
{
+ using namespace vf::basics::constant;
+
resetStepMeans = (int)step;
for (int level = minInitLevel; level <= maxInitLevel; level++) {
@@ -147,10 +151,10 @@ void AverageValuesSimulationObserver::resetDataMeans(real step)
//////////////////////////////////////////////////////////////////////////
// compute average values
//////////////////////////////////////////////////////////////////////////
- (*av)(AvVx, ix1, ix2, ix3) = 0.0;
- (*av)(AvVy, ix1, ix2, ix3) = 0.0;
- (*av)(AvVz, ix1, ix2, ix3) = 0.0;
- (*av)(AvP, ix1, ix2, ix3) = 0.0;
+ (*av)(AvVx, ix1, ix2, ix3) = c0o1;
+ (*av)(AvVy, ix1, ix2, ix3) = c0o1;
+ (*av)(AvVz, ix1, ix2, ix3) = c0o1;
+ (*av)(AvP, ix1, ix2, ix3) = c0o1;
//////////////////////////////////////////////////////////////////////////
}
}
@@ -333,6 +337,7 @@ void AverageValuesSimulationObserver::addData(const SPtr<Block3D> block)
void AverageValuesSimulationObserver::calculateAverageValues(real timeStep)
{
using namespace D3Q27System;
+ using namespace vf::basics::constant;
// Funktionszeiger
calcMacros = NULL;
@@ -391,50 +396,50 @@ void AverageValuesSimulationObserver::calculateAverageValues(real timeStep)
// mean velocity
(*av)(AvVx, ix1, ix2, ix3) =
((*av)(AvVx, ix1, ix2, ix3) * timeStepAfterResetMeans + vx) /
- (timeStepAfterResetMeans + 1.0);
+ (timeStepAfterResetMeans + c1o1);
(*av)(AvVy, ix1, ix2, ix3) =
((*av)(AvVy, ix1, ix2, ix3) * timeStepAfterResetMeans + vy) /
- (timeStepAfterResetMeans + 1.0);
+ (timeStepAfterResetMeans + c1o1);
(*av)(AvVz, ix1, ix2, ix3) =
((*av)(AvVz, ix1, ix2, ix3) * timeStepAfterResetMeans + vz) /
- (timeStepAfterResetMeans + 1.0);
+ (timeStepAfterResetMeans + c1o1);
// rms
(*av)(AvVxx, ix1, ix2, ix3) =
((vx - (*av)(AvVx, ix1, ix2, ix3)) * (vx - (*av)(AvVx, ix1, ix2, ix3)) +
(*av)(AvVxx, ix1, ix2, ix3) * timeStepAfterResetRMS) /
- (timeStepAfterResetRMS + 1.0);
+ (timeStepAfterResetRMS + c1o1);
(*av)(AvVyy, ix1, ix2, ix3) =
((vy - (*av)(AvVy, ix1, ix2, ix3)) * (vy - (*av)(AvVy, ix1, ix2, ix3)) +
(*av)(AvVyy, ix1, ix2, ix3) * timeStepAfterResetRMS) /
- (timeStepAfterResetRMS + 1.0);
+ (timeStepAfterResetRMS + c1o1);
(*av)(AvVzz, ix1, ix2, ix3) =
((vz - (*av)(AvVz, ix1, ix2, ix3)) * (vz - (*av)(AvVz, ix1, ix2, ix3)) +
(*av)(AvVzz, ix1, ix2, ix3) * timeStepAfterResetRMS) /
- (timeStepAfterResetRMS + 1.0);
+ (timeStepAfterResetRMS + c1o1);
// cross-correlations
(*av)(AvVxy, ix1, ix2, ix3) =
((vx - (*av)(AvVx, ix1, ix2, ix3)) * (vy - (*av)(AvVy, ix1, ix2, ix3)) +
(*av)(AvVxy, ix1, ix2, ix3) * timeStepAfterResetRMS) /
- (timeStepAfterResetRMS + 1.0);
+ (timeStepAfterResetRMS + c1o1);
(*av)(AvVxz, ix1, ix2, ix3) =
((vx - (*av)(AvVx, ix1, ix2, ix3)) * (vz - (*av)(AvVz, ix1, ix2, ix3)) +
(*av)(AvVxz, ix1, ix2, ix3) * timeStepAfterResetRMS) /
- (timeStepAfterResetRMS + 1.0);
+ (timeStepAfterResetRMS + c1o1);
(*av)(AvVyz, ix1, ix2, ix3) =
((vy - (*av)(AvVy, ix1, ix2, ix3)) * (vz - (*av)(AvVz, ix1, ix2, ix3)) +
(*av)(AvVyz, ix1, ix2, ix3) * timeStepAfterResetRMS) /
- (timeStepAfterResetRMS + 1.0);
+ (timeStepAfterResetRMS + c1o1);
// mean and rms press
(*av)(AvP, ix1, ix2, ix3) =
((*av)(AvP, ix1, ix2, ix3) * timeStepAfterResetMeans + press) /
- (timeStepAfterResetMeans + 1.0);
+ (timeStepAfterResetMeans + c1o1);
(*av)(AvPrms, ix1, ix2, ix3) =
((press - (*av)(AvP, ix1, ix2, ix3)) * (press - (*av)(AvP, ix1, ix2, ix3)) +
(*av)(AvPrms, ix1, ix2, ix3) * timeStepAfterResetRMS) /
- (timeStepAfterResetRMS + 1.0);
+ (timeStepAfterResetRMS + c1o1);
//////////////////////////////////////////////////////////////////////////
}
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateForcesSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateForcesSimulationObserver.cpp
index 8610c5df9e4b56496c3dc3ba1c25fabfd355f294..3b501424f1426b99a341906fbc3d3b638204ce20 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateForcesSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateForcesSimulationObserver.cpp
@@ -103,15 +103,17 @@ void CalculateForcesSimulationObserver::collectData(real step)
//////////////////////////////////////////////////////////////////////////
void CalculateForcesSimulationObserver::calculateForces()
{
- forceX1global = 0.0;
- forceX2global = 0.0;
- forceX3global = 0.0;
+ using namespace vf::basics::constant;
+
+ forceX1global = c0o1;
+ forceX2global = c0o1;
+ forceX3global = c0o1;
for (SPtr<D3Q27Interactor> interactor : interactors) {
for (BcNodeIndicesMap::value_type t : interactor->getBcNodeIndicesMap()) {
- real forceX1 = 0.0;
- real forceX2 = 0.0;
- real forceX3 = 0.0;
+ real forceX1 = c0o1;
+ real forceX2 = c0o1;
+ real forceX3 = c0o1;
SPtr<Block3D> block = t.first;
std::set<std::vector<int>> &transNodeIndicesSet = t.second;
@@ -172,9 +174,9 @@ void CalculateForcesSimulationObserver::calculateForces()
rvalues = comm->gather(values);
if (comm->getProcessID() == comm->getRoot()) {
- forceX1global = 0.0;
- forceX2global = 0.0;
- forceX3global = 0.0;
+ forceX1global = c0o1;
+ forceX2global = c0o1;
+ forceX3global = c0o1;
for (int i = 0; i < (int)rvalues.size(); i += 3) {
forceX1global += rvalues[i];
@@ -187,7 +189,9 @@ void CalculateForcesSimulationObserver::calculateForces()
UbTupleDouble3 CalculateForcesSimulationObserver::getForces(int x1, int x2, int x3, SPtr<DistributionArray3D> distributions,
SPtr<BoundaryConditions> bc)
{
- UbTupleDouble3 force(0.0, 0.0, 0.0);
+ using namespace vf::basics::constant;
+
+ UbTupleDouble3 force(c0o1, c0o1, c0o1);
if (bc) {
// references to tuple "force"
@@ -217,14 +221,16 @@ UbTupleDouble3 CalculateForcesSimulationObserver::getForces(int x1, int x2, int
//////////////////////////////////////////////////////////////////////////
void CalculateForcesSimulationObserver::calculateCoefficients()
{
+ using namespace vf::basics::constant;
+
real F1 = forceX1global;
real F2 = forceX2global;
real F3 = forceX3global;
// return 2*F/(rho*v*v*a);
- C1 = 2.0 * F1 / (v * v * a);
- C2 = 2.0 * F2 / (v * v * a);
- C3 = 2.0 * F3 / (v * v * a);
+ C1 = c2o1 * F1 / (v * v * a);
+ C2 = c2o1 * F2 / (v * v * a);
+ C3 = c2o1 * F3 / (v * v * a);
}
//////////////////////////////////////////////////////////////////////////
void CalculateForcesSimulationObserver::addInteractor(SPtr<D3Q27Interactor> interactor) { interactors.push_back(interactor); }
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateTorqueSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateTorqueSimulationObserver.cpp
index 768fbbb26241edfe5771bf056b6b83be21b02312..09409ccf11dfb38c5b615f05a667701444d6f208 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateTorqueSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/CalculateTorqueSimulationObserver.cpp
@@ -82,9 +82,11 @@ void CalculateTorqueSimulationObserver::collectData( real step )
//////////////////////////////////////////////////////////////////////////
void CalculateTorqueSimulationObserver::calculateForces()
{
- torqueX1global = 0.0;
- torqueX2global = 0.0;
- torqueX3global = 0.0;
+ using namespace vf::basics::constant;
+
+ torqueX1global = c0o1;
+ torqueX2global = c0o1;
+ torqueX3global = c0o1;
for(SPtr<D3Q27Interactor> interactor : interactors)
{
@@ -94,9 +96,9 @@ void CalculateTorqueSimulationObserver::calculateForces()
for(BcNodeIndicesMap::value_type t : interactor->getBcNodeIndicesMap())
{
- real torqueX1 = 0.0;
- real torqueX2 = 0.0;
- real torqueX3 = 0.0;
+ real torqueX1 = c0o1;
+ real torqueX2 = c0o1;
+ real torqueX3 = c0o1;
SPtr<Block3D> block = t.first;
std::set< std::vector<int> >& transNodeIndicesSet = t.second;
@@ -167,9 +169,9 @@ void CalculateTorqueSimulationObserver::calculateForces()
rvalues = comm->gather(values);
if (comm->getProcessID() == comm->getRoot())
{
- torqueX1global = 0.0;
- torqueX2global = 0.0;
- torqueX3global = 0.0;
+ torqueX1global = c0o1;
+ torqueX2global = c0o1;
+ torqueX3global = c0o1;
for (int i = 0; i < (int)rvalues.size(); i+=3)
{
@@ -182,7 +184,9 @@ void CalculateTorqueSimulationObserver::calculateForces()
//////////////////////////////////////////////////////////////////////////
UbTupleDouble3 CalculateTorqueSimulationObserver::getForces(int x1, int x2, int x3, SPtr<DistributionArray3D> distributions, SPtr<BoundaryConditions> bc)
{
- UbTupleDouble3 force(0.0,0.0,0.0);
+ using namespace vf::basics::constant;
+
+ UbTupleDouble3 force(c0o1,c0o1,c0o1);
if(bc)
{
@@ -217,7 +221,7 @@ UbTupleDouble3 CalculateTorqueSimulationObserver::getForces(int x1, int x2, int
UbTupleDouble3 CalculateTorqueSimulationObserver::getForcesFromMoments(int x1, int x2, int x3, SPtr<ILBMKernel> kernel, SPtr<DistributionArray3D> distributions, SPtr<BoundaryConditions> bc, real nx, real ny, real nz)
{
using namespace vf::basics::constant;
- UbTupleDouble3 force(0.0, 0.0, 0.0);
+ UbTupleDouble3 force(c0o1, c0o1, c0o1);
if (bc) {
@@ -252,7 +256,7 @@ UbTupleDouble3 CalculateTorqueSimulationObserver::getForcesFromMoments(int x1, i
UbTupleDouble3 CalculateTorqueSimulationObserver::getForcesFromStressTensor(int x1, int x2, int x3, SPtr<ILBMKernel> kernel, SPtr<DistributionArray3D> distributions, SPtr<BoundaryConditions> bc, real nx, real ny, real nz)
{
using namespace vf::basics::constant;
- UbTupleDouble3 force(0.0, 0.0, 0.0);
+ UbTupleDouble3 force(c0o1, c0o1, c0o1);
if (bc) {
real f[D3Q27System::ENDF + 1];
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/ForceCalculator.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/ForceCalculator.cpp
index 9a39ce11ed15e939e9fc32eaeb15d541675387aa..25e0e30a5e0c0d5789284df71b6f550da6856006 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/ForceCalculator.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/ForceCalculator.cpp
@@ -20,9 +20,11 @@ ForceCalculator::~ForceCalculator() = default;
Vector3D ForceCalculator::getForces(int x1, int x2, int x3, SPtr<DistributionArray3D> distributions,
SPtr<BoundaryConditions> bc, const Vector3D &boundaryVelocity) const
{
- real forceX1 = 0;
- real forceX2 = 0;
- real forceX3 = 0;
+ using namespace vf::basics::constant;
+
+ real forceX1 = c0o1;
+ real forceX2 = c0o1;
+ real forceX3 = c0o1;
if (bc) {
for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++) {
if (bc->hasNoSlipBoundaryFlag(fdir) || bc->hasVelocityBoundaryFlag(fdir)) {
@@ -31,7 +33,7 @@ Vector3D ForceCalculator::getForces(int x1, int x2, int x3, SPtr<DistributionArr
const real fnbr = distributions->getDistributionInvForDirection(
x1 + D3Q27System::DX1[invDir], x2 + D3Q27System::DX2[invDir], x3 + D3Q27System::DX3[invDir], fdir);
- real correction[3] = { 0.0, 0.0, 0.0 };
+ real correction[3] = { c0o1, c0o1, c0o1 };
if (bc->hasVelocityBoundaryFlag(fdir)) {
const real forceTerm = f - fnbr;
correction[0] = forceTerm * boundaryVelocity[0];
@@ -54,15 +56,17 @@ Vector3D ForceCalculator::getForces(int x1, int x2, int x3, SPtr<DistributionArr
void ForceCalculator::calculateForces(std::vector<SPtr<D3Q27Interactor>> interactors)
{
- forceX1global = 0.0;
- forceX2global = 0.0;
- forceX3global = 0.0;
+ using namespace vf::basics::constant;
+
+ forceX1global = c0o1;
+ forceX2global = c0o1;
+ forceX3global = c0o1;
for (const auto &interactor : interactors) {
for (const auto &t : interactor->getBcNodeIndicesMap()) {
- real forceX1 = 0.0;
- real forceX2 = 0.0;
- real forceX3 = 0.0;
+ real forceX1 = c0o1;
+ real forceX2 = c0o1;
+ real forceX3 = c0o1;
SPtr<Block3D> block = t.first;
SPtr<ILBMKernel> kernel = block->getKernel();
@@ -104,6 +108,8 @@ void ForceCalculator::calculateForces(std::vector<SPtr<D3Q27Interactor>> interac
void ForceCalculator::gatherGlobalForces()
{
+ using namespace vf::basics::constant;
+
std::vector<real>
values; // intel compiler 17 dasn't support this { forceX1global , forceX2global, forceX3global };
values.push_back(forceX1global);
@@ -112,9 +118,9 @@ void ForceCalculator::gatherGlobalForces()
std::vector<real> rvalues = comm->gather(values);
if (comm->isRoot()) {
- forceX1global = 0.0;
- forceX2global = 0.0;
- forceX3global = 0.0;
+ forceX1global = c0o1;
+ forceX2global = c0o1;
+ forceX3global = c0o1;
for (int i = 0; i < (int)rvalues.size(); i += 3) {
forceX1global += rvalues[i];
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp
index 7eabcd2849f2fca11cb057357492fa1062c46dce..89938deb91f5c2a13cf04d073a1589fbe65ba6c2 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp
@@ -36,6 +36,8 @@ IntegrateValuesHelper::~IntegrateValuesHelper() = default;
//////////////////////////////////////////////////////////////////////////
void IntegrateValuesHelper::init(int level)
{
+ using namespace vf::basics::constant;
+
root = comm->isRoot();
real orgX1, orgX2, orgX3;
@@ -49,8 +51,8 @@ void IntegrateValuesHelper::init(int level)
maxInitLevel = level;
}
- real numSolids = 0.0;
- real numFluids = 0.0;
+ real numSolids = c0o1;
+ real numFluids = c0o1;
for (int level_it = minInitLevel; level_it <= maxInitLevel; level_it++) {
std::vector<SPtr<Block3D>> blockVector;
grid->getBlocks(level_it, gridRank, blockVector);
@@ -101,8 +103,8 @@ void IntegrateValuesHelper::init(int level)
rvalues = comm->gather(values);
if (root) {
- numberOfSolidNodes = 0.0;
- numberOfFluidsNodes = 0.0;
+ numberOfSolidNodes = c0o1;
+ numberOfFluidsNodes = c0o1;
int rsize = (int)rvalues.size();
int vsize = (int)values.size();
for (int i = 0; i < rsize; i += vsize) {
@@ -230,21 +232,23 @@ void IntegrateValuesHelper::calculateMQ()
//////////////////////////////////////////////////////////////////////////
void IntegrateValuesHelper::clearData()
{
- sRho = 0.0;
- sVx1 = 0.0;
- sVx2 = 0.0;
- sVx3 = 0.0;
- sCellVolume = 0.0;
+ using namespace vf::basics::constant;
+
+ sRho = c0o1;
+ sVx1 = c0o1;
+ sVx2 = c0o1;
+ sVx3 = c0o1;
+ sCellVolume = c0o1;
// sVm = 0.0;
// sPress = 0.0;
// numberOfFluidsNodes = 0.0;
- sAvVx1 = 0.0;
- sAvVx2 = 0.0;
- sAvVx3 = 0.0;
- sTSx1 = 0.0;
- sTSx2 = 0.0;
- sTSx3 = 0.0;
- sTSx1x3 = 0.0;
+ sAvVx1 = c0o1;
+ sAvVx2 = c0o1;
+ sAvVx3 = c0o1;
+ sTSx1 = c0o1;
+ sTSx2 = c0o1;
+ sTSx3 = c0o1;
+ sTSx1x3 = c0o1;
}
//////////////////////////////////////////////////////////////////////////
real IntegrateValuesHelper::getNumberOfFluidsNodes() { return this->numberOfFluidsNodes; }
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/PressureDifferenceSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/PressureDifferenceSimulationObserver.cpp
index 9b3c63f407b9fac00de6177a369fec2cb3e74a82..ac786be093d161a4da1450d5f7ffca8e4c7e5358 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/PressureDifferenceSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/PressureDifferenceSimulationObserver.cpp
@@ -22,6 +22,8 @@ PressureDifferenceSimulationObserver::PressureDifferenceSimulationObserver(SPtr<
: SimulationObserver(grid, s), path(path), h1(h1), h2(h2), comm(comm)
{
+ using namespace vf::basics::constant;
+
if (comm->getProcessID() == comm->getRoot()) {
std::ofstream ostr;
std::string fname = path;
@@ -64,7 +66,7 @@ PressureDifferenceSimulationObserver::PressureDifferenceSimulationObserver(SPtr<
ostr << std::endl;
ostr.close();
- factor1 = (1.0 / 3.0) * rhoReal * (uReal / uLB) * (uReal / uLB);
+ factor1 = (c1o1 / c3o1) * rhoReal * (uReal / uLB) * (uReal / uLB);
factor2 = rhoReal * (uReal / uLB) * (uReal / uLB);
}
}
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/QCriterionSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/QCriterionSimulationObserver.cpp
index 010d9ff664e22519ceb169c549ecc05307655ed4..e1181e9df12ea737779c543b7cba5011c3027442 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/QCriterionSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/QCriterionSimulationObserver.cpp
@@ -89,6 +89,8 @@ void QCriterionSimulationObserver::clearData()
//////////////////////////////////////////////////////////////////////////
void QCriterionSimulationObserver::addData(const SPtr<Block3D> block)
{
+ using namespace vf::basics::constant;
+
UbTupleDouble3 org = grid->getBlockWorldCoordinates(block);
// UbTupleDouble3 blockLengths = grid->getBlockLengths(block);
UbTupleDouble3 nodeOffset = grid->getNodeOffset(block);
@@ -147,16 +149,16 @@ void QCriterionSimulationObserver::addData(const SPtr<Block3D> block)
getNeighborVelocities(0, 0, 1, ix1, ix2, ix3, block, vT, vB);
//////////////////////////////////
// derivatives
- real duxdy = (vN[xdir] - vS[xdir]) * 0.5;
- real duydx = (vE[ydir] - vW[ydir]) * 0.5;
- real duxdz = (vT[xdir] - vB[xdir]) * 0.5;
- real duzdx = (vE[zdir] - vW[zdir]) * 0.5;
- real duydz = (vT[ydir] - vB[ydir]) * 0.5;
- real duzdy = (vN[zdir] - vS[zdir]) * 0.5;
+ real duxdy = (vN[xdir] - vS[xdir]) * c1o2;
+ real duydx = (vE[ydir] - vW[ydir]) * c1o2;
+ real duxdz = (vT[xdir] - vB[xdir]) * c1o2;
+ real duzdx = (vE[zdir] - vW[zdir]) * c1o2;
+ real duydz = (vT[ydir] - vB[ydir]) * c1o2;
+ real duzdy = (vN[zdir] - vS[zdir]) * c1o2;
- real duxdx = (vE[xdir] - vW[xdir]) * 0.5;
- real duydy = (vN[ydir] - vS[ydir]) * 0.5;
- real duzdz = (vT[zdir] - vB[zdir]) * 0.5;
+ real duxdx = (vE[xdir] - vW[xdir]) * c1o2;
+ real duydy = (vN[ydir] - vS[ydir]) * c1o2;
+ real duzdz = (vT[zdir] - vB[zdir]) * c1o2;
real scaleFactor =
(real)(1
@@ -203,7 +205,9 @@ void QCriterionSimulationObserver::addData(const SPtr<Block3D> block)
void QCriterionSimulationObserver::getNeighborVelocities(int offx, int offy, int offz, int ix1, int ix2, int ix3,
const SPtr<Block3D> block, real *vE, real *vW)
{
- SPtr<ILBMKernel> kernel = block->getKernel();
+ using namespace vf::basics::constant;
+
+ SPtr<ILBMKernel> kernel = block->getKernel();
SPtr<BCArray3D> bcArray = kernel->getBCSet()->getBCArray();
SPtr<DistributionArray3D> distributions = kernel->getDataSet()->getFdistributions();
@@ -306,9 +310,9 @@ void QCriterionSimulationObserver::getNeighborVelocities(int offx, int offy, int
computeVelocity(fW2, vW2, compressible);
computeVelocity(f0, v0, compressible);
// second order non-symetric interpolation
- vW[0] = v0[0] * 1.5 - vW[0] + 0.5 * vW2[0];
- vW[1] = v0[1] * 1.5 - vW[1] + 0.5 * vW2[1];
- vW[2] = v0[2] * 1.5 - vW[2] + 0.5 * vW2[2];
+ vW[0] = v0[0] * c3o2 - vW[0] + c1o2 * vW2[0];
+ vW[1] = v0[1] * c3o2 - vW[1] + c1o2 * vW2[1];
+ vW[2] = v0[2] * c3o2 - vW[2] + c1o2 * vW2[2];
// throw UbException(UB_EXARGS,"Parallel or Non-Uniform Simulation -- not yet implemented");
} else {
SPtr<ILBMKernel> kernelW = blockNeighW->getKernel();
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/ShearStressSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/ShearStressSimulationObserver.cpp
index 92c8f5f60344019cff472851104b86e5838302a8..0a754b6d21864c53bba10b3698e9314d3b7fdb8e 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/ShearStressSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/ShearStressSimulationObserver.cpp
@@ -16,6 +16,8 @@ ShearStressSimulationObserver::ShearStressSimulationObserver(SPtr<Grid3D> grid,
SPtr<UbScheduler> s, SPtr<UbScheduler> rs)
: SimulationObserver(grid, s), Resetscheduler(rs), path(path), writer(writer)
{
+ using namespace vf::basics::constant;
+
std::shared_ptr<vf::mpi::Communicator> comm = vf::mpi::Communicator::getInstance();
normals.push_back(0);
normals.push_back(0);
@@ -30,7 +32,7 @@ ShearStressSimulationObserver::ShearStressSimulationObserver(SPtr<Grid3D> grid,
for (SPtr<Block3D> block : blockVector[level]) {
UbTupleInt3 nx = grid->getBlockNX();
SPtr<ShearStressValuesArray3D> shearStressValues = SPtr<ShearStressValuesArray3D>(
- new ShearStressValuesArray3D(14, val<1>(nx) + 1, val<2>(nx) + 1, val<3>(nx) + 1, 0.0));
+ new ShearStressValuesArray3D(14, val<1>(nx) + 1, val<2>(nx) + 1, val<3>(nx) + 1, c0o1));
block->getKernel()->getDataSet()->setShearStressValues(shearStressValues);
}
}
@@ -126,6 +128,7 @@ void ShearStressSimulationObserver::calculateShearStress(real timeStep)
{
using namespace vf::lbm::dir;
using namespace D3Q27System;
+ using namespace vf::basics::constant;
real f[27];
real vx, vy, vz, sxx, syy, szz, sxy, syz, sxz;
@@ -182,23 +185,23 @@ void ShearStressSimulationObserver::calculateShearStress(real timeStep)
vz = ((((f[DIR_PPP] - f[DIR_MMM]) + (f[DIR_PMP] - f[DIR_MPM])) + ((f[DIR_MPP] - f[DIR_PMM]) + (f[DIR_MMP] - f[DIR_PPM]))) +
(((f[DIR_0MP] - f[DIR_0PM]) + (f[DIR_0PP] - f[DIR_0MM])) + ((f[DIR_M0P] - f[DIR_P0M]) + (f[DIR_P0P] - f[DIR_M0M]))) + (f[DIR_00P] - f[DIR_00M]));
- sxy = 3.0 * collFactor / (collFactor - 1.0) *
+ sxy = c3o1 * collFactor / (collFactor - c1o1) *
(((f[DIR_PPP] + f[DIR_MMM]) - (f[DIR_PMP] + f[DIR_MPM])) + (-(f[DIR_PMM] + f[DIR_MPP]) + (f[DIR_MMP] + f[DIR_PPM])) +
(((f[DIR_PP0] + f[DIR_MM0]) - (f[DIR_PM0] + f[DIR_MP0]))) - vx * vy);
- sxz = 3.0 * collFactor / (collFactor - 1.0) *
+ sxz = c3o1 * collFactor / (collFactor - c1o1) *
(((f[DIR_PPP] + f[DIR_MMM]) + (f[DIR_PMP] + f[DIR_MPM])) + (-(f[DIR_PMM] + f[DIR_MPP]) - (f[DIR_MMP] + f[DIR_PPM])) +
((f[DIR_P0P] + f[DIR_M0M]) - (f[DIR_P0M] + f[DIR_M0P])) - vx * vz);
- syz = 3.0 * collFactor / (collFactor - 1.0) *
+ syz = c3o1 * collFactor / (collFactor - c1o1) *
(((f[DIR_PPP] + f[DIR_MMM]) - (f[DIR_PMP] + f[DIR_MPM])) + ((f[DIR_PMM] + f[DIR_MPP]) - (f[DIR_MMP] + f[DIR_PPM])) +
(-(f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM])) - vy * vz);
- real dxxMyy = 3.0 / 2.0 * collFactor / (collFactor - 1.0) *
+ real dxxMyy = c3o1 / c2o1 * collFactor / (collFactor - c1o1) *
(((f[DIR_P0P] + f[DIR_M0M]) + (f[DIR_P0M] + f[DIR_M0P])) - ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM])) +
((f[DIR_P00] + f[DIR_M00]) - (f[DIR_0P0] + f[DIR_0M0])) - vx * vx + vy * vy);
- real dxxMzz = 3.0 / 2.0 * collFactor / (collFactor - 1.0) *
+ real dxxMzz = c3o1 / c2o1 * collFactor / (collFactor - c1o1) *
((((f[DIR_PP0] + f[DIR_MM0]) + (f[DIR_PM0] + f[DIR_MP0])) - ((f[DIR_0PM] + f[DIR_0MP]) + (f[DIR_0PP] + f[DIR_0MM]))) +
((f[DIR_P00] + f[DIR_M00]) - (f[DIR_00P] + f[DIR_00M])) - vx * vx + vz * vz);
@@ -206,25 +209,25 @@ void ShearStressSimulationObserver::calculateShearStress(real timeStep)
// f[NW]))-((f[TE] + f[BW])+(f[BE]+ f[TW])))
// +((f[N] + f[S])-(f[T] + f[B])) -vy*vy +vz*vz);
- sxx = (dxxMyy + dxxMzz) / 3.0; // weil dxxPyyPzz=0
+ sxx = (dxxMyy + dxxMzz) / c3o1; // weil dxxPyyPzz=0
- syy = (dxxMzz - 2 * dxxMyy) / 3.0;
+ syy = (dxxMzz - c2o1 * dxxMyy) / c3o1;
- szz = (dxxMyy - 2 * dxxMzz) / 3.0;
+ szz = (dxxMyy - c2o1 * dxxMzz) / c3o1;
//////////////////////////////////////////////////////////////////////////
// compute average values
//////////////////////////////////////////////////////////////////////////
- (*ssv)(AvVx, ix1, ix2, ix3) = ((*ssv)(AvVx, ix1, ix2, ix3) * timeStep + vx) / (timeStep + 1.0);
- (*ssv)(AvVy, ix1, ix2, ix3) = ((*ssv)(AvVy, ix1, ix2, ix3) * timeStep + vy) / (timeStep + 1.0);
- (*ssv)(AvVz, ix1, ix2, ix3) = ((*ssv)(AvVz, ix1, ix2, ix3) * timeStep + vz) / (timeStep + 1.0);
-
- (*ssv)(AvSxx, ix1, ix2, ix3) = ((*ssv)(AvSxx, ix1, ix2, ix3) * timeStep + sxx) / (timeStep + 1.0);
- (*ssv)(AvSyy, ix1, ix2, ix3) = ((*ssv)(AvSyy, ix1, ix2, ix3) * timeStep + syy) / (timeStep + 1.0);
- (*ssv)(AvSzz, ix1, ix2, ix3) = ((*ssv)(AvSzz, ix1, ix2, ix3) * timeStep + szz) / (timeStep + 1.0);
- (*ssv)(AvSxy, ix1, ix2, ix3) = ((*ssv)(AvSxy, ix1, ix2, ix3) * timeStep + sxy) / (timeStep + 1.0);
- (*ssv)(AvSyz, ix1, ix2, ix3) = ((*ssv)(AvSyz, ix1, ix2, ix3) * timeStep + syz) / (timeStep + 1.0);
- (*ssv)(AvSxz, ix1, ix2, ix3) = ((*ssv)(AvSxz, ix1, ix2, ix3) * timeStep + sxz) / (timeStep + 1.0);
+ (*ssv)(AvVx, ix1, ix2, ix3) = ((*ssv)(AvVx, ix1, ix2, ix3) * timeStep + vx) / (timeStep + c1o1);
+ (*ssv)(AvVy, ix1, ix2, ix3) = ((*ssv)(AvVy, ix1, ix2, ix3) * timeStep + vy) / (timeStep + c1o1);
+ (*ssv)(AvVz, ix1, ix2, ix3) = ((*ssv)(AvVz, ix1, ix2, ix3) * timeStep + vz) / (timeStep + c1o1);
+
+ (*ssv)(AvSxx, ix1, ix2, ix3) = ((*ssv)(AvSxx, ix1, ix2, ix3) * timeStep + sxx) / (timeStep + c1o1);
+ (*ssv)(AvSyy, ix1, ix2, ix3) = ((*ssv)(AvSyy, ix1, ix2, ix3) * timeStep + syy) / (timeStep + c1o1);
+ (*ssv)(AvSzz, ix1, ix2, ix3) = ((*ssv)(AvSzz, ix1, ix2, ix3) * timeStep + szz) / (timeStep + c1o1);
+ (*ssv)(AvSxy, ix1, ix2, ix3) = ((*ssv)(AvSxy, ix1, ix2, ix3) * timeStep + sxy) / (timeStep + c1o1);
+ (*ssv)(AvSyz, ix1, ix2, ix3) = ((*ssv)(AvSyz, ix1, ix2, ix3) * timeStep + syz) / (timeStep + c1o1);
+ (*ssv)(AvSxz, ix1, ix2, ix3) = ((*ssv)(AvSxz, ix1, ix2, ix3) * timeStep + sxz) / (timeStep + c1o1);
}
}
}
@@ -233,6 +236,8 @@ void ShearStressSimulationObserver::calculateShearStress(real timeStep)
//////////////////////////////////////////////////////////////////////////
void ShearStressSimulationObserver::addData()
{
+ using namespace vf::basics::constant;
+
// Diese Daten werden geschrieben:
datanames.resize(0);
datanames.emplace_back("y^plus");
@@ -322,16 +327,16 @@ void ShearStressSimulationObserver::addData()
real nvty = vty / normVt;
real nvtz = vtz / normVt;
- real sx = 0.5 * ((*ssv)(AvSxx, ix1, ix2, ix3) * nvtx + (*ssv)(AvSxy, ix1, ix2, ix3) * nvty +
+ real sx = c1o2 * ((*ssv)(AvSxx, ix1, ix2, ix3) * nvtx + (*ssv)(AvSxy, ix1, ix2, ix3) * nvty +
(*ssv)(AvSxz, ix1, ix2, ix3) * nvtz);
- real sy = 0.5 * ((*ssv)(AvSxy, ix1, ix2, ix3) * nvtx + (*ssv)(AvSyy, ix1, ix2, ix3) * nvty +
+ real sy = c1o2 * ((*ssv)(AvSxy, ix1, ix2, ix3) * nvtx + (*ssv)(AvSyy, ix1, ix2, ix3) * nvty +
(*ssv)(AvSyz, ix1, ix2, ix3) * nvtz);
- real sz = 0.5 * ((*ssv)(AvSxz, ix1, ix2, ix3) * nvtx + (*ssv)(AvSyz, ix1, ix2, ix3) * nvty +
+ real sz = c1o2 * ((*ssv)(AvSxz, ix1, ix2, ix3) * nvtx + (*ssv)(AvSyz, ix1, ix2, ix3) * nvty +
(*ssv)(AvSzz, ix1, ix2, ix3) * nvtz);
real sabs = sqrt(sx * sx + sy * sy + sz * sz);
- real viscosity = (1.0 / 3.0) * (1.0 / collFactor - 0.5);
- real rho = 1.0;
+ real viscosity = (c1o1 / c3o1) * (c1o1 / collFactor - c1o2);
+ real rho = c1o1;
real utau = sqrt(viscosity / rho * sabs);
// double q=(*av)(ix1,ix2,ix3,normalq) ;
@@ -355,6 +360,8 @@ void ShearStressSimulationObserver::reset(real step)
//////////////////////////////////////////////////////////////////////////
void ShearStressSimulationObserver::resetData(real /*step*/)
{
+ using namespace vf::basics::constant;
+
for (int level = minInitLevel; level <= maxInitLevel; level++) {
for (const auto &block : blockVector[level]) {
if (block) {
@@ -383,16 +390,16 @@ void ShearStressSimulationObserver::resetData(real /*step*/)
//////////////////////////////////////////////////////////////////////////
// compute average values
//////////////////////////////////////////////////////////////////////////
- (*ssv)(AvVx, ix1, ix2, ix3) = 0.0;
- (*ssv)(AvVy, ix1, ix2, ix3) = 0.0;
- (*ssv)(AvVz, ix1, ix2, ix3) = 0.0;
-
- (*ssv)(AvSxx, ix1, ix2, ix3) = 0.0;
- (*ssv)(AvSyy, ix1, ix2, ix3) = 0.0;
- (*ssv)(AvSzz, ix1, ix2, ix3) = 0.0;
- (*ssv)(AvSxy, ix1, ix2, ix3) = 0.0;
- (*ssv)(AvSyz, ix1, ix2, ix3) = 0.0;
- (*ssv)(AvSxz, ix1, ix2, ix3) = 0.0;
+ (*ssv)(AvVx, ix1, ix2, ix3) = c0o1;
+ (*ssv)(AvVy, ix1, ix2, ix3) = c0o1;
+ (*ssv)(AvVz, ix1, ix2, ix3) = c0o1;
+
+ (*ssv)(AvSxx, ix1, ix2, ix3) = c0o1;
+ (*ssv)(AvSyy, ix1, ix2, ix3) = c0o1;
+ (*ssv)(AvSzz, ix1, ix2, ix3) = c0o1;
+ (*ssv)(AvSxy, ix1, ix2, ix3) = c0o1;
+ (*ssv)(AvSyz, ix1, ix2, ix3) = c0o1;
+ (*ssv)(AvSxz, ix1, ix2, ix3) = c0o1;
//////////////////////////////////////////////////////////////////////////
}
}
@@ -409,10 +416,11 @@ void ShearStressSimulationObserver::findPlane(int ix1, int ix2, int ix3, SPtr<Gr
real &B, real &C, real &D, real &ii)
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
- real x1plane = 0.0, y1plane = 0.0, z1plane = 0.0;
- real x2plane = 0.0, y2plane = 0.0, z2plane = 0.0;
- real x3plane = 0.0, y3plane = 0.0, z3plane = 0.0;
+ real x1plane = c0o1, y1plane = c0o1, z1plane = c0o1;
+ real x2plane = c0o1, y2plane = c0o1, z2plane = c0o1;
+ real x3plane = c0o1, y3plane = c0o1, z3plane = c0o1;
SPtr<BoundaryConditions> bcPtr;
real dx = grid->getDeltaX(block);
SPtr<ILBMKernel> kernel = block->getKernel();
@@ -839,6 +847,7 @@ bool ShearStressSimulationObserver::checkUndefindedNodes(SPtr<BCArray3D> bcArray
void ShearStressSimulationObserver::initDistance()
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
for (const auto &interactor : interactors) {
// typedef std::map<SPtr<Block3D>, std::set< std::vector<int> > > TransNodeIndicesMap;
@@ -906,7 +915,7 @@ void ShearStressSimulationObserver::initDistance()
continue;
//////get normal and distance//////
- real A, B, C, D, ii = 0.0;
+ real A, B, C, D, ii = c0o1;
findPlane(ix1, ix2, ix3, grid, block, A, B, C, D, ii);
Vector3D pointplane1 = grid->getNodeCoordinates(block, ix1, ix2, ix3);
real ix1ph = pointplane1[0];
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/TimeAveragedValuesSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/TimeAveragedValuesSimulationObserver.cpp
index ebd65f625600a1c68f48d00c33a79976ea6d1a5a..e71e0a1f211ef7703010f24ec5bac99b78233110 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/TimeAveragedValuesSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/TimeAveragedValuesSimulationObserver.cpp
@@ -354,6 +354,8 @@ void TimeAveragedValuesSimulationObserver::addData(const SPtr<Block3D> block)
//////////////////////////////////////////////////////////////////////////
void TimeAveragedValuesSimulationObserver::calculateAverageValues(real timeSteps)
{
+ using namespace vf::basics::constant;
+
for (int level = minInitLevel; level <= maxInitLevel; level++) {
int i;
const int block_size = (int) blockVector[level].size();
@@ -433,25 +435,25 @@ void TimeAveragedValuesSimulationObserver::calculateAverageValues(real timeSteps
if ((options & Triplecorrelations) == Triplecorrelations) {
// mean triple-correlations
(*at)(Vxxx, ix1, ix2, ix3) =
- (*at)(Vxxx, ix1, ix2, ix3) / timeSteps - 3.0 * uxx * ux + 2.0 * ux * ux * ux;
+ (*at)(Vxxx, ix1, ix2, ix3) / timeSteps - c3o1 * uxx * ux + c2o1 * ux * ux * ux;
(*at)(Vxxy, ix1, ix2, ix3) = (*at)(Vxxy, ix1, ix2, ix3) / timeSteps -
- 2.0 * uxy * ux - uxx * uy + 2.0 * ux * ux * uy;
+ c2o1 * uxy * ux - uxx * uy + c2o1 * ux * ux * uy;
(*at)(Vxxz, ix1, ix2, ix3) = (*at)(Vxxz, ix1, ix2, ix3) / timeSteps -
- 2.0 * uxz * ux - uxx * uz + 2.0 * ux * ux * uz;
+ c2o1 * uxz * ux - uxx * uz + c2o1 * ux * ux * uz;
(*at)(Vyyy, ix1, ix2, ix3) =
- (*at)(Vyyy, ix1, ix2, ix3) / timeSteps - 3.0 * uyy * uy + 2.0 * uy * uy * uy;
+ (*at)(Vyyy, ix1, ix2, ix3) / timeSteps - c3o1 * uyy * uy + c2o1 * uy * uy * uy;
(*at)(Vyyx, ix1, ix2, ix3) = (*at)(Vyyx, ix1, ix2, ix3) / timeSteps -
- 2.0 * uxy * uy - uyy * ux + 2.0 * uy * uy * ux;
+ c2o1 * uxy * uy - uyy * ux + c2o1 * uy * uy * ux;
(*at)(Vyyz, ix1, ix2, ix3) = (*at)(Vyyz, ix1, ix2, ix3) / timeSteps -
- 2.0 * uyz * uy - uyy * uz + 2.0 * uy * uy * uz;
+ c2o1 * uyz * uy - uyy * uz + c2o1 * uy * uy * uz;
(*at)(Vzzz, ix1, ix2, ix3) =
- (*at)(Vzzz, ix1, ix2, ix3) / timeSteps - 3.0 * uzz * uz + 2.0 * uz * uz * uz;
+ (*at)(Vzzz, ix1, ix2, ix3) / timeSteps - c3o1 * uzz * uz + c2o1 * uz * uz * uz;
(*at)(Vzzx, ix1, ix2, ix3) = (*at)(Vzzx, ix1, ix2, ix3) / timeSteps -
- 2.0 * uxz * uz - uzz * ux + 2.0 * uz * uz * ux;
+ c2o1 * uxz * uz - uzz * ux + c2o1 * uz * uz * ux;
(*at)(Vzzy, ix1, ix2, ix3) = (*at)(Vzzy, ix1, ix2, ix3) / timeSteps -
- 2.0 * uyz * uz - uzz * uy + 2.0 * uz * uz * uy;
+ c2o1 * uyz * uz - uzz * uy + c2o1 * uz * uz * uy;
(*at)(Vxyz, ix1, ix2, ix3) = (*at)(Vxyz, ix1, ix2, ix3) / timeSteps - uxy * uz -
- uxz * uy - uyz * ux + 2.0 * ux * uy * uz;
+ uxz * uy - uyz * ux + c2o1 * ux * uy * uz;
}
//////////////////////////////////////////////////////////////////////////
}
@@ -577,6 +579,7 @@ void TimeAveragedValuesSimulationObserver::calculateSubtotal(real step)
void TimeAveragedValuesSimulationObserver::planarAverage(real step)
{
std::ofstream ostr;
+ using namespace vf::basics::constant;
if (root) {
int istep = int(step);
@@ -642,7 +645,7 @@ void TimeAveragedValuesSimulationObserver::planarAverage(real step)
if (root) {
real numberOfFluidsNodes = intValHelp.getNumberOfFluidsNodes();
if (numberOfFluidsNodes > 0) {
- ostr << j + 0.5 * dx << std::setprecision(15);
+ ostr << j + c1o2 * dx << std::setprecision(15);
// mean density
if ((options & Density) == Density) {
diff --git a/src/cpu/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp b/src/cpu/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp
index dbfb8907dacdad96849812a1cf5b01ccb52e1483..0c0f84b27f55d194b155688cdcae29c3e8397d83 100644
--- a/src/cpu/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp
+++ b/src/cpu/VirtualFluidsCore/Utilities/ChangeRandomQs.hpp
@@ -27,10 +27,10 @@ namespace Utilities
if (bc->hasNoSlipBoundaryFlag(fdir))
{
const int invDir = D3Q27System::INVDIR[fdir];
- real q = (real) bc->getQ(invDir);
+ float q = bc->getQ(invDir);
//double r = (double)UbRandom::rand(-50, 50);
- real r = (real)UbRandom::rand(-10, 10);
- real q_temp = q + q/r;
+ float r = (float)UbRandom::rand(-10, 10);
+ float q_temp = q + q/r;
if (q_temp < 0.0)
{
q_temp = 0.0001f;
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
index 2e24a2e26d1709b5c1af33d2210f269ca59f2e40..68ebae1d664d33832a6e760daf0835abbb87c75c 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
@@ -43,10 +43,12 @@
InitDistributionsBlockVisitor::InitDistributionsBlockVisitor() : Block3DVisitor(0, D3Q27System::MAXLEVEL)
{
- this->setVx1(0.0);
- this->setVx2(0.0);
- this->setVx3(0.0);
- this->setRho(0.0);
+ using namespace vf::basics::constant;
+
+ this->setVx1(c0o1);
+ this->setVx2(c0o1);
+ this->setVx3(c0o1);
+ this->setRho(c0o1);
}
//////////////////////////////////////////////////////////////////////////
void InitDistributionsBlockVisitor::setVx1(const mu::Parser &parser)
@@ -125,6 +127,7 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
{
using namespace D3Q27System;
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
if(!block) UB_THROW( UbException(UB_EXARGS,"block is not exist") );
@@ -214,38 +217,38 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
real vx2Minusx3 = muVx2.Eval();
real vx3Minusx3 = muVx3.Eval();
- real ax=(vx1Plusx1-vx1Minusx1)/(2.0*deltaX)*dx;
- real bx=(vx2Plusx1-vx2Minusx1)/(2.0*deltaX)*dx;
- real cx=(vx3Plusx1-vx3Minusx1)/(2.0*deltaX)*dx;
+ real ax=(vx1Plusx1-vx1Minusx1)/(c2o1*deltaX)*dx;
+ real bx=(vx2Plusx1-vx2Minusx1)/(c2o1*deltaX)*dx;
+ real cx=(vx3Plusx1-vx3Minusx1)/(c2o1*deltaX)*dx;
- real ay=(vx1Plusx2-vx1Minusx2)/(2.0*deltaX)*dx;
- real by=(vx2Plusx2-vx2Minusx2)/(2.0*deltaX)*dx;
- real cy=(vx3Plusx2-vx3Minusx2)/(2.0*deltaX)*dx;
+ real ay=(vx1Plusx2-vx1Minusx2)/(c2o1*deltaX)*dx;
+ real by=(vx2Plusx2-vx2Minusx2)/(c2o1*deltaX)*dx;
+ real cy=(vx3Plusx2-vx3Minusx2)/(c2o1*deltaX)*dx;
- real az=(vx1Plusx3-vx1Minusx3)/(2.0*deltaX)*dx;
- real bz=(vx2Plusx3-vx2Minusx3)/(2.0*deltaX)*dx;
- real cz=(vx3Plusx3-vx3Minusx3)/(2.0*deltaX)*dx;
- real eps_new=1.0;
- real op = 1.;
+ real az=(vx1Plusx3-vx1Minusx3)/(c2o1*deltaX)*dx;
+ real bz=(vx2Plusx3-vx2Minusx3)/(c2o1*deltaX)*dx;
+ real cz=(vx3Plusx3-vx3Minusx3)/(c2o1*deltaX)*dx;
+ real eps_new=c1o1;
+ real op = c1o1;
real feq[27];
calcFeqsFct(feq,rho,vx1,vx2,vx3);
- real f_E = eps_new *((5.*ax*o + 5.*by*o + 5.*cz*o - 8.*ax*op + 4.*by*op + 4.*cz*op)/(54.*o*op));
- real f_N = f_E + eps_new *((2.*(ax - by))/(9.*o));
- real f_T = f_E + eps_new *((2.*(ax - cz))/(9.*o));
- real f_NE = eps_new *(-(5.*cz*o + 3.*(ay + bx)*op - 2.*cz*op + ax*(5.*o + op) + by*(5.*o + op))/(54.*o*op));
- real f_SE = f_NE + eps_new *(( ay + bx )/(9.*o));
- real f_TE = eps_new *(-(5.*cz*o + by*(5.*o - 2.*op) + 3.*(az + cx)*op + cz*op + ax*(5.*o + op))/(54.*o*op));
- real f_BE = f_TE + eps_new *(( az + cx )/(9.*o));
- real f_TN = eps_new *(-(5.*ax*o + 5.*by*o + 5.*cz*o - 2.*ax*op + by*op + 3.*bz*op + 3.*cy*op + cz*op)/(54.*o*op));
- real f_BN = f_TN + eps_new *(( bz + cy )/(9.*o));
- real f_ZERO = eps_new *((5.*(ax + by + cz))/(9.*op));
- real f_TNE = eps_new *(-(ay + az + bx + bz + cx + cy)/(72.*o));
- real f_TSW = - eps_new *((ay + bx)/(36.*o)) - f_TNE;
- real f_TSE = - eps_new *((az + cx)/(36.*o)) - f_TNE;
- real f_TNW = - eps_new *((bz + cy)/(36.*o)) - f_TNE;
+ real f_E = eps_new *((c5o1*ax*o + c5o1*by*o + c5o1*cz*o - c8o1*ax*op + c4o1*by*op + c4o1*cz*op)/(c54o1*o*op));
+ real f_N = f_E + eps_new *((c2o1*(ax - by))/(c9o1*o));
+ real f_T = f_E + eps_new *((c2o1*(ax - cz))/(c9o1*o));
+ real f_NE = eps_new *(-(c5o1*cz*o + c3o1*(ay + bx)*op - c2o1*cz*op + ax*(c5o1*o + op) + by*(c5o1*o + op))/(c54o1*o*op));
+ real f_SE = f_NE + eps_new *(( ay + bx )/(c9o1*o));
+ real f_TE = eps_new *(-(c5o1*cz*o + by*(c5o1*o - c2o1*op) + c3o1*(az + cx)*op + cz*op + ax*(c5o1*o + op))/(c54o1*o*op));
+ real f_BE = f_TE + eps_new *(( az + cx )/(c9o1*o));
+ real f_TN = eps_new *(-(c5o1*ax*o + c5o1*by*o + c5o1*cz*o - c2o1*ax*op + by*op + c3o1*bz*op + c3o1*cy*op + cz*op)/(c54o1*o*op));
+ real f_BN = f_TN + eps_new *(( bz + cy )/(c9o1*o));
+ real f_ZERO = eps_new *((c5o1*(ax + by + cz))/(c9o1*op));
+ real f_TNE = eps_new *(-(ay + az + bx + bz + cx + cy)/(c72o1*o));
+ real f_TSW = - eps_new *((ay + bx)/(c36o1*o)) - f_TNE;
+ real f_TSE = - eps_new *((az + cx)/(c36o1*o)) - f_TNE;
+ real f_TNW = - eps_new *((bz + cy)/(c36o1*o)) - f_TNE;
f[DIR_P00] = f_E + feq[DIR_P00];
@@ -298,7 +301,9 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
//////////////////////////////////////////////////////////////////////////
void InitDistributionsBlockVisitor::checkFunction(mu::Parser fct)
{
- real x1 = 1.0, x2 = 1.0, x3 = 1.0;
+ using namespace vf::basics::constant;
+
+ real x1 = c1o1, x2 = c1o1, x3 = c1o1;
fct.DefineVar("x1", &x1);
fct.DefineVar("x2", &x2);
fct.DefineVar("x3", &x3);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp
index 660363e22e7c315a40df596aa95137f5589fff72..21d06dffa9f1fb09a1b5d7baa43d72393eeda47b 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsWithInterpolationGridVisitor.cpp
@@ -163,6 +163,8 @@ void InitDistributionsWithInterpolationGridVisitor::copyRemoteBlock(SPtr<Block3D
void InitDistributionsWithInterpolationGridVisitor::interpolateLocalBlockCoarseToFine(SPtr<Block3D> oldBlock,
SPtr<Block3D> newBlock)
{
+ using namespace vf::basics::constant;
+
D3Q27ICell icellC;
D3Q27ICell icellF;
real xoff, yoff, zoff;
@@ -225,14 +227,14 @@ void InitDistributionsWithInterpolationGridVisitor::interpolateLocalBlockCoarseT
iProcessor->interpolateCoarseToFine(icellC, icellF, xoff, yoff, zoff);
} else {
for (int i = 0; i < 27; i++) {
- icellF.BSW[i] = 0.0;
- icellF.BSE[i] = 0.0;
- icellF.BNW[i] = 0.0;
- icellF.BNE[i] = 0.0;
- icellF.TSW[i] = 0.0;
- icellF.TSE[i] = 0.0;
- icellF.TNW[i] = 0.0;
- icellF.TNE[i] = 0.0;
+ icellF.BSW[i] = c0o1;
+ icellF.BSE[i] = c0o1;
+ icellF.BNW[i] = c0o1;
+ icellF.BNE[i] = c0o1;
+ icellF.TSW[i] = c0o1;
+ icellF.TSE[i] = c0o1;
+ icellF.TNW[i] = c0o1;
+ icellF.TNE[i] = c0o1;
}
// std::string err = "For "+oldBlock->toString()+
// " x1="+UbSystem::toString(val<1>(oldGridIndexMin))+
@@ -254,6 +256,8 @@ void InitDistributionsWithInterpolationGridVisitor::interpolateLocalBlockCoarseT
void InitDistributionsWithInterpolationGridVisitor::interpolateRemoteBlockCoarseToFine(SPtr<Block3D> oldBlock,
SPtr<Block3D> newBlock)
{
+ using namespace vf::basics::constant;
+
int newGridRank = newGrid->getRank();
int oldBlockRank = oldBlock->getRank();
int newBlockRank = newBlock->getRank();
@@ -363,14 +367,14 @@ void InitDistributionsWithInterpolationGridVisitor::interpolateRemoteBlockCoarse
iProcessor->interpolateCoarseToFine(icellC, icellF, xoff, yoff, zoff);
} else {
for (int i = 0; i < 27; i++) {
- icellF.BSW[i] = 0.0;
- icellF.BSE[i] = 0.0;
- icellF.BNW[i] = 0.0;
- icellF.BNE[i] = 0.0;
- icellF.TSW[i] = 0.0;
- icellF.TSE[i] = 0.0;
- icellF.TNW[i] = 0.0;
- icellF.TNE[i] = 0.0;
+ icellF.BSW[i] = c0o1;
+ icellF.BSE[i] = c0o1;
+ icellF.BNW[i] = c0o1;
+ icellF.BNE[i] = c0o1;
+ icellF.TSW[i] = c0o1;
+ icellF.TSE[i] = c0o1;
+ icellF.TNW[i] = c0o1;
+ icellF.TNE[i] = c0o1;
}
// std::string err = "For "+oldBlock->toString()+
// " x1="+UbSystem::toString(val<1>(oldGridIndexMin))+
@@ -393,6 +397,8 @@ void InitDistributionsWithInterpolationGridVisitor::interpolateRemoteBlockCoarse
void InitDistributionsWithInterpolationGridVisitor::interpolateLocalBlockFineToCoarse(SPtr<Block3D> oldBlock,
SPtr<Block3D> newBlock)
{
+ using namespace vf::basics::constant;
+
real icellC[27];
D3Q27ICell icellF;
real xoff, yoff, zoff;
@@ -456,14 +462,14 @@ void InitDistributionsWithInterpolationGridVisitor::interpolateLocalBlockFineToC
iProcessor->interpolateFineToCoarse(icellF, icellC, xoff, yoff, zoff);
} else {
for (int i = 0; i < 27; i++) {
- icellF.BSW[i] = 0.0;
- icellF.BSE[i] = 0.0;
- icellF.BNW[i] = 0.0;
- icellF.BNE[i] = 0.0;
- icellF.TSW[i] = 0.0;
- icellF.TSE[i] = 0.0;
- icellF.TNW[i] = 0.0;
- icellF.TNE[i] = 0.0;
+ icellF.BSW[i] = c0o1;
+ icellF.BSE[i] = c0o1;
+ icellF.BNW[i] = c0o1;
+ icellF.BNE[i] = c0o1;
+ icellF.TSW[i] = c0o1;
+ icellF.TSE[i] = c0o1;
+ icellF.TNW[i] = c0o1;
+ icellF.TNE[i] = c0o1;
}
// std::string err = "For "+oldBlock->toString()+
// " x1="+UbSystem::toString(val<1>(oldGridIndexMin))+
@@ -485,6 +491,8 @@ void InitDistributionsWithInterpolationGridVisitor::interpolateLocalBlockFineToC
void InitDistributionsWithInterpolationGridVisitor::interpolateRemoteBlockFineToCoarse(SPtr<Block3D> oldBlock,
SPtr<Block3D> newBlock)
{
+ using namespace vf::basics::constant;
+
int newGridRank = newGrid->getRank();
int oldBlockRank = oldBlock->getRank();
int newBlockRank = newBlock->getRank();
@@ -594,14 +602,14 @@ void InitDistributionsWithInterpolationGridVisitor::interpolateRemoteBlockFineTo
iProcessor->interpolateFineToCoarse(icellF, icellC, xoff, yoff, zoff);
} else {
for (int i = 0; i < 27; i++) {
- icellF.BSW[i] = 0.0;
- icellF.BSE[i] = 0.0;
- icellF.BNW[i] = 0.0;
- icellF.BNE[i] = 0.0;
- icellF.TSW[i] = 0.0;
- icellF.TSE[i] = 0.0;
- icellF.TNW[i] = 0.0;
- icellF.TNE[i] = 0.0;
+ icellF.BSW[i] = c0o1;
+ icellF.BSE[i] = c0o1;
+ icellF.BNW[i] = c0o1;
+ icellF.BNE[i] = c0o1;
+ icellF.TSW[i] = c0o1;
+ icellF.TSE[i] = c0o1;
+ icellF.TNW[i] = c0o1;
+ icellF.TNE[i] = c0o1;
}
// std::string err = "For "+oldBlock->toString()+
// " x1="+UbSystem::toString(val<1>(oldGridIndexMin))+
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
index be0c694bc733ff3b9ebd808d533757f98eefe73c..79c354e5f28724279447a4751ba81c3b6eba0958 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitThixotropyBlockVisitor.cpp
@@ -1,4 +1,4 @@
-//=======================================================================================
+ //=======================================================================================
// ____ ____ __ ______ __________ __ __ __ __
// \ \ | | | | | _ \ |___ ___| | | | | / \ | |
// \ \ | | | | | |_) | | | | | | | / \ | |
@@ -44,6 +44,8 @@
InitThixotropyBlockVisitor::InitThixotropyBlockVisitor()
: Block3DVisitor(0, D3Q27System::MAXLEVEL)
{
+ using namespace vf::basics::constant;
+
//this->setVx1(0.0);
//this->setVx2(0.0);
//this->setVx3(0.0);
@@ -52,7 +54,7 @@ InitThixotropyBlockVisitor::InitThixotropyBlockVisitor()
//this->setf2(0.0);
//this->setf3(0.0);
//this->setConcentration(0.0);
- this->setLambda(0.0);
+ this->setLambda(c0o1);
}
//////////////////////////////////////////////////////////////////////////
//InitThixotropyBlockVisitor::InitThixotropyBlockVisitor(LBMReal lambda /*LBMReal nu, LBMReal D, LBMReal rho, LBMReal vx1, LBMReal vx2, LBMReal vx3, LBMReal c, LBMReal f1, LBMReal f2, LBMReal f3*/)
@@ -218,6 +220,7 @@ void InitThixotropyBlockVisitor::setLambda(real lambda)
void InitThixotropyBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
{
using namespace D3Q27System;
+ using namespace vf::basics::constant;
if(!block) UB_THROW( UbException(UB_EXARGS,"block is not exist") );
@@ -283,7 +286,7 @@ void InitThixotropyBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
real lambda = muLambda.Eval();
- calcFeqsFct(h,lambda,0.0,0.0,0.0);
+ calcFeqsFct(h,lambda,c0o1,c0o1,c0o1);
distributions->setDistribution(h, ix1, ix2, ix3);
distributions->setDistributionInv(h, ix1, ix2, ix3);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
index fed80000c562b8aafdefd83ff2791781a8907f7a..1111429442829adca999bf6c587819419c8bd9fb 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SetKernelBlockVisitor.cpp
@@ -124,6 +124,8 @@ void SetKernelBlockVisitor::throwExceptionIfNotEnoughMemory(const SPtr<Grid3D> &
real SetKernelBlockVisitor::getRequiredPhysicalMemory(const SPtr<Grid3D> &grid) const
{
+ using namespace vf::basics::constant;
+
unsigned long long numberOfNodesPerBlockWithGhostLayer;
auto numberOfBlocks = (unsigned long long)grid->getNumberOfBlocks();
auto blockNx = grid->getBlockNX();
@@ -133,7 +135,7 @@ real SetKernelBlockVisitor::getRequiredPhysicalMemory(const SPtr<Grid3D> &grid)
(val<2>(blockNx) + ghostLayer) * (val<3>(blockNx) + ghostLayer);
auto needMemAll =
- real(numberOfNodesPerBlockWithGhostLayer * (27 * sizeof(real) + sizeof(int) + sizeof(float) * 4));
+ real(numberOfNodesPerBlockWithGhostLayer * (c27o1 * sizeof(real) + sizeof(int) + sizeof(float) * c4o1));
return needMemAll / real(numberOfProcesses);
}
diff --git a/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
index 4cdcfb80bb4aa3119f1e2ca4dfcd19c2abf381dd..520c7cfa30f582438b188643d8f913fd6ce2bd0b 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/SpongeLayerBlockVisitor.cpp
@@ -24,6 +24,7 @@ SpongeLayerBlockVisitor::~SpongeLayerBlockVisitor() = default;
void SpongeLayerBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
if (!boundingBox) {
UB_THROW(UbException(UB_EXARGS, "The bounding box isn't set!"));
@@ -78,22 +79,22 @@ void SpongeLayerBlockVisitor::visit(SPtr<Grid3D> grid, SPtr<Block3D> block)
int ibX1 = block->getX1();
int ibMax = val<1>(ixMax) - val<1>(ixMin) + 1;
real index = (real)(ibX1 - val<1>(ixMin) + 1);
- newCollFactor = oldCollFactor - (oldCollFactor - 1.0) / (real)(ibMax)*index;
+ newCollFactor = oldCollFactor - (oldCollFactor - c1o1) / (real)(ibMax)*index;
} else if (dir == DIR_M00) {
int ibX1 = block->getX1();
int ibMax = val<1>(ixMax) - val<1>(ixMin) + 1;
real index = (real)(ibX1 - val<1>(ixMin) + 1);
- newCollFactor = (oldCollFactor - 1.0) / (real)(ibMax)*index;
+ newCollFactor = (oldCollFactor - c1o1) / (real)(ibMax)*index;
} else if (dir == DIR_00P) {
int ibX3 = block->getX3();
int ibMax = val<3>(ixMax) - val<3>(ixMin) + 1;
real index = (real)(ibX3 - val<3>(ixMin) + 1);
- newCollFactor = oldCollFactor - (oldCollFactor - 1.0) / (real)(ibMax)*index;
+ newCollFactor = oldCollFactor - (oldCollFactor - c1o1) / (real)(ibMax)*index;
} else if (dir == DIR_00M) {
int ibX3 = block->getX3();
int ibMax = val<3>(ixMax) - val<3>(ixMin) + 1;
real index = (real)(ibX3 - val<3>(ixMin) + 1);
- newCollFactor = (oldCollFactor - 1.0) / (real)(ibMax)*index;
+ newCollFactor = (oldCollFactor - c1o1) / (real)(ibMax)*index;
} else
UB_THROW(UbException(UB_EXARGS, "Problem: no orthogonal sponge layer!"));