diff --git a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh index 071888dfddb0381fe66e3ae4b540fe40dadb856d..43772a00aa3ba2641617d0c631e4446081b56676 100644 --- a/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh +++ b/dumux/discretization/staggered/freeflow/staggeredgeometryhelper.hh @@ -52,7 +52,7 @@ struct PairData std::bitset<upwindSchemeOrder> hasParallelNeighbor; std::array<GridIndexType, upwindSchemeOrder> parallelDofs; - std::array<Scalar, upwindSchemeOrder+1> parallelDistances; // TODO store only two distances, not three. + std::array<Scalar, upwindSchemeOrder> parallelCellWidths; bool hasNormalNeighbor = false; std::pair<GridIndexType, GridIndexType> normalPair; SmallLocalIndexType localNormalFaceIdx; @@ -424,17 +424,14 @@ private: { if (facetIsNormal_(intersection.indexInInside(), parallelLocalIdx)) { - // Store the parallel dimension of self cell in the direction of the axis auto parallelAxisIdx = directionIndex(intersection); - pairData_[numPairParallelIdx].parallelDistances[0] = setParallelPairDistances_(element_, parallelAxisIdx); - if (intersection.neighbor()) { // If the normal intersection has a neighboring cell, go in and store the parallel information. const auto& outerElement = intersection.outside(); pairData_[numPairParallelIdx].hasParallelNeighbor.set(0, true); pairData_[numPairParallelIdx].parallelDofs[0] = gridView_.indexSet().subIndex(outerElement, parallelLocalIdx, codimIntersection); - pairData_[numPairParallelIdx].parallelDistances[1] = setParallelPairDistances_(outerElement, parallelAxisIdx); + pairData_[numPairParallelIdx].parallelCellWidths[0] = setParallelPairCellWidths_(outerElement, parallelAxisIdx); } else // No parallel neighbor available { @@ -457,7 +454,7 @@ private: void fillParallelPairData_(std::true_type) { // set basic global positions and stencil size definitions - const auto numParallelFaces = pairData_[0].parallelDistances.size(); + const auto numParallelFaces = pairData_[0].parallelCellWidths.size(); const auto& elementCenter = element_.geometry().center(); // get the parallel Dofs @@ -472,12 +469,11 @@ private: { auto parallelAxisIdx = directionIndex(intersection); auto localNormalIntersectionIndex = intersection.indexInInside(); - parallelElementStack.push(element_); + auto e = element_; bool keepStacking = (parallelElementStack.size() < numParallelFaces); while(keepStacking) { - auto e = parallelElementStack.top(); for(const auto& normalIntersection : intersections(gridView_, e)) { if( normalIntersection.indexInInside() == localNormalIntersectionIndex ) @@ -489,20 +485,18 @@ private: } else { - keepStacking = false; + keepStacking = false; } } } + e = parallelElementStack.top(); } while(!parallelElementStack.empty()) { - if(parallelElementStack.size() > 1) - { - pairData_[numPairParallelIdx].hasParallelNeighbor.set(parallelElementStack.size()-2, true); - pairData_[numPairParallelIdx].parallelDofs[parallelElementStack.size()-2] = gridView_.indexSet().subIndex(parallelElementStack.top(), parallelLocalIdx, codimIntersection); - } - pairData_[numPairParallelIdx].parallelDistances[parallelElementStack.size()-1] = setParallelPairDistances_(parallelElementStack.top(), parallelAxisIdx); + pairData_[numPairParallelIdx].hasParallelNeighbor.set(parallelElementStack.size()-1, true); + pairData_[numPairParallelIdx].parallelDofs[parallelElementStack.size()-1] = gridView_.indexSet().subIndex(parallelElementStack.top(), parallelLocalIdx, codimIntersection); + pairData_[numPairParallelIdx].parallelCellWidths[parallelElementStack.size()-1] = setParallelPairCellWidths_(parallelElementStack.top(), parallelAxisIdx); parallelElementStack.pop(); } @@ -516,10 +510,6 @@ private: const auto virtualFirstParallelFaceDofPos = intersection_.geometry().center() + distance; pairData_[numPairParallelIdx].virtualFirstParallelFaceDofPos = std::move(virtualFirstParallelFaceDofPos); - - // The distance is saved doubled because with scvf.cellCenteredSelfToFirstParallelDistance - // an average between parallelDistances[0] and parallelDistances[1] will be computed - pairData_[numPairParallelIdx].parallelDistances[0] = std::move(distance.two_norm() * 2); } numPairParallelIdx++; } @@ -578,9 +568,8 @@ private: } //! Sets the information about the parallel distances - Scalar setParallelPairDistances_(const Element& element, const int parallelAxisIdx) const + Scalar setParallelPairCellWidths_(const Element& element, const int parallelAxisIdx) const { - Scalar distance = 0; std::vector<GlobalPosition> faces; faces.reserve(numfacets); for (const auto& intElement : intersections(gridView_, element)) @@ -590,27 +579,17 @@ private: switch(parallelAxisIdx) { case 0: - { - distance = (faces[1] - faces[0]).two_norm(); - break; - } + return (faces[1] - faces[0]).two_norm(); case 1: - { - distance = (faces[3] - faces[2]).two_norm(); - break; - } + return (faces[3] - faces[2]).two_norm(); case 2: { assert(dim == 3); - distance = (faces[5] - faces[4]).two_norm(); - break; + return (faces[5] - faces[4]).two_norm(); } default: - { DUNE_THROW(Dune::InvalidStateException, "Something went terribly wrong"); - } } - return distance; } Intersection intersection_; //!< The intersection of interest diff --git a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh index 9f163725c141730cb9a1906fa099adbddbb668d5..4842324a2b85c2ef415d5a9c9021d9ba811900e6 100644 --- a/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh +++ b/dumux/discretization/staggered/freeflow/subcontrolvolumeface.hh @@ -252,6 +252,20 @@ public: return isGhostFace_; } + //! Returns the length of the face in a certain direction (adaptation of area() for 3d) + Scalar faceLength(const int localSubFaceIdx) const + { + if (dimworld == 3) + { + if (localSubFaceIdx < 2) + return (corner(1) - corner(0)).two_norm(); + else + return (corner(2) - corner(0)).two_norm(); + } + else + return (corner(1) - corner(0)).two_norm(); + } + /*! * \brief Check if the face has a parallel neighbor * @@ -333,8 +347,13 @@ public: */ Scalar cellCenteredParallelDistance(const int localSubFaceIdx, const int parallelDegreeIdx) const { - return (pairData(localSubFaceIdx).parallelDistances[parallelDegreeIdx] + - pairData(localSubFaceIdx).parallelDistances[parallelDegreeIdx+1]) * 0.5; + if (parallelDegreeIdx == 0) + return (faceLength(localSubFaceIdx) + pairData(localSubFaceIdx).parallelCellWidths[0]) * 0.5; + else + { + assert((parallelDegreeIdx == 1) && "Only the width of the first two parallel cells (indicies 0 and 1) is stored for each scvf."); + return (pairData(localSubFaceIdx).parallelCellWidths[0] + pairData(localSubFaceIdx).parallelCellWidths[1]) * 0.5; + } } diff --git a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh index 39ae0b74f5c856bf8352fcd5d316e51b0499a879..544eeac6c4cad4a9681bc9abae0d23b650428993 100644 --- a/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh +++ b/dumux/freeflow/navierstokes/staggered/staggeredupwindfluxvariables.hh @@ -550,7 +550,7 @@ private: distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0); distances[1] = ownScvf.cellCenteredParallelDistance(oppositeSubFaceIdx, 0); if(ownScvf.hasParallelNeighbor(localSubFaceIdx, 0)) - distances[2] = ownScvf.pairData(localSubFaceIdx).parallelDistances[1]; + distances[2] = ownScvf.pairData(localSubFaceIdx).parallelCellWidths[0]; else distances[2] = ownScvf.area() / 2.0; } @@ -558,7 +558,7 @@ private: { distances[0] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 0); distances[1] = ownScvf.cellCenteredParallelDistance(localSubFaceIdx, 1); - distances[2] = ownScvf.pairData(localSubFaceIdx).parallelDistances[1]; + distances[2] = ownScvf.pairData(localSubFaceIdx).parallelCellWidths[0]; } return distances; diff --git a/test/freeflow/navierstokes/kovasznay/CMakeLists.txt b/test/freeflow/navierstokes/kovasznay/CMakeLists.txt index 3b2f353413d0744c1d8f05116eb4fbadeee4c1e0..6d46e035ac713adaf9aab4fa0845e77ae98ae265 100644 --- a/test/freeflow/navierstokes/kovasznay/CMakeLists.txt +++ b/test/freeflow/navierstokes/kovasznay/CMakeLists.txt @@ -12,7 +12,7 @@ dune_add_test(NAME test_ff_navierstokes_kovasznay dune_add_test(NAME test_ff_navierstokes_kovasznay_higherorder SOURCES main.cc - LABELS freeflow + LABELS freeflow navierstokes CMAKE_GUARD HAVE_UMFPACK COMPILE_DEFINITIONS UPWINDSCHEMEORDER=2 COMMAND ${CMAKE_SOURCE_DIR}/bin/testing/runtest.py