From e4ce1e4086ca94437a54f064294af20b7dc36a76 Mon Sep 17 00:00:00 2001
From: Ned Coltman <edward.coltman@iws.uni-stuttgart.de>
Date: Fri, 9 Oct 2020 09:42:56 +0200
Subject: [PATCH] [rans][walldistance] Rewrite wall distance calculation
* Include non- normal axis distances
* Cleanup wall and flow normal axis variables
* Include scvf corners as possible wall locations
---
.../freeflow/subcontrolvolumeface.hh | 2 +
dumux/freeflow/rans/problem.hh | 96 ++++++++++++-------
dumux/freeflow/rans/twoeq/kepsilon/problem.hh | 58 ++++++-----
dumux/freeflow/rans/volumevariables.hh | 8 +-
dumux/freeflow/rans/zeroeq/problem.hh | 23 ++---
dumux/freeflow/rans/zeroeq/volumevariables.hh | 4 +-
test/freeflow/rans/problem.hh | 2 +-
7 files changed, 107 insertions(+), 86 deletions(-)
diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
index 9805aa626a..6a7d6a76af 100644
--- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
+++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh
@@ -142,6 +142,8 @@ class FreeFlowStaggeredSubControlVolumeFace
public:
using GlobalPosition = typename T::GlobalPosition;
+ static constexpr int numCornersPerFace = 2 * (dim - 1);
+
//! State the traits public and thus export all types
using Traits = T;
diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh
index 1b5721d002..c46d30aa35 100644
--- a/dumux/freeflow/rans/problem.hh
+++ b/dumux/freeflow/rans/problem.hh
@@ -76,6 +76,7 @@ class RANSProblemBase : public NavierStokesProblem<TypeTag>
using GlobalPosition = typename SubControlVolumeFace::GlobalPosition;
static constexpr auto dim = GridView::dimension;
+ static constexpr int numCorners = SubControlVolumeFace::numCornersPerFace;
using DimVector = GlobalPosition;
using DimMatrix = Dune::FieldMatrix<Scalar, dim, dim>;
@@ -111,15 +112,19 @@ public:
velocityGradients_.resize(this->gridGeometry().elementMapper().size(), DimMatrix(0.0));
stressTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
vorticityTensorScalarProduct_.resize(this->gridGeometry().elementMapper().size(), 0.0);
- flowNormalAxis_.resize(this->gridGeometry().elementMapper().size(), 0);
- wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), 1);
+ flowNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedFlowNormalAxis_);
+ wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedWallNormalAxis_);
kinematicViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0);
sandGrainRoughness_.resize(this->gridGeometry().elementMapper().size(), 0.0);
- // retrieve all wall intersections and corresponding elements
- std::vector<unsigned int> wallElements;
- std::vector<GlobalPosition> wallPositions;
- std::vector<unsigned int> wallNormalAxisTemp;
+ // store the element indicies for all elements with an intersection on the wall
+ std::vector<unsigned int> wallElementIndicies;
+
+ // for each wall element, store the location of the face center and each corner.
+ std::vector<std::array<GlobalPosition, numCorners+1>> wallPositions;
+
+ // for each wall element, store the faces normal axis
+ std::vector<unsigned int> wallFaceNormalAxis;
const auto gridView = this->gridGeometry().gridView();
auto fvGeometry = localView(this->gridGeometry());
@@ -135,46 +140,51 @@ public:
if (asImp_().isOnWall(scvf))
{
- wallElements.push_back(this->gridGeometry().elementMapper().index(element));
- wallPositions.push_back(scvf.center());
- wallNormalAxisTemp.push_back(scvf.directionIndex());
+ // element has an scvf on the wall, store element index
+ wallElementIndicies.push_back(this->gridGeometry().elementMapper().index(element));
+
+ // store the location of the wall adjacent face's center and all corners
+ std::array<GlobalPosition, numCorners+1> wallElementPosition;
+ wallElementPosition[0] = scvf.center();
+ for (int i = 1; i <= numCorners; i++)
+ wallElementPosition[i] = scvf.corner(i-1);
+ wallPositions.push_back(wallElementPosition);
+
+ // Store the wall adjacent face's normal direction
+ wallFaceNormalAxis.push_back(scvf.directionIndex());
}
}
}
+ // output the number of wall adjacent faces. Check that this is non-zero.
std::cout << "NumWallIntersections=" << wallPositions.size() << std::endl;
if (wallPositions.size() == 0)
DUNE_THROW(Dune::InvalidStateException,
"No wall intersections have been found. Make sure that the isOnWall(globalPos) is working properly.");
- // search for shortest distance to wall for each element
+ // search for shortest distance to the wall for each element
for (const auto& element : elements(gridView))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
+
+ // Store the cell center position for each element
cellCenter_[elementIdx] = element.geometry().center();
+
for (unsigned int i = 0; i < wallPositions.size(); ++i)
{
- static const int problemWallNormalAxis
- = getParamFromGroup<int>(this->paramGroup(), "RANS.WallNormalAxis", -1);
- int searchAxis = problemWallNormalAxis;
-
- // search along wall normal axis of the intersection
- if (problemWallNormalAxis < 0 || problemWallNormalAxis >= dim)
- {
- searchAxis = wallNormalAxisTemp[i];
- }
+ // Find the minimum distance from the cell center to the wall face (center and corners)
+ std::array<Scalar,numCorners+1> cellToWallDistances;
+ for (unsigned int j = 0; j < wallPositions[i].size(); j++)
+ cellToWallDistances[j] = (cellCenter(elementIdx) - wallPositions[i][j]).two_norm();
+ Scalar distanceToWall = *std::min_element(cellToWallDistances.begin(), cellToWallDistances.end());
- GlobalPosition global = element.geometry().center();
- global -= wallPositions[i];
- // second and argument ensures to use only aligned elements
- if (abs(global[searchAxis]) < wallDistance_[elementIdx]
- && abs(global[searchAxis]) < global.two_norm() + 1e-8
- && abs(global[searchAxis]) > global.two_norm() - 1e-8)
+ if (distanceToWall < wallDistance_[elementIdx])
{
- wallDistance_[elementIdx] = abs(global[searchAxis]);
- wallElementIdx_[elementIdx] = wallElements[i];
- wallNormalAxis_[elementIdx] = searchAxis;
- sandGrainRoughness_[elementIdx] = asImp_().sandGrainRoughnessAtPos(wallPositions[i]);
+ wallDistance_[elementIdx] = distanceToWall;
+ wallElementIdx_[elementIdx] = wallElementIndicies[i];
+ if ( !(hasParam("RANS.WallNormalAxis")) )
+ wallNormalAxis_[elementIdx] = wallFaceNormalAxis[i];
+ sandGrainRoughness_[elementIdx] = asImp_().sandGrainRoughnessAtPos(wallPositions[i][0]);
}
}
}
@@ -232,9 +242,6 @@ public:
DUNE_THROW(Dune::InvalidStateException,
"You have to call updateStaticWallProperties() once before you call updateDynamicWallProperties().");
- static const int flowNormalAxis
- = getParamFromGroup<int>(this->paramGroup(), "RANS.FlowNormalAxis", -1);
-
// re-initialize min and max values
velocityMaximum_.assign(this->gridGeometry().elementMapper().size(), DimVector(1e-16));
velocityMinimum_.assign(this->gridGeometry().elementMapper().size(), DimVector(std::numeric_limits<Scalar>::max()));
@@ -529,19 +536,38 @@ public:
}
public:
+ int wallNormalAxis(const int elementIdx) const
+ { return wallNormalAxis_[elementIdx]; }
+
+ int flowNormalAxis(const int elementIdx) const
+ { return flowNormalAxis_[elementIdx]; }
+
+ unsigned int wallElementIndex(const int elementIdx) const
+ { return wallElementIdx_[elementIdx]; }
+
+ Scalar wallDistance(const int elementIdx) const
+ { return wallDistance_[elementIdx]; }
+
+
bool calledUpdateStaticWallProperties = false;
- std::vector<unsigned int> wallElementIdx_;
+
+ const int fixedFlowNormalAxis_ = getParam<int>("RANS.FlowNormalAxis", 0);
+ const int fixedWallNormalAxis_ = getParam<int>("RANS.WallNormalAxis", 1);
+
+ std::vector<unsigned int> wallNormalAxis_;
+ std::vector<unsigned int> flowNormalAxis_;
std::vector<Scalar> wallDistance_;
+ std::vector<unsigned int> wallElementIdx_;
std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIdx_;
std::vector<GlobalPosition> cellCenter_;
std::vector<DimVector> velocity_;
std::vector<DimVector> velocityMaximum_;
std::vector<DimVector> velocityMinimum_;
std::vector<DimMatrix> velocityGradients_;
+
std::vector<Scalar> stressTensorScalarProduct_;
std::vector<Scalar> vorticityTensorScalarProduct_;
- std::vector<unsigned int> wallNormalAxis_;
- std::vector<unsigned int> flowNormalAxis_;
+
std::vector<Scalar> kinematicViscosity_;
std::vector<Scalar> sandGrainRoughness_;
diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
index 1efb32c5b4..7b0382947f 100644
--- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
+++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh
@@ -144,17 +144,17 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
- unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementIdx];
- unsigned int neighborIdx0 = asImp_().neighborIdx_[elementIdx][wallNormalAxis][0];
- unsigned int neighborIdx1 = asImp_().neighborIdx_[elementIdx][wallNormalAxis][1];
+ unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
+ unsigned int neighborIndex0 = asImp_().neighborIndex(elementIdx, wallNormalAxis, 0);
+ unsigned int neighborIndex1 = asImp_().neighborIndex(elementIdx, wallNormalAxis, 1);
numElementsInNearWallRegion = inNearWallRegion(elementIdx)
? numElementsInNearWallRegion + 1
: numElementsInNearWallRegion + 0;
- if ((!inNearWallRegion(elementIdx) && (inNearWallRegion(neighborIdx0) || inNearWallRegion(neighborIdx1)))
- || (!inNearWallRegion(elementIdx) && elementIdx == asImp_().wallElementIdx_[elementIdx])
- || (inNearWallRegion(elementIdx) && (asImp_().wallElementIdx_[neighborIdx0] != asImp_().wallElementIdx_[neighborIdx1])))
+ if ((!inNearWallRegion(elementIdx) && (inNearWallRegion(neighborIndex0) || inNearWallRegion(neighborIndex1)))
+ || (!inNearWallRegion(elementIdx) && elementIdx == asImp_().wallElementIndex(elementIdx))
+ || (inNearWallRegion(elementIdx) && (asImp_().wallElementIndex(neighborIndex0) != asImp_().wallElementIndex(neighborIndex1))))
{
- matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]] = elementIdx;
+ matchingPointIdx_[asImp_().wallElementIndex(elementIdx)] = elementIdx;
}
}
std::cout << "numElementsInNearWallRegion: " << numElementsInNearWallRegion << std::endl;
@@ -174,7 +174,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
- unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]];
+ unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIndex(elementIdx)];
Scalar scalingFactor = storedDynamicEddyViscosity_[matchingPointIdx]
/ zeroEqDynamicEddyViscosity_[matchingPointIdx];
@@ -187,7 +187,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
- unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]];
+ unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIndex(elementIdx)];
if (isMatchingPoint(elementIdx))
{
zeroEqDynamicEddyViscosity_[matchingPointIdx] = storedDynamicEddyViscosity_[matchingPointIdx];
@@ -201,7 +201,7 @@ public:
*/
bool inNearWallRegion(unsigned int elementIdx) const
{
- unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx];
+ unsigned int wallElementIdx = asImp_().wallElementIndex(elementIdx);
unsigned int matchingPointIdx = matchingPointIdx_[wallElementIdx];
return wallElementIdx == matchingPointIdx ? yPlusNominal(elementIdx) < yPlusThreshold_
: yPlus(elementIdx) < yPlusThreshold_;
@@ -211,23 +211,23 @@ public:
* \brief Returns if an element is the matching point
*/
bool isMatchingPoint(unsigned int elementIdx) const
- { return matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]] == elementIdx; }
+ { return matchingPointIdx_[asImp_().wallElementIndex(elementIdx)] == elementIdx; }
/*!
* \brief Returns the \f$ y^+ \f$ value at an element center
*/
const Scalar yPlus(unsigned int elementIdx) const
{
- return asImp_().wallDistance_[elementIdx] * uStar(elementIdx)
- / asImp_().kinematicViscosity_[elementIdx];
+ return asImp_().wallDistance(elementIdx) * uStar(elementIdx)
+ / asImp_().kinematicViscosity(elementIdx);
}
/*!
* \brief Returns the nominal \f$ y^+ \f$ value at an element center
*/
const Scalar yPlusNominal(unsigned int elementIdx) const
{
- return asImp_().wallDistance_[elementIdx] * uStarNominal(elementIdx)
- / asImp_().kinematicViscosity_[elementIdx];
+ return asImp_().wallDistance(elementIdx) * uStarNominal(elementIdx)
+ / asImp_().kinematicViscosity(elementIdx);
}
/*!
@@ -244,14 +244,14 @@ public:
Scalar mixingLength = 0.0;
if (yPlusValue > 0.0)
{
- mixingLength = asImp_().karmanConstant() * asImp_().wallDistance_[elementIdx]
+ mixingLength = asImp_().karmanConstant() * asImp_().wallDistance(elementIdx)
* (1.0 - exp(-yPlusValue / 26.0 ))
/ sqrt(1.0 - exp(-0.26 * yPlusValue));
}
- unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementIdx];
- unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementIdx];
- Scalar velocityGradient = asImp_().velocityGradients_[elementIdx][flowNormalAxis][wallNormalAxis];
+ unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
+ unsigned int flowNormalAxis = asImp_().flowNormalAxis(elementIdx);
+ Scalar velocityGradient = asImp_().velocityGradient(elementIdx, flowNormalAxis, wallNormalAxis);
return mixingLength * mixingLength * abs(velocityGradient) * storedDensity_[elementIdx];
}
@@ -260,11 +260,11 @@ public:
{
using std::abs;
using std::sqrt;
- unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx];
- unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementIdx];
- unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementIdx];
- return sqrt(asImp_().kinematicViscosity_[wallElementIdx]
- * abs(asImp_().velocityGradients_[wallElementIdx][flowNormalAxis][wallNormalAxis]));
+ unsigned int wallElementIdx = asImp_().wallElementIndex(elementIdx);
+ unsigned int wallNormalAxis = asImp_().wallNormalAxis(elementIdx);
+ unsigned int flowNormalAxis = asImp_().flowNormalAxis(elementIdx);
+ return sqrt(asImp_().kinematicViscosity(wallElementIdx)
+ * abs(asImp_().velocityGradient(wallElementIdx, flowNormalAxis, wallNormalAxis)));
}
//! \brief Returns the nominal wall shear stress velocity (accounts for poor approximation of viscous sublayer)
@@ -272,9 +272,8 @@ public:
{
using std::pow;
using std::sqrt;
- unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]];
- return pow(cMu(), 0.25)
- * sqrt(storedTurbulentKineticEnergy_[matchingPointIdx]);
+ unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIndex(elementIdx)];
+ return pow(cMu(), 0.25) * sqrt(storedTurbulentKineticEnergy_[matchingPointIdx]);
}
/*!
@@ -283,7 +282,7 @@ public:
const Scalar dissipationWallFunction(unsigned int elementIdx) const
{
return uStarNominal(elementIdx) * uStarNominal(elementIdx) * uStarNominal(elementIdx)
- / asImp_().karmanConstant() / asImp_().wallDistance_[elementIdx];
+ / asImp_().karmanConstant() / asImp_().wallDistance(elementIdx);
}
/*!
@@ -291,8 +290,7 @@ public:
*/
const Scalar turbulentKineticEnergyWallFunction(unsigned int elementIdx) const
{
- unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx];
- unsigned int matchingPointIdx = matchingPointIdx_[wallElementIdx];
+ unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIndex(elementIdx)];
return storedTurbulentKineticEnergy_[matchingPointIdx];
}
diff --git a/dumux/freeflow/rans/volumevariables.hh b/dumux/freeflow/rans/volumevariables.hh
index e1ed768072..cecf13e0ff 100644
--- a/dumux/freeflow/rans/volumevariables.hh
+++ b/dumux/freeflow/rans/volumevariables.hh
@@ -93,15 +93,14 @@ public:
// calculate characteristic properties of the turbulent flow
elementIdx_ = problem.gridGeometry().elementMapper().index(element);
- wallElementIdx_ = problem.wallElementIdx_[elementIdx_];
- wallDistance_ = problem.wallDistance_[elementIdx_];
+ wallDistance_ = problem.wallDistance(elementIdx_);
velocity_ = problem.velocity_[elementIdx_];
velocityMaximum_ = problem.velocityMaximum_[wallElementIdx_];
velocityGradients_ = problem.velocityGradients_[elementIdx_];
- const auto flowNormalAxis = problem.flowNormalAxis_[elementIdx_];
- const auto wallNormalAxis = problem.wallNormalAxis_[elementIdx_];
uStar_ = sqrt(problem.kinematicViscosity_[wallElementIdx_]
* abs(problem.velocityGradients_[wallElementIdx_][flowNormalAxis][wallNormalAxis]));
+ const auto flowNormalAxis = problem.flowNormalAxis(elementIdx_);
+ const auto wallNormalAxis = problem.wallNormalAxis(elementIdx_);
uStar_ = max(uStar_, 1e-10); // zero values lead to numerical problems in some turbulence models
yPlus_ = wallDistance_ * uStar_ / problem.kinematicViscosity_[elementIdx_];
uPlus_ = velocity_[flowNormalAxis] / uStar_;
@@ -258,7 +257,6 @@ protected:
DimVector velocityMaximum_;
DimMatrix velocityGradients_;
std::size_t elementIdx_;
- std::size_t wallElementIdx_;
Scalar wallDistance_;
Scalar karmanConstant_;
Scalar uStar_;
diff --git a/dumux/freeflow/rans/zeroeq/problem.hh b/dumux/freeflow/rans/zeroeq/problem.hh
index 86e09ad0a7..b6479fbf26 100644
--- a/dumux/freeflow/rans/zeroeq/problem.hh
+++ b/dumux/freeflow/rans/zeroeq/problem.hh
@@ -134,14 +134,14 @@ public:
}
if (ksPlus > 2000.)
{
- std::cout << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << this->cellCenter_[this->wallElementIdx_[elementIdx]]
+ std::cout << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
<< " is not in the valid range (ksPlus < 2000),"
<< " for high ksPlus values the roughness function reaches a turning point."<< std::endl;
DUNE_THROW(Dune::InvalidStateException, "Unphysical roughness behavior.");
}
else if (ksPlus > 0.0 && ksPlus < 4.535 && !printedRangeWarning)
{
- Dune::dinfo << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << this->cellCenter_[this->wallElementIdx_[elementIdx]]
+ Dune::dinfo << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << asImp_().cellCenter(asImp_().wallElementIndex(elementIdx))
<< " is not in the valid range (ksPlus > 4.535) and now set to 0.0"<< std::endl;
ksPlus = 0.0;
printedRangeWarning = true;
@@ -186,10 +186,10 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
- unsigned int wallElementIdx = this->wallElementIdx_[elementIdx];
- Scalar wallDistance = this->wallDistance_[elementIdx] + additionalRoughnessLength_[elementIdx];
- unsigned int flowNormalAxis = this->flowNormalAxis_[elementIdx];
- unsigned int wallNormalAxis = this->wallNormalAxis_[elementIdx];
+ Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
+ unsigned int flowNormalAxis = this->flowNormalAxis(elementIdx);
+ unsigned int wallNormalAxis = this->wallNormalAxis(elementIdx);
+
Scalar omegaAbs = abs(this->velocityGradients_[elementIdx][flowNormalAxis][wallNormalAxis]
- this->velocityGradients_[elementIdx][wallNormalAxis][flowNormalAxis]);
@@ -211,8 +211,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
- unsigned int wallElementIdx = this->wallElementIdx_[elementIdx];
- Scalar wallDistance = this->wallDistance_[elementIdx] + additionalRoughnessLength_[elementIdx];
+ Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
Scalar maxVelocityNorm = 0.0;
Scalar minVelocityNorm = 0.0;
@@ -238,8 +237,7 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
- unsigned int wallElementIdx = this->wallElementIdx_[elementIdx];
- Scalar wallDistance = this->wallDistance_[elementIdx] + additionalRoughnessLength_[elementIdx];
+ Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
// checks if sign switches, by multiplication
Scalar check = kinematicEddyViscosityDifference[wallElementIdx] * kinematicEddyViscosityDifference[elementIdx];
@@ -254,11 +252,10 @@ public:
for (const auto& element : elements(this->gridGeometry().gridView()))
{
unsigned int elementIdx = this->gridGeometry().elementMapper().index(element);
- unsigned int wallElementIdx = this->wallElementIdx_[elementIdx];
- Scalar wallDistance = this->wallDistance_[elementIdx] + additionalRoughnessLength_[elementIdx];
+ Scalar effectiveWallDistance = asImp_().wallDistance(elementIdx) + additionalRoughnessLength_[elementIdx];
kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityInner[elementIdx];
- if (wallDistance >= switchingPosition[wallElementIdx])
+ if (effectiveWallDistance >= switchingPosition[asImp_().wallElementIndex(elementIdx)])
{
kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityOuter[elementIdx];
}
diff --git a/dumux/freeflow/rans/zeroeq/volumevariables.hh b/dumux/freeflow/rans/zeroeq/volumevariables.hh
index e8b3560be8..6a444925fd 100644
--- a/dumux/freeflow/rans/zeroeq/volumevariables.hh
+++ b/dumux/freeflow/rans/zeroeq/volumevariables.hh
@@ -108,8 +108,8 @@ public:
using std::exp;
using std::sqrt;
Scalar kinematicEddyViscosity = 0.0;
- unsigned int flowNormalAxis = problem.flowNormalAxis_[RANSParentType::elementIdx()];
- unsigned int wallNormalAxis = problem.wallNormalAxis_[RANSParentType::elementIdx()];
+ unsigned int flowNormalAxis = problem.flowNormalAxis(RANSParentType::elementIdx());
+ unsigned int wallNormalAxis = problem.wallNormalAxis(RANSParentType::elementIdx());
Scalar velGrad = abs(RANSParentType::velocityGradients()[flowNormalAxis][wallNormalAxis]);
if (modelName.compare("none") == 0)
diff --git a/test/freeflow/rans/problem.hh b/test/freeflow/rans/problem.hh
index 32fd52316f..926bf9386c 100644
--- a/test/freeflow/rans/problem.hh
+++ b/test/freeflow/rans/problem.hh
@@ -429,7 +429,7 @@ private:
{
static_assert(ModelTraits::turbulenceModel() == TurbulenceModel::komega, "Only valid for Komega");
// For the komega model we set a fixed value for the dissipation
- const auto wallDistance = ParentType::wallDistance_[elementIdx];
+ const auto wallDistance = ParentType::wallDistance(elementIdx);
using std::pow;
values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx]
/ (ParentType::betaOmega() * wallDistance * wallDistance);
--
GitLab