Commit e4ce1e40 authored by Ned Coltman's avatar Ned Coltman Committed by Melanie Lipp

[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
parent fb2a953c
......@@ -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;
......
......@@ -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_;
......
......@@ -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];
}
......
......@@ -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_;
......
......@@ -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];
}
......
......@@ -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)
......
......@@ -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);
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment