From 95c9f2ee5722fc0596876d86503bd80c3b25f8a3 Mon Sep 17 00:00:00 2001
From: alena <akaranchuk@list.ru>
Date: Mon, 17 Jul 2023 12:08:46 +0200
Subject: [PATCH] Change digit constants to char
---
.../BoundaryConditions/BCStrategy.cpp | 6 +-
.../NonReflectingInflowBCStrategy.cpp | 112 ++++++------
...lectingOutflowWithRelaxationBCStrategy.cpp | 110 ++++++------
.../VelocityWithDensityBCStrategy.cpp | 4 +-
...onsDoubleGhostLayerFullVectorConnector.cpp | 15 +-
.../ThreeDistributionsFullVectorConnector.cpp | 15 +-
.../Interactors/D3Q27Interactor.cpp | 92 +++++-----
.../D3Q27TriFaceMeshInteractor.cpp | 159 +++++++++---------
.../VirtualFluidsCore/LBM/BGKLBMKernel.cpp | 2 +-
.../LBM/InitDensityLBMKernel.cpp | 2 +-
.../CompressibleOffsetInterpolator.cpp | 6 +-
.../CompressibleOffsetMomentsInterpolator.cpp | 16 +-
.../IncompressibleOffsetInterpolator.cpp | 90 +++++-----
.../IntegrateValuesHelper.cpp | 4 +-
...teBoundaryConditionsSimulationObserver.cpp | 16 +-
.../InitDistributionsBlockVisitor.cpp | 2 +-
16 files changed, 336 insertions(+), 315 deletions(-)
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCStrategy.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCStrategy.cpp
index 3331f2453..5458fa09c 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/BCStrategy.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/BCStrategy.cpp
@@ -54,18 +54,20 @@ void BCStrategy::setBcPointer(SPtr<BoundaryConditions> bcPtr) { this->bcPtr = bc
//////////////////////////////////////////////////////////////////////////
void BCStrategy::setCompressible(bool c)
{
+ using namespace vf::basics::constant;
+
compressible = c;
if (this->compressible) {
calcFeqsForDirFct = &D3Q27System::getCompFeqForDirection;
calcMacrosFct = &D3Q27System::calcCompMacroscopicValues;
calcFeqFct = &D3Q27System::calcCompFeq;
- compressibleFactor = 1.0;
+ compressibleFactor = c1o1;
} else {
calcFeqsForDirFct = &D3Q27System::getIncompFeqForDirection;
calcMacrosFct = &D3Q27System::calcIncompMacroscopicValues;
calcFeqFct = &D3Q27System::calcIncompFeq;
- compressibleFactor = 0.0;
+ compressibleFactor = c0o1;
}
}
//////////////////////////////////////////////////////////////////////////
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingInflowBCStrategy.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingInflowBCStrategy.cpp
index 7f5eaf840..b1cf37e4b 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingInflowBCStrategy.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingInflowBCStrategy.cpp
@@ -98,14 +98,14 @@ void NonReflectingInflowBCStrategy::applyBC()
LBMReal rho, vx1, vx2, vx3;
calcMacrosFct(f, rho, vx1, vx2, vx3);
//vx1 = 0.;
- LBMReal BCVeloWeight = 0.5;
+ LBMReal BCVeloWeight = c1o2;
// LBMReal velocity = 0.004814077025232405;
// LBMReal velocity = 0.00057735;
//LBMReal velocity = 0.04;
// LBMReal velocity = 0.01;
// LBMReal velocity = 1./112.;
// LBMReal velocity = 1./126.;
- LBMReal velocity = 1./200.;
+ LBMReal velocity = c1o100/2;
// LBMReal velocity = 0.005;
//LBMReal delf =(-velocity+vx1)*0.5 ;
LBMReal delf;
@@ -114,15 +114,15 @@ void NonReflectingInflowBCStrategy::applyBC()
case DIR_P00:
delf = (-velocity + vx1) * BCVeloWeight;
// delf = (-velocity ) * BCVeloWeight;
- f[DIR_P00] = ftemp[DIR_P00] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_P00] - delf* WEIGTH[DIR_P00];
- f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PP0]- delf* WEIGTH[DIR_PP0];
- f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PM0]- delf* WEIGTH[DIR_PM0];
- f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_P0P]- delf* WEIGTH[DIR_P0P];
- f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_P0M]- delf* WEIGTH[DIR_P0M];
- f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
- f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
- f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
- f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
+ f[DIR_P00] = ftemp[DIR_P00] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_P00] - delf* WEIGTH[DIR_P00];
+ f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PP0]- delf* WEIGTH[DIR_PP0];
+ f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PM0]- delf* WEIGTH[DIR_PM0];
+ f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_P0P]- delf* WEIGTH[DIR_P0P];
+ f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_P0M]- delf* WEIGTH[DIR_P0M];
+ f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
+ f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
+ f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
+ f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
//f[DIR_P00] = (ftemp[DIR_P00] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_P00]) *
// (1 - BCVeloWeight) +
// (ftemp[DIR_M00] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_M00] +
@@ -181,23 +181,23 @@ void NonReflectingInflowBCStrategy::applyBC()
break;
case DIR_M00:
delf = (-velocity - vx1) * BCVeloWeight;
- f[DIR_M00] = ftemp[DIR_M00] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_M00] -
+ f[DIR_M00] = ftemp[DIR_M00] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_M00] -
delf * WEIGTH[DIR_M00];
- f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MP0] -
+ f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MP0] -
delf * WEIGTH[DIR_MP0];
- f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MM0] -
+ f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MM0] -
delf * WEIGTH[DIR_MM0];
- f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_M0P] -
+ f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_M0P] -
delf * WEIGTH[DIR_M0P];
- f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_M0M] -
+ f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_M0M] -
delf * WEIGTH[DIR_M0M];
- f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MPP] -
+ f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MPP] -
delf * WEIGTH[DIR_MPP];
- f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MMP] -
+ f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MMP] -
delf * WEIGTH[DIR_MMP];
- f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MPM] -
+ f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MPM] -
delf * WEIGTH[DIR_MPM];
- f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MMM] -
+ f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MMM] -
delf * WEIGTH[DIR_MMM];
distributions->setDistributionInvForDirection(f[DIR_M00], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], DIR_P00);
@@ -212,23 +212,23 @@ void NonReflectingInflowBCStrategy::applyBC()
break;
case DIR_0P0:
delf = (-velocity + vx2) * BCVeloWeight;
- f[DIR_0P0] = ftemp[DIR_0P0] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_0P0] -
+ f[DIR_0P0] = ftemp[DIR_0P0] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_0P0] -
delf * WEIGTH[DIR_0P0];
- f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_PP0] -
+ f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_PP0] -
delf * WEIGTH[DIR_PP0];
- f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_MP0] -
+ f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_MP0] -
delf * WEIGTH[DIR_MP0];
- f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_0PP] -
+ f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_0PP] -
delf * WEIGTH[DIR_0PP];
- f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_0PM] -
+ f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_0PM] -
delf * WEIGTH[DIR_0PM];
- f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_PPP] -
+ f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_PPP] -
delf * WEIGTH[DIR_PPP];
- f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_MPP] -
+ f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_MPP] -
delf * WEIGTH[DIR_MPP];
- f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_PPM] -
+ f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_PPM] -
delf * WEIGTH[DIR_PPM];
- f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_MPM] -
+ f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_MPM] -
delf * WEIGTH[DIR_MPM];
distributions->setDistributionInvForDirection(f[DIR_0P0], x1 + DX1[DIR_0M0], x2 + DX2[DIR_0M0], x3 + DX3[DIR_0M0], DIR_0M0);
@@ -243,23 +243,23 @@ void NonReflectingInflowBCStrategy::applyBC()
break;
case DIR_0M0:
delf = (-velocity - vx2) * BCVeloWeight;
- f[DIR_0M0] = ftemp[DIR_0M0] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_0M0] -
+ f[DIR_0M0] = ftemp[DIR_0M0] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_0M0] -
delf * WEIGTH[DIR_0M0];
- f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_PM0] -
+ f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_PM0] -
delf * WEIGTH[DIR_PM0];
- f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_MM0] -
+ f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_MM0] -
delf * WEIGTH[DIR_MM0];
- f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_0MP] -
+ f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_0MP] -
delf * WEIGTH[DIR_0MP];
- f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_0MM] -
+ f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_0MM] -
delf * WEIGTH[DIR_0MM];
- f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_PMP] -
+ f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_PMP] -
delf * WEIGTH[DIR_PMP];
- f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_MMP] -
+ f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_MMP] -
delf * WEIGTH[DIR_MMP];
- f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_PMM] -
+ f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_PMM] -
delf * WEIGTH[DIR_PMM];
- f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_MMM] -
+ f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_MMM] -
delf * WEIGTH[DIR_MMM];
distributions->setDistributionInvForDirection(f[DIR_0M0], x1 + DX1[DIR_0P0], x2 + DX2[DIR_0P0], x3 + DX3[DIR_0P0], DIR_0P0);
@@ -274,23 +274,23 @@ void NonReflectingInflowBCStrategy::applyBC()
break;
case DIR_00P:
delf = (-velocity + vx3) * BCVeloWeight;
- f[DIR_00P] = ftemp[DIR_00P] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_00P] -
+ f[DIR_00P] = ftemp[DIR_00P] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_00P] -
delf * WEIGTH[DIR_00P];
- f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_P0P] -
+ f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_P0P] -
delf * WEIGTH[DIR_P0P];
- f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_M0P] -
+ f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_M0P] -
delf * WEIGTH[DIR_M0P];
- f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_0PP] -
+ f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_0PP] -
delf * WEIGTH[DIR_0PP];
- f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_0MP] -
+ f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_0MP] -
delf * WEIGTH[DIR_0MP];
- f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_PPP] -
+ f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_PPP] -
delf * WEIGTH[DIR_PPP];
- f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_MPP] -
+ f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_MPP] -
delf * WEIGTH[DIR_MPP];
- f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_PMP] -
+ f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_PMP] -
delf * WEIGTH[DIR_PMP];
- f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_MMP] -
+ f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_MMP] -
delf * WEIGTH[DIR_MMP];
distributions->setDistributionInvForDirection(f[DIR_00P], x1 + DX1[DIR_00M], x2 + DX2[DIR_00M], x3 + DX3[DIR_00M], DIR_00M);
@@ -305,23 +305,23 @@ void NonReflectingInflowBCStrategy::applyBC()
break;
case DIR_00M:
delf = (-velocity - vx3) * BCVeloWeight;
- f[DIR_00M] = ftemp[DIR_00M] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_00M] -
+ f[DIR_00M] = ftemp[DIR_00M] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_00M] -
delf * WEIGTH[DIR_00M];
- f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_P0M] -
+ f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_P0M] -
delf * WEIGTH[DIR_P0M];
- f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_M0M] -
+ f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_M0M] -
delf * WEIGTH[DIR_M0M];
- f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_0PM] -
+ f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_0PM] -
delf * WEIGTH[DIR_0PM];
- f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_0MM] -
+ f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_0MM] -
delf * WEIGTH[DIR_0MM];
- f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_PPM] -
+ f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_PPM] -
delf * WEIGTH[DIR_PPM];
- f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_MPM] -
+ f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_MPM] -
delf * WEIGTH[DIR_MPM];
- f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_PMM] -
+ f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_PMM] -
delf * WEIGTH[DIR_PMM];
- f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_MMM] -
+ f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_MMM] -
delf * WEIGTH[DIR_MMM];
distributions->setDistributionInvForDirection(f[DIR_00M], x1 + DX1[DIR_00P], x2 + DX2[DIR_00P], x3 + DX3[DIR_00P], DIR_00P);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowWithRelaxationBCStrategy.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowWithRelaxationBCStrategy.cpp
index 3e10421e5..f0fb4ccae 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowWithRelaxationBCStrategy.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/NonReflectingOutflowWithRelaxationBCStrategy.cpp
@@ -98,18 +98,18 @@ void NonReflectingOutflowWithRelaxationBCStrategy::applyBC()
LBMReal rho, vx1, vx2, vx3;
calcMacrosFct(f, rho, vx1, vx2, vx3);
- LBMReal delf = rho*0.01;
+ LBMReal delf = rho* c1o100;
switch (direction) {
case DIR_P00:
- f[DIR_P00] = ftemp[DIR_P00] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_P00] - delf* WEIGTH[DIR_P00];
- f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PP0]- delf* WEIGTH[DIR_PP0];
- f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PM0]- delf* WEIGTH[DIR_PM0];
- f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_P0P]- delf* WEIGTH[DIR_P0P];
- f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_P0M]- delf* WEIGTH[DIR_P0M];
- f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
- f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
- f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
- f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 + vx1) + (1.0 - one_over_sqrt3 - vx1) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
+ f[DIR_P00] = ftemp[DIR_P00] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_P00] - delf* WEIGTH[DIR_P00];
+ f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PP0]- delf* WEIGTH[DIR_PP0];
+ f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PM0]- delf* WEIGTH[DIR_PM0];
+ f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_P0P]- delf* WEIGTH[DIR_P0P];
+ f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_P0M]- delf* WEIGTH[DIR_P0M];
+ f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
+ f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
+ f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
+ f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 + vx1) + (c1o1 - one_over_sqrt3 - vx1) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
distributions->setDistributionInvForDirection(f[DIR_P00], x1 + DX1[DIR_M00], x2 + DX2[DIR_M00], x3 + DX3[DIR_M00], DIR_M00);
distributions->setDistributionInvForDirection(f[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
@@ -122,15 +122,15 @@ void NonReflectingOutflowWithRelaxationBCStrategy::applyBC()
distributions->setDistributionInvForDirection(f[DIR_PMM], x1 + DX1[DIR_MPP], x2 + DX2[DIR_MPP], x3 + DX3[DIR_MPP], DIR_MPP);
break;
case DIR_M00:
- f[DIR_M00] = ftemp[DIR_M00] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_M00]- delf* WEIGTH[DIR_M00];
- f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MP0]- delf* WEIGTH[DIR_MP0];
- f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MM0]- delf* WEIGTH[DIR_MM0];
- f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_M0P]- delf* WEIGTH[DIR_M0P];
- f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_M0M]- delf* WEIGTH[DIR_M0M];
- f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MPP]- delf* WEIGTH[DIR_MPP];
- f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MMP]- delf* WEIGTH[DIR_MMP];
- f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MPM]- delf* WEIGTH[DIR_MPM];
- f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx1) + (1.0 - one_over_sqrt3 + vx1) * f[DIR_MMM]- delf* WEIGTH[DIR_MMM];
+ f[DIR_M00] = ftemp[DIR_M00] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_M00]- delf* WEIGTH[DIR_M00];
+ f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MP0]- delf* WEIGTH[DIR_MP0];
+ f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MM0]- delf* WEIGTH[DIR_MM0];
+ f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_M0P]- delf* WEIGTH[DIR_M0P];
+ f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_M0M]- delf* WEIGTH[DIR_M0M];
+ f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MPP]- delf* WEIGTH[DIR_MPP];
+ f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MMP]- delf* WEIGTH[DIR_MMP];
+ f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MPM]- delf* WEIGTH[DIR_MPM];
+ f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx1) + (c1o1 - one_over_sqrt3 + vx1) * f[DIR_MMM]- delf* WEIGTH[DIR_MMM];
distributions->setDistributionInvForDirection(f[DIR_M00], x1 + DX1[DIR_P00], x2 + DX2[DIR_P00], x3 + DX3[DIR_P00], DIR_P00);
distributions->setDistributionInvForDirection(f[DIR_MP0], x1 + DX1[DIR_PM0], x2 + DX2[DIR_PM0], x3 + DX3[DIR_PM0], DIR_PM0);
@@ -143,15 +143,15 @@ void NonReflectingOutflowWithRelaxationBCStrategy::applyBC()
distributions->setDistributionInvForDirection(f[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
break;
case DIR_0P0:
- f[DIR_0P0] = ftemp[DIR_0P0] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_0P0]- delf* WEIGTH[DIR_0P0];
- f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_PP0]- delf* WEIGTH[DIR_PP0];
- f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_MP0]- delf* WEIGTH[DIR_MP0];
- f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_0PP]- delf* WEIGTH[DIR_0PP];
- f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_0PM]- delf* WEIGTH[DIR_0PM];
- f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
- f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_MPP]- delf* WEIGTH[DIR_MPP];
- f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
- f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 + vx2) + (1.0 - one_over_sqrt3 - vx2) * f[DIR_MPM]- delf* WEIGTH[DIR_MPM];
+ f[DIR_0P0] = ftemp[DIR_0P0] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_0P0]- delf* WEIGTH[DIR_0P0];
+ f[DIR_PP0] = ftemp[DIR_PP0] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_PP0]- delf* WEIGTH[DIR_PP0];
+ f[DIR_MP0] = ftemp[DIR_MP0] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_MP0]- delf* WEIGTH[DIR_MP0];
+ f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_0PP]- delf* WEIGTH[DIR_0PP];
+ f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_0PM]- delf* WEIGTH[DIR_0PM];
+ f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
+ f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_MPP]- delf* WEIGTH[DIR_MPP];
+ f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
+ f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 + vx2) + (c1o1 - one_over_sqrt3 - vx2) * f[DIR_MPM]- delf* WEIGTH[DIR_MPM];
distributions->setDistributionInvForDirection(f[DIR_0P0], x1 + DX1[DIR_0M0], x2 + DX2[DIR_0M0], x3 + DX3[DIR_0M0], DIR_0M0);
distributions->setDistributionInvForDirection(f[DIR_PP0], x1 + DX1[DIR_MM0], x2 + DX2[DIR_MM0], x3 + DX3[DIR_MM0], DIR_MM0);
@@ -164,15 +164,15 @@ void NonReflectingOutflowWithRelaxationBCStrategy::applyBC()
distributions->setDistributionInvForDirection(f[DIR_MPM], x1 + DX1[DIR_PMP], x2 + DX2[DIR_PMP], x3 + DX3[DIR_PMP], DIR_PMP);
break;
case DIR_0M0:
- f[DIR_0M0] = ftemp[DIR_0M0] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_0M0]- delf* WEIGTH[DIR_0M0];
- f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_PM0]- delf* WEIGTH[DIR_PM0];
- f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_MM0]- delf* WEIGTH[DIR_MM0];
- f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_0MP]- delf* WEIGTH[DIR_0MP];
- f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_0MM]- delf* WEIGTH[DIR_0MM];
- f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
- f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_MMP]- delf* WEIGTH[DIR_MMP];
- f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
- f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx2) + (1.0 - one_over_sqrt3 + vx2) * f[DIR_MMM]- delf* WEIGTH[DIR_MMM];
+ f[DIR_0M0] = ftemp[DIR_0M0] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_0M0]- delf* WEIGTH[DIR_0M0];
+ f[DIR_PM0] = ftemp[DIR_PM0] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_PM0]- delf* WEIGTH[DIR_PM0];
+ f[DIR_MM0] = ftemp[DIR_MM0] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_MM0]- delf* WEIGTH[DIR_MM0];
+ f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_0MP]- delf* WEIGTH[DIR_0MP];
+ f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_0MM]- delf* WEIGTH[DIR_0MM];
+ f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
+ f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_MMP]- delf* WEIGTH[DIR_MMP];
+ f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
+ f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx2) + (c1o1 - one_over_sqrt3 + vx2) * f[DIR_MMM]- delf* WEIGTH[DIR_MMM];
distributions->setDistributionInvForDirection(f[DIR_0M0], x1 + DX1[DIR_0P0], x2 + DX2[DIR_0P0], x3 + DX3[DIR_0P0], DIR_0P0);
distributions->setDistributionInvForDirection(f[DIR_PM0], x1 + DX1[DIR_MP0], x2 + DX2[DIR_MP0], x3 + DX3[DIR_MP0], DIR_MP0);
@@ -185,15 +185,15 @@ void NonReflectingOutflowWithRelaxationBCStrategy::applyBC()
distributions->setDistributionInvForDirection(f[DIR_MMM], x1 + DX1[DIR_PPP], x2 + DX2[DIR_PPP], x3 + DX3[DIR_PPP], DIR_PPP);
break;
case DIR_00P:
- f[DIR_00P] = ftemp[DIR_00P] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_00P]- delf* WEIGTH[DIR_00P];
- f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_P0P]- delf* WEIGTH[DIR_P0P];
- f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_M0P]- delf* WEIGTH[DIR_M0P];
- f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_0PP]- delf* WEIGTH[DIR_0PP];
- f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_0MP]- delf* WEIGTH[DIR_0MP];
- f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
- f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_MPP]- delf* WEIGTH[DIR_MPP];
- f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
- f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 + vx3) + (1.0 - one_over_sqrt3 - vx3) * f[DIR_MMP]- delf* WEIGTH[DIR_MMP];
+ f[DIR_00P] = ftemp[DIR_00P] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_00P]- delf* WEIGTH[DIR_00P];
+ f[DIR_P0P] = ftemp[DIR_P0P] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_P0P]- delf* WEIGTH[DIR_P0P];
+ f[DIR_M0P] = ftemp[DIR_M0P] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_M0P]- delf* WEIGTH[DIR_M0P];
+ f[DIR_0PP] = ftemp[DIR_0PP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_0PP]- delf* WEIGTH[DIR_0PP];
+ f[DIR_0MP] = ftemp[DIR_0MP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_0MP]- delf* WEIGTH[DIR_0MP];
+ f[DIR_PPP] = ftemp[DIR_PPP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_PPP]- delf* WEIGTH[DIR_PPP];
+ f[DIR_MPP] = ftemp[DIR_MPP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_MPP]- delf* WEIGTH[DIR_MPP];
+ f[DIR_PMP] = ftemp[DIR_PMP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_PMP]- delf* WEIGTH[DIR_PMP];
+ f[DIR_MMP] = ftemp[DIR_MMP] * (one_over_sqrt3 + vx3) + (c1o1 - one_over_sqrt3 - vx3) * f[DIR_MMP]- delf* WEIGTH[DIR_MMP];
distributions->setDistributionInvForDirection(f[DIR_00P], x1 + DX1[DIR_00M], x2 + DX2[DIR_00M], x3 + DX3[DIR_00M], DIR_00M);
distributions->setDistributionInvForDirection(f[DIR_P0P], x1 + DX1[DIR_M0M], x2 + DX2[DIR_M0M], x3 + DX3[DIR_M0M], DIR_M0M);
@@ -206,15 +206,15 @@ void NonReflectingOutflowWithRelaxationBCStrategy::applyBC()
distributions->setDistributionInvForDirection(f[DIR_MMP], x1 + DX1[DIR_PPM], x2 + DX2[DIR_PPM], x3 + DX3[DIR_PPM], DIR_PPM);
break;
case DIR_00M:
- f[DIR_00M] = ftemp[DIR_00M] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_00M]- delf* WEIGTH[DIR_00M];
- f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_P0M]- delf* WEIGTH[DIR_P0M];
- f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_M0M]- delf* WEIGTH[DIR_M0M];
- f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_0PM]- delf* WEIGTH[DIR_0PM];
- f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_0MM]- delf* WEIGTH[DIR_0MM];
- f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
- f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_MPM]- delf* WEIGTH[DIR_MPM];
- f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
- f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx3) + (1.0 - one_over_sqrt3 + vx3) * f[DIR_MMM]- delf* WEIGTH[DIR_MMM];
+ f[DIR_00M] = ftemp[DIR_00M] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_00M]- delf* WEIGTH[DIR_00M];
+ f[DIR_P0M] = ftemp[DIR_P0M] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_P0M]- delf* WEIGTH[DIR_P0M];
+ f[DIR_M0M] = ftemp[DIR_M0M] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_M0M]- delf* WEIGTH[DIR_M0M];
+ f[DIR_0PM] = ftemp[DIR_0PM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_0PM]- delf* WEIGTH[DIR_0PM];
+ f[DIR_0MM] = ftemp[DIR_0MM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_0MM]- delf* WEIGTH[DIR_0MM];
+ f[DIR_PPM] = ftemp[DIR_PPM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_PPM]- delf* WEIGTH[DIR_PPM];
+ f[DIR_MPM] = ftemp[DIR_MPM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_MPM]- delf* WEIGTH[DIR_MPM];
+ f[DIR_PMM] = ftemp[DIR_PMM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_PMM]- delf* WEIGTH[DIR_PMM];
+ f[DIR_MMM] = ftemp[DIR_MMM] * (one_over_sqrt3 - vx3) + (c1o1 - one_over_sqrt3 + vx3) * f[DIR_MMM]- delf* WEIGTH[DIR_MMM];
distributions->setDistributionInvForDirection(f[DIR_00M], x1 + DX1[DIR_00P], x2 + DX2[DIR_00P], x3 + DX3[DIR_00P], DIR_00P);
distributions->setDistributionInvForDirection(f[DIR_P0M], x1 + DX1[DIR_M0P], x2 + DX2[DIR_M0P], x3 + DX3[DIR_M0P], DIR_M0P);
diff --git a/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCStrategy.cpp b/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCStrategy.cpp
index 1b579877a..0278ee095 100644
--- a/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCStrategy.cpp
+++ b/src/cpu/VirtualFluidsCore/BoundaryConditions/VelocityWithDensityBCStrategy.cpp
@@ -55,6 +55,8 @@ void VelocityWithDensityBCStrategy::addDistributions(SPtr<DistributionArray3D> d
//////////////////////////////////////////////////////////////////////////
void VelocityWithDensityBCStrategy::applyBC()
{
+ using namespace vf::basics::constant;
+
//velocity bc for non reflecting pressure bc
real f[D3Q27System::ENDF+1];
//real feq[D3Q27System::ENDF+1];
@@ -63,7 +65,7 @@ void VelocityWithDensityBCStrategy::applyBC()
calcMacrosFct(f, drho, vx1, vx2, vx3);
//calcFeqFct(feq, drho, vx1, vx2, vx3);
- rho = 1.0+drho*compressibleFactor;
+ rho = c1o1+drho*compressibleFactor;
for (int fdir = D3Q27System::FSTARTDIR; fdir <= D3Q27System::FENDDIR; fdir++)
{
diff --git a/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp
index 8334b93d2..a94da62ef 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsDoubleGhostLayerFullVectorConnector.cpp
@@ -51,6 +51,7 @@ ThreeDistributionsDoubleGhostLayerFullVectorConnector::ThreeDistributionsDoubleG
void ThreeDistributionsDoubleGhostLayerFullVectorConnector::init()
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
FullVectorConnector::init();
@@ -64,26 +65,26 @@ void ThreeDistributionsDoubleGhostLayerFullVectorConnector::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:
@@ -92,7 +93,7 @@ void ThreeDistributionsDoubleGhostLayerFullVectorConnector::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/ThreeDistributionsFullVectorConnector.cpp b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsFullVectorConnector.cpp
index 1b4f243ee..65315b4f3 100644
--- a/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsFullVectorConnector.cpp
+++ b/src/cpu/VirtualFluidsCore/Connectors/ThreeDistributionsFullVectorConnector.cpp
@@ -51,6 +51,7 @@ ThreeDistributionsFullVectorConnector::ThreeDistributionsFullVectorConnector(SPt
void ThreeDistributionsFullVectorConnector::init()
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
FullVectorConnector::init();
@@ -63,26 +64,26 @@ void ThreeDistributionsFullVectorConnector::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:
@@ -91,7 +92,7 @@ void ThreeDistributionsFullVectorConnector::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/Interactors/D3Q27Interactor.cpp b/src/cpu/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
index 02a94f573..d9c605a58 100644
--- a/src/cpu/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
+++ b/src/cpu/VirtualFluidsCore/Interactors/D3Q27Interactor.cpp
@@ -88,80 +88,81 @@ D3Q27Interactor::~D3Q27Interactor() = default;
void D3Q27Interactor::initRayVectors()
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
int fdir;
real c1oS2 = vf::basics::constant::one_over_sqrt2;
real c1oS3 = vf::basics::constant::one_over_sqrt3;
fdir = DIR_P00;
- rayX1[fdir] = 1.0;
- rayX2[fdir] = 0.0;
- rayX3[fdir] = 0.0;
+ rayX1[fdir] = c1o1;
+ rayX2[fdir] = c0o1;
+ rayX3[fdir] = c0o1;
fdir = DIR_M00;
- rayX1[fdir] = -1.0;
- rayX2[fdir] = 0.0;
- rayX3[fdir] = 0.0;
+ rayX1[fdir] = -c1o1;
+ rayX2[fdir] = c0o1;
+ rayX3[fdir] = c0o1;
fdir = DIR_0P0;
- rayX1[fdir] = 0.0;
- rayX2[fdir] = 1.0;
- rayX3[fdir] = 0.0;
+ rayX1[fdir] = c0o1;
+ rayX2[fdir] = c1o1;
+ rayX3[fdir] = c0o1;
fdir = DIR_0M0;
- rayX1[fdir] = 0.0;
- rayX2[fdir] = -1.0;
- rayX3[fdir] = 0.0;
+ rayX1[fdir] = c0o1;
+ rayX2[fdir] = -c1o1;
+ rayX3[fdir] = c0o1;
fdir = DIR_00P;
- rayX1[fdir] = 0.0;
- rayX2[fdir] = 0.0;
- rayX3[fdir] = 1.0;
+ rayX1[fdir] = c0o1;
+ rayX2[fdir] = c0o1;
+ rayX3[fdir] = c1o1;
fdir = DIR_00M;
- rayX1[fdir] = 0.0;
- rayX2[fdir] = 0.0;
- rayX3[fdir] = -1.0;
+ rayX1[fdir] = c0o1;
+ rayX2[fdir] = c0o1;
+ rayX3[fdir] = -c1o1;
fdir = DIR_PP0;
rayX1[fdir] = c1oS2;
rayX2[fdir] = c1oS2;
- rayX3[fdir] = 0.0;
+ rayX3[fdir] = c0o1;
fdir = DIR_MM0;
rayX1[fdir] = -c1oS2;
rayX2[fdir] = -c1oS2;
- rayX3[fdir] = 0.0;
+ rayX3[fdir] = c0o1;
fdir = DIR_PM0;
rayX1[fdir] = c1oS2;
rayX2[fdir] = -c1oS2;
- rayX3[fdir] = 0.0;
+ rayX3[fdir] = c0o1;
fdir = DIR_MP0;
rayX1[fdir] = -c1oS2;
rayX2[fdir] = c1oS2;
- rayX3[fdir] = 0.0;
+ rayX3[fdir] = c0o1;
fdir = DIR_P0P;
rayX1[fdir] = c1oS2;
- rayX2[fdir] = 0.0;
+ rayX2[fdir] = c0o1;
rayX3[fdir] = c1oS2;
fdir = DIR_M0M;
rayX1[fdir] = -c1oS2;
- rayX2[fdir] = 0.0;
+ rayX2[fdir] = c0o1;
rayX3[fdir] = -c1oS2;
fdir = DIR_P0M;
rayX1[fdir] = c1oS2;
- rayX2[fdir] = 0.0;
+ rayX2[fdir] = c0o1;
rayX3[fdir] = -c1oS2;
fdir = DIR_M0P;
rayX1[fdir] = -c1oS2;
- rayX2[fdir] = 0.0;
+ rayX2[fdir] = c0o1;
rayX3[fdir] = c1oS2;
fdir = DIR_0PP;
- rayX1[fdir] = 0.0;
+ rayX1[fdir] = c0o1;
rayX2[fdir] = c1oS2;
rayX3[fdir] = c1oS2;
fdir = DIR_0MM;
- rayX1[fdir] = 0.0;
+ rayX1[fdir] = c0o1;
rayX2[fdir] = -c1oS2;
rayX3[fdir] = -c1oS2;
fdir = DIR_0PM;
- rayX1[fdir] = 0.0;
+ rayX1[fdir] = c0o1;
rayX2[fdir] = c1oS2;
rayX3[fdir] = -c1oS2;
fdir = DIR_0MP;
- rayX1[fdir] = 0.0;
+ rayX1[fdir] = c0o1;
rayX2[fdir] = -c1oS2;
rayX3[fdir] = c1oS2;
@@ -285,6 +286,7 @@ void D3Q27Interactor::updateInteractor(const real ×tep)
bool D3Q27Interactor::setDifferencesToGbObject3D(const SPtr<Block3D> block)
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
if (!block)
return false;
@@ -297,7 +299,7 @@ bool D3Q27Interactor::setDifferencesToGbObject3D(const SPtr<Block3D> block)
solidNodeIndicesMap[block] = set<UbTupleInt3>();
set<UbTupleInt3> &solidNodeIndices = solidNodeIndicesMap[block];
- real timestep = 0;
+ real timestep = c0o1;
bool oneEntryGotBC = false;
bool gotQs = false;
SPtr<BoundaryConditions> bc;
@@ -396,11 +398,11 @@ bool D3Q27Interactor::setDifferencesToGbObject3D(const SPtr<Block3D> block)
rayX2[fdir], rayX3[fdir]);
q /= distNeigh[fdir];
- // assert(UbMath::lessEqual(q, 1.0));
+ // assert(UbMath::lessEqual(q, c1o1));
- if (UbMath::inClosedInterval(q, 1.0, 1.0))
- q = 1.0;
- if (UbMath::greater(q, 0.0) && UbMath::lessEqual(q, 1.0)) {
+ if (UbMath::inClosedInterval(q, c1o1, c1o1))
+ q = c1o1;
+ if (UbMath::greater(q, c0o1) && UbMath::lessEqual(q, c1o1)) {
//#pragma omp critical (BC_CHANGE)
{
bc = bcArray->getBC(ix1, ix2, ix3);
@@ -410,9 +412,9 @@ bool D3Q27Interactor::setDifferencesToGbObject3D(const SPtr<Block3D> block)
}
if (bc->hasNoSlipBoundary()) {
- bc->setBoundaryVelocityX1(0.0);
- bc->setBoundaryVelocityX2(0.0);
- bc->setBoundaryVelocityX3(0.0);
+ bc->setBoundaryVelocityX1(c0o1);
+ bc->setBoundaryVelocityX2(c0o1);
+ bc->setBoundaryVelocityX3(c0o1);
}
for (int index = (int)BCs.size() - 1; index >= 0; --index)
@@ -499,7 +501,7 @@ bool D3Q27Interactor::setDifferencesToGbObject3D(const SPtr<Block3D> block)
GbLine3D *clippedLine = this->geoObject3D->createClippedLine3D(pointA, pointB);
if (clippedLine) {
- real q = 0.0;
+ real q = c0o1;
if (!this->isInverseSolid()) // A is outside
{
real distanceAB = pointA.getDistance(&pointB); // pointA to B
@@ -518,16 +520,16 @@ bool D3Q27Interactor::setDifferencesToGbObject3D(const SPtr<Block3D> block)
pointB.getX1Coordinate(), pointB.getX2Coordinate(),
pointB.getX3Coordinate(), pointIsOnBoundary) &&
pointIsOnBoundary) {
- // A is definitely inside, B is exactly on ObjectBoundary => q = 1.0
- q = 1.0;
+ // A is definitely inside, B is exactly on ObjectBoundary => q = c1o1
+ q = c1o1;
} else {
- q = 0.0;
+ q = c0o1;
}
}
- if (UbMath::inClosedInterval(q, 1.0, 1.0))
- q = 1.0;
- if (UbMath::lessEqual(q, 1.0) && UbMath::greater(q, 0.0)) {
+ if (UbMath::inClosedInterval(q, c1o1, c1o1))
+ q = c1o1;
+ if (UbMath::lessEqual(q, c1o1) && UbMath::greater(q, c0o1)) {
{
bc = bcArray->getBC(ix1, ix2, ix3);
if (!bc) {
diff --git a/src/cpu/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp b/src/cpu/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp
index c30a4cc8c..db41bba85 100644
--- a/src/cpu/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp
+++ b/src/cpu/VirtualFluidsCore/Interactors/D3Q27TriFaceMeshInteractor.cpp
@@ -123,6 +123,7 @@ bool D3Q27TriFaceMeshInteractor::setDifferencesToGbObject3D(const SPtr<Block3D>
void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
UBLOGML(logDEBUG1, "\nLBMTriFaceMeshInteractor - setQs start ");
if (!this->grid.lock())
@@ -168,7 +169,7 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
// grobe Blocklaengen
SPtr<CoordinateTransformation3D> trafo = grid.lock()->getCoordinateTransformator();
double cblockDeltaX1, cblockDeltaX2, cblockDeltaX3, delta;
- cblockDeltaX1 = cblockDeltaX2 = cblockDeltaX3 = delta = 1.0 / (double)(1 << coarsestInitLevel);
+ cblockDeltaX1 = cblockDeltaX2 = cblockDeltaX3 = delta = c1o1 / (double)(1 << coarsestInitLevel);
if (trafo) {
cblockDeltaX1 = trafo->getX1CoordinateScaling() * delta;
cblockDeltaX2 = trafo->getX2CoordinateScaling() * delta;
@@ -208,12 +209,12 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
// blockDeltaCalculator->getMinX1Delta(level) -> ein geinlinter wert wird geholt
// TODO: set 5.0 as variable parameter in constructor, default 2.0
- deltaMinX1[level] = (float)(5.0 * nodeDeltaX1); // kein minus da unten -deltaMin
- deltaMinX2[level] = (float)(5.0 * nodeDeltaX2);
- deltaMinX3[level] = (float)(5.0 * nodeDeltaX3);
- deltaMaxX1[level] = (float)(5.0 * nodeDeltaX1);
- deltaMaxX2[level] = (float)(5.0 * nodeDeltaX2);
- deltaMaxX3[level] = (float)(5.0 * nodeDeltaX3);
+ deltaMinX1[level] = (float)(c5o1 * nodeDeltaX1); // kein minus da unten -deltaMin
+ deltaMinX2[level] = (float)(c5o1 * nodeDeltaX2);
+ deltaMinX3[level] = (float)(c5o1 * nodeDeltaX3);
+ deltaMaxX1[level] = (float)(c5o1 * nodeDeltaX1);
+ deltaMaxX2[level] = (float)(c5o1 * nodeDeltaX2);
+ deltaMaxX3[level] = (float)(c5o1 * nodeDeltaX3);
}
//////////////////////////////////////////////////////////////////////////
@@ -237,7 +238,7 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
std::vector<GbTriFaceMesh3D::Vertex> &nodes = *mesh->getNodes();
std::map<SPtr<Block3D>, std::set<UbTupleInt3>> tmpSolidNodesFromOtherInteractors;
- int onePercent = UbMath::integerRounding(triangles.size() * 0.01);
+ int onePercent = UbMath::integerRounding(triangles.size() * c1o100);
if (onePercent == 0)
onePercent = 1;
UbTimer setQTimer;
@@ -330,13 +331,13 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
blockMaxX[1] = (float)(val<2>(coords) + val<2>(deltas) + deltaMaxX2[level]);
blockMaxX[2] = (float)(val<3>(coords) + val<3>(deltas) + deltaMaxX3[level]);
- boxCenter[0] = (float)(0.5 * (blockMaxX[0] + blockMinX[0]));
- boxCenter[1] = (float)(0.5 * (blockMaxX[1] + blockMinX[1]));
- boxCenter[2] = (float)(0.5 * (blockMaxX[2] + blockMinX[2]));
+ boxCenter[0] = (float)(c1o2 * (blockMaxX[0] + blockMinX[0]));
+ boxCenter[1] = (float)(c1o2 * (blockMaxX[1] + blockMinX[1]));
+ boxCenter[2] = (float)(c1o2 * (blockMaxX[2] + blockMinX[2]));
- halfBoxSize[0] = (float)(0.5 * (blockMaxX[0] - blockMinX[0]));
- halfBoxSize[1] = (float)(0.5 * (blockMaxX[1] - blockMinX[1]));
- halfBoxSize[2] = (float)(0.5 * (blockMaxX[2] - blockMinX[2]));
+ halfBoxSize[0] = (float)(c1o2 * (blockMaxX[0] - blockMinX[0]));
+ halfBoxSize[1] = (float)(c1o2 * (blockMaxX[1] - blockMinX[1]));
+ halfBoxSize[2] = (float)(c1o2 * (blockMaxX[2] - blockMinX[2]));
// wenn dreieck "vergroesserten cube" nicht schneidet/beruehrt -> keine BC moeglich -> continue
if (!GbMeshTools3D::triBoxOverlap(boxCenter, halfBoxSize, triPoints)) {
@@ -445,11 +446,11 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
a = e1x1 * px1 + e1x2 * px2 + e1x3 * px3;
if (fabs(a) < 1.E-10)
continue;
- f = 1.0 / a;
+ f = c1o1 / a;
// u = f * ( s dot p)
u = f * (sx1 * px1 + sx2 * px2 + sx3 * px3);
- if (u < -1.E-10 || u > 1.0 + 1.E-10)
+ if (u < -1.E-10 || u > c1o1 + 1.E-10)
continue;
// q = s x e1
@@ -459,7 +460,7 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
// v = f*(e2 dot q)
v = f * (this->rayX1[fdir] * qx1 + this->rayX2[fdir] * qx2 + this->rayX3[fdir] * qx3);
- if (v < -1.E-10 || (u + v) > 1.0 + 1.E-10)
+ if (v < -1.E-10 || (u + v) > c1o1 + 1.E-10)
continue;
// t = f * (e2 dot q)
@@ -492,7 +493,7 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
(z1 - internX3) * (z1 - internX3));
q_ehsan /= nodeDeltaToNeigh[level][fdir];
q = q_ehsan;
- if (UbMath::greater(q, 1.0) || UbMath::lessEqual(q, 0.0))
+ if (UbMath::greater(q, c1o1) || UbMath::lessEqual(q, 0.0))
continue;
// gefundenes q auf gueltigkeit pruefen
@@ -507,9 +508,9 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
continue;
}
- if (UbMath::inClosedInterval(q, 1.0, 1.0))
- q = 1.0;
- if (UbMath::greater(q, 0.0) && UbMath::lessEqual(q, 1.0)) {
+ if (UbMath::inClosedInterval(q, c1o1, c1o1))
+ q = c1o1;
+ if (UbMath::greater(q, 0.0) && UbMath::lessEqual(q, c1o1)) {
gotQs = blockGotBCs = true;
bc = bcMatrix->getBC(ix1, ix2, ix3);
@@ -591,6 +592,7 @@ void D3Q27TriFaceMeshInteractor::setQs(const real &timeStep)
void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
{
using namespace vf::lbm::dir;
+ using namespace vf::basics::constant;
UBLOGML(logDEBUG1, "\nLBMTriFaceMeshInteractor - initInteractor start ");
if (!this->grid.lock())
@@ -642,7 +644,7 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
// grobe Blocklaengen
SPtr<CoordinateTransformation3D> trafo = grid.lock()->getCoordinateTransformator();
double cblockDeltaX1, cblockDeltaX2, cblockDeltaX3, delta;
- cblockDeltaX1 = cblockDeltaX2 = cblockDeltaX3 = delta = 1.0 / (double)(1 << coarsestInitLevel);
+ cblockDeltaX1 = cblockDeltaX2 = cblockDeltaX3 = delta = c1o1 / (double)(1 << coarsestInitLevel);
if (trafo) {
cblockDeltaX1 = trafo->getX1CoordinateScaling() * delta;
cblockDeltaX2 = trafo->getX2CoordinateScaling() * delta;
@@ -735,13 +737,13 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
// notwendige variablen initialisieren (u.a. blockDeltas des groben levels)
float triPoints[3][3];
- real vx1 = 0.0, vx2 = 0.0, vx3 = 0.0;
+ real vx1 = c0o1, vx2 = c0o1, vx3 = c0o1;
unsigned counterTriBoxOverlap = 0, counterAABBTriFace = 0, counterHalfspace = 0, counterBilligOBB = 0;
std::vector<GbTriFaceMesh3D::TriFace> &triangles = *mesh->getTriangles();
std::vector<GbTriFaceMesh3D::Vertex> &nodes = *mesh->getNodes();
std::map<SPtr<Block3D>, std::set<std::vector<int>>> tmpSolidNodesFromOtherInteractors;
- int onePercent = UbMath::integerRounding(triangles.size() * 0.01);
+ int onePercent = UbMath::integerRounding(triangles.size() * c1o100);
if (onePercent == 0)
onePercent = 1;
UbTimer setQTimer;
@@ -845,13 +847,13 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
blockMaxX[1] = (float)(val<2>(coords) + val<2>(deltas) + deltaMaxX2[level]);
blockMaxX[2] = (float)(val<3>(coords) + val<3>(deltas) + deltaMaxX3[level]);
- boxCenter[0] = (float)(0.5 * (blockMaxX[0] + blockMinX[0]));
- boxCenter[1] = (float)(0.5 * (blockMaxX[1] + blockMinX[1]));
- boxCenter[2] = (float)(0.5 * (blockMaxX[2] + blockMinX[2]));
+ boxCenter[0] = (float)(c1o2 * (blockMaxX[0] + blockMinX[0]));
+ boxCenter[1] = (float)(c1o2 * (blockMaxX[1] + blockMinX[1]));
+ boxCenter[2] = (float)(c1o2 * (blockMaxX[2] + blockMinX[2]));
- halfBoxSize[0] = (float)(0.5 * (blockMaxX[0] - blockMinX[0]));
- halfBoxSize[1] = (float)(0.5 * (blockMaxX[1] - blockMinX[1]));
- halfBoxSize[2] = (float)(0.5 * (blockMaxX[2] - blockMinX[2]));
+ halfBoxSize[0] = (float)(c1o2 * (blockMaxX[0] - blockMinX[0]));
+ halfBoxSize[1] = (float)(c1o2 * (blockMaxX[1] - blockMinX[1]));
+ halfBoxSize[2] = (float)(c1o2 * (blockMaxX[2] - blockMinX[2]));
// wenn dreieck "vergroesserten cube" nicht schneidet/beruehrt -> keine BC moeglich -> continue
if (!GbMeshTools3D::triBoxOverlap(boxCenter, halfBoxSize, triPoints)) {
@@ -887,9 +889,9 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
double qEinflussDelta = 1.1 * sqrt(nodeDx1 * nodeDx1 + nodeDx2 * nodeDx2 + nodeDx3 * nodeDx3);
for (int ix3 = indexMinX3; ix3 < indexMaxX3; ix3++) {
- internX3 = val<3>(coords) + nodeDx3 * ix3 - 0.5 * nodeDx3;
+ internX3 = val<3>(coords) + nodeDx3 * ix3 - c1o2 * nodeDx3;
for (int ix2 = indexMinX2; ix2 < indexMaxX2; ix2++) {
- internX2 = val<2>(coords) + nodeDx2 * ix2 - 0.5 * nodeDx2;
+ internX2 = val<2>(coords) + nodeDx2 * ix2 - c1o2 * nodeDx2;
for (int ix1 = indexMinX1; ix1 < indexMaxX1; ix1++) {
// int blx1 =block->getX1();
@@ -930,7 +932,7 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
}
}
- internX1 = val<1>(coords) + nodeDx1 * ix1 - 0.5 * nodeDx1;
+ internX1 = val<1>(coords) + nodeDx1 * ix1 - c1o2 * nodeDx1;
//////////////////////////////////////////////////////////////////////////
// Punkt in AABB von Dreieck?
@@ -996,11 +998,11 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
a = e1x1 * px1 + e1x2 * px2 + e1x3 * px3;
if (fabs(a) < 1.E-10)
continue;
- f = 1.0 / a;
+ f = c1o1 / a;
// u = f * ( s dot p)
u = f * (sx1 * px1 + sx2 * px2 + sx3 * px3);
- if (u < -1.E-10 || u > 1.0 + 1.E-10)
+ if (u < -1.E-10 || u > c1o1 + 1.E-10)
continue;
// q = s x e1
@@ -1010,7 +1012,7 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
// v = f*(e2 dot q)
v = f * (this->rayX1[fdir] * qx1 + this->rayX2[fdir] * qx2 + this->rayX3[fdir] * qx3);
- if (v < -1.E-10 || (u + v) > 1.0 + 1.E-10)
+ if (v < -1.E-10 || (u + v) > c1o1 + 1.E-10)
continue;
// t = f * (e2 dot q)
@@ -1029,9 +1031,9 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
continue;
}
- if (UbMath::inClosedInterval(q, 1.0, 1.0))
- q = 1.0;
- if (UbMath::greater(q, 0.0) && UbMath::lessEqual(q, 1.0)) {
+ if (UbMath::inClosedInterval(q, c1o1, c1o1))
+ q = c1o1;
+ if (UbMath::greater(q, 0.0) && UbMath::lessEqual(q, c1o1)) {
// if( !solidFromOtherInteractor ) //--> Knoten schon solid-->BC setzen
// ueberfluessig SG changed to if( solidFromOtherInteractor ) //--> Knoten schon
// solid-->BC setzen ueberfluessig
@@ -1211,9 +1213,9 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
// int ride=0;
// }
if (flagField(bx1, bx2, bx3) == UNDEF_FLAG) {
- if (mesh->isPointInGbObject3D(val<1>(coords) + bx1 * nodeDeltaX1 - 0.5 * nodeDeltaX1,
- val<2>(coords) + bx2 * nodeDeltaX2 - 0.5 * nodeDeltaX2,
- val<3>(coords) + bx3 * nodeDeltaX3 - 0.5 * nodeDeltaX3)) {
+ if (mesh->isPointInGbObject3D(val<1>(coords) + bx1 * nodeDeltaX1 - c1o2 * nodeDeltaX1,
+ val<2>(coords) + bx2 * nodeDeltaX2 - c1o2 * nodeDeltaX2,
+ val<3>(coords) + bx3 * nodeDeltaX3 - c1o2 * nodeDeltaX3)) {
(this->*gridFill)(flagField, bx1, bx2, bx3, SOLID_FLAG);
} else {
(this->*gridFill)(flagField, bx1, bx2, bx3, FLUID_FLAG);
@@ -1309,6 +1311,8 @@ void D3Q27TriFaceMeshInteractor::initInteractor2(const real &timeStep)
//////////////////////////////////////////////////////////////////////////
void D3Q27TriFaceMeshInteractor::refineBlockGridToLevel(int level, real startDistance, real stopDistance)
{
+ using namespace vf::basics::constant;
+
UBLOG(logDEBUG1, "D3Q27TriFaceMeshInteractor::refineBlockGridToLevel - start");
// ToDo: evtl checken, ob man noch einen HalbraumCheck für StopDistance einbaut
@@ -1365,27 +1369,27 @@ void D3Q27TriFaceMeshInteractor::refineBlockGridToLevel(int level, real startDis
// Boundingbox um massgebliches dx erweitern -> zur Bestimmung der zu testenden Bloecke
if (triangle.nx > 1.0E-8) {
- minX1 = triangle.getMinX(nodes) + 1.05 * triangle.nx * startDistance;
- maxX1 = triangle.getMaxX(nodes) + 1.05 * triangle.nx * stopDistance;
+ minX1 = triangle.getMinX(nodes) + c21o20 * triangle.nx * startDistance;
+ maxX1 = triangle.getMaxX(nodes) + c21o20 * triangle.nx * stopDistance;
} else {
- minX1 = triangle.getMinX(nodes) + 1.05 * triangle.nx * stopDistance;
- maxX1 = triangle.getMaxX(nodes) + 1.05 * triangle.nx * startDistance;
+ minX1 = triangle.getMinX(nodes) + c21o20 * triangle.nx * stopDistance;
+ maxX1 = triangle.getMaxX(nodes) + c21o20 * triangle.nx * startDistance;
}
if (triangle.ny > 1.0E-8) {
- minX2 = triangle.getMinY(nodes) + 1.05 * triangle.ny * startDistance;
- maxX2 = triangle.getMaxY(nodes) + 1.05 * triangle.ny * stopDistance;
+ minX2 = triangle.getMinY(nodes) + c21o20 * triangle.ny * startDistance;
+ maxX2 = triangle.getMaxY(nodes) + c21o20 * triangle.ny * stopDistance;
} else {
- minX2 = triangle.getMinY(nodes) + 1.05 * triangle.ny * stopDistance;
- maxX2 = triangle.getMaxY(nodes) + 1.05 * triangle.ny * startDistance;
+ minX2 = triangle.getMinY(nodes) + c21o20 * triangle.ny * stopDistance;
+ maxX2 = triangle.getMaxY(nodes) + c21o20 * triangle.ny * startDistance;
}
if (triangle.nz > 1.0E-8) {
- minX3 = triangle.getMinZ(nodes) + 1.05 * triangle.nz * startDistance;
- maxX3 = triangle.getMaxZ(nodes) + 1.05 * triangle.nz * stopDistance;
+ minX3 = triangle.getMinZ(nodes) + c21o20 * triangle.nz * startDistance;
+ maxX3 = triangle.getMaxZ(nodes) + c21o20 * triangle.nz * stopDistance;
} else {
- minX3 = triangle.getMinZ(nodes) + 1.05 * triangle.nz * stopDistance;
- maxX3 = triangle.getMaxZ(nodes) + 1.05 * triangle.nz * startDistance;
+ minX3 = triangle.getMinZ(nodes) + c21o20 * triangle.nz * stopDistance;
+ maxX3 = triangle.getMaxZ(nodes) + c21o20 * triangle.nz * startDistance;
}
int flag = 0;
@@ -1435,15 +1439,15 @@ void D3Q27TriFaceMeshInteractor::refineBlockGridToLevel(int level, real startDis
// blockseite ermitteln (skalarprodukt dreiecks-normale, vector (midTri->midCub) )
// je nachdem muss für den massgeblichen block start oder stopdistance verwendet werden
// liegt block auf pos seite -> stopdistance ansonsten startdistance
- double skalarprod = triangle.nx * (0.5 * (x1a + x1b) - triangle.getX1Centroid(nodes)) +
- triangle.ny * (0.5 * (x2a + x2b) - triangle.getX2Centroid(nodes)) +
- triangle.nz * (0.5 * (x3a + x3b) - triangle.getX3Centroid(nodes));
+ double skalarprod = triangle.nx * (c1o2 * (x1a + x1b) - triangle.getX1Centroid(nodes)) +
+ triangle.ny * (c1o2 * (x2a + x2b) - triangle.getX2Centroid(nodes)) +
+ triangle.nz * (c1o2 * (x3a + x3b) - triangle.getX3Centroid(nodes));
- double blockdelta = 1.05 * stopDistance;
+ double blockdelta = c21o20 * stopDistance;
if (skalarprod < 1.E-8)
- blockdelta = -1.05 * startDistance; // startDistance<0!!
+ blockdelta = -c21o20 * startDistance; // startDistance<0!!
else if (fabs(skalarprod) < 1.E-8)
- blockdelta = 1.05 * UbMath::max(-startDistance, stopDistance);
+ blockdelta = c21o20 * UbMath::max(-startDistance, stopDistance);
// block anpassen
blockMinX[0] = (float)(val<1>(coords) - blockdelta);
@@ -1454,13 +1458,13 @@ void D3Q27TriFaceMeshInteractor::refineBlockGridToLevel(int level, real startDis
blockMaxX[1] = (float)(val<2>(coords) + val<2>(deltas) + blockdelta);
blockMaxX[2] = (float)(val<3>(coords) + val<3>(deltas) + blockdelta);
- boxCenter[0] = (float)(0.5 * (blockMaxX[0] + blockMinX[0]));
- boxCenter[1] = (float)(0.5 * (blockMaxX[1] + blockMinX[1]));
- boxCenter[2] = (float)(0.5 * (blockMaxX[2] + blockMinX[2]));
+ boxCenter[0] = (float)(c1o2 * (blockMaxX[0] + blockMinX[0]));
+ boxCenter[1] = (float)(c1o2 * (blockMaxX[1] + blockMinX[1]));
+ boxCenter[2] = (float)(c1o2 * (blockMaxX[2] + blockMinX[2]));
- halfBoxSize[0] = (float)(0.5 * (blockMaxX[0] - blockMinX[0]));
- halfBoxSize[1] = (float)(0.5 * (blockMaxX[1] - blockMinX[1]));
- halfBoxSize[2] = (float)(0.5 * (blockMaxX[2] - blockMinX[2]));
+ halfBoxSize[0] = (float)(c1o2 * (blockMaxX[0] - blockMinX[0]));
+ halfBoxSize[1] = (float)(c1o2 * (blockMaxX[1] - blockMinX[1]));
+ halfBoxSize[2] = (float)(c1o2 * (blockMaxX[2] - blockMinX[2]));
GbTriFaceMesh3D::Vertex &v1_ = nodes[triangle.v1];
GbTriFaceMesh3D::Vertex &v2_ = nodes[triangle.v2];
@@ -1558,10 +1562,11 @@ UbTupleDouble3 D3Q27TriFaceMeshInteractor::getForces()
//}
////return getForcesTriangle();
// this->calculateForces();
+ using namespace vf::basics::constant;
- real forceX1 = 0.0;
- real forceX2 = 0.0;
- real forceX3 = 0.0;
+ real forceX1 = c0o1;
+ real forceX2 = c0o1;
+ real forceX3 = c0o1;
// double area = 0.0;
@@ -1580,9 +1585,11 @@ UbTupleDouble3 D3Q27TriFaceMeshInteractor::getForces()
//////////////////////////////////////////////////////////////////////////
UbTupleDouble3 D3Q27TriFaceMeshInteractor::getForcesTriangle()
{
- real forceX1 = 0.0;
- real forceX2 = 0.0;
- real forceX3 = 0.0;
+ using namespace vf::basics::constant;
+
+ real forceX1 = c0o1;
+ real forceX2 = c0o1;
+ real forceX3 = c0o1;
// D3Q19BlockGrid& grid.lock() = dynamic_cast<D3Q19BlockGrid&>(*this->grid.lock());
//// CoordinateTransformation3D *trafo = this->grid.lock()->getTransformation();
@@ -1829,6 +1836,8 @@ string D3Q27TriFaceMeshInteractor::toString()
//////////////////////////////////////////////////////////////////////////
void D3Q27TriFaceMeshInteractor::reinitWithStoredQs(const real & /*timeStep*/)
{
+ using namespace vf::basics::constant;
+
// alle solid Bloecke wieder solid setzen
std::vector<SPtr<Block3D>> &solidBlocks = this->getSolidBlockSet();
for (size_t i = 0; i < solidBlocks.size(); i++) {
@@ -1879,9 +1888,9 @@ void D3Q27TriFaceMeshInteractor::reinitWithStoredQs(const real & /*timeStep*/)
// es handelt sich un ein statisches Objekt und beim Propeller gibt
// es Schwierigkeiten an den Flügelspitzen, dass kann daher kommen,
// dass dort zuviel bc-flaggs sind und mit Geschwindigkeit ergibt dies ziemlich grosse Werte
- bc->setBoundaryVelocityX1(0.0);
- bc->setBoundaryVelocityX2(0.0);
- bc->setBoundaryVelocityX3(0.0);
+ bc->setBoundaryVelocityX1(c0o1);
+ bc->setBoundaryVelocityX2(c0o1);
+ bc->setBoundaryVelocityX3(c0o1);
// TODO: HACK GEHOERT NICHT HIERHIER!!! - end
bool gotQs = false;
diff --git a/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp
index b50103ff8..0145e633d 100644
--- a/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/BGKLBMKernel.cpp
@@ -138,7 +138,7 @@ void BGKLBMKernel::calculate(int step)
vx3 = f[DIR_00P] - f[DIR_00M] + f[DIR_P0P] - f[DIR_M0M] - f[DIR_P0M] + f[DIR_M0P] + f[DIR_0PP] - f[DIR_0MM] - f[DIR_0PM] + f[DIR_0MP] + f[DIR_PPP] +
f[DIR_MMP] + f[DIR_PMP] + f[DIR_MPP] - f[DIR_PPM] - f[DIR_MMM] - f[DIR_PMM] - f[DIR_MPM];
- real cu_sq = 1.5 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
+ real cu_sq = c3o2 * (vx1 * vx1 + vx2 * vx2 + vx3 * vx3);
feq[DIR_000] = c8o27 * (drho - cu_sq);
feq[DIR_P00] = c2o27 * (drho + c3o1 * (vx1) + c9o2 * (vx1) * (vx1)-cu_sq);
diff --git a/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp b/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
index bf68b4d43..574014992 100644
--- a/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/InitDensityLBMKernel.cpp
@@ -957,7 +957,7 @@ void InitDensityLBMKernel::calculate(int /*step*/)
//vx2 = vx2+(vvy-vx2);
//vx3 = vx3+(vvz-vx3);
- real cu_sq = 1.5*(vx1*vx1+vx2*vx2+vx3*vx3);
+ real cu_sq = c3o2*(vx1*vx1+vx2*vx2+vx3*vx3);
feq[DIR_000] = c8o27*(drho-cu_sq);
feq[DIR_P00] = c2o27*(drho+c3o1*(vx1)+c9o2*(vx1)*(vx1)-cu_sq);
diff --git a/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetInterpolator.cpp b/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetInterpolator.cpp
index 509115d26..a5450b3b4 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetInterpolator.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetInterpolator.cpp
@@ -375,7 +375,7 @@ void CompressibleOffsetInterpolator::calcInterpolatedCoefficiets(const D3Q27ICel
x_BE = c1o4*eps_new*(-((c2o1*axx - c3o1*axz - c2o1*bxy - c6o1*cxx + cxz))/(c54o1*o));
x_TN = c1o4*eps_new*(-((-c4o1*axx + bxy + c3o1*bxz + c3o1*cxy + cxz))/(c54o1*o));
x_BN = c1o4*eps_new*(-((-c4o1*axx + bxy - c3o1*bxz - c3o1*cxy + cxz))/(c54o1*o));
- x_ZERO = 0.;
+ x_ZERO = c0o1;
x_TNE = c1o4*eps_new*(-((axy + axz + c2o1*bxx + bxz + c2o1*cxx + cxy))/(c72o1*o));
x_TSW = c1o4*eps_new*(((-axy + axz - c2o1*bxx + bxz + c2o1*cxx + cxy))/(c72o1*o));
x_TSE = c1o4*eps_new*(((axy - axz + c2o1*bxx + bxz - c2o1*cxx + cxy))/(c72o1*o));
@@ -390,7 +390,7 @@ void CompressibleOffsetInterpolator::calcInterpolatedCoefficiets(const D3Q27ICel
y_BE = c1o4*eps_new*(-((axy - c3o1*ayz - c4o1*byy - c3o1*cxy + cyz))/(c54o1*o));
y_TN = c1o4*eps_new*(-((-c2o1*axy + c2o1*byy + c3o1*byz + c6o1*cyy + cyz))/(c54o1*o));
y_BN = c1o4*eps_new*(-((-c2o1*axy + c2o1*byy - c3o1*byz - c6o1*cyy + cyz))/(c54o1*o));
- y_ZERO = 0.;
+ y_ZERO = c0o1;
y_TNE = c1o4*eps_new*(-((c2o1*ayy + ayz + bxy + byz + cxy + c2o1*cyy))/(c72o1*o));
y_TSW = c1o4*eps_new*(((-c2o1*ayy + ayz - bxy + byz + cxy + c2o1*cyy))/(c72o1*o));
y_TSE = c1o4*eps_new*(((c2o1*ayy - ayz + bxy + byz - cxy + c2o1*cyy))/(c72o1*o));
@@ -405,7 +405,7 @@ void CompressibleOffsetInterpolator::calcInterpolatedCoefficiets(const D3Q27ICel
z_BE = c1o4*eps_new*(-((axz - c6o1*azz - c2o1*byz - c3o1*cxz + c2o1*czz))/(c54o1*o));
z_TN = c1o4*eps_new*(-((-c2o1*axz + byz + c6o1*bzz + c3o1*cyz + c2o1*czz))/(c54o1*o));
z_BN = c1o4*eps_new*(-((-c2o1*axz + byz - c6o1*bzz - c3o1*cyz + c2o1*czz))/(c54o1*o));
- z_ZERO = 0.;
+ z_ZERO = c0o1;
z_TNE = c1o4*eps_new*(-((ayz + c2o1*azz + bxz + c2o1*bzz + cxz + cyz))/(c72o1*o));
z_TSW = c1o4*eps_new*(((-ayz + c2o1*azz - bxz + c2o1*bzz + cxz + cyz))/(c72o1*o));
z_TSE = c1o4*eps_new*(((ayz - c2o1*azz + bxz + c2o1*bzz - cxz + cyz))/(c72o1*o));
diff --git a/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetMomentsInterpolator.cpp b/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetMomentsInterpolator.cpp
index e7d891df0..9003b1206 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetMomentsInterpolator.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/Interpolation/CompressibleOffsetMomentsInterpolator.cpp
@@ -46,14 +46,14 @@ void CompressibleOffsetMomentsInterpolator::interpolateCoarseToFine(D3Q27ICell&
vf::lbm::Coefficients coefficients;
calculateCoefficients(coefficients, icellC, omegaC, xoff, yoff, zoff);
- vf::lbm::interpolateCF(icellF.BSW, omegaF, vf::basics::constant::c1o2, coefficients, -0.25, -0.25, -0.25);
- vf::lbm::interpolateCF(icellF.BNE, omegaF, vf::basics::constant::c1o2, coefficients, 0.25, 0.25, -0.25);
- vf::lbm::interpolateCF(icellF.TNW, omegaF, vf::basics::constant::c1o2, coefficients, -0.25, 0.25, 0.25);
- vf::lbm::interpolateCF(icellF.TSE, omegaF, vf::basics::constant::c1o2, coefficients, 0.25, -0.25, 0.25);
- vf::lbm::interpolateCF(icellF.BNW, omegaF, vf::basics::constant::c1o2, coefficients, -0.25, 0.25, -0.25);
- vf::lbm::interpolateCF(icellF.BSE, omegaF, vf::basics::constant::c1o2, coefficients, 0.25, -0.25, -0.25);
- vf::lbm::interpolateCF(icellF.TSW, omegaF, vf::basics::constant::c1o2, coefficients, -0.25, -0.25, 0.25);
- vf::lbm::interpolateCF(icellF.TNE, omegaF, vf::basics::constant::c1o2, coefficients, 0.25, 0.25, 0.25);
+ vf::lbm::interpolateCF(icellF.BSW, omegaF, vf::basics::constant::c1o2, coefficients, -c1o4, -c1o4, -c1o4);
+ vf::lbm::interpolateCF(icellF.BNE, omegaF, vf::basics::constant::c1o2, coefficients, c1o4, c1o4, -c1o4);
+ vf::lbm::interpolateCF(icellF.TNW, omegaF, vf::basics::constant::c1o2, coefficients, -c1o4, c1o4, c1o4);
+ vf::lbm::interpolateCF(icellF.TSE, omegaF, vf::basics::constant::c1o2, coefficients, c1o4, -c1o4, c1o4);
+ vf::lbm::interpolateCF(icellF.BNW, omegaF, vf::basics::constant::c1o2, coefficients, -c1o4, c1o4, -c1o4);
+ vf::lbm::interpolateCF(icellF.BSE, omegaF, vf::basics::constant::c1o2, coefficients, c1o4, -c1o4, -c1o4);
+ vf::lbm::interpolateCF(icellF.TSW, omegaF, vf::basics::constant::c1o2, coefficients, -c1o4, -c1o4, c1o4);
+ vf::lbm::interpolateCF(icellF.TNE, omegaF, vf::basics::constant::c1o2, coefficients, c1o4, c1o4, c1o4);
}
void CompressibleOffsetMomentsInterpolator::interpolateFineToCoarse(D3Q27ICell& icellF, real* icellC, real xoff, real yoff, real zoff)
diff --git a/src/cpu/VirtualFluidsCore/LBM/Interpolation/IncompressibleOffsetInterpolator.cpp b/src/cpu/VirtualFluidsCore/LBM/Interpolation/IncompressibleOffsetInterpolator.cpp
index 0cd441255..a00c26619 100644
--- a/src/cpu/VirtualFluidsCore/LBM/Interpolation/IncompressibleOffsetInterpolator.cpp
+++ b/src/cpu/VirtualFluidsCore/LBM/Interpolation/IncompressibleOffsetInterpolator.cpp
@@ -435,56 +435,56 @@ void IncompressibleOffsetInterpolator::calcInterpolatedCoefficiets(const D3Q27IC
f_BE = eps_new*(-(ax - c3o1*az - c2o1*by - c3o1*cx + cz+c2o1*kxxMyyAverage-kxxMzzAverage-c3o1*kxzAverage)/(c54o1*o));
f_TN = eps_new*(-(-c2o1*ax + by + c3o1*bz + c3o1*cy + cz-kxxMyyAverage-kxxMzzAverage+c3o1*kyzAverage)/(c54o1*o));
f_BN = eps_new*(-(-c2o1*ax + by - c3o1*bz - c3o1*cy + cz-kxxMyyAverage-kxxMzzAverage-c3o1*kyzAverage)/(c54o1*o));
- f_ZERO = 0.;
+ f_ZERO = c0o1;
f_TNE = eps_new*(-(ay + az + bx + bz + cx + cy+kxyAverage+kxzAverage+kyzAverage)/(c72o1*o));
f_TSW = eps_new*((-ay + az - bx + bz + cx + cy-kxyAverage+kxzAverage+kyzAverage)/(c72o1*o));
f_TSE = eps_new*((ay - az + bx + bz - cx + cy+kxyAverage-kxzAverage+kyzAverage)/(c72o1*o));
f_TNW = eps_new*((ay + az + bx - bz + cx - cy+kxyAverage+kxzAverage-kyzAverage)/(c72o1*o));
- x_E = 0.25*eps_new*((c2o1*(-c4o1*axx + bxy + cxz))/(c27o1*o));
- x_N = 0.25*eps_new*((c2o1*(c2o1*axx - c2o1*bxy + cxz))/(c27o1*o));
- x_T = 0.25*eps_new*((c2o1*(c2o1*axx + bxy - c2o1*cxz))/(c27o1*o));
- x_NE = 0.25*eps_new*(-((c2o1*axx + c3o1*axy + c6o1*bxx + bxy - c2o1*cxz))/(c54o1*o));
- x_SE = 0.25*eps_new*(-((c2o1*axx - c3o1*axy - c6o1*bxx + bxy - c2o1*cxz))/(c54o1*o));
- x_TE = 0.25*eps_new*(-((c2o1*axx + c3o1*axz - c2o1*bxy + c6o1*cxx + cxz))/(c54o1*o));
- x_BE = 0.25*eps_new*(-((c2o1*axx - c3o1*axz - c2o1*bxy - c6o1*cxx + cxz))/(c54o1*o));
- x_TN = 0.25*eps_new*(-((-c4o1*axx + bxy + c3o1*bxz + c3o1*cxy + cxz))/(c54o1*o));
- x_BN = 0.25*eps_new*(-((-c4o1*axx + bxy - c3o1*bxz - c3o1*cxy + cxz))/(c54o1*o));
- x_ZERO = 0.;
- x_TNE = 0.25*eps_new*(-((axy + axz + c2o1*bxx + bxz + c2o1*cxx + cxy))/(c72o1*o));
- x_TSW = 0.25*eps_new*(((-axy + axz - c2o1*bxx + bxz + c2o1*cxx + cxy))/(c72o1*o));
- x_TSE = 0.25*eps_new*(((axy - axz + c2o1*bxx + bxz - c2o1*cxx + cxy))/(c72o1*o));
- x_TNW = 0.25*eps_new*(((axy + axz + c2o1*bxx - bxz + c2o1*cxx - cxy))/(c72o1*o));
-
- y_E = 0.25*eps_new*(c2o1*(-c2o1*axy + c2o1*byy + cyz))/(c27o1*o);
- y_N = 0.25*eps_new*(c2o1*(axy - c4o1*byy + cyz))/(c27o1*o);
- y_T = 0.25*eps_new*(c2o1*(axy + c2o1*byy - c2o1*cyz))/(c27o1*o);
- y_NE = 0.25*eps_new*(-((axy + c6o1*ayy + c3o1*bxy + c2o1*byy - c2o1*cyz))/(c54o1*o));
- y_SE = 0.25*eps_new*(-((axy - c6o1*ayy - c3o1*bxy + c2o1*byy - c2o1*cyz))/(c54o1*o));
- y_TE = 0.25*eps_new*(-((axy + c3o1*ayz - c4o1*byy + c3o1*cxy + cyz))/(c54o1*o));
- y_BE = 0.25*eps_new*(-((axy - c3o1*ayz - c4o1*byy - c3o1*cxy + cyz))/(c54o1*o));
- y_TN = 0.25*eps_new*(-((-c2o1*axy + c2o1*byy + c3o1*byz + c6o1*cyy + cyz))/(c54o1*o));
- y_BN = 0.25*eps_new*(-((-c2o1*axy + c2o1*byy - c3o1*byz - c6o1*cyy + cyz))/(c54o1*o));
- y_ZERO = 0.;
- y_TNE = 0.25*eps_new*(-((c2o1*ayy + ayz + bxy + byz + cxy + c2o1*cyy))/(c72o1*o));
- y_TSW = 0.25*eps_new*(((-c2o1*ayy + ayz - bxy + byz + cxy + c2o1*cyy))/(c72o1*o));
- y_TSE = 0.25*eps_new*(((c2o1*ayy - ayz + bxy + byz - cxy + c2o1*cyy))/(c72o1*o));
- y_TNW = 0.25*eps_new*(((c2o1*ayy + ayz + bxy - byz + cxy - c2o1*cyy))/(c72o1*o));
-
- z_E = 0.25*eps_new*((c2o1*(-c2o1*axz + byz + c2o1*czz))/(c27o1*o));
- z_N = 0.25*eps_new*((c2o1*(axz - c2o1*byz + c2o1*czz))/(c27o1*o));
- z_T = 0.25*eps_new*((c2o1*(axz + byz - c4o1*czz))/(c27o1*o));
- z_NE = 0.25*eps_new*(-((axz + c3o1*ayz + c3o1*bxz + byz - c4o1*czz))/(c54o1*o));
- z_SE = 0.25*eps_new*(-((axz - c3o1*ayz - c3o1*bxz + byz - c4o1*czz))/(c54o1*o));
- z_TE = 0.25*eps_new*(-((axz + c6o1*azz - c2o1*byz + c3o1*cxz + c2o1*czz))/(c54o1*o));
- z_BE = 0.25*eps_new*(-((axz - c6o1*azz - c2o1*byz - c3o1*cxz + c2o1*czz))/(c54o1*o));
- z_TN = 0.25*eps_new*(-((-c2o1*axz + byz + c6o1*bzz + c3o1*cyz + c2o1*czz))/(c54o1*o));
- z_BN = 0.25*eps_new*(-((-c2o1*axz + byz - c6o1*bzz - c3o1*cyz + c2o1*czz))/(c54o1*o));
- z_ZERO = 0.;
- z_TNE = 0.25*eps_new*(-((ayz + c2o1*azz + bxz + c2o1*bzz + cxz + cyz))/(c72o1*o));
- z_TSW = 0.25*eps_new*(((-ayz + c2o1*azz - bxz + c2o1*bzz + cxz + cyz))/(c72o1*o));
- z_TSE = 0.25*eps_new*(((ayz - c2o1*azz + bxz + c2o1*bzz - cxz + cyz))/(c72o1*o));
- z_TNW = 0.25*eps_new*(((ayz + c2o1*azz + bxz - c2o1*bzz + cxz - cyz))/(c72o1*o));
+ x_E = c1o4*eps_new*((c2o1*(-c4o1*axx + bxy + cxz))/(c27o1*o));
+ x_N = c1o4*eps_new*((c2o1*(c2o1*axx - c2o1*bxy + cxz))/(c27o1*o));
+ x_T = c1o4*eps_new*((c2o1*(c2o1*axx + bxy - c2o1*cxz))/(c27o1*o));
+ x_NE = c1o4*eps_new*(-((c2o1*axx + c3o1*axy + c6o1*bxx + bxy - c2o1*cxz))/(c54o1*o));
+ x_SE = c1o4*eps_new*(-((c2o1*axx - c3o1*axy - c6o1*bxx + bxy - c2o1*cxz))/(c54o1*o));
+ x_TE = c1o4*eps_new*(-((c2o1*axx + c3o1*axz - c2o1*bxy + c6o1*cxx + cxz))/(c54o1*o));
+ x_BE = c1o4*eps_new*(-((c2o1*axx - c3o1*axz - c2o1*bxy - c6o1*cxx + cxz))/(c54o1*o));
+ x_TN = c1o4*eps_new*(-((-c4o1*axx + bxy + c3o1*bxz + c3o1*cxy + cxz))/(c54o1*o));
+ x_BN = c1o4*eps_new*(-((-c4o1*axx + bxy - c3o1*bxz - c3o1*cxy + cxz))/(c54o1*o));
+ x_ZERO = c0o1;
+ x_TNE = c1o4*eps_new*(-((axy + axz + c2o1*bxx + bxz + c2o1*cxx + cxy))/(c72o1*o));
+ x_TSW = c1o4*eps_new*(((-axy + axz - c2o1*bxx + bxz + c2o1*cxx + cxy))/(c72o1*o));
+ x_TSE = c1o4*eps_new*(((axy - axz + c2o1*bxx + bxz - c2o1*cxx + cxy))/(c72o1*o));
+ x_TNW = c1o4*eps_new*(((axy + axz + c2o1*bxx - bxz + c2o1*cxx - cxy))/(c72o1*o));
+
+ y_E = c1o4*eps_new*(c2o1*(-c2o1*axy + c2o1*byy + cyz))/(c27o1*o);
+ y_N = c1o4*eps_new*(c2o1*(axy - c4o1*byy + cyz))/(c27o1*o);
+ y_T = c1o4*eps_new*(c2o1*(axy + c2o1*byy - c2o1*cyz))/(c27o1*o);
+ y_NE = c1o4*eps_new*(-((axy + c6o1*ayy + c3o1*bxy + c2o1*byy - c2o1*cyz))/(c54o1*o));
+ y_SE = c1o4*eps_new*(-((axy - c6o1*ayy - c3o1*bxy + c2o1*byy - c2o1*cyz))/(c54o1*o));
+ y_TE = c1o4*eps_new*(-((axy + c3o1*ayz - c4o1*byy + c3o1*cxy + cyz))/(c54o1*o));
+ y_BE = c1o4*eps_new*(-((axy - c3o1*ayz - c4o1*byy - c3o1*cxy + cyz))/(c54o1*o));
+ y_TN = c1o4*eps_new*(-((-c2o1*axy + c2o1*byy + c3o1*byz + c6o1*cyy + cyz))/(c54o1*o));
+ y_BN = c1o4*eps_new*(-((-c2o1*axy + c2o1*byy - c3o1*byz - c6o1*cyy + cyz))/(c54o1*o));
+ y_ZERO = c0o1;
+ y_TNE = c1o4*eps_new*(-((c2o1*ayy + ayz + bxy + byz + cxy + c2o1*cyy))/(c72o1*o));
+ y_TSW = c1o4*eps_new*(((-c2o1*ayy + ayz - bxy + byz + cxy + c2o1*cyy))/(c72o1*o));
+ y_TSE = c1o4*eps_new*(((c2o1*ayy - ayz + bxy + byz - cxy + c2o1*cyy))/(c72o1*o));
+ y_TNW = c1o4*eps_new*(((c2o1*ayy + ayz + bxy - byz + cxy - c2o1*cyy))/(c72o1*o));
+
+ z_E = c1o4*eps_new*((c2o1*(-c2o1*axz + byz + c2o1*czz))/(c27o1*o));
+ z_N = c1o4*eps_new*((c2o1*(axz - c2o1*byz + c2o1*czz))/(c27o1*o));
+ z_T = c1o4*eps_new*((c2o1*(axz + byz - c4o1*czz))/(c27o1*o));
+ z_NE = c1o4*eps_new*(-((axz + c3o1*ayz + c3o1*bxz + byz - c4o1*czz))/(c54o1*o));
+ z_SE = c1o4*eps_new*(-((axz - c3o1*ayz - c3o1*bxz + byz - c4o1*czz))/(c54o1*o));
+ z_TE = c1o4*eps_new*(-((axz + c6o1*azz - c2o1*byz + c3o1*cxz + c2o1*czz))/(c54o1*o));
+ z_BE = c1o4*eps_new*(-((axz - c6o1*azz - c2o1*byz - c3o1*cxz + c2o1*czz))/(c54o1*o));
+ z_TN = c1o4*eps_new*(-((-c2o1*axz + byz + c6o1*bzz + c3o1*cyz + c2o1*czz))/(c54o1*o));
+ z_BN = c1o4*eps_new*(-((-c2o1*axz + byz - c6o1*bzz - c3o1*cyz + c2o1*czz))/(c54o1*o));
+ z_ZERO = c0o1;
+ z_TNE = c1o4*eps_new*(-((ayz + c2o1*azz + bxz + c2o1*bzz + cxz + cyz))/(c72o1*o));
+ z_TSW = c1o4*eps_new*(((-ayz + c2o1*azz - bxz + c2o1*bzz + cxz + cyz))/(c72o1*o));
+ z_TSE = c1o4*eps_new*(((ayz - c2o1*azz + bxz + c2o1*bzz - cxz + cyz))/(c72o1*o));
+ z_TNW = c1o4*eps_new*(((ayz + c2o1*azz + bxz - c2o1*bzz + cxz - cyz))/(c72o1*o));
xy_E = c1o16*eps_new *(( c2o1*cxyz)/(c27o1*o));
xy_N = c1o16*eps_new *(( c2o1*cxyz)/(c27o1*o));
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp
index 89938deb9..720cb2809 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/IntegrateValuesHelper.cpp
@@ -174,6 +174,8 @@ void IntegrateValuesHelper::calculateAV()
//////////////////////////////////////////////////////////////////////////
void IntegrateValuesHelper::calculateMQ()
{
+ using namespace vf::basics::constant;
+
real f[D3Q27System::ENDF + 1];
real vx1, vx2, vx3, rho;
clearData();
@@ -186,7 +188,7 @@ void IntegrateValuesHelper::calculateMQ()
for (CalcNodes cn : cnodes) {
SPtr<ILBMKernel> kernel = cn.block->getKernel();
- real dx = 1.0 / (real)(1 << cn.block->getLevel());
+ real dx = c1o1 / (real)(1 << cn.block->getLevel());
real cellVolume = dx * dx * dx;
if (kernel->getCompressible()) {
diff --git a/src/cpu/VirtualFluidsCore/SimulationObservers/WriteBoundaryConditionsSimulationObserver.cpp b/src/cpu/VirtualFluidsCore/SimulationObservers/WriteBoundaryConditionsSimulationObserver.cpp
index 9d09db9e2..8bbe4cc95 100644
--- a/src/cpu/VirtualFluidsCore/SimulationObservers/WriteBoundaryConditionsSimulationObserver.cpp
+++ b/src/cpu/VirtualFluidsCore/SimulationObservers/WriteBoundaryConditionsSimulationObserver.cpp
@@ -132,6 +132,8 @@ void WriteBoundaryConditionsSimulationObserver::clearData()
//////////////////////////////////////////////////////////////////////////
void WriteBoundaryConditionsSimulationObserver::addDataGeo(SPtr<Block3D> block)
{
+ using namespace vf::basics::constant;
+
UbTupleDouble3 org = grid->getBlockWorldCoordinates(block);
UbTupleDouble3 nodeOffset = grid->getNodeOffset(block);
real dx = grid->getDeltaX(block);
@@ -178,22 +180,22 @@ void WriteBoundaryConditionsSimulationObserver::addDataGeo(SPtr<Block3D> block)
float(val<3>(org) - val<3>(nodeOffset) + ix3 * dx)));
if (!bcArray->hasBC(ix1, ix2, ix3)) {
- data[0].push_back(0.0);
+ data[0].push_back(c0o1);
} else if (bcArray->getBC(ix1, ix2, ix3)->hasNoSlipBoundary())
- data[0].push_back(1.0);
+ data[0].push_back(c1o1);
else if (bcArray->getBC(ix1, ix2, ix3)->hasVelocityBoundary())
- data[0].push_back(2.0);
+ data[0].push_back(c2o1);
else if (bcArray->getBC(ix1, ix2, ix3)->hasDensityBoundary())
- data[0].push_back(3.0);
+ data[0].push_back(c3o1);
else if (bcArray->getBC(ix1, ix2, ix3)->hasSlipBoundary())
- data[0].push_back(4.0);
+ data[0].push_back(c4o1);
// else
// data[0].push_back(5.0);
if (bcArray->isSolid(ix1, ix2, ix3)) {
- data[1].push_back(1.0);
+ data[1].push_back(c1o1);
} else {
- data[1].push_back(0.0);
+ data[1].push_back(c0o1);
}
data[2].push_back(level);
diff --git a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
index 68ebae1d6..0cc445b30 100644
--- a/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
+++ b/src/cpu/VirtualFluidsCore/Visitors/InitDistributionsBlockVisitor.cpp
@@ -182,7 +182,7 @@ void InitDistributionsBlockVisitor::visit(const SPtr<Grid3D> grid, SPtr<Block3D>
rho = muRho.Eval();
//x-derivative
- real deltaX=dx*0.5;
+ real deltaX=dx*c1o2;
x1 = coords[0]+deltaX;
real vx1Plusx1 = muVx1.Eval();
real vx2Plusx1 = muVx2.Eval();
--
GitLab