diff --git a/dumux/freeflow/rans/oneeq/problem.hh b/dumux/freeflow/rans/oneeq/problem.hh index d5a7df3d4b237a95ee879a8b116b6f400a1f5d55..6788b7506b083b4a6c6b570378d094365bf32381 100644 --- a/dumux/freeflow/rans/oneeq/problem.hh +++ b/dumux/freeflow/rans/oneeq/problem.hh @@ -96,7 +96,7 @@ public: for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); @@ -107,26 +107,26 @@ public: PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars); auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(priVars)); // NOTE: first update the turbulence quantities - storedViscosityTilde_[elementID] = elemSol[0][Indices::viscosityTildeIdx]; + storedViscosityTilde_[elementIdx] = elemSol[0][Indices::viscosityTildeIdx]; // NOTE: then update the volVars VolumeVariables volVars; volVars.update(elemSol, asImp_(), element, scv); - storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity(); + storedDynamicEddyViscosity_[elementIdx] = volVars.calculateEddyViscosity(); } } // calculate cell-center-averaged velocity gradients, maximum, and minimum values for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); for (unsigned int dimIdx = 0; dimIdx < Grid::dimension; ++dimIdx) { - storedViscosityTildeGradient_[elementID][dimIdx] - = (storedViscosityTilde_[ParentType::neighborID_[elementID][dimIdx][1]] - - storedViscosityTilde_[ParentType::neighborID_[elementID][dimIdx][0]]) - / (ParentType::cellCenter_[ParentType::neighborID_[elementID][dimIdx][1]][dimIdx] - - ParentType::cellCenter_[ParentType::neighborID_[elementID][dimIdx][0]][dimIdx]); + storedViscosityTildeGradient_[elementIdx][dimIdx] + = (storedViscosityTilde_[ParentType::neighborIdx_[elementIdx][dimIdx][1]] + - storedViscosityTilde_[ParentType::neighborIdx_[elementIdx][dimIdx][0]]) + / (ParentType::cellCenter_[ParentType::neighborIdx_[elementIdx][dimIdx][1]][dimIdx] + - ParentType::cellCenter_[ParentType::neighborIdx_[elementIdx][dimIdx][0]][dimIdx]); } auto fvGeometry = localView(this->fvGridGeometry()); @@ -139,13 +139,13 @@ public: // face Value Scalar dirichletViscosityTilde = asImp_().dirichlet(element, scvf)[Indices::viscosityTildeIdx]; - unsigned int neighborID = ParentType::neighborID_[elementID][normDim][0]; - if (scvf.center()[normDim] < ParentType::cellCenter_[elementID][normDim]) - neighborID = ParentType::neighborID_[elementID][normDim][1]; + unsigned int neighborIdx = ParentType::neighborIdx_[elementIdx][normDim][0]; + if (scvf.center()[normDim] < ParentType::cellCenter_[elementIdx][normDim]) + neighborIdx = ParentType::neighborIdx_[elementIdx][normDim][1]; - storedViscosityTildeGradient_[elementID][normDim] - = (storedViscosityTilde_[neighborID] - dirichletViscosityTilde) - / (ParentType::cellCenter_[neighborID][normDim] - scvf.center()[normDim]); + storedViscosityTildeGradient_[elementIdx][normDim] + = (storedViscosityTilde_[neighborIdx] - dirichletViscosityTilde) + / (ParentType::cellCenter_[neighborIdx][normDim] - scvf.center()[normDim]); } } } diff --git a/dumux/freeflow/rans/oneeq/volumevariables.hh b/dumux/freeflow/rans/oneeq/volumevariables.hh index 9cc779ade7e2846dda8e6ed2fd6948030b933961..112d23e94f0080ceadcc16b99d737259eb133ce8 100644 --- a/dumux/freeflow/rans/oneeq/volumevariables.hh +++ b/dumux/freeflow/rans/oneeq/volumevariables.hh @@ -86,12 +86,12 @@ public: { RANSParentType::updateRANSProperties(elemSol, problem, element, scv); viscosityTilde_ = elemSol[0][Indices::viscosityTildeIdx]; - storedViscosityTilde_ = problem.storedViscosityTilde_[RANSParentType::elementID()]; - storedViscosityTildeGradient_ = problem.storedViscosityTildeGradient_[RANSParentType::elementID()]; - stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()]; - vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct_[RANSParentType::elementID()]; + storedViscosityTilde_ = problem.storedViscosityTilde_[RANSParentType::elementIdx()]; + storedViscosityTildeGradient_ = problem.storedViscosityTildeGradient_[RANSParentType::elementIdx()]; + stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()]; + vorticityTensorScalarProduct_ = problem.vorticityTensorScalarProduct_[RANSParentType::elementIdx()]; if (problem.useStoredEddyViscosity_) - RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementID()]); + RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]); else RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity()); RANSParentType::calculateEddyDiffusivity(problem); diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh index 255a3163ef6169c12c47ae29c4b02240abb4113e..450a4d387f3a33705fe8d7ced4fab6073b46978c 100644 --- a/dumux/freeflow/rans/problem.hh +++ b/dumux/freeflow/rans/problem.hh @@ -92,9 +92,9 @@ public: calledUpdateStaticWallProperties = true; // update size and initial values of the global vectors - wallElementID_.resize(this->fvGridGeometry().elementMapper().size()); + wallElementIdx_.resize(this->fvGridGeometry().elementMapper().size()); wallDistance_.resize(this->fvGridGeometry().elementMapper().size(), std::numeric_limits<Scalar>::max()); - neighborID_.resize(this->fvGridGeometry().elementMapper().size()); + neighborIdx_.resize(this->fvGridGeometry().elementMapper().size()); cellCenter_.resize(this->fvGridGeometry().elementMapper().size(), GlobalPosition(0.0)); velocity_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(0.0)); velocityMaximum_.resize(this->fvGridGeometry().elementMapper().size(), DimVector(0.0)); @@ -141,8 +141,8 @@ public: // search for shortest distance to wall for each element for (const auto& element : elements(gridView)) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - cellCenter_[elementID] = element.geometry().center(); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + cellCenter_[elementIdx] = element.geometry().center(); for (unsigned int i = 0; i < wallPositions.size(); ++i) { static const int problemWallNormalAxis @@ -158,26 +158,26 @@ public: GlobalPosition global = element.geometry().center(); global -= wallPositions[i]; // second and argument ensures to use only aligned elements - if (abs(global[searchAxis]) < wallDistance_[elementID] + if (abs(global[searchAxis]) < wallDistance_[elementIdx] && abs(global[searchAxis]) < global.two_norm() + 1e-8 && abs(global[searchAxis]) > global.two_norm() - 1e-8) { - wallDistance_[elementID] = abs(global[searchAxis]); - wallElementID_[elementID] = wallElements[i]; - wallNormalAxis_[elementID] = searchAxis; - sandGrainRoughness_[elementID] = asImp_().sandGrainRoughnessAtPos(wallPositions[i]); + wallDistance_[elementIdx] = abs(global[searchAxis]); + wallElementIdx_[elementIdx] = wallElements[i]; + wallNormalAxis_[elementIdx] = searchAxis; + sandGrainRoughness_[elementIdx] = asImp_().sandGrainRoughnessAtPos(wallPositions[i]); } } } - // search for neighbor IDs + // search for neighbor Idxs for (const auto& element : elements(gridView)) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) { - neighborID_[elementID][dimIdx][0] = elementID; - neighborID_[elementID][dimIdx][1] = elementID; + neighborIdx_[elementIdx][dimIdx][0] = elementIdx; + neighborIdx_[elementIdx][dimIdx][1] = elementIdx; } for (const auto& intersection : intersections(gridView, element)) @@ -185,18 +185,18 @@ public: if (intersection.boundary()) continue; - unsigned int neighborID = this->fvGridGeometry().elementMapper().index(intersection.outside()); + unsigned int neighborIdx = this->fvGridGeometry().elementMapper().index(intersection.outside()); for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) { - if (abs(cellCenter_[elementID][dimIdx] - cellCenter_[neighborID][dimIdx]) > 1e-8) + if (abs(cellCenter_[elementIdx][dimIdx] - cellCenter_[neighborIdx][dimIdx]) > 1e-8) { - if (cellCenter_[elementID][dimIdx] > cellCenter_[neighborID][dimIdx]) + if (cellCenter_[elementIdx][dimIdx] > cellCenter_[neighborIdx][dimIdx]) { - neighborID_[elementID][dimIdx][0] = neighborID; + neighborIdx_[elementIdx][dimIdx][0] = neighborIdx; } - if (cellCenter_[elementID][dimIdx] < cellCenter_[neighborID][dimIdx]) + if (cellCenter_[elementIdx][dimIdx] < cellCenter_[neighborIdx][dimIdx]) { - neighborID_[elementID][dimIdx][1] = neighborID; + neighborIdx_[elementIdx][dimIdx][1] = neighborIdx; } } } @@ -235,7 +235,7 @@ public: { auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); // calculate velocities DimVector velocityTemp(0.0); @@ -246,44 +246,44 @@ public: velocityTemp[scvf.directionIndex()] += numericalSolutionFace; } for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) - velocity_[elementID][dimIdx] = velocityTemp[dimIdx] * 0.5; // faces are equidistant to cell center + velocity_[elementIdx][dimIdx] = velocityTemp[dimIdx] * 0.5; // faces are equidistant to cell center } // calculate cell-center-averaged velocity gradients, maximum, and minimum values for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int wallElementID = wallElementID_[elementID]; + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + unsigned int wallElementIdx = wallElementIdx_[elementIdx]; Scalar maxVelocity = 0.0; for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) { for (unsigned int velIdx = 0; velIdx < dim; ++velIdx) { - velocityGradients_[elementID][velIdx][dimIdx] - = (velocity_[neighborID_[elementID][dimIdx][1]][velIdx] - - velocity_[neighborID_[elementID][dimIdx][0]][velIdx]) - / (cellCenter_[neighborID_[elementID][dimIdx][1]][dimIdx] - - cellCenter_[neighborID_[elementID][dimIdx][0]][dimIdx]); + velocityGradients_[elementIdx][velIdx][dimIdx] + = (velocity_[neighborIdx_[elementIdx][dimIdx][1]][velIdx] + - velocity_[neighborIdx_[elementIdx][dimIdx][0]][velIdx]) + / (cellCenter_[neighborIdx_[elementIdx][dimIdx][1]][dimIdx] + - cellCenter_[neighborIdx_[elementIdx][dimIdx][0]][dimIdx]); } - if (abs(velocity_[elementID][dimIdx]) > abs(velocityMaximum_[wallElementID][dimIdx])) + if (abs(velocity_[elementIdx][dimIdx]) > abs(velocityMaximum_[wallElementIdx][dimIdx])) { - velocityMaximum_[wallElementID][dimIdx] = velocity_[elementID][dimIdx]; + velocityMaximum_[wallElementIdx][dimIdx] = velocity_[elementIdx][dimIdx]; } - if (abs(velocity_[elementID][dimIdx]) < abs(velocityMinimum_[wallElementID][dimIdx])) + if (abs(velocity_[elementIdx][dimIdx]) < abs(velocityMinimum_[wallElementIdx][dimIdx])) { - velocityMinimum_[wallElementID][dimIdx] = velocity_[elementID][dimIdx]; + velocityMinimum_[wallElementIdx][dimIdx] = velocity_[elementIdx][dimIdx]; } if (0 <= flowNormalAxis && flowNormalAxis < dim) { - flowNormalAxis_[elementID] = flowNormalAxis; + flowNormalAxis_[elementIdx] = flowNormalAxis; } - else if (abs(maxVelocity) < abs(velocity_[elementID][dimIdx])) + else if (abs(maxVelocity) < abs(velocity_[elementIdx][dimIdx])) { - maxVelocity = abs(velocity_[elementID][dimIdx]); - flowNormalAxis_[elementID] = dimIdx; + maxVelocity = abs(velocity_[elementIdx][dimIdx]); + flowNormalAxis_[elementIdx] = dimIdx; } } @@ -302,13 +302,13 @@ public: Scalar dirichletVelocity = asImp_().dirichlet(element, scvf)[Indices::velocity(velIdx)]; - unsigned int neighborID = neighborID_[elementID][scvfNormDim][0]; - if (scvf.center()[scvfNormDim] < cellCenter_[elementID][scvfNormDim]) - neighborID = neighborID_[elementID][scvfNormDim][1]; + unsigned int neighborIdx = neighborIdx_[elementIdx][scvfNormDim][0]; + if (scvf.center()[scvfNormDim] < cellCenter_[elementIdx][scvfNormDim]) + neighborIdx = neighborIdx_[elementIdx][scvfNormDim][1]; - velocityGradients_[elementID][velIdx][scvfNormDim] - = (velocity_[neighborID][velIdx] - dirichletVelocity) - / (cellCenter_[neighborID][scvfNormDim] - scvf.center()[scvfNormDim]); + velocityGradients_[elementIdx][velIdx][scvfNormDim] + = (velocity_[neighborIdx][velIdx] - dirichletVelocity) + / (cellCenter_[neighborIdx][scvfNormDim] - scvf.center()[scvfNormDim]); } } @@ -327,14 +327,14 @@ public: unsigned int normalNormDim = normalFace.directionIndex(); if (normalFace.boundary() && (asImp_().boundaryTypes(element, normalFace).isBJS(Indices::velocity(velIdx)))) { - unsigned int neighborID = neighborID_[elementID][normalNormDim][0]; - if (normalFace.center()[normalNormDim] < cellCenter_[elementID][normalNormDim]) - neighborID = neighborID_[elementID][normalNormDim][1]; - - bjsVelocityAverage[normalNormDim] += ParentType::bjsVelocity(scvf, normalFace, localSubFaceIdx, velocity_[elementID][velIdx]); - if (bjsNumFaces[normalNormDim] > 0 && neighborID != bjsNeighbor[normalNormDim]) - DUNE_THROW(Dune::InvalidStateException, "Two different neighborID should not occur"); - bjsNeighbor[normalNormDim] = neighborID; + unsigned int neighborIdx = neighborIdx_[elementIdx][normalNormDim][0]; + if (normalFace.center()[normalNormDim] < cellCenter_[elementIdx][normalNormDim]) + neighborIdx = neighborIdx_[elementIdx][normalNormDim][1]; + + bjsVelocityAverage[normalNormDim] += ParentType::bjsVelocity(scvf, normalFace, localSubFaceIdx, velocity_[elementIdx][velIdx]); + if (bjsNumFaces[normalNormDim] > 0 && neighborIdx != bjsNeighbor[normalNormDim]) + DUNE_THROW(Dune::InvalidStateException, "Two different neighborIdx should not occur"); + bjsNeighbor[normalNormDim] = neighborIdx; normalNormCoordinate[normalNormDim] = normalFace.center()[normalNormDim]; bjsNumFaces[normalNormDim]++; } @@ -344,12 +344,12 @@ public: if (bjsNumFaces[dirIdx] == 0) continue; - unsigned int neighborID = bjsNeighbor[dirIdx]; + unsigned int neighborIdx = bjsNeighbor[dirIdx]; bjsVelocityAverage[dirIdx] /= bjsNumFaces[dirIdx]; - velocityGradients_[elementID][velIdx][dirIdx] - = (velocity_[neighborID][velIdx] - bjsVelocityAverage[dirIdx]) - / (cellCenter_[neighborID][dirIdx] - normalNormCoordinate[dirIdx]); + velocityGradients_[elementIdx][velIdx][dirIdx] + = (velocity_[neighborIdx][velIdx] - bjsVelocityAverage[dirIdx]) + / (cellCenter_[neighborIdx][dirIdx] - normalNormCoordinate[dirIdx]); } } @@ -358,23 +358,23 @@ public: // calculate or call all secondary variables for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); Dune::FieldMatrix<Scalar, GridView::dimension, GridView::dimension> stressTensor(0.0); for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) { for (unsigned int velIdx = 0; velIdx < dim; ++velIdx) { - stressTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementID][dimIdx][velIdx] - + 0.5 * velocityGradients_[elementID][velIdx][dimIdx]; + stressTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementIdx][dimIdx][velIdx] + + 0.5 * velocityGradients_[elementIdx][velIdx][dimIdx]; } } - stressTensorScalarProduct_[elementID] = 0.0; + stressTensorScalarProduct_[elementIdx] = 0.0; for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) { for (unsigned int velIdx = 0; velIdx < dim; ++velIdx) { - stressTensorScalarProduct_[elementID] += stressTensor[dimIdx][velIdx] * stressTensor[dimIdx][velIdx]; + stressTensorScalarProduct_[elementIdx] += stressTensor[dimIdx][velIdx] * stressTensor[dimIdx][velIdx]; } } @@ -383,16 +383,16 @@ public: { for (unsigned int velIdx = 0; velIdx < dim; ++velIdx) { - vorticityTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementID][dimIdx][velIdx] - - 0.5 * velocityGradients_[elementID][velIdx][dimIdx]; + vorticityTensor[dimIdx][velIdx] = 0.5 * velocityGradients_[elementIdx][dimIdx][velIdx] + - 0.5 * velocityGradients_[elementIdx][velIdx][dimIdx]; } } - vorticityTensorScalarProduct_[elementID] = 0.0; + vorticityTensorScalarProduct_[elementIdx] = 0.0; for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) { for (unsigned int velIdx = 0; velIdx < dim; ++velIdx) { - vorticityTensorScalarProduct_[elementID] += vorticityTensor[dimIdx][velIdx] * vorticityTensor[dimIdx][velIdx]; + vorticityTensorScalarProduct_[elementIdx] += vorticityTensor[dimIdx][velIdx] * vorticityTensor[dimIdx][velIdx]; } } @@ -409,7 +409,7 @@ public: VolumeVariables volVars; volVars.update(elemSol, asImp_(), element, scv); - kinematicViscosity_[elementID] = volVars.viscosity() / volVars.density(); + kinematicViscosity_[elementIdx] = volVars.viscosity() / volVars.density(); } } } @@ -466,9 +466,9 @@ public: public: bool calledUpdateStaticWallProperties = false; - std::vector<unsigned int> wallElementID_; + std::vector<unsigned int> wallElementIdx_; std::vector<Scalar> wallDistance_; - std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborID_; + std::vector<std::array<std::array<unsigned int, 2>, dim>> neighborIdx_; std::vector<GlobalPosition> cellCenter_; std::vector<DimVector> velocity_; std::vector<DimVector> velocityMaximum_; diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh index c09cbc9cf9b1f711e45f99f5d647be8dde5b43c8..8333634b667667ae275829d1ca840dab697b7ef8 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh @@ -98,7 +98,7 @@ public: ParentType::updateStaticWallProperties(); // update size and initial values of the global vectors - matchingPointID_.resize(this->fvGridGeometry().elementMapper().size(), 0); + matchingPointIdx_.resize(this->fvGridGeometry().elementMapper().size(), 0); storedDensity_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); storedDissipation_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); storedTurbulentKineticEnergy_.resize(this->fvGridGeometry().elementMapper().size(), 0.0); @@ -118,7 +118,7 @@ public: // update the stored eddy viscosities for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); @@ -129,13 +129,13 @@ public: PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars); auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(priVars)); // NOTE: first update the turbulence quantities - storedDissipation_[elementID] = elemSol[0][Indices::dissipationEqIdx]; - storedTurbulentKineticEnergy_[elementID] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; + storedDissipation_[elementIdx] = elemSol[0][Indices::dissipationEqIdx]; + storedTurbulentKineticEnergy_[elementIdx] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; // NOTE: then update the volVars VolumeVariables volVars; volVars.update(elemSol, asImp_(), element, scv); - storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity(); - storedDensity_[elementID] = volVars.density(); + storedDynamicEddyViscosity_[elementIdx] = volVars.calculateEddyViscosity(); + storedDensity_[elementIdx] = volVars.density(); } } @@ -143,18 +143,18 @@ public: unsigned int numElementsInNearWallRegion = 0; for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementID]; - unsigned int neighborID0 = asImp_().neighborID_[elementID][wallNormalAxis][0]; - unsigned int neighborID1 = asImp_().neighborID_[elementID][wallNormalAxis][1]; - numElementsInNearWallRegion = inNearWallRegion(elementID) + unsigned int elementIdx = this->fvGridGeometry().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]; + numElementsInNearWallRegion = inNearWallRegion(elementIdx) ? numElementsInNearWallRegion + 1 : numElementsInNearWallRegion + 0; - if ((!inNearWallRegion(elementID) && (inNearWallRegion(neighborID0) || inNearWallRegion(neighborID1))) - || (!inNearWallRegion(elementID) && elementID == asImp_().wallElementID_[elementID]) - || (inNearWallRegion(elementID) && (asImp_().wallElementID_[neighborID0] != asImp_().wallElementID_[neighborID1]))) + if ((!inNearWallRegion(elementIdx) && (inNearWallRegion(neighborIdx0) || inNearWallRegion(neighborIdx1))) + || (!inNearWallRegion(elementIdx) && elementIdx == asImp_().wallElementIdx_[elementIdx]) + || (inNearWallRegion(elementIdx) && (asImp_().wallElementIdx_[neighborIdx0] != asImp_().wallElementIdx_[neighborIdx1]))) { - matchingPointID_[asImp_().wallElementID_[elementID]] = elementID; + matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]] = elementIdx; } } std::cout << "numElementsInNearWallRegion: " << numElementsInNearWallRegion << std::endl; @@ -162,8 +162,8 @@ public: // calculate the potential zeroeq eddy viscosities for two-layer model for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - zeroEqDynamicEddyViscosity_[elementID] = zeroEqEddyViscosityModel(elementID); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + zeroEqDynamicEddyViscosity_[elementIdx] = zeroEqEddyViscosityModel(elementIdx); } // then make them match at the matching point @@ -173,24 +173,24 @@ public: { for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int matchingPointID = matchingPointID_[asImp_().wallElementID_[elementID]]; + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]]; - Scalar scalingFactor = storedDynamicEddyViscosity_[matchingPointID] - / zeroEqDynamicEddyViscosity_[matchingPointID]; - if (!isMatchingPoint(elementID) + Scalar scalingFactor = storedDynamicEddyViscosity_[matchingPointIdx] + / zeroEqDynamicEddyViscosity_[matchingPointIdx]; + if (!isMatchingPoint(elementIdx) && !std::isnan(scalingFactor) && !std::isinf(scalingFactor)) { - zeroEqDynamicEddyViscosity_[elementID] *= scalingFactor; + zeroEqDynamicEddyViscosity_[elementIdx] *= scalingFactor; } } for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int matchingPointID = matchingPointID_[asImp_().wallElementID_[elementID]]; - if (isMatchingPoint(elementID)) + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]]; + if (isMatchingPoint(elementIdx)) { - zeroEqDynamicEddyViscosity_[matchingPointID] = storedDynamicEddyViscosity_[matchingPointID]; + zeroEqDynamicEddyViscosity_[matchingPointIdx] = storedDynamicEddyViscosity_[matchingPointIdx]; } } } @@ -199,109 +199,109 @@ public: /*! * \brief Returns if an element is located in the near-wall region */ - const bool inNearWallRegion(unsigned int elementID) const + const bool inNearWallRegion(unsigned int elementIdx) const { - unsigned int wallElementID = asImp_().wallElementID_[elementID]; - unsigned int matchingPointID = matchingPointID_[wallElementID]; - return wallElementID == matchingPointID ? yPlusNominal(elementID) < yPlusThreshold_ - : yPlus(elementID) < yPlusThreshold_; + unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx]; + unsigned int matchingPointIdx = matchingPointIdx_[wallElementIdx]; + return wallElementIdx == matchingPointIdx ? yPlusNominal(elementIdx) < yPlusThreshold_ + : yPlus(elementIdx) < yPlusThreshold_; } /*! * \brief Returns if an element is the matching point */ - const bool isMatchingPoint(unsigned int elementID) const - { return matchingPointID_[asImp_().wallElementID_[elementID]] == elementID; } + const bool isMatchingPoint(unsigned int elementIdx) const + { return matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]] == elementIdx; } /*! * \brief Returns the \f$ y^+ \f$ value at an element center */ - const Scalar yPlus(unsigned int elementID) const + const Scalar yPlus(unsigned int elementIdx) const { - return asImp_().wallDistance_[elementID] * uStar(elementID) - / asImp_().kinematicViscosity_[elementID]; + 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 elementID) const + const Scalar yPlusNominal(unsigned int elementIdx) const { - return asImp_().wallDistance_[elementID] * uStarNominal(elementID) - / asImp_().kinematicViscosity_[elementID]; + return asImp_().wallDistance_[elementIdx] * uStarNominal(elementIdx) + / asImp_().kinematicViscosity_[elementIdx]; } /*! * \brief Returns the kinematic eddy viscosity of a 0-Eq. model */ - const Scalar zeroEqEddyViscosityModel(unsigned int elementID) const + const Scalar zeroEqEddyViscosityModel(unsigned int elementIdx) const { using std::abs; using std::exp; using std::sqrt; // use VanDriest's model - Scalar yPlusValue = yPlus(elementID); + Scalar yPlusValue = yPlus(elementIdx); Scalar mixingLength = 0.0; if (yPlusValue > 0.0) { - mixingLength = asImp_().karmanConstant() * asImp_().wallDistance_[elementID] + mixingLength = asImp_().karmanConstant() * asImp_().wallDistance_[elementIdx] * (1.0 - exp(-yPlusValue / 26.0 )) / sqrt(1.0 - exp(-0.26 * yPlusValue)); } - unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementID]; - unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementID]; - Scalar velocityGradient = asImp_().velocityGradients_[elementID][flowNormalAxis][wallNormalAxis]; - return mixingLength * mixingLength * abs(velocityGradient) * storedDensity_[elementID]; + unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementIdx]; + unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementIdx]; + Scalar velocityGradient = asImp_().velocityGradients_[elementIdx][flowNormalAxis][wallNormalAxis]; + return mixingLength * mixingLength * abs(velocityGradient) * storedDensity_[elementIdx]; } //! \brief Returns the wall shear stress velocity - const Scalar uStar(unsigned int elementID) const + const Scalar uStar(unsigned int elementIdx) const { using std::abs; using std::sqrt; - unsigned int wallElementID = asImp_().wallElementID_[elementID]; - unsigned int wallNormalAxis = asImp_().wallNormalAxis_[elementID]; - unsigned int flowNormalAxis = asImp_().flowNormalAxis_[elementID]; - return sqrt(asImp_().kinematicViscosity_[wallElementID] - * abs(asImp_().velocityGradients_[wallElementID][flowNormalAxis][wallNormalAxis])); + 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])); } //! \brief Returns the nominal wall shear stress velocity (accounts for poor approximation of viscous sublayer) - const Scalar uStarNominal(unsigned int elementID) const + const Scalar uStarNominal(unsigned int elementIdx) const { using std::pow; using std::sqrt; - unsigned int matchingPointID = matchingPointID_[asImp_().wallElementID_[elementID]]; + unsigned int matchingPointIdx = matchingPointIdx_[asImp_().wallElementIdx_[elementIdx]]; return pow(cMu(), 0.25) - * sqrt(storedTurbulentKineticEnergy_[matchingPointID]); + * sqrt(storedTurbulentKineticEnergy_[matchingPointIdx]); } /*! * \brief Returns the dissipation calculated from the wall function consideration */ - const Scalar dissipationWallFunction(unsigned int elementID) const + const Scalar dissipationWallFunction(unsigned int elementIdx) const { - return uStarNominal(elementID) * uStarNominal(elementID) * uStarNominal(elementID) - / asImp_().karmanConstant() / asImp_().wallDistance_[elementID]; + return uStarNominal(elementIdx) * uStarNominal(elementIdx) * uStarNominal(elementIdx) + / asImp_().karmanConstant() / asImp_().wallDistance_[elementIdx]; } /*! * \brief Returns the turbulentKineticEnergy calculated from the wall function consideration */ - const Scalar turbulentKineticEnergyWallFunction(unsigned int elementID) const + const Scalar turbulentKineticEnergyWallFunction(unsigned int elementIdx) const { - unsigned int wallElementID = asImp_().wallElementID_[elementID]; - unsigned int matchingPointID = matchingPointID_[wallElementID]; - return storedTurbulentKineticEnergy_[matchingPointID]; + unsigned int wallElementIdx = asImp_().wallElementIdx_[elementIdx]; + unsigned int matchingPointIdx = matchingPointIdx_[wallElementIdx]; + return storedTurbulentKineticEnergy_[matchingPointIdx]; } //! \brief Returns the nominal wall shear stress (accounts for poor approximation of viscous sublayer) - const Scalar tangentialMomentumWallFunction(unsigned int elementID, Scalar velocity) const + const Scalar tangentialMomentumWallFunction(unsigned int elementIdx, Scalar velocity) const { using std::log; - Scalar velocityNominal = uStarNominal(elementID) * (1.0 / asImp_().karmanConstant() * log(yPlusNominal(elementID)) + 5.0); - return uStarNominal(elementID) * uStarNominal(elementID) + Scalar velocityNominal = uStarNominal(elementIdx) * (1.0 / asImp_().karmanConstant() * log(yPlusNominal(elementIdx)) + 5.0); + return uStarNominal(elementIdx) * uStarNominal(elementIdx) * velocity / velocityNominal; } @@ -310,11 +310,11 @@ public: const SubControlVolumeFace& localSubFace, const int& eqIdx) const { - unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = asImp_().fvGridGeometry().elementMapper().index(element); auto bcTypes = asImp_().boundaryTypes(element, localSubFace); return asImp_().isOnWall(localSubFace.center()) && bcTypes.isDirichlet(eqIdx) - && isMatchingPoint(elementID); + && isMatchingPoint(elementIdx); } //! \brief Returns an additional wall function momentum flux (only needed for RANS models) @@ -325,8 +325,8 @@ public: const SubControlVolumeFace& scvf, const SubControlVolumeFace& localSubFace) const { - unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element); - return FacePrimaryVariables(asImp_().tangentialMomentumWallFunction(elementID, elemFaceVars[scvf].velocitySelf()) + unsigned int elementIdx = asImp_().fvGridGeometry().elementMapper().index(element); + return FacePrimaryVariables(asImp_().tangentialMomentumWallFunction(elementIdx, elemFaceVars[scvf].velocitySelf()) * elemVolVars[scvf.insideScvIdx()].density()); } @@ -382,7 +382,7 @@ public: { using std::log; auto wallFunctionFlux = CellCenterPrimaryVariables(0.0); - unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = asImp_().fvGridGeometry().elementMapper().index(element); // component mass fluxes for (int compIdx = 0; compIdx < ModelTraits::numComponents(); ++compIdx) @@ -399,9 +399,9 @@ public: - elemVolVars[scvf.insideScvIdx()].moleFraction(compIdx)) * elemVolVars[scvf.insideScvIdx()].molarDensity() * moleToMassConversionFactor - * uStarNominal(elementID) + * uStarNominal(elementIdx) / asImp_().turbulentSchmidtNumber() - / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementID) * 9.793) + / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementIdx) * 9.793) + pFunction(schmidtNumber, asImp_().turbulentSchmidtNumber())); } @@ -423,7 +423,7 @@ public: { using std::log; auto wallFunctionFlux = CellCenterPrimaryVariables(0.0); - unsigned int elementID = asImp_().fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = asImp_().fvGridGeometry().elementMapper().index(element); // energy fluxes Scalar prandtlNumber = elemVolVars[scvf.insideScvIdx()].kinematicViscosity() * elemVolVars[scvf.insideScvIdx()].density() @@ -434,9 +434,9 @@ public: - elemVolVars[scvf.insideScvIdx()].temperature()) * elemVolVars[scvf.insideScvIdx()].density() * elemVolVars[scvf.insideScvIdx()].heatCapacity() - * uStarNominal(elementID) + * uStarNominal(elementIdx) / asImp_().turbulentPrandtlNumber() - / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementID) * 9.793) + / (1. / asImp_().karmanConstant() * log(yPlusNominal(elementIdx) * 9.793) + pFunction(prandtlNumber, asImp_().turbulentPrandtlNumber())); return wallFunctionFlux; @@ -459,7 +459,7 @@ public: } public: - std::vector<unsigned int> matchingPointID_; + std::vector<unsigned int> matchingPointIdx_; std::vector<Scalar> storedDensity_; std::vector<Scalar> storedDissipation_; std::vector<Scalar> storedTurbulentKineticEnergy_; diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh index 86bdc00b9db397a84d5cf79b28814d98a3e5879e..f17c61dc41c2eec40431e4641e72971f8771db0b 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/localresidual.hh @@ -152,7 +152,7 @@ protected: elemVolVars, elemFaceVars, elemBcTypes, elemFluxVarsCache); for (auto&& scvf : scvfs(fvGeometry)) { - unsigned int elementID = problem.fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = problem.fvGridGeometry().elementMapper().index(element); const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); const auto& insideVolVars = elemVolVars[insideScv]; @@ -160,14 +160,14 @@ protected: if (insideVolVars.inNearWallRegion()) { residual[Indices::turbulentKineticEnergyEqIdx - cellCenterOffset] - = insideVolVars.turbulentKineticEnergy() - problem.turbulentKineticEnergyWallFunction(elementID); + = insideVolVars.turbulentKineticEnergy() - problem.turbulentKineticEnergyWallFunction(elementIdx); } // fixed value for the dissipation if (insideVolVars.inNearWallRegion() || insideVolVars.isMatchingPoint()) { residual[Indices::dissipationEqIdx - cellCenterOffset] - = insideVolVars.dissipation() - problem.dissipationWallFunction(elementID); + = insideVolVars.dissipation() - problem.dissipationWallFunction(elementIdx); } } } diff --git a/dumux/freeflow/rans/twoeq/kepsilon/volumevariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/volumevariables.hh index 329a43a8e0a2fa5e064bd1bf45747aa00730e99b..0c9ccacc18118d314c3880168530fea1f44f48fd 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/volumevariables.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/volumevariables.hh @@ -83,25 +83,25 @@ public: const SubControlVolume& scv) { RANSParentType::updateRANSProperties(elemSol, problem, element, scv); - isMatchingPoint_ = problem.isMatchingPoint(RANSParentType::elementID()); - inNearWallRegion_ = problem.inNearWallRegion(RANSParentType::elementID()); + isMatchingPoint_ = problem.isMatchingPoint(RANSParentType::elementIdx()); + inNearWallRegion_ = problem.inNearWallRegion(RANSParentType::elementIdx()); turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx]; dissipation_ = elemSol[0][Indices::dissipationIdx]; - storedDissipation_ = problem.storedDissipation_[RANSParentType::elementID()]; - storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()]; - stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()]; - const Scalar uStarNominal = problem.uStarNominal(RANSParentType::elementID()); - const auto flowNormalAxis = problem.flowNormalAxis_[RANSParentType::elementID()]; - yPlusNominal_ = RANSParentType::wallDistance() * uStarNominal / problem.kinematicViscosity_[RANSParentType::elementID()]; + storedDissipation_ = problem.storedDissipation_[RANSParentType::elementIdx()]; + storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementIdx()]; + stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()]; + const Scalar uStarNominal = problem.uStarNominal(RANSParentType::elementIdx()); + const auto flowNormalAxis = problem.flowNormalAxis_[RANSParentType::elementIdx()]; + yPlusNominal_ = RANSParentType::wallDistance() * uStarNominal / problem.kinematicViscosity_[RANSParentType::elementIdx()]; uPlusNominal_ = RANSParentType::velocity()[flowNormalAxis] / uStarNominal; cMu_ = problem.cMu(); if (problem.useStoredEddyViscosity_) - RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementID()]); + RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]); else RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity()); if (inNearWallRegion_) { - RANSParentType::setDynamicEddyViscosity_(problem.zeroEqDynamicEddyViscosity_[RANSParentType::elementID()]); + RANSParentType::setDynamicEddyViscosity_(problem.zeroEqDynamicEddyViscosity_[RANSParentType::elementIdx()]); } RANSParentType::calculateEddyDiffusivity(problem); RANSParentType::calculateEddyThermalConductivity(problem); diff --git a/dumux/freeflow/rans/twoeq/komega/problem.hh b/dumux/freeflow/rans/twoeq/komega/problem.hh index 81a3a5a9065f1e9303cb71ed44d98c09d67f6209..05c6d95d000a56035c7ed5ade3ab5379e68ff7a0 100644 --- a/dumux/freeflow/rans/twoeq/komega/problem.hh +++ b/dumux/freeflow/rans/twoeq/komega/problem.hh @@ -101,7 +101,7 @@ public: for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); @@ -112,30 +112,30 @@ public: PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars); auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(priVars)); // NOTE: first update the turbulence quantities - storedDissipation_[elementID] = elemSol[0][Indices::dissipationEqIdx]; - storedTurbulentKineticEnergy_[elementID] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; + storedDissipation_[elementIdx] = elemSol[0][Indices::dissipationEqIdx]; + storedTurbulentKineticEnergy_[elementIdx] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; // NOTE: then update the volVars VolumeVariables volVars; volVars.update(elemSol, asImp_(), element, scv); - storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity(*this); + storedDynamicEddyViscosity_[elementIdx] = volVars.calculateEddyViscosity(*this); } } // calculate cell-centered gradients for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) { - unsigned backwardNeighbor = ParentType::neighborID_[elementID][dimIdx][0]; - unsigned forwardNeighbor = ParentType::neighborID_[elementID][dimIdx][1]; - storedTurbulentKineticEnergyGradient_[elementID][dimIdx] + unsigned backwardNeighbor = ParentType::neighborIdx_[elementIdx][dimIdx][0]; + unsigned forwardNeighbor = ParentType::neighborIdx_[elementIdx][dimIdx][1]; + storedTurbulentKineticEnergyGradient_[elementIdx][dimIdx] = (storedTurbulentKineticEnergy_[forwardNeighbor] - storedTurbulentKineticEnergy_[backwardNeighbor]) / (ParentType::cellCenter_[forwardNeighbor][dimIdx] - ParentType::cellCenter_[backwardNeighbor][dimIdx]); - storedDissipationGradient_[elementID][dimIdx] + storedDissipationGradient_[elementIdx][dimIdx] = (storedDissipation_[forwardNeighbor] - storedDissipation_[backwardNeighbor]) / (ParentType::cellCenter_[forwardNeighbor][dimIdx] diff --git a/dumux/freeflow/rans/twoeq/komega/volumevariables.hh b/dumux/freeflow/rans/twoeq/komega/volumevariables.hh index 15eed5071e2a10e8cf5e321c19e4d6f3f8496b94..fc83edec1dc9417333658b020b1398394cb44191 100644 --- a/dumux/freeflow/rans/twoeq/komega/volumevariables.hh +++ b/dumux/freeflow/rans/twoeq/komega/volumevariables.hh @@ -87,13 +87,13 @@ public: betaOmega_ = problem.betaOmega(); turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx]; dissipation_ = elemSol[0][Indices::dissipationIdx]; - storedDissipation_ = problem.storedDissipation_[RANSParentType::elementID()]; - storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()]; - storedDissipationGradient_ = problem.storedDissipationGradient_[RANSParentType::elementID()]; - storedTurbulentKineticEnergyGradient_ = problem.storedTurbulentKineticEnergyGradient_[RANSParentType::elementID()]; - stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()]; + storedDissipation_ = problem.storedDissipation_[RANSParentType::elementIdx()]; + storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementIdx()]; + storedDissipationGradient_ = problem.storedDissipationGradient_[RANSParentType::elementIdx()]; + storedTurbulentKineticEnergyGradient_ = problem.storedTurbulentKineticEnergyGradient_[RANSParentType::elementIdx()]; + stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()]; if (problem.useStoredEddyViscosity_) - RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementID()]); + RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]); else RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity(problem)); RANSParentType::calculateEddyDiffusivity(problem); diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh index dedbcdfa42a782816ac3f148d8b96c709b3c9190..742b6e790ccd69081d348555c92fd80ba2f2e9b9 100644 --- a/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh +++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/problem.hh @@ -92,7 +92,7 @@ public: for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); @@ -103,12 +103,12 @@ public: PrimaryVariables priVars = makePriVarsFromCellCenterPriVars<PrimaryVariables>(cellCenterPriVars); auto elemSol = elementSolution<typename FVGridGeometry::LocalView>(std::move(priVars)); // NOTE: first update the turbulence quantities - storedDissipationTilde_[elementID] = elemSol[0][Indices::dissipationEqIdx]; - storedTurbulentKineticEnergy_[elementID] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; + storedDissipationTilde_[elementIdx] = elemSol[0][Indices::dissipationEqIdx]; + storedTurbulentKineticEnergy_[elementIdx] = elemSol[0][Indices::turbulentKineticEnergyEqIdx]; // NOTE: then update the volVars VolumeVariables volVars; volVars.update(elemSol, asImp_(), element, scv); - storedDynamicEddyViscosity_[elementID] = volVars.calculateEddyViscosity(); + storedDynamicEddyViscosity_[elementIdx] = volVars.calculateEddyViscosity(); } } } diff --git a/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh b/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh index a79bc6d4051cbf04ac33a8376e1d96a306e3fde7..c691a231025b698cc55d0e7d93f8cf4aee058dc5 100644 --- a/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh +++ b/dumux/freeflow/rans/twoeq/lowrekepsilon/volumevariables.hh @@ -85,11 +85,11 @@ public: RANSParentType::updateRANSProperties(elemSol, problem, element, scv); turbulentKineticEnergy_ = elemSol[0][Indices::turbulentKineticEnergyIdx]; dissipationTilde_ = elemSol[0][Indices::dissipationIdx]; - storedDissipationTilde_ = problem.storedDissipationTilde_[RANSParentType::elementID()]; - storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementID()]; - stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementID()]; + storedDissipationTilde_ = problem.storedDissipationTilde_[RANSParentType::elementIdx()]; + storedTurbulentKineticEnergy_ = problem.storedTurbulentKineticEnergy_[RANSParentType::elementIdx()]; + stressTensorScalarProduct_ = problem.stressTensorScalarProduct_[RANSParentType::elementIdx()]; if (problem.useStoredEddyViscosity_) - RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementID()]); + RANSParentType::setDynamicEddyViscosity_(problem.storedDynamicEddyViscosity_[RANSParentType::elementIdx()]); else RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity()); RANSParentType::calculateEddyDiffusivity(problem); diff --git a/dumux/freeflow/rans/volumevariables.hh b/dumux/freeflow/rans/volumevariables.hh index 1461d93d73a4ece4399b6141b927a3edd61f8d7d..77c1c86b7698729ea8a23610bd6e3204ad3a5ff0 100644 --- a/dumux/freeflow/rans/volumevariables.hh +++ b/dumux/freeflow/rans/volumevariables.hh @@ -92,27 +92,27 @@ public: using std::sqrt; // calculate characteristic properties of the turbulent flow - elementID_ = problem.fvGridGeometry().elementMapper().index(element); - wallElementID_ = problem.wallElementID_[elementID_]; - wallDistance_ = problem.wallDistance_[elementID_]; - velocity_ = problem.velocity_[elementID_]; - velocityMaximum_ = problem.velocityMaximum_[wallElementID_]; - velocityGradients_ = problem.velocityGradients_[elementID_]; - const auto flowNormalAxis = problem.flowNormalAxis_[elementID_]; - const auto wallNormalAxis = problem.wallNormalAxis_[elementID_]; - uStar_ = sqrt(problem.kinematicViscosity_[wallElementID_] - * abs(problem.velocityGradients_[wallElementID_][flowNormalAxis][wallNormalAxis])); + elementIdx_ = problem.fvGridGeometry().elementMapper().index(element); + wallElementIdx_ = problem.wallElementIdx_[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])); uStar_ = max(uStar_, 1e-10); // zero values lead to numerical problems in some turbulence models - yPlus_ = wallDistance_ * uStar_ / problem.kinematicViscosity_[elementID_]; + yPlus_ = wallDistance_ * uStar_ / problem.kinematicViscosity_[elementIdx_]; uPlus_ = velocity_[flowNormalAxis] / uStar_; karmanConstant_ = problem.karmanConstant(); } /*! - * \brief Return the element ID of the control volume. + * \brief Return the element Idx of the control volume. */ - unsigned int elementID() const - { return elementID_; } + unsigned int elementIdx() const + { return elementIdx_; } /*! * \brief Return the velocity vector \f$\mathrm{[m/s]}\f$ at the control volume center. @@ -262,8 +262,8 @@ protected: DimVector velocity_; DimVector velocityMaximum_; DimMatrix velocityGradients_; - std::size_t elementID_; - std::size_t wallElementID_; + 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 4db036eb3310a333030ab7c9e846d5561f0ce961..8d8eedd023a31546d5915d5eb2ae4c23e6f367b6 100644 --- a/dumux/freeflow/rans/zeroeq/problem.hh +++ b/dumux/freeflow/rans/zeroeq/problem.hh @@ -108,7 +108,7 @@ public: bool printedRangeWarning = false; for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); auto fvGeometry = localView(this->fvGridGeometry()); fvGeometry.bindElement(element); @@ -127,26 +127,26 @@ public: VolumeVariables volVars; volVars.update(elemSol, asImp_(), element, scv); - Scalar ksPlus = this->sandGrainRoughness_[elementID] * volVars.uStar() / volVars.kinematicViscosity(); + Scalar ksPlus = this->sandGrainRoughness_[elementIdx] * volVars.uStar() / volVars.kinematicViscosity(); if (ksPlus > 0 && eddyViscosityModel_.compare("baldwinLomax") == 0) { DUNE_THROW(Dune::NotImplemented, "Roughness is not implemented for the Baldwin-Lomax model."); } if (ksPlus > 2000.) { - std::cout << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << this->cellCenter_[this->wallElementID_[elementID]] + std::cout << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << this->cellCenter_[this->wallElementIdx_[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->wallElementID_[elementID]] + Dune::dinfo << "info: equivalent sand grain roughness ks+=" << ksPlus << " at " << this->cellCenter_[this->wallElementIdx_[elementIdx]] << " is not in the valid range (ksPlus > 4.535) and now set to 0.0"<< std::endl; ksPlus = 0.0; printedRangeWarning = true; } - additionalRoughnessLength_[elementID] = 0.9 / (volVars.uStar() / volVars.kinematicViscosity()) + additionalRoughnessLength_[elementIdx] = 0.9 / (volVars.uStar() / volVars.kinematicViscosity()) * (sqrt(ksPlus) - ksPlus * exp(-ksPlus / 6.0)); } } @@ -185,82 +185,82 @@ public: // (1) calculate inner viscosity and Klebanoff function for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int wallElementID = this->wallElementID_[elementID]; - Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID]; - unsigned int flowNormalAxis = this->flowNormalAxis_[elementID]; - unsigned int wallNormalAxis = this->wallNormalAxis_[elementID]; - - Scalar omegaAbs = abs(this->velocityGradients_[elementID][flowNormalAxis][wallNormalAxis] - - this->velocityGradients_[elementID][wallNormalAxis][flowNormalAxis]); - Scalar uStar = sqrt(this->kinematicViscosity_[wallElementID] - * abs(this->velocityGradients_[wallElementID][flowNormalAxis][wallNormalAxis])); - Scalar yPlus = wallDistance * uStar / this->kinematicViscosity_[elementID]; + unsigned int elementIdx = this->fvGridGeometry().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 omegaAbs = abs(this->velocityGradients_[elementIdx][flowNormalAxis][wallNormalAxis] + - this->velocityGradients_[elementIdx][wallNormalAxis][flowNormalAxis]); + Scalar uStar = sqrt(this->kinematicViscosity_[wallElementIdx] + * abs(this->velocityGradients_[wallElementIdx][flowNormalAxis][wallNormalAxis])); + Scalar yPlus = wallDistance * uStar / this->kinematicViscosity_[elementIdx]; Scalar mixingLength = this->karmanConstant() * wallDistance * (1.0 - exp(-yPlus / aPlus)); - kinematicEddyViscosityInner[elementID] = mixingLength * mixingLength * omegaAbs; + kinematicEddyViscosityInner[elementIdx] = mixingLength * mixingLength * omegaAbs; Scalar f = wallDistance * omegaAbs * (1.0 - exp(-yPlus / aPlus)); - if (f > storedFMax[wallElementID]) + if (f > storedFMax[wallElementIdx]) { - storedFMax[wallElementID] = f; - storedYFMax[wallElementID] = wallDistance; + storedFMax[wallElementIdx] = f; + storedYFMax[wallElementIdx] = wallDistance; } } // (2) calculate outer viscosity for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int wallElementID = this->wallElementID_[elementID]; - Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID]; + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + unsigned int wallElementIdx = this->wallElementIdx_[elementIdx]; + Scalar wallDistance = this->wallDistance_[elementIdx] + additionalRoughnessLength_[elementIdx]; Scalar maxVelocityNorm = 0.0; Scalar minVelocityNorm = 0.0; for (unsigned dimIdx = 0; dimIdx < dim; ++dimIdx) { - maxVelocityNorm += this->velocityMaximum_[wallElementID][dimIdx] - * this->velocityMaximum_[wallElementID][dimIdx]; - minVelocityNorm += this->velocityMinimum_[wallElementID][dimIdx] - * this->velocityMinimum_[wallElementID][dimIdx]; + maxVelocityNorm += this->velocityMaximum_[wallElementIdx][dimIdx] + * this->velocityMaximum_[wallElementIdx][dimIdx]; + minVelocityNorm += this->velocityMinimum_[wallElementIdx][dimIdx] + * this->velocityMinimum_[wallElementIdx][dimIdx]; } Scalar deltaU = sqrt(maxVelocityNorm) - sqrt(minVelocityNorm); - Scalar yFMax = storedYFMax[wallElementID]; - Scalar fMax = storedFMax[wallElementID]; + Scalar yFMax = storedYFMax[wallElementIdx]; + Scalar fMax = storedFMax[wallElementIdx]; Scalar fWake = min(yFMax * fMax, cWake * yFMax * deltaU * deltaU / fMax); Scalar fKleb = 1.0 / (1.0 + 5.5 * pow(cKleb * wallDistance / yFMax, 6.0)); - kinematicEddyViscosityOuter[elementID] = k * cCP * fWake * fKleb; + kinematicEddyViscosityOuter[elementIdx] = k * cCP * fWake * fKleb; - kinematicEddyViscosityDifference[elementID] - = kinematicEddyViscosityInner[elementID] - kinematicEddyViscosityOuter[elementID]; + kinematicEddyViscosityDifference[elementIdx] + = kinematicEddyViscosityInner[elementIdx] - kinematicEddyViscosityOuter[elementIdx]; } // (3) switching point for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int wallElementID = this->wallElementID_[elementID]; - Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID]; + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + unsigned int wallElementIdx = this->wallElementIdx_[elementIdx]; + Scalar wallDistance = this->wallDistance_[elementIdx] + additionalRoughnessLength_[elementIdx]; // checks if sign switches, by multiplication - Scalar check = kinematicEddyViscosityDifference[wallElementID] * kinematicEddyViscosityDifference[elementID]; + Scalar check = kinematicEddyViscosityDifference[wallElementIdx] * kinematicEddyViscosityDifference[elementIdx]; if (check < 0 // means sign has switched - && switchingPosition[wallElementID] > wallDistance) + && switchingPosition[wallElementIdx] > wallDistance) { - switchingPosition[wallElementID] = wallDistance; + switchingPosition[wallElementIdx] = wallDistance; } } // (4) finally determine eddy viscosity for (const auto& element : elements(this->fvGridGeometry().gridView())) { - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - unsigned int wallElementID = this->wallElementID_[elementID]; - Scalar wallDistance = this->wallDistance_[elementID] + additionalRoughnessLength_[elementID]; + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + unsigned int wallElementIdx = this->wallElementIdx_[elementIdx]; + Scalar wallDistance = this->wallDistance_[elementIdx] + additionalRoughnessLength_[elementIdx]; - kinematicEddyViscosity_[elementID] = kinematicEddyViscosityInner[elementID]; - if (wallDistance >= switchingPosition[wallElementID]) + kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityInner[elementIdx]; + if (wallDistance >= switchingPosition[wallElementIdx]) { - kinematicEddyViscosity_[elementID] = kinematicEddyViscosityOuter[elementID]; + kinematicEddyViscosity_[elementIdx] = kinematicEddyViscosityOuter[elementIdx]; } } } diff --git a/dumux/freeflow/rans/zeroeq/volumevariables.hh b/dumux/freeflow/rans/zeroeq/volumevariables.hh index b2d0e2764b2fb2516d5d20d54d2aa7c70b3bc20e..2b4cdced6bc7f8440f4909da73a67d65840ab199 100644 --- a/dumux/freeflow/rans/zeroeq/volumevariables.hh +++ b/dumux/freeflow/rans/zeroeq/volumevariables.hh @@ -83,7 +83,7 @@ public: const SubControlVolume& scv) { RANSParentType::updateRANSProperties(elemSol, problem, element, scv); - additionalRoughnessLength_ = problem.additionalRoughnessLength_[RANSParentType::elementID()]; + additionalRoughnessLength_ = problem.additionalRoughnessLength_[RANSParentType::elementIdx()]; yPlusRough_ = wallDistanceRough() * RANSParentType::uStar() / RANSParentType::kinematicViscosity(); RANSParentType::setDynamicEddyViscosity_(calculateEddyViscosity(elemSol, problem, element, scv, problem.eddyViscosityModel_)); RANSParentType::calculateEddyDiffusivity(problem); @@ -110,8 +110,8 @@ public: using std::exp; using std::sqrt; Scalar kinematicEddyViscosity = 0.0; - unsigned int flowNormalAxis = problem.flowNormalAxis_[RANSParentType::elementID()]; - unsigned int wallNormalAxis = problem.wallNormalAxis_[RANSParentType::elementID()]; + 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) @@ -132,7 +132,7 @@ public: } else if (modelName.compare("baldwinLomax") == 0) { - kinematicEddyViscosity = problem.kinematicEddyViscosity_[RANSParentType::elementID()]; + kinematicEddyViscosity = problem.kinematicEddyViscosity_[RANSParentType::elementIdx()]; } else { diff --git a/test/freeflow/rans/pipelauferproblem.hh b/test/freeflow/rans/pipelauferproblem.hh index 915d94fd5ed90c7b2c744c3661d1d68ac01ebffa..b3ea81f4235e44cb1738d20006b53ec595d4e648 100644 --- a/test/freeflow/rans/pipelauferproblem.hh +++ b/test/freeflow/rans/pipelauferproblem.hh @@ -308,9 +308,9 @@ public: PrimaryVariables values(initialAtPos(globalPos)); #if KOMEGA using std::pow; - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - const auto wallDistance = ParentType::wallDistance_[elementID]; - values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementID] + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + const auto wallDistance = ParentType::wallDistance_[elementIdx]; + values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx] / (ParentType::betaOmega() * pow(wallDistance, 2)); #endif return values; diff --git a/test/freeflow/ransnc/flatplatetestproblem.hh b/test/freeflow/ransnc/flatplatetestproblem.hh index 723868298f52de7cab857c6f6a591eaf9286d0cb..39bffa6f7c5839adfc08a512588adbeb6c15d605 100644 --- a/test/freeflow/ransnc/flatplatetestproblem.hh +++ b/test/freeflow/ransnc/flatplatetestproblem.hh @@ -339,9 +339,9 @@ public: PrimaryVariables values(initialAtPos(globalPos)); #if KOMEGA using std::pow; - unsigned int elementID = this->fvGridGeometry().elementMapper().index(element); - const auto wallDistance = ParentType::wallDistance_[elementID]; - values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementID] + unsigned int elementIdx = this->fvGridGeometry().elementMapper().index(element); + const auto wallDistance = ParentType::wallDistance_[elementIdx]; + values[Indices::dissipationEqIdx] = 6.0 * ParentType::kinematicViscosity_[elementIdx] / (ParentType::betaOmega() * pow(wallDistance, 2)); #endif return values;