diff --git a/dumux/common/deprecated.hh b/dumux/common/deprecated.hh index 76e3b726529ebc1e4ae51eef511bb36d952462e2..905938784b50dd0794231504a6d5c6a78d36887c 100644 --- a/dumux/common/deprecated.hh +++ b/dumux/common/deprecated.hh @@ -54,6 +54,20 @@ template void update(Mapper& mapper) { mapper.update(); }; +template +using HasIsOnWallDetector = decltype(std::declval().isOnWallAtPos(std::declval())); + +template +static constexpr bool hasIsOnWall() +{ return Dune::Std::is_detected::value; } + +template +using HasWallBCDetector = decltype(std::declval().hasWall()); + +template +static constexpr bool hasHasWallBC() +{ return Dune::Std::is_detected::value; } + } // end namespace Deprecated #endif diff --git a/dumux/discretization/walldistance.hh b/dumux/discretization/walldistance.hh index 5fd9fdd523beb5f8775af5a257677faafc81aca9..8ed975881f132fa39874a8d6fb8e383d03afd7c0 100644 --- a/dumux/discretization/walldistance.hh +++ b/dumux/discretization/walldistance.hh @@ -51,6 +51,8 @@ class WallDistance using GridView = typename GridGeometry::GridView; using GridIndexType = typename IndexTraits::GridIndex; using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using FVElementGeometry = typename GridGeometry::LocalView; + using Element = typename GridGeometry::GridView::template Codim<0>::Entity; using Scalar = typename GridView::Grid::ctype; using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; @@ -131,7 +133,7 @@ public: */ template WallDistance(std::shared_ptr gridGeometry, LocationTag tag) - : WallDistance(gridGeometry, tag, [](const SubControlVolumeFace& scvf) { return true; }) {} + : WallDistance(gridGeometry, tag, [](const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) { return true; }) {} //! caller has to make sure the lifetime of grid geometry exceeds the lifetime of wall distance template @@ -194,7 +196,7 @@ private: for (const auto& scvf : scvfs(fvGeometry)) { - if (scvf.boundary() && considerFace(scvf)) + if (scvf.boundary() && considerFace(fvGeometry, scvf)) { const auto& geo = scvf.geometry(); CornerStorage corners; diff --git a/dumux/freeflow/rans/boundarytypes.hh b/dumux/freeflow/rans/boundarytypes.hh new file mode 100644 index 0000000000000000000000000000000000000000..5d715030600ba68ef980c4a464e9c0909bd6d39f --- /dev/null +++ b/dumux/freeflow/rans/boundarytypes.hh @@ -0,0 +1,106 @@ +// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- +// vi: set et ts=4 sw=4 sts=4: +/***************************************************************************** + * See the file COPYING for full copying permissions. * + * * + * This program is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * This program is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have received a copy of the GNU General Public License * + * along with this program. If not, see . * + *****************************************************************************/ +/*! + * \file + * \ingroup RANSModel + * \copydoc Dumux::RANSBoundaryTypes + */ +#ifndef FREEFLOW_RANS_BOUNDARY_TYPES_HH +#define FREEFLOW_RANS_BOUNDARY_TYPES_HH + +#include +#include + +namespace Dumux { + +/*! + * \ingroup RANSModel + * \brief Class to specify the type of a boundary condition for the RANS extension to the Navier-Stokes model. + */ +template +class RANSBoundaryTypes : public NavierStokesBoundaryTypes +{ + using ParentType = NavierStokesBoundaryTypes; + static constexpr auto dimWorld = ModelTraits::dim(); + using Indices = typename ModelTraits::Indices; + static_assert(dimWorld > 1, "Wall conditions cannot be set for 1D domains."); +public: + RANSBoundaryTypes() + { + for (int eqIdx=0; eqIdx < numEq; ++eqIdx) + resetEq(eqIdx); + } + + void setWall() + { + isWall_ = true; + for (int eqIdx=0; eqIdx < numEq; ++eqIdx) + { + if constexpr (dimWorld == 3) + { + if ((eqIdx == Indices::velocityXIdx) + || (eqIdx == Indices::velocityYIdx) + || (eqIdx == Indices::velocityZIdx)) + BoundaryTypes::boundaryInfo_[eqIdx].isDirichlet = true; + } + else if constexpr (dimWorld == 2) + { + if ((eqIdx == Indices::velocityXIdx) + || (eqIdx == Indices::velocityYIdx)) + BoundaryTypes::boundaryInfo_[eqIdx].isDirichlet = true; + } + else + DUNE_THROW(Dune::NotImplemented, "1D Turbulence models are not supported"); + + if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 1) + { + if (eqIdx == Indices::viscosityTildeIdx) + BoundaryTypes::boundaryInfo_[eqIdx].isDirichlet = true; + } + else if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 2) + { + if (eqIdx == Indices::turbulentKineticEnergyEqIdx || eqIdx == Indices::dissipationEqIdx) + BoundaryTypes::boundaryInfo_[eqIdx].isDirichlet = true; + } + } + } + + /*! + * \brief Returns true if some equation is used to specify a + * Wall boundary condition. + */ + bool hasWall() const + { return isWall_; } + + /*! + * \brief Reset the boundary types for one equation. + */ + void resetEq(const int eqIdx) + { + ParentType::resetEq(eqIdx); + isWall_ = false; + } + +protected: + bool isWall_; +}; + +} // end namespace Dumux + +#endif diff --git a/dumux/freeflow/rans/problem.hh b/dumux/freeflow/rans/problem.hh index 884bda4e5d7233de04ae3fbcd2b1d1be99700f10..9e08345c41b2e4283ebb33773a741ddb44230a49 100644 --- a/dumux/freeflow/rans/problem.hh +++ b/dumux/freeflow/rans/problem.hh @@ -28,12 +28,13 @@ #include #include +#include #include #include -#include #include +#include +#include #include - #include "model.hh" namespace Dumux { @@ -94,6 +95,7 @@ class RANSProblemBase : public NavierStokesProblem }; public: + /*! * \brief The constructor * \param gridGeometry The finite volume grid geometry @@ -101,145 +103,41 @@ public: */ RANSProblemBase(std::shared_ptr gridGeometry, const std::string& paramGroup = "") : ParentType(gridGeometry, paramGroup) - { } - - /*! - * \brief Update the static (solution independent) relations to the walls - * - * This function determines all element with a wall intersection, - * the wall distances and the relation to the neighboring elements. - */ - void updateStaticWallProperties() { - using std::abs; - std::cout << "Update static wall properties. "; - calledUpdateStaticWallProperties = true; + if ( !(hasParamInGroup(this->paramGroup(), "RANS.IsFlatWallBounded"))) + { + std::cout << "The parameter \"Rans.IsFlatWallBounded\" is not specified. \n" + << " -- Based on the grid and the boundary conditions specified by the user," + << " this parameter is set to be "<< std::boolalpha << isFlatWallBounded() << "\n"; + } // update size and initial values of the global vectors - wallElementIdx_.resize(this->gridGeometry().elementMapper().size()); wallDistance_.resize(this->gridGeometry().elementMapper().size(), std::numeric_limits::max()); neighborIdx_.resize(this->gridGeometry().elementMapper().size()); - cellCenter_.resize(this->gridGeometry().elementMapper().size(), GlobalPosition(0.0)); velocity_.resize(this->gridGeometry().elementMapper().size(), DimVector(0.0)); 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); flowDirectionAxis_.resize(this->gridGeometry().elementMapper().size(), fixedFlowDirectionAxis_); - wallNormalAxis_.resize(this->gridGeometry().elementMapper().size(), fixedWallNormalAxis_); storedViscosity_.resize(this->gridGeometry().elementMapper().size(), 0.0); storedDensity_.resize(this->gridGeometry().elementMapper().size(), 0.0); + } - if ( !(hasParamInGroup(this->paramGroup(), "RANS.IsFlatWallBounded"))) - { - std::cout << "The parameter \"Rans.IsFlatWallBounded\" is not specified. \n" - << " -- Based on the grid and the isOnWallAtPos function specified by the user," - << " this parameter is set to be "<< std::boolalpha << isFlatWallBounded() << "\n"; - } - - std::vector wallElements; - - const auto gridView = this->gridGeometry().gridView(); - - auto fvGeometry = localView(this->gridGeometry()); - for (const auto& element : elements(gridView)) - { - fvGeometry.bindElement(element); - for (const auto& scvf : scvfs(fvGeometry)) - { - // only search for walls at a global boundary - if (!scvf.boundary()) - continue; - - if (asImp_().isOnWall(scvf)) - { - WallElementInformation wallElementInformation; - - // store the location of the wall adjacent face's center and all corners - wallElementInformation.wallFaceCenter = scvf.center(); - for (int i = 0; i < numCorners; i++) - wallElementInformation.wallFaceCorners[i] = scvf.corner(i); - - // Store the wall element index and face's normal direction (used only with isFlatWallBounded on) - wallElementInformation.wallElementIdx = this->gridGeometry().elementMapper().index(element); - wallElementInformation.wallFaceNormalAxis = scvf.directionIndex(); - - wallElements.push_back(wallElementInformation); - } - } - } - // output the number of wall adjacent faces. Check that this is non-zero. - std::cout << "NumWallIntersections=" << wallElements.size() << std::endl; - if (wallElements.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 the wall for each element - for (const auto& element : elements(gridView)) - { - // Store the cell center position for each element - unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); - cellCenter_[elementIdx] = element.geometry().center(); - - for (unsigned int i = 0; i < wallElements.size(); ++i) - { - // Find the minimum distance from the cell center to the wall face (center and corners) - std::array cellToWallDistances; - for (unsigned int j = 0; j < numCorners; j++) - cellToWallDistances[j] = (cellCenter(elementIdx) - wallElements[i].wallFaceCorners[j]).two_norm(); - cellToWallDistances[numCorners] = (cellCenter(elementIdx) - wallElements[i].wallFaceCenter).two_norm(); - Scalar distanceToWall = *std::min_element(cellToWallDistances.begin(), cellToWallDistances.end()); - - if (distanceToWall < wallDistance_[elementIdx]) - { - wallDistance_[elementIdx] = distanceToWall; - // If isFlatWallBounded, the corresonding wall element is stored for each element - if (isFlatWallBounded()) - { - wallElementIdx_[elementIdx] = wallElements[i].wallElementIdx; - if ( !(hasParam("RANS.WallNormalAxis")) ) - wallNormalAxis_[elementIdx] = wallElements[i].wallFaceNormalAxis; - } - } - } - } - - // search for neighbor Idxs - for (const auto& element : elements(gridView)) - { - unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); - for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) - { - neighborIdx_[elementIdx][dimIdx][0] = elementIdx; - neighborIdx_[elementIdx][dimIdx][1] = elementIdx; - } - - for (const auto& intersection : intersections(gridView, element)) - { - if (intersection.boundary()) - continue; - - unsigned int neighborIdx = this->gridGeometry().elementMapper().index(intersection.outside()); - for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) - { - if (abs(cellCenter(elementIdx)[dimIdx] - cellCenter(neighborIdx)[dimIdx]) > 1e-8) - { - if (cellCenter(elementIdx)[dimIdx] > cellCenter(neighborIdx)[dimIdx]) - neighborIdx_[elementIdx][dimIdx][0] = neighborIdx; + /*! + * \brief Update the static (solution independent) relations to the walls and neighbors + */ + void updateStaticWallProperties() + { + std::cout << "Update static wall properties. "; + calledUpdateStaticWallProperties = true; - if (cellCenter(elementIdx)[dimIdx] < cellCenter(neighborIdx)[dimIdx]) - neighborIdx_[elementIdx][dimIdx][1] = neighborIdx; - } - } - } - } + checkForWalls_(); + findWallDistances_(); + findNeighborIndices_(); } /*! - * \brief Update the dynamic (solution dependent) relations to the walls - * - * The basic function calcuates the cell-centered velocities and - * the respective gradients. - * Further, the kinematic viscosity at the wall is stored. + * \brief Update the dynamic (solution dependent) turbulence parameters * * \param curSol The solution vector. */ @@ -298,23 +196,13 @@ public: * * \param scvf The sub control volume face. */ + [[deprecated("The isOnWall and IsOnWallAtPos functions will be removed after release 3.5. " + "Please use the Rans specific boundarytypes, and mark wall boundaries with the setWall() function.")]] bool isOnWall(const SubControlVolumeFace& scvf) const { return asImp_().isOnWallAtPos(scvf.center()); } - /*! - * \brief Returns whether a given point is on a wall - * - * \param globalPos The position in global coordinates. - */ - bool isOnWallAtPos(const GlobalPosition &globalPos) const - { - // Throw an exception if no walls are implemented - DUNE_THROW(Dune::InvalidStateException, - "The problem does not provide an isOnWall() method."); - } - bool isFlatWallBounded() const { static const bool hasAlignedWalls = hasAlignedWalls_(); @@ -358,24 +246,30 @@ public: int wallNormalAxis(const int elementIdx) const { if (!isFlatWallBounded()) - DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, models requiring a wallNormalAxis can only be used for flat wall bounded flows. " - << "\n If your geometry is a flat channel, please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); + DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, models requiring a wallNormalAxis " + << "can only be used for flat wall bounded flows. " + << "\n If your geometry is a flat channel, " + << "please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); return wallNormalAxis_[elementIdx]; } int flowDirectionAxis(const int elementIdx) const { if (!isFlatWallBounded()) - DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, models requiring a flowDirectionAxis can only be used for flat wall bounded flows. " - << "\n If your geometry is a flat channel, please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); + DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, models requiring a flowDirectionAxis " + << "can only be used for flat wall bounded flows. " + << "\n If your geometry is a flat channel, " + << "please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); return flowDirectionAxis_[elementIdx]; } unsigned int wallElementIndex(const int elementIdx) const { if (!isFlatWallBounded()) - DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, models requiring a wallElementIndex can only be used for flat wall bounded flows. " - << "\n If your geometry is a flat channel, please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); + DUNE_THROW(Dune::NotImplemented, "\n Due to grid/geometric concerns, models requiring a wallElementIndex " + << "can only be used for flat wall bounded flows. " + << "\n If your geometry is a flat channel, " + << "please set the runtime parameter RANS.IsFlatWallBounded to true. \n"); return wallElementIdx_[elementIdx]; } @@ -384,7 +278,10 @@ public: { return wallDistance_[elementIdx]; } GlobalPosition cellCenter(const int elementIdx) const - { return cellCenter_[elementIdx]; } + { + const auto& element = this->gridGeometry().element(elementIdx); + return element.geometry().center(); + } unsigned int neighborIndex(const int elementIdx, const int dimIdx, const int sideIdx) const { return neighborIdx_[elementIdx][dimIdx][sideIdx];} @@ -444,9 +341,19 @@ private: fvGeometry.bindElement(element); for (const auto& scvf : scvfs(fvGeometry)) { - // only search for walls at a global boundary - if (!scvf.boundary() && asImp_().isOnWall(scvf)) - wallFaceAxis.push_back(scvf.directionIndex()); + // Remove this check after release 3.5. IsOnWall Interface is deprecated + if constexpr (Deprecated::hasIsOnWall()) + { + // Remove this part + if (!scvf.boundary() && asImp_().isOnWall(scvf)) // only search for walls at a global boundary + wallFaceAxis.push_back(scvf.directionIndex()); + } + else + { + // Keep this part + if (!scvf.boundary() && asImp_().boundaryTypes(element, scvf).hasWall()) // only search for walls at a global boundary + wallFaceAxis.push_back(scvf.directionIndex()); + } } } @@ -454,6 +361,122 @@ private: return std::all_of(wallFaceAxis.begin(), wallFaceAxis.end(), [firstDir=wallFaceAxis[0]](auto dir){ return (dir == firstDir);} ) ; } + void checkForWalls_() + { + for (const auto& element : elements(this->gridGeometry().gridView())) + { + auto fvGeometry = localView(this->gridGeometry()); + fvGeometry.bindElement(element); + for (auto&& scvf : scvfs(fvGeometry)) + { + if constexpr (Deprecated::hasHasWallBC()) + { + if (asImp_().boundaryTypes(element, scvf).hasWall()) + return; + } + else + noSetWallCompilerWarning_(); + } + } + + // If reached, no walls were found, throw exception. Remove check after 3.5 + if constexpr (!Deprecated::hasIsOnWall()) + DUNE_THROW(Dune::InvalidStateException, "No walls are are specified with the setWall() function"); + } + + /*! + * \brief Use the boundary search algorithm to find the shortest distance to a wall for each element + * + * Also store the wall element's index, and its direction in the case of flat wall bounded problems + */ + void findWallDistances_() + { + // Remove this check after release 3.5. IsOnWall Interface is deprecated + if constexpr (Deprecated::hasIsOnWall()) + { + WallDistance wallInformation(this->gridGeometry(), WallDistance::atElementCenters, + [this] (const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) + { return asImp_().isOnWall(scvf); }); + wallDistance_ = wallInformation.wallDistance(); + storeWallElementAndDirectionIndex_(wallInformation.wallData()); + } + else + { + WallDistance wallInformation(this->gridGeometry(), WallDistance::atElementCenters, + [this] (const FVElementGeometry& fvGeometry, const SubControlVolumeFace& scvf) + { return asImp_().boundaryTypes(fvGeometry.element(), scvf).hasWall(); }); + wallDistance_ = wallInformation.wallDistance(); + storeWallElementAndDirectionIndex_(wallInformation.wallData()); + } + } + + template + void storeWallElementAndDirectionIndex_(const WallData& wallData) + { + // The wall Direction Index is used for flat quadrilateral channel problems only + if (!(GridGeometry::discMethod == DiscretizationMethod::staggered)) + DUNE_THROW(Dune::NotImplemented, "The wall direction Index can only be calculated for quadrilateral structured grids"); + + // If isFlatWallBounded, the corresonding wall element is stored for each element + if (isFlatWallBounded()) + { + wallNormalAxis_.resize(wallData.size()); + wallElementIdx_.resize(wallData.size()); + + for (const auto& element : elements(this->gridGeometry().gridView())) + { + unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); + wallElementIdx_[elementIdx] = wallData[elementIdx].eIdx; + if ( ! (hasParam("RANS.WallNormalAxis")) ) + { + GlobalPosition wallOuterNormal = wallData[elementIdx].scvfOuterNormal; + if constexpr (dim == 2) // 2D + wallNormalAxis_[elementIdx] = (wallOuterNormal[0] == 1) ? 0 : 1; + else // 3D + wallNormalAxis_[elementIdx] = (wallOuterNormal[0] == 1) ? 0 : ((wallOuterNormal[1] == 1) ? 1 : 2); + } + else + wallNormalAxis_[elementIdx] = fixedWallNormalAxis_; + } + } + } + + /*! + * \brief Store all direct neighbor indicies for each element + */ + void findNeighborIndices_() + { + // search for neighbor Idxs + for (const auto& element : elements(this->gridGeometry().gridView())) + { + unsigned int elementIdx = this->gridGeometry().elementMapper().index(element); + for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) + { + neighborIdx_[elementIdx][dimIdx][0] = elementIdx; + neighborIdx_[elementIdx][dimIdx][1] = elementIdx; + } + + for (const auto& intersection : intersections(this->gridGeometry().gridView(), element)) + { + if (intersection.boundary()) + continue; + + unsigned int neighborIdx = this->gridGeometry().elementMapper().index(intersection.outside()); + for (unsigned int dimIdx = 0; dimIdx < dim; ++dimIdx) + { + if (abs(cellCenter(elementIdx)[dimIdx] - cellCenter(neighborIdx)[dimIdx]) > 1e-8) + { + if (cellCenter(elementIdx)[dimIdx] > cellCenter(neighborIdx)[dimIdx]) + neighborIdx_[elementIdx][dimIdx][0] = neighborIdx; + + if (cellCenter(elementIdx)[dimIdx] < cellCenter(neighborIdx)[dimIdx]) + neighborIdx_[elementIdx][dimIdx][1] = neighborIdx; + } + } + } + } + } + void calculateCCVelocities_(const SolutionVector& curSol) { auto fvGeometry = localView(this->gridGeometry()); @@ -708,6 +731,11 @@ private: } } + [[deprecated("The isOnWall and IsOnWallAtPos functions will be removed after release 3.5. " + "Please use the Rans specific boundarytypes. " + "Mark wall boundaries in the rans problems with the setWall() function.")]] + void noSetWallCompilerWarning_(){} + const int fixedFlowDirectionAxis_ = getParam("RANS.FlowDirectionAxis", 0); const int fixedWallNormalAxis_ = getParam("RANS.WallNormalAxis", 1); @@ -715,7 +743,6 @@ private: std::vector flowDirectionAxis_; std::vector wallDistance_; std::vector wallElementIdx_; - std::vector cellCenter_; std::vector, dim>> neighborIdx_; std::vector velocity_; diff --git a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh index 8936d1054ccfca40550b7150683e8e29ef02e506..c80cd4c99572fda0ad8e99e12666ea4cda3d84a7 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/problem.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/problem.hh @@ -27,6 +27,7 @@ #include #include +#include #include #include #include @@ -63,6 +64,7 @@ class RANSProblemImpl : public RANSProblemBa using FVElementGeometry = typename GetPropType::LocalView; using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; + using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; using VolumeVariables = GetPropType; using SolutionVector = GetPropType; @@ -309,9 +311,22 @@ public: { unsigned int elementIdx = asImp_().gridGeometry().elementMapper().index(element); auto bcTypes = asImp_().boundaryTypes(element, localSubFace); - return asImp_().isOnWall(localSubFace) - && bcTypes.isDirichlet(eqIdx) - && isMatchingPoint(elementIdx); + + // Remove this check after release 3.5. IsOnWall Interface is deprecated + if constexpr (Deprecated::hasIsOnWall()) + { + // Remove this part + return asImp_().isOnWall(localSubFace) + && bcTypes.isDirichlet(eqIdx) + && isMatchingPoint(elementIdx); + } + else + { + // Keep this part + return bcTypes.hasWall() + && bcTypes.isDirichlet(eqIdx) + && isMatchingPoint(elementIdx); + } } //! \brief Returns an additional wall function momentum flux (only needed for RANS models) diff --git a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh index 0eca6c4ecc9ab239a1cd91c942382e21ecd41326..85fc3c92692170de99a7cf19fcbced5e9f02b35c 100644 --- a/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh +++ b/dumux/freeflow/rans/twoeq/kepsilon/staggered/fluxvariables.hh @@ -26,6 +26,7 @@ #include #include +#include #include #include #include @@ -67,6 +68,7 @@ class KEpsilonFluxVariablesImpl; using FVElementGeometry = typename GridGeometry::LocalView; using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; + using GlobalPosition = typename SubControlVolumeFace::GlobalPosition; using Extrusion = Extrusion_t; using GridView = typename GridGeometry::GridView; using Problem = GetPropType; @@ -146,20 +148,46 @@ public: } const auto bcTypes = problem.boundaryTypes(element, scvf); - if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::turbulentKineticEnergyEqIdx) - || bcTypes.isSymmetry() - || problem.isOnWall(scvf)))) + + // Remove this check after release 3.5. IsOnWall Interface is deprecated + if constexpr (Deprecated::hasIsOnWall()) { - if (!(insideVolVars.isMatchingPoint() && outsideVolVars.isMatchingPoint()) - || !(insideVolVars.isMatchingPoint() && outsideVolVars.inNearWallRegion()) - || !(insideVolVars.inNearWallRegion() && outsideVolVars.isMatchingPoint())) + noSetWallCompilerWarning_(); + // Remove this part + if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::turbulentKineticEnergyEqIdx) + || bcTypes.isSymmetry() + || problem.isOnWall(scvf)))) { - flux[turbulentKineticEnergyEqIdx] - += coeff_k / distance - * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) - * Extrusion::area(scvf); + if (!(insideVolVars.isMatchingPoint() && outsideVolVars.isMatchingPoint()) + || !(insideVolVars.isMatchingPoint() && outsideVolVars.inNearWallRegion()) + || !(insideVolVars.inNearWallRegion() && outsideVolVars.isMatchingPoint())) + { + flux[turbulentKineticEnergyEqIdx] + += coeff_k / distance + * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) + * Extrusion::area(scvf); + } } } + else + { + // Keep this part + if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::turbulentKineticEnergyEqIdx) + || bcTypes.isSymmetry() + || bcTypes.hasWall()))) + { + if (!(insideVolVars.isMatchingPoint() && outsideVolVars.isMatchingPoint()) + || !(insideVolVars.isMatchingPoint() && outsideVolVars.inNearWallRegion()) + || !(insideVolVars.inNearWallRegion() && outsideVolVars.isMatchingPoint())) + { + flux[turbulentKineticEnergyEqIdx] + += coeff_k / distance + * (insideVolVars.turbulentKineticEnergy() - outsideVolVars.turbulentKineticEnergy()) + * Extrusion::area(scvf); + } + } + } + if (!(scvf.boundary() && (bcTypes.isOutflow(Indices::dissipationEqIdx) || bcTypes.isSymmetry()))) { @@ -189,6 +217,13 @@ public: + 2.0 / ModelTraits::dim() * insideVolVars.density() * insideVolVars.turbulentKineticEnergy() * Extrusion::area(scvf) * scvf.directionSign() * insideVolVars.extrusionFactor(); } + +private: + + [[deprecated("The isOnWall and IsOnWallAtPos functions will be removed after release 3.5. " + "Please use the Rans specific boundarytypes. " + "Mark wall boundaries in the rans problems with the setWall() function.")]] + void noSetWallCompilerWarning_(){} }; } // end namespace diff --git a/test/discretization/test_walldistance.cc b/test/discretization/test_walldistance.cc index 8b30ecb5c8c80be3f7198631bacc6a55f2fdfb31..afd9152c0227c654c9be8927d49382d6d3264c54 100644 --- a/test/discretization/test_walldistance.cc +++ b/test/discretization/test_walldistance.cc @@ -75,7 +75,7 @@ void testElement(const GridGeometry& gridGeometry, const std::string& outputName const auto bottom = gridGeometry.bBoxMin()[GridGeometry::Grid::dimension-1]; // Do not consider scvfs at top or bottom of domain. - const auto considerScvf = [&](const auto& scvf){ + const auto considerScvf = [&](const auto& fvGeometry, const auto& scvf){ return scvf.ipGlobal()[GridGeometry::Grid::dimension-1] < top - 1e-6 && scvf.ipGlobal()[GridGeometry::Grid::dimension-1] > bottom + 1e-6; }; diff --git a/test/freeflow/rans/problem.hh b/test/freeflow/rans/problem.hh index 9a419eec20c3488033f6b0184f8daa7e47a4ae46..73ad3d0a6d61e9e81825aac035ac2dca96812fa9 100644 --- a/test/freeflow/rans/problem.hh +++ b/test/freeflow/rans/problem.hh @@ -32,7 +32,7 @@ #include #include -#include +#include #include #include #include @@ -51,16 +51,17 @@ class PipeLauferProblem : public RANSProblem { using ParentType = RANSProblem; - using BoundaryTypes = Dumux::NavierStokesBoundaryTypes::numEq()>; using FluidSystem = GetPropType; using FluidState = GetPropType; using GridGeometry = GetPropType; - using Indices = typename GetPropType::Indices; using PrimaryVariables = GetPropType; using NumEqVector = Dumux::NumEqVector; using Scalar = GetPropType; + static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld; using ModelTraits = GetPropType; + using Indices = typename GetPropType::Indices; + using BoundaryTypes = Dumux::RANSBoundaryTypes; using Element = typename GridGeometry::GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using FVElementGeometry = typename GetPropType::LocalView; @@ -69,7 +70,6 @@ class PipeLauferProblem : public RANSProblem using TimeLoopPtr = std::shared_ptr>; - static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld; public: PipeLauferProblem(std::shared_ptr gridGeometry) @@ -112,12 +112,6 @@ public: */ // \{ - bool isOnWallAtPos(const GlobalPosition &globalPos) const - { - return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_ - || globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; - } - /*! * \brief Returns the temperature [K] within the domain for the isothermal model. */ @@ -146,11 +140,13 @@ public: values.setDirichlet(Indices::pressureIdx); else { - // walls and inflow values.setDirichlet(Indices::velocityXIdx); values.setDirichlet(Indices::velocityYIdx); } + if (isLowerWall_(globalPos) || isUpperWall_(globalPos)) // All walls + values.setWall(); + #if NONISOTHERMAL if(isOutlet_(globalPos)) values.setOutflow(Indices::energyEqIdx); @@ -191,7 +187,7 @@ public: PrimaryVariables values(initialAtPos(globalPos)); #if NONISOTHERMAL - values[Indices::temperatureIdx] = isOnWallAtPos(globalPos) ? wallTemperature_ : inletTemperature_; + values[Indices::temperatureIdx] = (isLowerWall_(globalPos) || isUpperWall_(globalPos)) ? wallTemperature_ : inletTemperature_; #endif return values; @@ -230,10 +226,10 @@ public: values[Indices::velocityXIdx] = time() > initializationTime_ ? inletVelocity_ : time() / initializationTime_ * inletVelocity_; - if (isOnWallAtPos(globalPos)) + if (isLowerWall_(globalPos) || isUpperWall_(globalPos)) values[Indices::velocityXIdx] = 0.0; #if NONISOTHERMAL - values[Indices::temperatureIdx] = isOnWallAtPos(globalPos) ? wallTemperature_ : inletTemperature_; + values[Indices::temperatureIdx] = (isLowerWall_(globalPos) || isUpperWall_(globalPos)) ? wallTemperature_ : inletTemperature_; #endif // turbulence model-specific initial conditions setInitialAtPos_(values, globalPos); @@ -255,6 +251,12 @@ private: bool isOutlet_(const GlobalPosition& globalPos) const { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; } + bool isLowerWall_(const GlobalPosition& globalPos) const + { return globalPos[1] < this->gridGeometry().bBoxMin()[1] + eps_; } + + bool isUpperWall_(const GlobalPosition& globalPos) const + { return globalPos[1] > this->gridGeometry().bBoxMax()[1] - eps_; } + //! Initial conditions for the komega, kepsilon and lowrekepsilon turbulence models void setInitialAtPos_([[maybe_unused]] PrimaryVariables& values, [[maybe_unused]] const GlobalPosition &globalPos) const @@ -264,7 +266,7 @@ private: else if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 1) // one equation models { values[Indices::viscosityTildeIdx] = viscosityTilde_; - if (isOnWallAtPos(globalPos)) + if (isLowerWall_(globalPos) || isUpperWall_(globalPos)) values[Indices::viscosityTildeIdx] = 0.0; } else // two equation models @@ -272,7 +274,7 @@ private: static_assert(numTurbulenceEq(ModelTraits::turbulenceModel()) == 2, "Only reached by 2eq models"); values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_; values[Indices::dissipationIdx] = dissipation_; - if (isOnWallAtPos(globalPos)) + if (isLowerWall_(globalPos) || isUpperWall_(globalPos)) { values[Indices::turbulentKineticEnergyIdx] = 0.0; values[Indices::dissipationIdx] = 0.0; @@ -330,7 +332,7 @@ private: { // For the komega model we set a fixed dissipation (omega) for all cells at the wall for (const auto& scvf : scvfs(fvGeometry)) - if (isOnWallAtPos(scvf.center()) && pvIdx == Indices::dissipationIdx) + if (this->boundaryTypes(element, scvf).hasWall() && pvIdx == Indices::dissipationIdx) return true; return false; } diff --git a/test/freeflow/ransnc/problem.hh b/test/freeflow/ransnc/problem.hh index 10facf02990de049ab1e780b36222310a723895a..48021d7ff63e500da457e779b8236035f46ee413 100644 --- a/test/freeflow/ransnc/problem.hh +++ b/test/freeflow/ransnc/problem.hh @@ -29,7 +29,7 @@ #include #include -#include +#include #include #include #include @@ -53,16 +53,17 @@ class FlatPlateNCTestProblem : public RANSProblem { using ParentType = RANSProblem; - using BoundaryTypes = Dumux::NavierStokesBoundaryTypes::numEq()>; using FluidSystem = GetPropType; using FluidState = GetPropType; using GridGeometry = GetPropType; - using Indices = typename GetPropType::Indices; using PrimaryVariables = GetPropType; using NumEqVector = Dumux::NumEqVector; using Scalar = GetPropType; + static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld; using ModelTraits = GetPropType; + using Indices = typename GetPropType::Indices; + using BoundaryTypes = Dumux::RANSBoundaryTypes; using Element = typename GridGeometry::GridView::template Codim<0>::Entity; using GlobalPosition = typename Element::Geometry::GlobalCoordinate; using FVElementGeometry = typename GetPropType::LocalView; @@ -71,7 +72,6 @@ class FlatPlateNCTestProblem : public RANSProblem using TimeLoopPtr = std::shared_ptr>; - static constexpr auto dimWorld = GridGeometry::GridView::dimensionworld; static constexpr auto transportEqIdx = Indices::conti0EqIdx + 1; static constexpr auto transportCompIdx = Indices::conti0EqIdx + 1; @@ -105,14 +105,6 @@ public: std::cout << std::endl; } - /*! - * \name Problem parameters - */ - // \{ - - bool isOnWallAtPos(const GlobalPosition& globalPos) const - { return globalPos[1] < eps_; } - /*! * \brief Returns the temperature within the domain in [K]. * @@ -157,10 +149,9 @@ public: values.setOutflow(Indices::energyEqIdx); #endif } - else if(isOnWallAtPos(globalPos)) + else if(isLowerWall_(globalPos)) { - values.setDirichlet(Indices::velocityXIdx); - values.setDirichlet(Indices::velocityYIdx); + values.setWall(); values.setNeumann(transportEqIdx); #if NONISOTHERMAL values.setDirichlet(Indices::temperatureIdx); @@ -205,11 +196,9 @@ public: { values[transportCompIdx] = (time() > 10.0) ? inletMoleFraction_ : 0.0; } + #if NONISOTHERMAL - if (time() > 10.0 && isOnWallAtPos(globalPos)) - { - values[Indices::temperatureIdx] = wallTemperature_; - } + values[Indices::temperatureIdx] = (isLowerWall_(globalPos) && time() > 10.0) ? wallTemperature_ : temperature(); #endif return values; @@ -250,7 +239,7 @@ public: #endif // block velocity profile values[Indices::velocityXIdx] = 0.0; - if (!isOnWallAtPos(globalPos)) + if (!isLowerWall_(globalPos)) values[Indices::velocityXIdx] = inletVelocity_; values[Indices::velocityYIdx] = 0.0; @@ -274,6 +263,9 @@ private: bool isOutlet_(const GlobalPosition& globalPos) const { return globalPos[0] > this->gridGeometry().bBoxMax()[0] - eps_; } + bool isLowerWall_(const GlobalPosition& globalPos) const + { return globalPos[1] < eps_; } + //! Initial conditions for the komega, kepsilon and lowrekepsilon turbulence models void setInitialAtPos_([[maybe_unused]] PrimaryVariables& values, [[maybe_unused]] const GlobalPosition &globalPos) const @@ -283,7 +275,7 @@ private: else if constexpr (numTurbulenceEq(ModelTraits::turbulenceModel()) == 1) // one equation models { values[Indices::viscosityTildeIdx] = viscosityTilde_; - if (isOnWallAtPos(globalPos)) + if (isLowerWall_(globalPos)) values[Indices::viscosityTildeIdx] = 0.0; } else // two equation models @@ -291,7 +283,7 @@ private: static_assert(numTurbulenceEq(ModelTraits::turbulenceModel()) == 2, "Only reached by 2eq models"); values[Indices::turbulentKineticEnergyIdx] = turbulentKineticEnergy_; values[Indices::dissipationIdx] = dissipation_; - if (isOnWallAtPos(globalPos)) + if (isLowerWall_(globalPos)) { values[Indices::turbulentKineticEnergyIdx] = 0.0; values[Indices::dissipationIdx] = 0.0; @@ -349,7 +341,7 @@ private: { // For the komega model we set a fixed dissipation (omega) for all cells at the wall for (const auto& scvf : scvfs(fvGeometry)) - if (isOnWallAtPos(scvf.center()) && pvIdx == Indices::dissipationIdx) + if (this->boundaryTypes(element, scvf).hasWall() && pvIdx == Indices::dissipationIdx) return true; return false; } diff --git a/test/io/vtk/test_vtk_staggeredfreeflowpvnames.cc b/test/io/vtk/test_vtk_staggeredfreeflowpvnames.cc index 2e9cfc15f3870eb0f1bd4f01eed95b42defcad4e..ee3cc5b3beaa25ad11deec68210231b8b05cff58 100644 --- a/test/io/vtk/test_vtk_staggeredfreeflowpvnames.cc +++ b/test/io/vtk/test_vtk_staggeredfreeflowpvnames.cc @@ -37,7 +37,7 @@ #include #include -#include +#include #include #include #include @@ -126,7 +126,8 @@ private: class MockProblem : public BaseProblem { using ParentType = BaseProblem; - using BoundaryTypes = Dumux::NavierStokesBoundaryTypes::numEq()>; + using BoundaryTypes = Dumux::RANSBoundaryTypes; + using PrimaryVariables = GetPropType; using Scalar = GetPropType; using Traits = GetPropType; Scalar eps_ = 1e-6; @@ -138,9 +139,15 @@ private: BoundaryTypes boundaryTypesAtPos(const GlobalPosition &globalPos) const { BoundaryTypes values; + if (globalPos[0] < eps_) + values.setWall(); return values; } + template + PrimaryVariables dirichletAtPos(const GlobalPosition &globalPos) const + { return PrimaryVariables(0.0); } + Scalar temperature() const { return 300; } @@ -160,9 +167,6 @@ private: void updateDynamicWallProperties(const U& u) { return ParentType::updateDynamicWallProperties(u); } - template - bool isOnWall(const Scvf& scvf) const - { return (scvf.center()[0] < eps_); } }; public: