| Line | Branch | Exec | Source |
|---|---|---|---|
| 1 | // -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- | ||
| 2 | // vi: set et ts=4 sw=4 sts=4: | ||
| 3 | // | ||
| 4 | // SPDX-FileCopyrightInfo: Copyright © DuMux Project contributors, see AUTHORS.md in root folder | ||
| 5 | // SPDX-License-Identifier: GPL-3.0-or-later | ||
| 6 | // | ||
| 7 | /*! | ||
| 8 | * \file | ||
| 9 | * \ingroup NavierStokesModel | ||
| 10 | * \copydoc Dumux::StaggeredUpwindHelper | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH | ||
| 13 | #define DUMUX_NAVIERSTOKES_STAGGERED_UPWINDVARIABLES_HH | ||
| 14 | |||
| 15 | #include <array> | ||
| 16 | #include <optional> | ||
| 17 | |||
| 18 | #include <dumux/common/math.hh> | ||
| 19 | #include <dumux/common/exceptions.hh> | ||
| 20 | #include <dumux/common/parameters.hh> | ||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | #include <dumux/common/typetraits/problem.hh> | ||
| 23 | |||
| 24 | #include <dumux/discretization/method.hh> | ||
| 25 | #include <dumux/freeflow/staggeredupwindmethods.hh> | ||
| 26 | #include "velocitygradients.hh" | ||
| 27 | |||
| 28 | namespace Dumux { | ||
| 29 | |||
| 30 | /*! | ||
| 31 | * \ingroup NavierStokesModel | ||
| 32 | * \brief The upwinding variables class for the Navier-Stokes model using the staggered grid discretization. | ||
| 33 | */ | ||
| 34 | template<class TypeTag, int upwindSchemeOrder> | ||
| 35 | class StaggeredUpwindHelper | ||
| 36 | { | ||
| 37 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 38 | |||
| 39 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 40 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 41 | using VolumeVariables = typename GridVolumeVariables::VolumeVariables; | ||
| 42 | |||
| 43 | using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache; | ||
| 44 | using UpwindScheme = typename GridFluxVariablesCache::UpwindScheme; | ||
| 45 | using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache; | ||
| 46 | |||
| 47 | using GridFaceVariables = typename GridVariables::GridFaceVariables; | ||
| 48 | using ElementFaceVariables = typename GridFaceVariables::LocalView; | ||
| 49 | using FaceVariables = typename GridFaceVariables::FaceVariables; | ||
| 50 | |||
| 51 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 52 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 53 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 54 | using GridView = typename GridGeometry::GridView; | ||
| 55 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 56 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 57 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 58 | using Indices = typename ModelTraits::Indices; | ||
| 59 | using SubControlVolume = typename FVElementGeometry::SubControlVolume; | ||
| 60 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 61 | using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 62 | using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>; | ||
| 63 | using BoundaryTypes = typename ProblemTraits<Problem>::BoundaryTypes; | ||
| 64 | using VelocityGradients = StaggeredVelocityGradients<Scalar, GridGeometry, BoundaryTypes, Indices>; | ||
| 65 | |||
| 66 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 67 | static constexpr bool useHigherOrder = upwindSchemeOrder > 1; | ||
| 68 | |||
| 69 | static_assert(upwindSchemeOrder <= 2, "Not implemented: Order higher than 2!"); | ||
| 70 | static_assert(upwindSchemeOrder <= 1 || GridFluxVariablesCache::cachingEnabled, | ||
| 71 | "Higher order upwind method requires caching enabled for the GridFluxVariablesCache!"); | ||
| 72 | static_assert(upwindSchemeOrder <= 1 || GridGeometry::cachingEnabled, | ||
| 73 | "Higher order upwind method requires caching enabled for the GridGeometry!"); | ||
| 74 | static_assert(upwindSchemeOrder <= 1 || GridFaceVariables::cachingEnabled, | ||
| 75 | "Higher order upwind method requires caching enabled for the GridFaceVariables!"); | ||
| 76 | static_assert(upwindSchemeOrder <= 1 || GridVolumeVariables::cachingEnabled, | ||
| 77 | "Higher order upwind method requires caching enabled for the GridGeometry!"); | ||
| 78 | |||
| 79 | public: | ||
| 80 | 405836084 | StaggeredUpwindHelper(const Element& element, | |
| 81 | const FVElementGeometry& fvGeometry, | ||
| 82 | const SubControlVolumeFace& scvf, | ||
| 83 | const ElementFaceVariables& elemFaceVars, | ||
| 84 | const ElementVolumeVariables& elemVolVars, | ||
| 85 | const UpwindScheme& upwindScheme) | ||
| 86 | : element_(element) | ||
| 87 | , fvGeometry_(fvGeometry) | ||
| 88 | , scvf_(scvf) | ||
| 89 | , elemFaceVars_(elemFaceVars) | ||
| 90 | 811672168 | , faceVars_(elemFaceVars[scvf]) | |
| 91 | , elemVolVars_(elemVolVars) | ||
| 92 | 811672168 | , upwindScheme_(upwindScheme) | |
| 93 | {} | ||
| 94 | |||
| 95 | /*! | ||
| 96 | * \brief Returns the momentum in the frontal direction. | ||
| 97 | * | ||
| 98 | * Checks if the model has higher order methods enabled and if the scvf in | ||
| 99 | * question is far enough from the boundary such that higher order methods can be employed. | ||
| 100 | * Then the corresponding set of momenta are collected and the prescribed | ||
| 101 | * upwinding method is used to calculate the momentum. | ||
| 102 | */ | ||
| 103 | 135725860 | FacePrimaryVariables computeUpwindFrontalMomentum(const bool selfIsUpstream) const | |
| 104 | { | ||
| 105 |
4/4✓ Branch 2 taken 68511040 times.
✓ Branch 3 taken 67214820 times.
✓ Branch 4 taken 46109799 times.
✓ Branch 5 taken 45514243 times.
|
271451720 | const auto density = elemVolVars_[scvf_.insideScvIdx()].density(); |
| 106 | |||
| 107 | // for higher order schemes do higher order upwind reconstruction | ||
| 108 | if constexpr (useHigherOrder) | ||
| 109 | { | ||
| 110 | // only second order is implemented so far | ||
| 111 |
2/2✓ Branch 0 taken 1925144 times.
✓ Branch 1 taken 22452 times.
|
1947596 | if (canDoFrontalSecondOrder_(selfIsUpstream)) |
| 112 | { | ||
| 113 | 3850288 | const auto distances = getFrontalDistances_(selfIsUpstream); | |
| 114 | 3850288 | const auto upwindMomenta = getFrontalSecondOrderUpwindMomenta_(density, selfIsUpstream); | |
| 115 | 1925144 | return upwindScheme_.tvd(upwindMomenta, distances, selfIsUpstream, upwindScheme_.tvdApproach()); | |
| 116 | } | ||
| 117 | } | ||
| 118 | |||
| 119 | // otherwise apply first order upwind scheme | ||
| 120 | 267601432 | const auto upwindMomenta = getFrontalFirstOrderUpwindMomenta_(density, selfIsUpstream); | |
| 121 | 535202864 | return upwindScheme_.upwind(upwindMomenta[0], upwindMomenta[1]); | |
| 122 | } | ||
| 123 | |||
| 124 | /*! | ||
| 125 | * \brief Returns the momentum in the lateral direction. | ||
| 126 | * | ||
| 127 | * Evaluates which face is upstream. | ||
| 128 | * Checks if the model has higher order methods enabled and if the scvf in | ||
| 129 | * question is far enough from the boundary such that higher order methods can be employed. | ||
| 130 | * Then the corresponding set of momenta are collected and the prescribed | ||
| 131 | * upwinding method is used to calculate the momentum. | ||
| 132 | */ | ||
| 133 | 270110224 | FacePrimaryVariables computeUpwindLateralMomentum(const bool selfIsUpstream, | |
| 134 | const SubControlVolumeFace& lateralFace, | ||
| 135 | const int localSubFaceIdx, | ||
| 136 | const std::optional<BoundaryTypes>& currentScvfBoundaryTypes, | ||
| 137 | const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const | ||
| 138 | { | ||
| 139 | // Check whether the own or the neighboring element is upstream. | ||
| 140 | // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof | ||
| 141 | // of interest is located. | ||
| 142 | 540220448 | const Scalar insideDensity = elemVolVars_[lateralFace.insideScvIdx()].density(); | |
| 143 |
2/2✓ Branch 2 taken 2216760 times.
✓ Branch 3 taken 1678432 times.
|
540220448 | const Scalar outsideDensity = elemVolVars_[lateralFace.outsideScvIdx()].density(); |
| 144 | |||
| 145 | // for higher order schemes do higher order upwind reconstruction | ||
| 146 | if constexpr (useHigherOrder) | ||
| 147 | { | ||
| 148 | // only second order is implemented so far | ||
| 149 | 3940427 | if (canDoLateralSecondOrder_(selfIsUpstream, localSubFaceIdx)) | |
| 150 | { | ||
| 151 | 3849957 | const auto upwindMomenta = getLateralSecondOrderUpwindMomenta_(insideDensity, outsideDensity, selfIsUpstream, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes); | |
| 152 | 3849957 | const auto distances = getLateralDistances_(localSubFaceIdx, selfIsUpstream); | |
| 153 | 3849957 | return upwindScheme_.tvd(upwindMomenta, distances, selfIsUpstream, upwindScheme_.tvdApproach()); | |
| 154 | } | ||
| 155 | } | ||
| 156 | |||
| 157 | // otherwise apply first order upwind scheme | ||
| 158 | 266260267 | const auto upwindMomenta = getLateralFirstOrderUpwindMomenta_(insideDensity, outsideDensity, selfIsUpstream, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes); | |
| 159 | 1065041068 | return upwindScheme_.upwind(upwindMomenta[0], upwindMomenta[1]); | |
| 160 | } | ||
| 161 | |||
| 162 | private: | ||
| 163 | /*! | ||
| 164 | * \brief Returns whether or not the face in question is far enough from the wall to handle higher order methods. | ||
| 165 | * | ||
| 166 | * Evaluates which face is upstream. | ||
| 167 | * If the face is upstream, and the scvf has a forward neighbor, higher order methods are possible. | ||
| 168 | * If the face is downstream, and the scvf has a backwards neighbor, higher order methods are possible. | ||
| 169 | * Otherwise, higher order methods are not possible. | ||
| 170 | */ | ||
| 171 | bool canDoFrontalSecondOrder_(bool selfIsUpstream) const | ||
| 172 | { | ||
| 173 | // Depending on selfIsUpstream I have to check if I have a forward or a backward neighbor to retrieve | ||
| 174 |
2/2✓ Branch 0 taken 1090576 times.
✓ Branch 1 taken 857020 times.
|
1947596 | return selfIsUpstream ? scvf_.hasForwardNeighbor(0) : scvf_.hasBackwardNeighbor(0); |
| 175 | } | ||
| 176 | |||
| 177 | /*! | ||
| 178 | * \brief Returns an array of the three momenta needed for second order upwinding methods. | ||
| 179 | */ | ||
| 180 | std::array<Scalar, 3> getFrontalSecondOrderUpwindMomenta_(const Scalar density, bool selfIsUpstream) const | ||
| 181 | { | ||
| 182 | static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled"); | ||
| 183 | 1925144 | const Scalar momentumSelf = faceVars_.velocitySelf() * density; | |
| 184 | 1925144 | const Scalar momentumOpposite = faceVars_.velocityOpposite() * density; | |
| 185 |
2/2✓ Branch 0 taken 1090576 times.
✓ Branch 1 taken 834568 times.
|
1925144 | if (selfIsUpstream) |
| 186 | 2181152 | return {momentumOpposite, momentumSelf, faceVars_.velocityForward(0)*density}; | |
| 187 | else | ||
| 188 | 1669136 | return {momentumSelf, momentumOpposite, faceVars_.velocityBackward(0)*density}; | |
| 189 | } | ||
| 190 | |||
| 191 | /*! | ||
| 192 | * \brief Returns an array of the two momenta needed for first order upwinding method | ||
| 193 | */ | ||
| 194 | std::array<Scalar, 2> getFrontalFirstOrderUpwindMomenta_(const Scalar density, bool selfIsUpstream) const | ||
| 195 | { | ||
| 196 | 133800716 | const Scalar momentumSelf = faceVars_.velocitySelf() * density; | |
| 197 | 133800716 | const Scalar momentumOpposite = faceVars_.velocityOpposite() * density; | |
| 198 |
2/2✓ Branch 0 taken 67442916 times.
✓ Branch 1 taken 66357800 times.
|
133800716 | if (selfIsUpstream) |
| 199 | return {momentumOpposite, momentumSelf}; | ||
| 200 | else | ||
| 201 | return {momentumSelf, momentumOpposite}; | ||
| 202 | } | ||
| 203 | |||
| 204 | /*! | ||
| 205 | * \brief Returns an array of distances needed for non-uniform higher order upwind schemes | ||
| 206 | * Depending on selfIsUpstream the downstream and the (up)upstream distances are saved. | ||
| 207 | * distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size} | ||
| 208 | */ | ||
| 209 | std::array<Scalar, 3> getFrontalDistances_(const bool selfIsUpstream) const | ||
| 210 | { | ||
| 211 | static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled"); | ||
| 212 |
2/2✓ Branch 0 taken 1090576 times.
✓ Branch 1 taken 834568 times.
|
1925144 | if (selfIsUpstream) |
| 213 | { | ||
| 214 | std::array<Scalar, 3> distances; | ||
| 215 | 2181152 | distances[0] = scvf_.selfToOppositeDistance(); | |
| 216 | 2181152 | distances[1] = scvf_.axisData().inAxisForwardDistances[0]; | |
| 217 | 2181152 | distances[2] = 0.5 * (scvf_.axisData().selfToOppositeDistance + scvf_.axisData().inAxisBackwardDistances[0]); | |
| 218 | 1090576 | return distances; | |
| 219 | } | ||
| 220 | else | ||
| 221 | { | ||
| 222 | std::array<Scalar, 3> distances; | ||
| 223 | 1669136 | distances[0] = scvf_.selfToOppositeDistance(); | |
| 224 | 1669136 | distances[1] = scvf_.axisData().inAxisBackwardDistances[0]; | |
| 225 | 1669136 | distances[2] = 0.5 * (scvf_.axisData().selfToOppositeDistance + scvf_.axisData().inAxisForwardDistances[0]); | |
| 226 | 834568 | return distances; | |
| 227 | } | ||
| 228 | } | ||
| 229 | |||
| 230 | /*! | ||
| 231 | * \brief Check if a second order approximation for the lateral part of the advective term can be used | ||
| 232 | * | ||
| 233 | * This helper function checks if the scvf of interest is not too near to the | ||
| 234 | * boundary so that a dof upstream with respect to the upstream dof is available. | ||
| 235 | */ | ||
| 236 | bool canDoLateralSecondOrder_(const bool selfIsUpstream, const int localSubFaceIdx) const | ||
| 237 | { | ||
| 238 | // If the self velocity is upstream, the downstream velocity can be assigned or retrieved | ||
| 239 | // from the boundary, even if there is no parallel neighbor. | ||
| 240 | // If the self velocity is downstream and there is no parallel neighbor I cannot use a second order approximation. | ||
| 241 |
6/6✓ Branch 0 taken 2216760 times.
✓ Branch 1 taken 1678432 times.
✓ Branch 2 taken 45235 times.
✓ Branch 3 taken 2171525 times.
✓ Branch 4 taken 45235 times.
✓ Branch 5 taken 2171525 times.
|
3895192 | return selfIsUpstream || scvf_.hasParallelNeighbor(localSubFaceIdx, 0); |
| 242 | } | ||
| 243 | |||
| 244 | /*! | ||
| 245 | * \brief Returns an array of momenta needed for higher order or calls a function to return an array for basic upwinding methods. | ||
| 246 | * TODO: In order to get a second order momentum upwind scheme for compressible flow the densities have to be evaluated | ||
| 247 | * at the same integration points / positions as the velocities. The currently implementation just takes the closest upwind density | ||
| 248 | * to compute the momentum as a crude approximation. | ||
| 249 | * | ||
| 250 | * ------------ | ||
| 251 | * | xxxx o | ||
| 252 | * | xxxx o | ||
| 253 | * | xxxx o | ||
| 254 | * -----------*bbbbbbbbbbb | ||
| 255 | * | yyyy o zzzz | | ||
| 256 | * | yyyy o zzzz | | ||
| 257 | * | yyyy o zzzz | | ||
| 258 | * ----------------------- | ||
| 259 | * If scvf_ is touching a corner, at which there is a Dirichlet condition for face b (half-control | ||
| 260 | * volumes x, y and z), the transported velocity at * is given. No upwinding or | ||
| 261 | * higher-order approximation for the velocity at * are required. This also means the transported | ||
| 262 | * velocity at * is the same for the half-control volumes y and z. | ||
| 263 | * | ||
| 264 | * ------------ | ||
| 265 | * | | | ||
| 266 | * | | | ||
| 267 | * | | | ||
| 268 | * ----------------------- | ||
| 269 | * | xxxx o wwww | | ||
| 270 | * | xxxx o wwww | | ||
| 271 | * | xxxx o wwww | | ||
| 272 | * ------+++++*~~~~~------ | ||
| 273 | * | yyyy o zzzz | | ||
| 274 | * | yyyy o zzzz | | ||
| 275 | * | yyyy o zzzz | | ||
| 276 | * ----------------------- | ||
| 277 | * If the flux over face + is calculated and a corner occurs a bit further away (here upper right), | ||
| 278 | * no special treatment of the corner geometry is provided. This means, the transported velocity at | ||
| 279 | * star (*) is obtained with an upwind scheme. In particularly, this also means | ||
| 280 | * that while x and y see the same velocity * for the flux over face x (continuity OK), z sees | ||
| 281 | * another velocity * for the flux over face ~ (continuity still OK, as only z and w have to use the | ||
| 282 | * same velocity at *). | ||
| 283 | */ | ||
| 284 | 3849957 | std::array<Scalar, 3> getLateralSecondOrderUpwindMomenta_(const Scalar insideDensity, | |
| 285 | const Scalar outsideDensity, | ||
| 286 | const bool selfIsUpstream, | ||
| 287 | const int localSubFaceIdx, | ||
| 288 | const std::optional<BoundaryTypes>& currentScvfBoundaryTypes, | ||
| 289 | const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const | ||
| 290 | { | ||
| 291 | static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled"); | ||
| 292 | |||
| 293 | // If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known, | ||
| 294 | // thus we always use this value for the computation of the transported momentum. | ||
| 295 |
10/12✓ Branch 0 taken 3828860 times.
✓ Branch 1 taken 21097 times.
✓ Branch 2 taken 3828860 times.
✓ Branch 3 taken 21097 times.
✓ Branch 4 taken 3823512 times.
✓ Branch 5 taken 5348 times.
✓ Branch 6 taken 3823512 times.
✓ Branch 7 taken 5348 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 3823512 times.
✗ Branch 10 not taken.
✓ Branch 11 taken 3823512 times.
|
7699914 | if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0) || scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) |
| 296 | { | ||
| 297 |
7/10✓ Branch 0 taken 21097 times.
✓ Branch 1 taken 5348 times.
✓ Branch 2 taken 21097 times.
✓ Branch 3 taken 5348 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 21097 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 21097 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 5348 times.
|
52890 | if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx)) |
| 298 | { | ||
| 299 | 5348 | const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx); | |
| 300 | 5348 | const Scalar boundaryMomentum = boundaryVelocity*insideDensity; | |
| 301 | |||
| 302 | 5348 | return {boundaryMomentum, boundaryMomentum, boundaryMomentum}; | |
| 303 | } | ||
| 304 |
2/4✓ Branch 0 taken 21097 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 21097 times.
✗ Branch 3 not taken.
|
42194 | else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0)) |
| 305 | { | ||
| 306 | 21097 | const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); | |
| 307 | 21097 | const Scalar boundaryMomentum = boundaryVelocity*insideDensity; | |
| 308 | |||
| 309 | 21097 | return {boundaryMomentum, boundaryMomentum, boundaryMomentum}; | |
| 310 | } | ||
| 311 | } | ||
| 312 | |||
| 313 |
2/2✓ Branch 0 taken 1653881 times.
✓ Branch 1 taken 2169631 times.
|
3823512 | if (selfIsUpstream) |
| 314 | { | ||
| 315 | std::array<Scalar, 3> momenta; | ||
| 316 |
2/2✓ Branch 0 taken 1186089 times.
✓ Branch 1 taken 467792 times.
|
1653881 | momenta[1] = faceVars_.velocitySelf()*insideDensity; |
| 317 |
4/4✓ Branch 0 taken 1186089 times.
✓ Branch 1 taken 467792 times.
✓ Branch 2 taken 1186089 times.
✓ Branch 3 taken 467792 times.
|
3307762 | momenta[0] = faceVars_.velocityParallel(localSubFaceIdx, 0)*insideDensity; |
| 318 | |||
| 319 | // The local index of the faces that is opposite to localSubFaceIdx | ||
| 320 |
2/2✓ Branch 0 taken 1186089 times.
✓ Branch 1 taken 467792 times.
|
1653881 | const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1; |
| 321 | |||
| 322 | // The "upstream-upstream" velocity is retrieved from the other parallel neighbor or from the boundary. | ||
| 323 |
4/4✓ Branch 0 taken 1616104 times.
✓ Branch 1 taken 37777 times.
✓ Branch 2 taken 1616104 times.
✓ Branch 3 taken 37777 times.
|
3307762 | if (scvf_.hasParallelNeighbor(oppositeSubFaceIdx, 0)) |
| 324 | 4848312 | momenta[2] = faceVars_.velocityParallel(oppositeSubFaceIdx, 0)*insideDensity; | |
| 325 | else | ||
| 326 | 37777 | momenta[2] = getParallelVelocityFromOppositeBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, oppositeSubFaceIdx)*insideDensity; | |
| 327 | 1653881 | return momenta; | |
| 328 | } | ||
| 329 | else | ||
| 330 | { | ||
| 331 | std::array<Scalar, 3> momenta; | ||
| 332 |
2/2✓ Branch 0 taken 2115753 times.
✓ Branch 1 taken 53878 times.
|
2169631 | momenta[0] = faceVars_.velocitySelf()*outsideDensity; |
| 333 |
4/4✓ Branch 0 taken 2115753 times.
✓ Branch 1 taken 53878 times.
✓ Branch 2 taken 2115753 times.
✓ Branch 3 taken 53878 times.
|
4339262 | momenta[1] = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity; |
| 334 | |||
| 335 | // If there is another parallel neighbor I can assign the "upstream-upstream" velocity, otherwise I retrieve it from the boundary. | ||
| 336 |
4/4✓ Branch 0 taken 2115753 times.
✓ Branch 1 taken 53878 times.
✓ Branch 2 taken 2115753 times.
✓ Branch 3 taken 53878 times.
|
4339262 | if (scvf_.hasParallelNeighbor(localSubFaceIdx, 1)) |
| 337 | 6347259 | momenta[2] = faceVars_.velocityParallel(localSubFaceIdx, 1)*outsideDensity; | |
| 338 | else | ||
| 339 | { | ||
| 340 |
3/6✓ Branch 1 taken 53878 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53878 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 53878 times.
✗ Branch 8 not taken.
|
161634 | const auto& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx); |
| 341 |
2/4✓ Branch 1 taken 53878 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53878 times.
✗ Branch 5 not taken.
|
107756 | const auto& elementParallel = fvGeometry_.gridGeometry().element(lateralFace.outsideScvIdx()); |
| 342 | 161634 | const auto& firstParallelScvf = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx()); | |
| 343 | 53878 | const auto& problem = elemVolVars_.gridVolVars().problem(); | |
| 344 | 53878 | const auto& boundaryTypes = problem.boundaryTypes(elementParallel, firstParallelScvf); | |
| 345 |
3/6✓ Branch 1 taken 53878 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 53878 times.
✗ Branch 5 not taken.
✓ Branch 7 taken 53878 times.
✗ Branch 8 not taken.
|
161634 | momenta[2] = getParallelVelocityFromOppositeBoundary_(elementParallel, firstParallelScvf, elemFaceVars_[firstParallelScvf], boundaryTypes, localSubFaceIdx)*outsideDensity; |
| 346 | } | ||
| 347 | 2169631 | return momenta; | |
| 348 | } | ||
| 349 | } | ||
| 350 | |||
| 351 | /*! | ||
| 352 | * \brief Returns an array of momenta needed for basic upwinding methods. | ||
| 353 | */ | ||
| 354 | 266260267 | std::array<Scalar, 2> getLateralFirstOrderUpwindMomenta_(const Scalar insideDensity, | |
| 355 | const Scalar outsideDensity, | ||
| 356 | const bool selfIsUpstream, | ||
| 357 | const int localSubFaceIdx, | ||
| 358 | const std::optional<BoundaryTypes>& currentScvfBoundaryTypes, | ||
| 359 | const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes) const | ||
| 360 | { | ||
| 361 | // If the lateral face lies on a boundary, we assume that the parallel velocity on the boundary is actually known, | ||
| 362 | // thus we always use this value for the computation of the transported momentum. | ||
| 363 |
4/10✓ Branch 0 taken 266260267 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 266260267 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 266260267 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 266260267 times.
✗ Branch 9 not taken.
✗ Branch 10 not taken.
|
532520534 | if ((scvf_.hasHalfParallelNeighbor(localSubFaceIdx) || scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) && dirichletParallelNeighbor_(localSubFaceIdx)) |
| 364 | { | ||
| 365 | ✗ | const Scalar boundaryVelocity = getParallelVelocityFromCorner_(localSubFaceIdx); | |
| 366 | ✗ | const Scalar boundaryMomentum = boundaryVelocity*outsideDensity; | |
| 367 | ✗ | return {boundaryMomentum, boundaryMomentum}; | |
| 368 | } | ||
| 369 |
4/4✓ Branch 0 taken 12564483 times.
✓ Branch 1 taken 253695784 times.
✓ Branch 2 taken 12564483 times.
✓ Branch 3 taken 253695784 times.
|
532520534 | else if (!scvf_.hasParallelNeighbor(localSubFaceIdx, 0)) |
| 370 | { | ||
| 371 | 12564483 | const Scalar boundaryVelocity = getParallelVelocityFromBoundary_(element_, scvf_, faceVars_, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); | |
| 372 | 12564483 | const Scalar boundaryMomentum = boundaryVelocity*insideDensity; | |
| 373 | 12564483 | return {boundaryMomentum, boundaryMomentum}; | |
| 374 | } | ||
| 375 | |||
| 376 |
2/2✓ Branch 0 taken 124632644 times.
✓ Branch 1 taken 129063140 times.
|
253695784 | const Scalar momentumParallel = faceVars_.velocityParallel(localSubFaceIdx, 0)*outsideDensity; |
| 377 | 253695784 | const Scalar momentumSelf = faceVars_.velocitySelf()*insideDensity; | |
| 378 |
2/2✓ Branch 0 taken 124632644 times.
✓ Branch 1 taken 129063140 times.
|
253695784 | if (selfIsUpstream) |
| 379 | 124632644 | return {momentumParallel, momentumSelf}; | |
| 380 | else | ||
| 381 | 129063140 | return {momentumSelf, momentumParallel}; | |
| 382 | } | ||
| 383 | |||
| 384 | /*! | ||
| 385 | * \brief Returns an array of distances needed for non-uniform higher order upwind schemes | ||
| 386 | * computes lateral distances {upstream to downstream distance, up-upstream to upstream distance, downstream staggered cell size} | ||
| 387 | */ | ||
| 388 | 3849957 | std::array<Scalar, 3> getLateralDistances_(const int localSubFaceIdx, const bool selfIsUpstream) const | |
| 389 | { | ||
| 390 | static_assert(useHigherOrder, "Should only be reached if higher order methods are enabled"); | ||
| 391 |
2/2✓ Branch 0 taken 1678432 times.
✓ Branch 1 taken 2171525 times.
|
3849957 | if (selfIsUpstream) |
| 392 | { | ||
| 393 | // The local index of the faces that is opposite to localSubFaceIdx | ||
| 394 |
2/2✓ Branch 0 taken 1203155 times.
✓ Branch 1 taken 475277 times.
|
1678432 | const int oppositeSubFaceIdx = localSubFaceIdx % 2 ? localSubFaceIdx - 1 : localSubFaceIdx + 1; |
| 395 | |||
| 396 | std::array<Scalar, 3> distances; | ||
| 397 |
4/4✓ Branch 0 taken 21097 times.
✓ Branch 1 taken 1657335 times.
✓ Branch 2 taken 21097 times.
✓ Branch 3 taken 1657335 times.
|
3356864 | distances[0] = scvf_.parallelDofsDistance(localSubFaceIdx, 0); |
| 398 |
4/4✓ Branch 0 taken 21097 times.
✓ Branch 1 taken 1657335 times.
✓ Branch 2 taken 21097 times.
✓ Branch 3 taken 1657335 times.
|
3356864 | distances[1] = scvf_.parallelDofsDistance(oppositeSubFaceIdx, 0); |
| 399 |
4/4✓ Branch 0 taken 21097 times.
✓ Branch 1 taken 1657335 times.
✓ Branch 2 taken 21097 times.
✓ Branch 3 taken 1657335 times.
|
3356864 | if (scvf_.hasParallelNeighbor(localSubFaceIdx, 0)) |
| 400 | distances[2] = scvf_.pairData(localSubFaceIdx).parallelCellWidths[0]; | ||
| 401 | else | ||
| 402 | 42194 | distances[2] = scvf_.area() / 2.0; | |
| 403 | 1678432 | return distances; | |
| 404 | } | ||
| 405 | else | ||
| 406 | { | ||
| 407 | std::array<Scalar, 3> distances; | ||
| 408 | 4343050 | distances[0] = scvf_.parallelDofsDistance(localSubFaceIdx, 0); | |
| 409 | 4343050 | distances[1] = scvf_.parallelDofsDistance(localSubFaceIdx, 1); | |
| 410 | 6514575 | distances[2] = scvf_.pairData(localSubFaceIdx).parallelCellWidths[0]; | |
| 411 | 2171525 | return distances; | |
| 412 | } | ||
| 413 | } | ||
| 414 | |||
| 415 | /*! | ||
| 416 | * \brief Return the outer parallel velocity for normal faces that are on the boundary and therefore have no neighbor. | ||
| 417 | * Calls the problem to retrieve a fixed value set on the boundary. | ||
| 418 | */ | ||
| 419 | 12677235 | Scalar getParallelVelocityFromBoundary_(const Element& element, | |
| 420 | const SubControlVolumeFace& scvf, | ||
| 421 | const FaceVariables& faceVars, | ||
| 422 | const std::optional<BoundaryTypes>& currentScvfBoundaryTypes, | ||
| 423 | const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes, | ||
| 424 | const int localSubFaceIdx) const | ||
| 425 | { | ||
| 426 | // If there is a Dirichlet condition for the pressure we assume zero gradient for the velocity, | ||
| 427 | // so the velocity at the boundary equal to that on the scvf. | ||
| 428 |
9/14✓ Branch 0 taken 12677235 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 12677235 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 12677235 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 12677235 times.
✗ Branch 7 not taken.
✓ Branch 8 taken 12677235 times.
✗ Branch 9 not taken.
✓ Branch 10 taken 8399711 times.
✓ Branch 11 taken 4277524 times.
✓ Branch 12 taken 8399711 times.
✓ Branch 13 taken 4277524 times.
|
25354470 | const bool useZeroGradient = lateralFaceBoundaryTypes && (lateralFaceBoundaryTypes->isSymmetry() || lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx)); |
| 429 | if (useZeroGradient) | ||
| 430 | 4277524 | return faceVars.velocitySelf(); | |
| 431 | |||
| 432 |
8/10✓ Branch 0 taken 8399711 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8399711 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 8287803 times.
✓ Branch 5 taken 111908 times.
✓ Branch 6 taken 8287803 times.
✓ Branch 7 taken 111908 times.
✓ Branch 8 taken 8287803 times.
✓ Branch 9 taken 111908 times.
|
16799422 | const bool lateralFaceHasBJS = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex())); |
| 433 | if (lateralFaceHasBJS) | ||
| 434 | 111908 | return VelocityGradients::beaversJosephVelocityAtLateralScvf(elemVolVars_.gridVolVars().problem(), element, fvGeometry_, scvf, faceVars, | |
| 435 | 111908 | currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); | |
| 436 | |||
| 437 |
4/8✓ Branch 0 taken 8287803 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8287803 times.
✗ Branch 3 not taken.
✗ Branch 4 not taken.
✓ Branch 5 taken 8287803 times.
✗ Branch 6 not taken.
✓ Branch 7 taken 8287803 times.
|
16575606 | const bool lateralFaceHasDirichletVelocity = lateralFaceBoundaryTypes && lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex())); |
| 438 | if (lateralFaceHasDirichletVelocity) | ||
| 439 | { | ||
| 440 | // ________________ | ||
| 441 | // ---------------o || frontal face of staggered half-control-volume | ||
| 442 | // | || # current scvf # ghostFace of interest, lies on boundary | ||
| 443 | // | || # x dof position | ||
| 444 | // | || x~~~~> vel.Self -- element boundaries | ||
| 445 | // | || # __ domain boundary | ||
| 446 | // | || # o position at which the boundary conditions will be evaluated | ||
| 447 | // ---------------# | ||
| 448 | |||
| 449 | 24863409 | const auto& lateralFace = fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localSubFaceIdx).localLateralFaceIdx); | |
| 450 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
24863409 | const auto ghostFace = makeStaggeredBoundaryFace(lateralFace, scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter); |
| 451 | 8287803 | const auto& problem = elemVolVars_.gridVolVars().problem(); | |
| 452 |
3/7✗ Branch 0 not taken.
✓ Branch 1 taken 8287803 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 8287803 times.
✗ Branch 4 not taken.
✓ Branch 5 taken 120800 times.
✗ Branch 6 not taken.
|
8287803 | return problem.dirichlet(element, ghostFace)[Indices::velocity(scvf.directionIndex())]; |
| 453 | } | ||
| 454 | |||
| 455 | // Neumann conditions are not well implemented | ||
| 456 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Something went wrong with the boundary conditions for the momentum equations at global position " | |
| 457 | << fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localSubFaceIdx).localLateralFaceIdx).center()); | ||
| 458 | } | ||
| 459 | |||
| 460 | /*! | ||
| 461 | * \brief Return a velocity value from a boundary for which the boundary conditions have to be checked. | ||
| 462 | */ | ||
| 463 | 91655 | Scalar getParallelVelocityFromOppositeBoundary_(const Element& element, | |
| 464 | const SubControlVolumeFace& scvf, | ||
| 465 | const FaceVariables& faceVars, | ||
| 466 | const std::optional<BoundaryTypes>& currentScvfBoundaryTypes, | ||
| 467 | const int localOppositeSubFaceIdx) const | ||
| 468 | { | ||
| 469 | // A ghost subface at the boundary is created, featuring the location of the sub face's center | ||
| 470 | 274965 | const auto& lateralOppositeScvf = fvGeometry_.scvf(scvf.insideScvIdx(), scvf.pairData(localOppositeSubFaceIdx).localLateralFaceIdx); | |
| 471 | 274965 | GlobalPosition center = scvf.pairData(localOppositeSubFaceIdx).lateralStaggeredFaceCenter + lateralOppositeScvf.center(); | |
| 472 | 91655 | center *= 0.5; | |
| 473 | |||
| 474 | // lateral face # lateralFace currently being assembled | ||
| 475 | // --------######## || frontal face of staggered half-control-volume | ||
| 476 | // | || | current scvf % lateralOppositeBoundaryFace of interest, lies on boundary | ||
| 477 | // | || | x dof position | ||
| 478 | // | || x~~~~> vel.Self -- element boundaries | ||
| 479 | // | || | __ domain boundary | ||
| 480 | // | || | o position at which the boundary conditions will be evaluated | ||
| 481 | // --------%%%c%%%o | ||
| 482 | // ________________ c center of opposite boundary face | ||
| 483 | |||
| 484 | |||
| 485 | // Get the boundary types of the lateral opposite boundary face | ||
| 486 | 91655 | const auto& problem = elemVolVars_.gridVolVars().problem(); | |
| 487 | 91655 | const auto lateralOppositeBoundaryFace = makeStaggeredBoundaryFace(lateralOppositeScvf, center); | |
| 488 | 91655 | const auto lateralOppositeFaceBoundaryTypes = problem.boundaryTypes(element, lateralOppositeBoundaryFace); | |
| 489 |
3/8✓ Branch 1 taken 91655 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 91655 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 91655 times.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 9 not taken.
|
274965 | return getParallelVelocityFromBoundary_(element, scvf, faceVars, currentScvfBoundaryTypes, lateralOppositeFaceBoundaryTypes, localOppositeSubFaceIdx); |
| 490 | } | ||
| 491 | |||
| 492 | /*! | ||
| 493 | * \brief Returns the boundary subcontrolvolumeface in a corner geometry. | ||
| 494 | * | ||
| 495 | * ------------ | ||
| 496 | * | xxxx o | ||
| 497 | * | xxxx o | ||
| 498 | * | xxxx o | ||
| 499 | * ------------bbbbbbbbbbb | ||
| 500 | * | yyyy o | | ||
| 501 | * | yyyy o | | ||
| 502 | * | yyyy o | | ||
| 503 | * ----------------------- | ||
| 504 | * | ||
| 505 | * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and | | ||
| 506 | * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one | ||
| 507 | * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding | ||
| 508 | * half-control volumes. In both cases, the returned boundaryScvf is the one marked by b. It needs to be the | ||
| 509 | * same boundaryScvf returned for the sake of flux continuity. | ||
| 510 | */ | ||
| 511 | 10696 | const SubControlVolumeFace& boundaryScvf_(const int localSubFaceIdx) const | |
| 512 | { | ||
| 513 |
2/4✓ Branch 0 taken 10696 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10696 times.
✗ Branch 3 not taken.
|
21392 | if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx)) |
| 514 | { | ||
| 515 | 42784 | return fvGeometry_.scvf(scvf_.outsideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx); | |
| 516 | } | ||
| 517 | ✗ | else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) | |
| 518 | { | ||
| 519 | /* | ||
| 520 | * ------------ | ||
| 521 | * | xxxxxxxx o | ||
| 522 | * | xxxxxxxx o | ||
| 523 | * | xxxxxxxx o | ||
| 524 | * lllllllllll-bbbbbbbbbbb | ||
| 525 | * | yyyyyyyy p | | ||
| 526 | * | yyyyyyyy p | | ||
| 527 | * | yyyyyyyy p | | ||
| 528 | * ----------------------- | ||
| 529 | * | ||
| 530 | * o: scvf_, l: lateralFace, p: parallelFace, b: returned scvf, x: scvf_ inside scv, y: lateralFace | ||
| 531 | * outside scv | ||
| 532 | */ | ||
| 533 | ✗ | const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx); | |
| 534 | ✗ | const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx()); | |
| 535 | |||
| 536 | ✗ | const auto& localLateralIdx = scvf_.pairData(localSubFaceIdx).localLateralFaceIdx; | |
| 537 | ✗ | const auto& localLateralOppositeIdx = (localLateralIdx % 2) ? (localLateralIdx - 1) : (localLateralIdx + 1); | |
| 538 | |||
| 539 | ✗ | return fvGeometry_.scvf(parallelFace.outsideScvIdx(), localLateralOppositeIdx); | |
| 540 | } | ||
| 541 | else | ||
| 542 | { | ||
| 543 | ✗ | DUNE_THROW(Dune::InvalidStateException, "The function boundaryScvf_ should only be called when hasHalfParallelNeighbor or hasCornerParallelNeighbor."); | |
| 544 | } | ||
| 545 | } | ||
| 546 | |||
| 547 | /*! | ||
| 548 | * \brief Gets the boundary element in a corner geometry. | ||
| 549 | * | ||
| 550 | * ------------ | ||
| 551 | * | xxxx o | ||
| 552 | * | xxxx o | ||
| 553 | * | xxxx o | ||
| 554 | * ----------------------- | ||
| 555 | * | yyyy o bbbbbbbb | | ||
| 556 | * | yyyy o bbbbbbbb | | ||
| 557 | * | yyyy o bbbbbbbb | | ||
| 558 | * ----------------------- | ||
| 559 | * | ||
| 560 | * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and | | ||
| 561 | * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one | ||
| 562 | * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding | ||
| 563 | * half-control volumes. In both cases, the returned boundaryElement is the one marked by b. It needs to be | ||
| 564 | * the same boundaryScvf returned for the sake of flux continuity. | ||
| 565 | */ | ||
| 566 | 10696 | Element boundaryElement_(const int localSubFaceIdx) const | |
| 567 | { | ||
| 568 |
2/4✓ Branch 0 taken 10696 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 10696 times.
✗ Branch 3 not taken.
|
21392 | if (scvf_.hasHalfParallelNeighbor(localSubFaceIdx)) |
| 569 | { | ||
| 570 | 21392 | return fvGeometry_.gridGeometry().element(scvf_.outsideScvIdx()); | |
| 571 | } | ||
| 572 | ✗ | else if (scvf_.hasCornerParallelNeighbor(localSubFaceIdx)) | |
| 573 | { | ||
| 574 | ✗ | const SubControlVolumeFace& lateralFace = fvGeometry_.scvf(scvf_.insideScvIdx(), scvf_.pairData(localSubFaceIdx).localLateralFaceIdx); | |
| 575 | ✗ | const SubControlVolumeFace& parallelFace = fvGeometry_.scvf(lateralFace.outsideScvIdx(), scvf_.localFaceIdx()); | |
| 576 | |||
| 577 | ✗ | return fvGeometry_.gridGeometry().element(parallelFace.outsideScvIdx()); | |
| 578 | } | ||
| 579 | else | ||
| 580 | { | ||
| 581 | ✗ | DUNE_THROW(Dune::InvalidStateException, "When entering boundaryElement_ scvf_ should have either hasHalfParallelNeighbor or hasCornerParallelNeighbor true. Not the case here."); | |
| 582 | } | ||
| 583 | } | ||
| 584 | |||
| 585 | /*! | ||
| 586 | * \brief Sets the bools hasDirichletCornerParallelNeighbor and hasDirichletHalfParallelNeighbor. | ||
| 587 | * | ||
| 588 | * ------------ | ||
| 589 | * | xxxx o | ||
| 590 | * | xxxx o | ||
| 591 | * | xxxx o | ||
| 592 | * ------------bbbbbbbbbbb | ||
| 593 | * | yyyy o | | ||
| 594 | * | yyyy o boundary | | ||
| 595 | * | yyyy o element | | ||
| 596 | * ----------------------- | ||
| 597 | * | ||
| 598 | * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and | | ||
| 599 | * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one | ||
| 600 | * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding | ||
| 601 | * half-control volumes. In both cases, we check if the face bbb, part of the edge of element boundaryElement, | ||
| 602 | * is a Dirichlet boundary. | ||
| 603 | */ | ||
| 604 | 5348 | bool dirichletParallelNeighbor_(const int localSubFaceIdx) const | |
| 605 | { | ||
| 606 | 5348 | const auto& problem = elemVolVars_.gridVolVars().problem(); | |
| 607 | 5348 | const Element& boundaryElement = boundaryElement_(localSubFaceIdx); | |
| 608 | 5348 | const SubControlVolumeFace& boundaryScvf = boundaryScvf_(localSubFaceIdx); | |
| 609 | |||
| 610 |
1/2✗ Branch 1 not taken.
✓ Branch 2 taken 5348 times.
|
5348 | return problem.boundaryTypes(boundaryElement, boundaryScvf).isDirichlet(Indices::velocity(scvf_.directionIndex())); |
| 611 | } | ||
| 612 | |||
| 613 | /*! | ||
| 614 | * \brief Gets the parallel velocity from a corner geometry. | ||
| 615 | * | ||
| 616 | * ------------ | ||
| 617 | * | xxxx o | ||
| 618 | * | xxxx o | ||
| 619 | * | xxxx o | ||
| 620 | * -----------*----------- | ||
| 621 | * | yyyy o | | ||
| 622 | * | yyyy o | | ||
| 623 | * | yyyy o | | ||
| 624 | * ----------------------- | ||
| 625 | * | ||
| 626 | * This function will be entered in such a corner geometry (there is no cell in the upper right, --- and | | ||
| 627 | * stand for the grid cells). The scvf_ will be one of the two ones denoted by o (upper one | ||
| 628 | * hasCornerParallelNeighbor, lower one hasHalfParallelNeighbor). x and y are the two possible corresponding | ||
| 629 | * half-control volumes. In both cases, the returned velocity is situated in the corner (*). | ||
| 630 | */ | ||
| 631 | 5348 | Scalar getParallelVelocityFromCorner_(const int localSubFaceIdx) const | |
| 632 | { | ||
| 633 | 5348 | const auto& problem = elemVolVars_.gridVolVars().problem(); | |
| 634 | 5348 | const Element& boundaryElement = boundaryElement_(localSubFaceIdx); | |
| 635 | 5348 | const SubControlVolumeFace& boundaryScvf = boundaryScvf_(localSubFaceIdx); | |
| 636 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
10696 | const auto ghostFace = makeStaggeredBoundaryFace(boundaryScvf, scvf_.pairData(localSubFaceIdx).lateralStaggeredFaceCenter); |
| 637 | |||
| 638 |
2/7✗ Branch 0 not taken.
✓ Branch 1 taken 5348 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 5348 times.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
|
5348 | return problem.dirichlet(boundaryElement, ghostFace)[Indices::velocity(scvf_.directionIndex())]; |
| 639 | } | ||
| 640 | |||
| 641 | const Element& element_; | ||
| 642 | const FVElementGeometry& fvGeometry_; | ||
| 643 | const SubControlVolumeFace& scvf_; | ||
| 644 | const ElementFaceVariables& elemFaceVars_; | ||
| 645 | const FaceVariables& faceVars_; | ||
| 646 | const ElementVolumeVariables& elemVolVars_; | ||
| 647 | const UpwindScheme& upwindScheme_; | ||
| 648 | }; | ||
| 649 | |||
| 650 | } // end namespace Dumux | ||
| 651 | |||
| 652 | #endif | ||
| 653 |