| 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::NavierStokesFluxVariablesImpl | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH | ||
| 13 | #define DUMUX_NAVIERSTOKES_STAGGERED_FLUXVARIABLES_HH | ||
| 14 | |||
| 15 | #include <array> | ||
| 16 | #include <optional> | ||
| 17 | #include <type_traits> | ||
| 18 | |||
| 19 | #include <dumux/common/math.hh> | ||
| 20 | #include <dumux/common/exceptions.hh> | ||
| 21 | #include <dumux/common/parameters.hh> | ||
| 22 | #include <dumux/common/properties.hh> | ||
| 23 | #include <dumux/common/typetraits/problem.hh> | ||
| 24 | |||
| 25 | #include <dumux/flux/fluxvariablesbase.hh> | ||
| 26 | #include <dumux/discretization/method.hh> | ||
| 27 | #include <dumux/discretization/extrusion.hh> | ||
| 28 | |||
| 29 | #include "staggeredupwindhelper.hh" | ||
| 30 | #include "velocitygradients.hh" | ||
| 31 | |||
| 32 | namespace Dumux { | ||
| 33 | |||
| 34 | /*! | ||
| 35 | * \ingroup NavierStokesModel | ||
| 36 | * \brief The flux variables class for the Navier-Stokes model using the staggered grid discretization. | ||
| 37 | */ | ||
| 38 | |||
| 39 | // forward declaration | ||
| 40 | template<class TypeTag, class DiscretizationMethod> | ||
| 41 | class NavierStokesFluxVariablesImpl; | ||
| 42 | |||
| 43 | template<class TypeTag> | ||
| 44 | class NavierStokesFluxVariablesImpl<TypeTag, DiscretizationMethods::Staggered> | ||
| 45 | : public FluxVariablesBase<GetPropType<TypeTag, Properties::Problem>, | ||
| 46 | typename GetPropType<TypeTag, Properties::GridGeometry>::LocalView, | ||
| 47 | typename GetPropType<TypeTag, Properties::GridVolumeVariables>::LocalView, | ||
| 48 | typename GetPropType<TypeTag, Properties::GridFluxVariablesCache>::LocalView> | ||
| 49 | { | ||
| 50 | using GridVariables = GetPropType<TypeTag, Properties::GridVariables>; | ||
| 51 | |||
| 52 | using GridVolumeVariables = typename GridVariables::GridVolumeVariables; | ||
| 53 | using ElementVolumeVariables = typename GridVolumeVariables::LocalView; | ||
| 54 | using VolumeVariables = typename GridVolumeVariables::VolumeVariables; | ||
| 55 | |||
| 56 | using GridFluxVariablesCache = typename GridVariables::GridFluxVariablesCache; | ||
| 57 | using FluxVariablesCache = typename GridFluxVariablesCache::FluxVariablesCache; | ||
| 58 | |||
| 59 | using GridFaceVariables = typename GridVariables::GridFaceVariables; | ||
| 60 | using ElementFaceVariables = typename GridFaceVariables::LocalView; | ||
| 61 | using FaceVariables = typename GridFaceVariables::FaceVariables; | ||
| 62 | |||
| 63 | using Problem = GetPropType<TypeTag, Properties::Problem>; | ||
| 64 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 65 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 66 | using GridView = typename GridGeometry::GridView; | ||
| 67 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 68 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 69 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 70 | using Indices = typename ModelTraits::Indices; | ||
| 71 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 72 | using FaceFrontalSubControlVolumeFace = typename GridGeometry::Traits::FaceFrontalSubControlVolumeFace; | ||
| 73 | using FaceLateralSubControlVolumeFace = typename GridGeometry::Traits::FaceLateralSubControlVolumeFace; | ||
| 74 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 75 | using CellCenterPrimaryVariables = GetPropType<TypeTag, Properties::CellCenterPrimaryVariables>; | ||
| 76 | using FacePrimaryVariables = GetPropType<TypeTag, Properties::FacePrimaryVariables>; | ||
| 77 | using BoundaryTypes = typename ProblemTraits<Problem>::BoundaryTypes; | ||
| 78 | using VelocityGradients = StaggeredVelocityGradients<Scalar, GridGeometry, BoundaryTypes, Indices>; | ||
| 79 | |||
| 80 | static constexpr bool normalizePressure = getPropValue<TypeTag, Properties::NormalizePressure>(); | ||
| 81 | |||
| 82 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 83 | |||
| 84 | static constexpr auto cellCenterIdx = GridGeometry::cellCenterIdx(); | ||
| 85 | static constexpr auto faceIdx = GridGeometry::faceIdx(); | ||
| 86 | |||
| 87 | static constexpr int upwindSchemeOrder = getPropValue<TypeTag, Properties::UpwindSchemeOrder>(); | ||
| 88 | static constexpr bool useHigherOrder = upwindSchemeOrder > 1; | ||
| 89 | |||
| 90 | public: | ||
| 91 | |||
| 92 | using HeatConductionType = GetPropType<TypeTag, Properties::HeatConductionType>; | ||
| 93 | |||
| 94 | |||
| 95 | /*! | ||
| 96 | * \brief Returns the advective flux over a sub control volume face. | ||
| 97 | * \param problem The object specifying the problem which ought to be simulated | ||
| 98 | * \param fvGeometry The element-centered finite volume geometry view | ||
| 99 | * \param elemVolVars All volume variables for the element | ||
| 100 | * \param elemFaceVars The face variables | ||
| 101 | * \param scvf The sub control volume face | ||
| 102 | * \param upwindTerm The uwind term (i.e. the advectively transported quantity) | ||
| 103 | */ | ||
| 104 | template<class UpwindTerm> | ||
| 105 | 1048814936 | static Scalar advectiveFluxForCellCenter(const Problem& problem, | |
| 106 | const FVElementGeometry& fvGeometry, | ||
| 107 | const ElementVolumeVariables& elemVolVars, | ||
| 108 | const ElementFaceVariables& elemFaceVars, | ||
| 109 | const SubControlVolumeFace &scvf, | ||
| 110 | UpwindTerm upwindTerm) | ||
| 111 | { | ||
| 112 |
2/2✓ Branch 0 taken 115 times.
✓ Branch 1 taken 542297701 times.
|
1048814936 | const Scalar velocity = elemFaceVars[scvf].velocitySelf(); |
| 113 |
2/2✓ Branch 0 taken 115 times.
✓ Branch 1 taken 542297701 times.
|
1048814936 | const bool insideIsUpstream = scvf.directionSign() == sign(velocity); |
| 114 |
5/8✓ Branch 0 taken 115 times.
✓ Branch 1 taken 542297701 times.
✓ Branch 3 taken 115 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 115 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 115 times.
✗ Branch 10 not taken.
|
1048814936 | static const Scalar upwindWeight = getParamFromGroup<Scalar>(problem.paramGroup(), "Flux.UpwindWeight"); |
| 115 | |||
| 116 | 2097629872 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 117 | 2097629872 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 118 | |||
| 119 |
2/2✓ Branch 0 taken 278244496 times.
✓ Branch 1 taken 264053320 times.
|
1048814936 | const auto& upstreamVolVars = insideIsUpstream ? insideVolVars : outsideVolVars; |
| 120 |
2/2✓ Branch 0 taken 278244496 times.
✓ Branch 1 taken 264053320 times.
|
1048814936 | const auto& downstreamVolVars = insideIsUpstream ? outsideVolVars : insideVolVars; |
| 121 | |||
| 122 | 1648262756 | const Scalar flux = (upwindWeight * upwindTerm(upstreamVolVars) + | |
| 123 | 1248630876 | (1.0 - upwindWeight) * upwindTerm(downstreamVolVars)) | |
| 124 | 1248630876 | * velocity * Extrusion::area(fvGeometry, scvf) * scvf.directionSign(); | |
| 125 | |||
| 126 | 1048814936 | return flux * extrusionFactor_(elemVolVars, scvf); | |
| 127 | } | ||
| 128 | |||
| 129 | /*! | ||
| 130 | * \brief Computes the flux for the cell center residual (mass balance). | ||
| 131 | * | ||
| 132 | * \verbatim | ||
| 133 | * scvf | ||
| 134 | * ---------------- | ||
| 135 | * | # current scvf # scvf over which fluxes are calculated | ||
| 136 | * | # | ||
| 137 | * | x #~~~~> vel.Self x dof position | ||
| 138 | * | # | ||
| 139 | * scvf | # -- element | ||
| 140 | * ---------------- | ||
| 141 | * scvf | ||
| 142 | * \endverbatim | ||
| 143 | */ | ||
| 144 | ✗ | CellCenterPrimaryVariables computeMassFlux(const Problem& problem, | |
| 145 | const Element& element, | ||
| 146 | const FVElementGeometry& fvGeometry, | ||
| 147 | const ElementVolumeVariables& elemVolVars, | ||
| 148 | const ElementFaceVariables& elemFaceVars, | ||
| 149 | const SubControlVolumeFace& scvf, | ||
| 150 | const FluxVariablesCache& fluxVarsCache) | ||
| 151 | { | ||
| 152 | // The advectively transported quantity (i.e density for a single-phase system). | ||
| 153 | 322239468 | auto upwindTerm = [](const auto& volVars) { return volVars.density(); }; | |
| 154 | |||
| 155 | // Call the generic flux function. | ||
| 156 |
1/8✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 138882 times.
✗ Branch 5 not taken.
✗ Branch 6 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
|
100386724 | CellCenterPrimaryVariables result(0.0); |
| 157 |
5/10✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 3 not taken.
✓ Branch 4 taken 138882 times.
✓ Branch 5 taken 140546 times.
✓ Branch 6 taken 41192 times.
✓ Branch 7 taken 97690 times.
✗ Branch 8 not taken.
✓ Branch 9 taken 8126 times.
✗ Branch 10 not taken.
|
100386724 | result[Indices::conti0EqIdx - ModelTraits::dim()] = advectiveFluxForCellCenter(problem, fvGeometry, elemVolVars, elemFaceVars, scvf, upwindTerm); |
| 158 | |||
| 159 | ✗ | return result; | |
| 160 | } | ||
| 161 | |||
| 162 | /*! | ||
| 163 | * \brief Returns the momentum flux over all staggered faces. | ||
| 164 | */ | ||
| 165 | 47252532 | FacePrimaryVariables computeMomentumFlux(const Problem& problem, | |
| 166 | const Element& element, | ||
| 167 | const SubControlVolumeFace& scvf, | ||
| 168 | const FVElementGeometry& fvGeometry, | ||
| 169 | const ElementVolumeVariables& elemVolVars, | ||
| 170 | const ElementFaceVariables& elemFaceVars, | ||
| 171 | const GridFluxVariablesCache& gridFluxVarsCache) | ||
| 172 | { | ||
| 173 | 47252532 | return computeFrontalMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache) + | |
| 174 | 47252532 | computeLateralMomentumFlux(problem, element, scvf, fvGeometry, elemVolVars, elemFaceVars, gridFluxVarsCache); | |
| 175 | } | ||
| 176 | |||
| 177 | /*! | ||
| 178 | * \brief Returns the frontal part of the momentum flux. | ||
| 179 | * This treats the flux over the staggered face at the center of an element, | ||
| 180 | * parallel to the current scvf where the velocity dof of interest lives. | ||
| 181 | * | ||
| 182 | * \verbatim | ||
| 183 | * scvf | ||
| 184 | * ---------======= == and # staggered half-control-volume | ||
| 185 | * | # | current scvf | ||
| 186 | * | # | # staggered face over which fluxes are calculated | ||
| 187 | * vel.Opp <~~| #~~> x~~~~> vel.Self | ||
| 188 | * | # | x dof position | ||
| 189 | * scvf | # | | ||
| 190 | * --------======== -- element | ||
| 191 | * scvf | ||
| 192 | * \endverbatim | ||
| 193 | */ | ||
| 194 | 145571324 | FacePrimaryVariables computeFrontalMomentumFlux(const Problem& problem, | |
| 195 | const Element& element, | ||
| 196 | const SubControlVolumeFace& scvf, | ||
| 197 | const FVElementGeometry& fvGeometry, | ||
| 198 | const ElementVolumeVariables& elemVolVars, | ||
| 199 | const ElementFaceVariables& elemFaceVars, | ||
| 200 | const GridFluxVariablesCache& gridFluxVarsCache) | ||
| 201 | { | ||
| 202 |
2/2✓ Branch 0 taken 135725860 times.
✓ Branch 1 taken 9845464 times.
|
145571324 | FacePrimaryVariables frontalFlux(0.0); |
| 203 | |||
| 204 | // The velocities of the dof at interest and the one of the opposite scvf. | ||
| 205 |
2/2✓ Branch 0 taken 135725860 times.
✓ Branch 1 taken 9845464 times.
|
145571324 | const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf(); |
| 206 |
2/2✓ Branch 0 taken 135725860 times.
✓ Branch 1 taken 9845464 times.
|
145571324 | const Scalar velocityOpposite = elemFaceVars[scvf].velocityOpposite(); |
| 207 |
2/2✓ Branch 0 taken 135725860 times.
✓ Branch 1 taken 9845464 times.
|
145571324 | const auto& faceVars = elemFaceVars[scvf]; |
| 208 | |||
| 209 | // Advective flux. | ||
| 210 |
2/2✓ Branch 0 taken 135725860 times.
✓ Branch 1 taken 9845464 times.
|
145571324 | if (problem.enableInertiaTerms()) |
| 211 | { | ||
| 212 | // Get the average velocity at the center of the element (i.e. the location of the staggered face). | ||
| 213 | 135725860 | const Scalar transportingVelocity = (velocitySelf + velocityOpposite) * 0.5; | |
| 214 | 135725860 | const bool selfIsUpstream = scvf.directionSign() != sign(transportingVelocity); | |
| 215 | |||
| 216 | 271451720 | StaggeredUpwindHelper<TypeTag, upwindSchemeOrder> upwindHelper(element, fvGeometry, scvf, elemFaceVars, elemVolVars, gridFluxVarsCache.staggeredUpwindMethods()); | |
| 217 | 271451720 | frontalFlux += upwindHelper.computeUpwindFrontalMomentum(selfIsUpstream) | |
| 218 | 135725860 | * transportingVelocity * -1.0 * scvf.directionSign(); | |
| 219 | } | ||
| 220 | |||
| 221 | // The volume variables within the current element. We only require those (and none of neighboring elements) | ||
| 222 | // because the fluxes are calculated over the staggered face at the center of the element. | ||
| 223 | 291142648 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 224 | |||
| 225 | // Diffusive flux. | ||
| 226 |
2/2✓ Branch 0 taken 51 times.
✓ Branch 1 taken 145571273 times.
|
145571324 | const Scalar velocityGrad_ii = VelocityGradients::velocityGradII(scvf, faceVars) * scvf.directionSign(); |
| 227 | |||
| 228 |
3/4✓ Branch 0 taken 51 times.
✓ Branch 1 taken 145571273 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
|
145571375 | static const bool enableUnsymmetrizedVelocityGradient |
| 229 |
2/4✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 51 times.
✗ Branch 5 not taken.
|
153 | = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false); |
| 230 |
1/2✓ Branch 0 taken 145571324 times.
✗ Branch 1 not taken.
|
145571324 | const Scalar factor = enableUnsymmetrizedVelocityGradient ? 1.0 : 2.0; |
| 231 | 291142648 | frontalFlux -= factor * insideVolVars.effectiveViscosity() * velocityGrad_ii; | |
| 232 | |||
| 233 | // The pressure term. | ||
| 234 | // If specified, the pressure can be normalized using the initial value on the scfv of interest. | ||
| 235 | // The scvf is used to normalize by the same value from the left and right side. | ||
| 236 | // Can potentially help to improve the condition number of the system matrix. | ||
| 237 | 145571324 | const Scalar pressure = normalizePressure ? | |
| 238 | 350308714 | insideVolVars.pressure() - problem.initial(scvf)[Indices::pressureIdx] | |
| 239 | : insideVolVars.pressure(); | ||
| 240 | |||
| 241 | // Account for the orientation of the staggered face's normal outer normal vector | ||
| 242 | // (pointing in opposite direction of the scvf's one). | ||
| 243 | 145571324 | frontalFlux += pressure * -1.0 * scvf.directionSign(); | |
| 244 | |||
| 245 | // Account for the staggered face's area. For rectangular elements, this equals the area of the scvf | ||
| 246 | // our velocity dof of interest lives on but with adjusted centroid | ||
| 247 | 291142648 | const auto& scv = fvGeometry.scv(scvf.insideScvIdx()); | |
| 248 | 291142648 | FaceFrontalSubControlVolumeFace frontalFace(scv.center(), scvf.area()); | |
| 249 | 582285296 | return frontalFlux * Extrusion::area(fvGeometry, frontalFace) * insideVolVars.extrusionFactor(); | |
| 250 | } | ||
| 251 | |||
| 252 | /*! | ||
| 253 | * \brief Returns the momentum flux over the staggered faces | ||
| 254 | * perpendicular to the scvf where the velocity dof of interest | ||
| 255 | * lives (coinciding with the element's scvfs). | ||
| 256 | * | ||
| 257 | * \verbatim | ||
| 258 | * scvf | ||
| 259 | * ---------####### || and # staggered half-control-volume | ||
| 260 | * | || | current scvf | ||
| 261 | * | || | # normal staggered sub faces over which fluxes are calculated | ||
| 262 | * | || x~~~~> vel.Self | ||
| 263 | * | || | x dof position | ||
| 264 | * scvf | || | | ||
| 265 | * --------######## -- element | ||
| 266 | * scvf | ||
| 267 | * \endverbatim | ||
| 268 | */ | ||
| 269 | 145571324 | FacePrimaryVariables computeLateralMomentumFlux(const Problem& problem, | |
| 270 | const Element& element, | ||
| 271 | const SubControlVolumeFace& scvf, | ||
| 272 | const FVElementGeometry& fvGeometry, | ||
| 273 | const ElementVolumeVariables& elemVolVars, | ||
| 274 | const ElementFaceVariables& elemFaceVars, | ||
| 275 | const GridFluxVariablesCache& gridFluxVarsCache) | ||
| 276 | { | ||
| 277 |
2/2✓ Branch 0 taken 2808642 times.
✓ Branch 1 taken 142762682 times.
|
145571324 | FacePrimaryVariables lateralFlux(0.0); |
| 278 |
2/2✓ Branch 0 taken 2808642 times.
✓ Branch 1 taken 142762682 times.
|
145571324 | const auto& faceVars = elemFaceVars[scvf]; |
| 279 |
2/2✓ Branch 0 taken 2808642 times.
✓ Branch 1 taken 142762682 times.
|
145571324 | const std::size_t numSubFaces = scvf.pairData().size(); |
| 280 | |||
| 281 | // If the current scvf is on a boundary, check if there is a Neumann BC for the stress in tangential direction. | ||
| 282 | // Create a boundaryTypes object (will be empty if not at a boundary). | ||
| 283 | 145571324 | std::optional<BoundaryTypes> currentScvfBoundaryTypes; | |
| 284 |
2/2✓ Branch 0 taken 2808642 times.
✓ Branch 1 taken 142762682 times.
|
145571324 | if (scvf.boundary()) |
| 285 | 2808642 | currentScvfBoundaryTypes.emplace(problem.boundaryTypes(element, scvf)); | |
| 286 | |||
| 287 | // Account for all sub faces. | ||
| 288 |
2/2✓ Branch 0 taken 291142648 times.
✓ Branch 1 taken 145571324 times.
|
436713972 | for (int localSubFaceIdx = 0; localSubFaceIdx < numSubFaces; ++localSubFaceIdx) |
| 289 | { | ||
| 290 |
2/2✓ Branch 0 taken 14651218 times.
✓ Branch 1 taken 276491430 times.
|
291142648 | const auto eIdx = scvf.insideScvIdx(); |
| 291 | // Get the face normal to the face the dof lives on. The staggered sub face coincides with half of this lateral face. | ||
| 292 |
4/4✓ Branch 0 taken 14651218 times.
✓ Branch 1 taken 276491430 times.
✓ Branch 2 taken 14651218 times.
✓ Branch 3 taken 276491430 times.
|
582285296 | const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); |
| 293 | |||
| 294 | // Create a boundaryTypes object (will be empty if not at a boundary). | ||
| 295 | 291142648 | std::optional<BoundaryTypes> lateralFaceBoundaryTypes; | |
| 296 | |||
| 297 | // Check if there is face/element parallel to our face of interest where the dof lives on. If there is no parallel neighbor, | ||
| 298 | // we are on a boundary where we have to check for boundary conditions. | ||
| 299 |
2/2✓ Branch 0 taken 14651218 times.
✓ Branch 1 taken 276491430 times.
|
291142648 | if (lateralFace.boundary()) |
| 300 | { | ||
| 301 | // Retrieve the boundary types that correspond to the center of the lateral scvf. As a convention, we always query | ||
| 302 | // the type of BCs at the center of the element's "actual" lateral scvf (not the face of the staggered control volume). | ||
| 303 | // The value of the BC will be evaluated at the center of the staggered face. | ||
| 304 | // --------T######V || frontal face of staggered half-control-volume | ||
| 305 | // | || | current scvf # lateral staggered face of interest (may lie on a boundary) | ||
| 306 | // | || | x dof position | ||
| 307 | // | || x~~~~> vel.Self -- element boundaries | ||
| 308 | // | || | T position at which the type of boundary conditions will be evaluated | ||
| 309 | // | || | (center of lateral scvf) | ||
| 310 | // ---------------- V position at which the value of the boundary conditions will be evaluated | ||
| 311 | // (center of the staggered lateral face) | ||
| 312 | 14651218 | lateralFaceBoundaryTypes.emplace(problem.boundaryTypes(element, lateralFace)); | |
| 313 | } | ||
| 314 | |||
| 315 | // If the current scvf is on a boundary and if there is a Neumann or Beavers-Joseph-(Saffmann) BC for the stress in tangential direction, | ||
| 316 | // assign this value for the lateral momentum flux. No further calculations are required. We assume that all lateral faces | ||
| 317 | // have the same type of BC (Neumann or Beavers-Joseph-(Saffmann)), but we sample the value at their actual positions. | ||
| 318 |
4/4✓ Branch 0 taken 5617284 times.
✓ Branch 1 taken 285525364 times.
✓ Branch 2 taken 5617284 times.
✓ Branch 3 taken 285525364 times.
|
582285296 | if (currentScvfBoundaryTypes) |
| 319 | { | ||
| 320 | // Handle Neumann BCs. | ||
| 321 |
2/4✓ Branch 0 taken 5617284 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 5617284 times.
✗ Branch 3 not taken.
|
11234568 | if (currentScvfBoundaryTypes->isNeumann(Indices::velocity(lateralFace.directionIndex()))) |
| 322 | { | ||
| 323 | // Get the location of the lateral staggered face's center. | ||
| 324 | ✗ | FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); | |
| 325 | ✗ | const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); | |
| 326 | ✗ | const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralStaggeredFaceCenter); | |
| 327 | ✗ | lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(lateralFace.directionIndex())] | |
| 328 | ✗ | * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf) * lateralFace.directionSign(); | |
| 329 | |||
| 330 | continue; | ||
| 331 | } | ||
| 332 | |||
| 333 | // Treat the edge case when our current scvf has a BJ condition while the lateral face has a Neumann BC. | ||
| 334 | // It is not clear how to evaluate the BJ condition here. | ||
| 335 | // For symmetry reasons, our own scvf should then have the same Neumann flux as the lateral face. | ||
| 336 | // TODO: We should clarify if this is the correct approach. | ||
| 337 |
10/10✓ Branch 0 taken 346528 times.
✓ Branch 1 taken 5270756 times.
✓ Branch 2 taken 346528 times.
✓ Branch 3 taken 5270756 times.
✓ Branch 4 taken 346528 times.
✓ Branch 5 taken 5270756 times.
✓ Branch 6 taken 18242 times.
✓ Branch 7 taken 328286 times.
✓ Branch 8 taken 18242 times.
✓ Branch 9 taken 328286 times.
|
16851852 | if (currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex())) && lateralFaceBoundaryTypes && |
| 338 |
2/4✓ Branch 0 taken 18242 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 18242 times.
✗ Branch 3 not taken.
|
36484 | lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex()))) |
| 339 | { | ||
| 340 | ✗ | FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); | |
| 341 | ✗ | const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); | |
| 342 | ✗ | const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter); | |
| 343 | ✗ | lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())] | |
| 344 | ✗ | * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf) * lateralFace.directionSign(); | |
| 345 | |||
| 346 | continue; | ||
| 347 | } | ||
| 348 | } | ||
| 349 | |||
| 350 | // Check if the lateral face (perpendicular to our current scvf) lies on a boundary. If yes, boundary conditions might need to be treated | ||
| 351 | // and further calculations can be skipped. | ||
| 352 |
2/2✓ Branch 0 taken 14651218 times.
✓ Branch 1 taken 276491430 times.
|
291142648 | if (lateralFace.boundary()) |
| 353 | { | ||
| 354 | // Check if we have a symmetry boundary condition. If yes, the tangential part of the momentum flux can be neglected | ||
| 355 | // and we may skip any further calculations for the given sub face. | ||
| 356 |
6/6✓ Branch 0 taken 13818570 times.
✓ Branch 1 taken 832648 times.
✓ Branch 2 taken 13818570 times.
✓ Branch 3 taken 832648 times.
✓ Branch 4 taken 13818570 times.
✓ Branch 5 taken 832648 times.
|
43953654 | if (lateralFaceBoundaryTypes->isSymmetry()) |
| 357 | continue; | ||
| 358 | |||
| 359 | // Handle Neumann boundary conditions. No further calculations are then required for the given sub face. | ||
| 360 |
2/4✓ Branch 0 taken 13818570 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 13818570 times.
✗ Branch 3 not taken.
|
27637140 | if (lateralFaceBoundaryTypes->isNeumann(Indices::velocity(scvf.directionIndex()))) |
| 361 | { | ||
| 362 | // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location | ||
| 363 | // the staggered faces's center. | ||
| 364 | ✗ | FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); | |
| 365 | ✗ | const auto& lateralStaggeredFaceCenter = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); | |
| 366 | ✗ | const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter); | |
| 367 | ✗ | lateralFlux += problem.neumann(element, fvGeometry, elemVolVars, elemFaceVars, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())] | |
| 368 | ✗ | * elemVolVars[lateralFace.insideScvIdx()].extrusionFactor() * Extrusion::area(fvGeometry, lateralScvf); | |
| 369 | |||
| 370 | continue; | ||
| 371 | } | ||
| 372 | |||
| 373 | // Handle wall-function fluxes (only required for RANS models) | ||
| 374 |
2/2✓ Branch 1 taken 3007923 times.
✓ Branch 2 taken 508848 times.
|
3516771 | if (incorporateWallFunction_(lateralFlux, problem, element, fvGeometry, scvf, elemVolVars, elemFaceVars, localSubFaceIdx)) |
| 375 | continue; | ||
| 376 | } | ||
| 377 | |||
| 378 | // Check the consistency of the boundary conditions, exactly one of the following must be set | ||
| 379 |
4/4✓ Branch 0 taken 13309722 times.
✓ Branch 1 taken 276491430 times.
✓ Branch 2 taken 13309722 times.
✓ Branch 3 taken 276491430 times.
|
579602304 | if (lateralFaceBoundaryTypes) |
| 380 | { | ||
| 381 | 13309722 | std::bitset<3> admittableBcTypes; | |
| 382 |
4/4✓ Branch 0 taken 8904604 times.
✓ Branch 1 taken 4405118 times.
✓ Branch 2 taken 8904604 times.
✓ Branch 3 taken 4405118 times.
|
26619444 | admittableBcTypes.set(0, lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx)); |
| 383 |
4/4✓ Branch 0 taken 4660362 times.
✓ Branch 1 taken 8649360 times.
✓ Branch 2 taken 4660362 times.
✓ Branch 3 taken 8649360 times.
|
31279806 | admittableBcTypes.set(1, lateralFaceBoundaryTypes->isDirichlet(Indices::velocity(scvf.directionIndex()))); |
| 384 |
6/6✓ Branch 0 taken 255244 times.
✓ Branch 1 taken 13054478 times.
✓ Branch 2 taken 255244 times.
✓ Branch 3 taken 13054478 times.
✓ Branch 4 taken 255244 times.
✓ Branch 5 taken 13054478 times.
|
39929166 | admittableBcTypes.set(2, lateralFaceBoundaryTypes->isBeaversJoseph(Indices::velocity(scvf.directionIndex()))); |
| 385 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 13309722 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 13309722 times.
|
26619444 | if (admittableBcTypes.count() != 1) |
| 386 | { | ||
| 387 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Invalid boundary conditions for lateral scvf " | |
| 388 | "for the momentum equations at global position " << lateralStaggeredFaceCenter_(scvf, localSubFaceIdx) | ||
| 389 | << ", current scvf global position " << scvf.center()); | ||
| 390 | } | ||
| 391 | } | ||
| 392 | |||
| 393 | // If none of the above boundary conditions apply for the given sub face, proceed to calculate the tangential momentum flux. | ||
| 394 |
2/2✓ Branch 0 taken 270110224 times.
✓ Branch 1 taken 19690928 times.
|
289801152 | if (problem.enableInertiaTerms()) |
| 395 |
1/2✓ Branch 1 taken 29634048 times.
✗ Branch 2 not taken.
|
270110224 | lateralFlux += computeAdvectivePartOfLateralMomentumFlux_(problem, fvGeometry, element, |
| 396 | scvf, elemVolVars, elemFaceVars, | ||
| 397 | gridFluxVarsCache, | ||
| 398 | currentScvfBoundaryTypes, lateralFaceBoundaryTypes, | ||
| 399 | localSubFaceIdx); | ||
| 400 | |||
| 401 |
1/2✓ Branch 1 taken 49324976 times.
✗ Branch 2 not taken.
|
289801152 | lateralFlux += computeDiffusivePartOfLateralMomentumFlux_(problem, fvGeometry, element, |
| 402 | scvf, elemVolVars, faceVars, | ||
| 403 | currentScvfBoundaryTypes, lateralFaceBoundaryTypes, | ||
| 404 | localSubFaceIdx); | ||
| 405 | } | ||
| 406 | 145571324 | return lateralFlux; | |
| 407 | } | ||
| 408 | |||
| 409 | |||
| 410 | /*! | ||
| 411 | * \brief Returns the momentum flux over an inflow or outflow boundary face. | ||
| 412 | * | ||
| 413 | * \verbatim | ||
| 414 | * scvf // | ||
| 415 | * ---------=======// == and # staggered half-control-volume | ||
| 416 | * | || #// current scvf | ||
| 417 | * | || #// # staggered boundary face over which fluxes are calculated | ||
| 418 | * | || x~~~~> vel.Self | ||
| 419 | * | || #// x dof position | ||
| 420 | * scvf | || #// | ||
| 421 | * --------========// -- element | ||
| 422 | * scvf // | ||
| 423 | * // boundary | ||
| 424 | * \endverbatim | ||
| 425 | */ | ||
| 426 | 2635378 | FacePrimaryVariables inflowOutflowBoundaryFlux(const Problem& problem, | |
| 427 | const SubControlVolumeFace& scvf, | ||
| 428 | const FVElementGeometry& fvGeometry, | ||
| 429 | const ElementVolumeVariables& elemVolVars, | ||
| 430 | const ElementFaceVariables& elemFaceVars) const | ||
| 431 | { | ||
| 432 | 2635378 | FacePrimaryVariables inOrOutflow(0.0); | |
| 433 | 2635378 | const auto& element = fvGeometry.element(); | |
| 434 | 5270756 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 435 | |||
| 436 | // Advective momentum flux. | ||
| 437 |
2/2✓ Branch 0 taken 2564614 times.
✓ Branch 1 taken 70764 times.
|
2635378 | if (problem.enableInertiaTerms()) |
| 438 | { | ||
| 439 | 2564614 | const Scalar velocitySelf = elemFaceVars[scvf].velocitySelf(); | |
| 440 | 5129228 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 441 |
4/4✓ Branch 0 taken 18215 times.
✓ Branch 1 taken 2546399 times.
✓ Branch 2 taken 18215 times.
✓ Branch 3 taken 2546399 times.
|
5129228 | const auto& upVolVars = (scvf.directionSign() == sign(velocitySelf)) ? |
| 442 | insideVolVars : outsideVolVars; | ||
| 443 | |||
| 444 | 5129228 | inOrOutflow += velocitySelf * velocitySelf * upVolVars.density(); | |
| 445 | } | ||
| 446 | |||
| 447 | // Apply a pressure at the boundary. | ||
| 448 | 2635378 | const Scalar boundaryPressure = normalizePressure | |
| 449 | 2709822 | ? (problem.dirichlet(element, scvf)[Indices::pressureIdx] - | |
| 450 | 2638354 | problem.initial(scvf)[Indices::pressureIdx]) | |
| 451 | : problem.dirichlet(element, scvf)[Indices::pressureIdx]; | ||
| 452 | 5270756 | inOrOutflow += boundaryPressure; | |
| 453 | |||
| 454 | // Account for the orientation of the face at the boundary, | ||
| 455 | 5270756 | return inOrOutflow * scvf.directionSign() * Extrusion::area(fvGeometry, scvf) * insideVolVars.extrusionFactor(); | |
| 456 | } | ||
| 457 | |||
| 458 | private: | ||
| 459 | |||
| 460 | /*! | ||
| 461 | * \brief Returns the advective momentum flux over the staggered face perpendicular to the scvf | ||
| 462 | * where the velocity dof of interest lives (coinciding with the element's scvfs). | ||
| 463 | * | ||
| 464 | * \verbatim | ||
| 465 | * ---------------- | ||
| 466 | * | | | ||
| 467 | * | transp. | | ||
| 468 | * | vel. |~~~~> vel.Parallel | ||
| 469 | * | ^ | | ||
| 470 | * | | | | ||
| 471 | * scvf ---------####### || and # staggered half-control-volume | ||
| 472 | * | || | current scvf | ||
| 473 | * | || | # normal staggered faces over which fluxes are calculated | ||
| 474 | * | || x~~~~> vel.Self | ||
| 475 | * | || | x dof position | ||
| 476 | * scvf | || | | ||
| 477 | * ---------####### -- elements | ||
| 478 | * scvf | ||
| 479 | * \endverbatim | ||
| 480 | */ | ||
| 481 | 270110224 | FacePrimaryVariables computeAdvectivePartOfLateralMomentumFlux_(const Problem& problem, | |
| 482 | const FVElementGeometry& fvGeometry, | ||
| 483 | const Element& element, | ||
| 484 | const SubControlVolumeFace& scvf, | ||
| 485 | const ElementVolumeVariables& elemVolVars, | ||
| 486 | const ElementFaceVariables& elemFaceVars, | ||
| 487 | const GridFluxVariablesCache& gridFluxVarsCache, | ||
| 488 | const std::optional<BoundaryTypes>& currentScvfBoundaryTypes, | ||
| 489 | const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes, | ||
| 490 | const int localSubFaceIdx) | ||
| 491 | { | ||
| 492 | 270110224 | const auto eIdx = scvf.insideScvIdx(); | |
| 493 | 540220448 | const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); | |
| 494 | |||
| 495 | // Get the transporting velocity, located at the scvf perpendicular to the current scvf where the dof | ||
| 496 | // of interest is located. | ||
| 497 | 540220448 | const Scalar transportingVelocity = [&]() | |
| 498 | { | ||
| 499 |
2/2✓ Branch 0 taken 264880736 times.
✓ Branch 1 taken 5229488 times.
|
270110224 | const auto& faceVars = elemFaceVars[scvf]; |
| 500 |
2/2✓ Branch 0 taken 264880736 times.
✓ Branch 1 taken 5229488 times.
|
270110224 | if (!scvf.boundary()) |
| 501 | 264880736 | return faceVars.velocityLateralInside(localSubFaceIdx); | |
| 502 | else | ||
| 503 | { | ||
| 504 | // Create a boundaryTypes object. Get the boundary conditions. We sample the type of BC at the center of the current scvf. | ||
| 505 | 5229488 | const auto bcTypes = problem.boundaryTypes(element, scvf); | |
| 506 | |||
| 507 |
1/2✓ Branch 0 taken 5229488 times.
✗ Branch 1 not taken.
|
5229488 | if (bcTypes.isDirichlet(Indices::velocity(lateralFace.directionIndex()))) |
| 508 | { | ||
| 509 | // Construct a temporary scvf which corresponds to the staggered sub face, featuring the location | ||
| 510 | // the staggered faces's center. | ||
| 511 | ✗ | const auto& lateralBoundaryFacePos = lateralStaggeredFaceCenter_(scvf, localSubFaceIdx); | |
| 512 | ✗ | const auto lateralBoundaryFace = makeStaggeredBoundaryFace(scvf, lateralBoundaryFacePos); | |
| 513 | ✗ | return problem.dirichlet(element, lateralBoundaryFace)[Indices::velocity(lateralFace.directionIndex())]; | |
| 514 | } | ||
| 515 |
4/4✓ Branch 0 taken 149608 times.
✓ Branch 1 taken 5079880 times.
✓ Branch 2 taken 149608 times.
✓ Branch 3 taken 5079880 times.
|
10458976 | else if (bcTypes.isBeaversJoseph(Indices::velocity(lateralFace.directionIndex()))) |
| 516 | { | ||
| 517 | 448824 | return VelocityGradients::beaversJosephVelocityAtCurrentScvf(problem, element, fvGeometry, scvf, faceVars, | |
| 518 | 149608 | currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); | |
| 519 | } | ||
| 520 | else | ||
| 521 | 5079880 | return faceVars.velocityLateralInside(localSubFaceIdx); | |
| 522 | } | ||
| 523 | 540220448 | }(); | |
| 524 | |||
| 525 | 270110224 | const bool selfIsUpstream = lateralFace.directionSign() == sign(transportingVelocity); | |
| 526 | 540220448 | StaggeredUpwindHelper<TypeTag, upwindSchemeOrder> upwindHelper(element, fvGeometry, scvf, elemFaceVars, elemVolVars, gridFluxVarsCache.staggeredUpwindMethods()); | |
| 527 | 270110224 | FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); | |
| 528 | return upwindHelper.computeUpwindLateralMomentum(selfIsUpstream, lateralFace, localSubFaceIdx, currentScvfBoundaryTypes, lateralFaceBoundaryTypes) | ||
| 529 | 270110224 | * transportingVelocity * lateralFace.directionSign() * Extrusion::area(fvGeometry, lateralScvf) * extrusionFactor_(elemVolVars, lateralFace); | |
| 530 | } | ||
| 531 | |||
| 532 | /*! | ||
| 533 | * \brief Returns the diffusive momentum flux over the staggered face perpendicular to the scvf | ||
| 534 | * where the velocity dof of interest lives (coinciding with the element's scvfs). | ||
| 535 | * | ||
| 536 | * \verbatim | ||
| 537 | * ---------------- | ||
| 538 | * | |vel. | ||
| 539 | * | in.norm. |Parallel | ||
| 540 | * | vel. |~~~~> | ||
| 541 | * | ^ | ^ out.norm.vel. | ||
| 542 | * | | | | | ||
| 543 | * scvf ---------#######::::::::: || and # staggered half-control-volume (own element) | ||
| 544 | * | || | curr. :: | ||
| 545 | * | || | scvf :: :: staggered half-control-volume (neighbor element) | ||
| 546 | * | || x~~~~> :: | ||
| 547 | * | || | vel. :: # lateral staggered faces over which fluxes are calculated | ||
| 548 | * scvf | || | Self :: | ||
| 549 | * ---------#######::::::::: x dof position | ||
| 550 | * scvf | ||
| 551 | * -- elements | ||
| 552 | * \endverbatim | ||
| 553 | */ | ||
| 554 | 289801152 | FacePrimaryVariables computeDiffusivePartOfLateralMomentumFlux_(const Problem& problem, | |
| 555 | const FVElementGeometry& fvGeometry, | ||
| 556 | const Element& element, | ||
| 557 | const SubControlVolumeFace& scvf, | ||
| 558 | const ElementVolumeVariables& elemVolVars, | ||
| 559 | const FaceVariables& faceVars, | ||
| 560 | const std::optional<BoundaryTypes>& currentScvfBoundaryTypes, | ||
| 561 | const std::optional<BoundaryTypes>& lateralFaceBoundaryTypes, | ||
| 562 | const int localSubFaceIdx) | ||
| 563 | { | ||
| 564 |
2/2✓ Branch 0 taken 51 times.
✓ Branch 1 taken 289801101 times.
|
289801152 | const auto eIdx = scvf.insideScvIdx(); |
| 565 |
4/4✓ Branch 0 taken 51 times.
✓ Branch 1 taken 289801101 times.
✓ Branch 2 taken 51 times.
✓ Branch 3 taken 289801101 times.
|
579602304 | const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); |
| 566 | |||
| 567 |
2/2✓ Branch 0 taken 51 times.
✓ Branch 1 taken 289801101 times.
|
289801152 | FacePrimaryVariables lateralDiffusiveFlux(0.0); |
| 568 | |||
| 569 |
3/4✓ Branch 0 taken 51 times.
✓ Branch 1 taken 289801101 times.
✓ Branch 3 taken 51 times.
✗ Branch 4 not taken.
|
289801203 | static const bool enableUnsymmetrizedVelocityGradient |
| 570 |
2/4✓ Branch 1 taken 51 times.
✗ Branch 2 not taken.
✓ Branch 4 taken 51 times.
✗ Branch 5 not taken.
|
153 | = getParamFromGroup<bool>(problem.paramGroup(), "FreeFlow.EnableUnsymmetrizedVelocityGradient", false); |
| 571 | |||
| 572 | // Get the volume variables of the own and the neighboring element. The neighboring | ||
| 573 | // element is adjacent to the staggered face normal to the current scvf | ||
| 574 | // where the dof of interest is located. | ||
| 575 | 579602304 | const auto& insideVolVars = elemVolVars[lateralFace.insideScvIdx()]; | |
| 576 | 579602304 | const auto& outsideVolVars = elemVolVars[lateralFace.outsideScvIdx()]; | |
| 577 | |||
| 578 | // Get the averaged viscosity at the staggered face normal to the current scvf. | ||
| 579 | 289801152 | const Scalar muAvg = lateralFace.boundary() | |
| 580 |
2/2✓ Branch 0 taken 13309722 times.
✓ Branch 1 taken 276491430 times.
|
289801152 | ? insideVolVars.effectiveViscosity() |
| 581 | 829474290 | : (insideVolVars.effectiveViscosity() + outsideVolVars.effectiveViscosity()) * 0.5; | |
| 582 | |||
| 583 | // Consider the shear stress caused by the gradient of the velocities normal to our face of interest. | ||
| 584 |
1/2✓ Branch 0 taken 289801152 times.
✗ Branch 1 not taken.
|
289801152 | if (!enableUnsymmetrizedVelocityGradient) |
| 585 | { | ||
| 586 | 306504960 | if (!scvf.boundary() || | |
| 587 |
8/10✓ Branch 0 taken 5567936 times.
✓ Branch 1 taken 284233216 times.
✓ Branch 2 taken 5567936 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 5567936 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 346528 times.
✓ Branch 7 taken 5221408 times.
✓ Branch 8 taken 346528 times.
✓ Branch 9 taken 234936 times.
|
295950552 | currentScvfBoundaryTypes->isDirichlet(Indices::velocity(lateralFace.directionIndex())) || |
| 588 |
3/4✓ Branch 0 taken 346528 times.
✓ Branch 1 taken 5221408 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 4986472 times.
|
10554408 | currentScvfBoundaryTypes->isBeaversJoseph(Indices::velocity(lateralFace.directionIndex()))) |
| 589 | { | ||
| 590 | 284579744 | const Scalar velocityGrad_ji = VelocityGradients::velocityGradJI(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); | |
| 591 | // Account for the orientation of the staggered normal face's outer normal vector. | ||
| 592 | 284579744 | lateralDiffusiveFlux -= muAvg * velocityGrad_ji * lateralFace.directionSign(); | |
| 593 | } | ||
| 594 | } | ||
| 595 | |||
| 596 | // Consider the shear stress caused by the gradient of the velocities parallel to our face of interest. | ||
| 597 | // If we have a Dirichlet condition for the pressure at the lateral face we assume to have a zero velocityGrad_ij velocity gradient | ||
| 598 | // so we can skip the computation. | ||
| 599 |
6/6✓ Branch 0 taken 13309722 times.
✓ Branch 1 taken 276491430 times.
✓ Branch 2 taken 8904604 times.
✓ Branch 3 taken 4405118 times.
✓ Branch 4 taken 8904604 times.
✓ Branch 5 taken 4405118 times.
|
289801152 | if (!lateralFace.boundary() || !lateralFaceBoundaryTypes->isDirichlet(Indices::pressureIdx)) |
| 600 | { | ||
| 601 | 285396034 | const Scalar velocityGrad_ij = VelocityGradients::velocityGradIJ(problem, element, fvGeometry, scvf, faceVars, currentScvfBoundaryTypes, lateralFaceBoundaryTypes, localSubFaceIdx); | |
| 602 | 285396034 | lateralDiffusiveFlux -= muAvg * velocityGrad_ij * lateralFace.directionSign(); | |
| 603 | } | ||
| 604 | |||
| 605 | // Account for the area of the staggered lateral face (0.5 of the coinciding scfv). | ||
| 606 | 289801152 | FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); | |
| 607 | 289801152 | return lateralDiffusiveFlux * Extrusion::area(fvGeometry, lateralScvf) * extrusionFactor_(elemVolVars, lateralFace); | |
| 608 | } | ||
| 609 | |||
| 610 | /*! | ||
| 611 | * \brief Get the location of the lateral staggered face's center. | ||
| 612 | * Only needed for boundary conditions if the current scvf or the lateral one is on a boundary. | ||
| 613 | * | ||
| 614 | * \verbatim | ||
| 615 | * --------#######o || frontal face of staggered half-control-volume | ||
| 616 | * | || | current scvf # lateral staggered face of interest (may lie on a boundary) | ||
| 617 | * | || | x dof position | ||
| 618 | * | || x~~~~> vel.Self -- element boundaries, current scvf may lie on a boundary | ||
| 619 | * | || | o position at which the boundary conditions will be evaluated | ||
| 620 | * | || | (lateralStaggeredFaceCenter) | ||
| 621 | * ---------------- | ||
| 622 | * \endverbatim | ||
| 623 | */ | ||
| 624 | ✗ | const GlobalPosition& lateralStaggeredFaceCenter_(const SubControlVolumeFace& scvf, const int localSubFaceIdx) const | |
| 625 | { | ||
| 626 | ✗ | return scvf.pairData(localSubFaceIdx).lateralStaggeredFaceCenter; | |
| 627 | }; | ||
| 628 | |||
| 629 | /*! | ||
| 630 | * \brief Get the location of the lateral staggered sub control volume face center. | ||
| 631 | * Only needed for boundary conditions if the current scvf or the lateral one is on a boundary. | ||
| 632 | * | ||
| 633 | * \verbatim | ||
| 634 | * --------###o#### || frontal face of staggered half-control-volume | ||
| 635 | * | || | current scvf # lateral staggered face of interest (may lie on a boundary) | ||
| 636 | * | || | x dof position | ||
| 637 | * | || x~~~~> vel.Self -- element boundaries, current scvf may lie on a boundary | ||
| 638 | * | || | o center of the lateral staggered scvf | ||
| 639 | * | || | | ||
| 640 | * ---------------- | ||
| 641 | * \endverbatim | ||
| 642 | */ | ||
| 643 | ✗ | GlobalPosition lateralStaggeredSCVFCenter_(const SubControlVolumeFace& lateralFace, | |
| 644 | const SubControlVolumeFace& currentFace, | ||
| 645 | const int localSubFaceIdx) const | ||
| 646 | { | ||
| 647 | ✗ | auto pos = lateralStaggeredFaceCenter_(currentFace, localSubFaceIdx) + lateralFace.center(); | |
| 648 | ✗ | pos *= 0.5; | |
| 649 | ✗ | return pos; | |
| 650 | } | ||
| 651 | |||
| 652 | //! helper function to get the averaged extrusion factor for a face | ||
| 653 | 1102718040 | static Scalar extrusionFactor_(const ElementVolumeVariables& elemVolVars, const SubControlVolumeFace& scvf) | |
| 654 | { | ||
| 655 | 2205436080 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 656 | 2205436080 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 657 |
1/2✓ Branch 0 taken 1102718040 times.
✗ Branch 1 not taken.
|
1102718040 | return harmonicMean(insideVolVars.extrusionFactor(), outsideVolVars.extrusionFactor()); |
| 658 | } | ||
| 659 | |||
| 660 | //! do nothing if no turbulence model is used | ||
| 661 | template<class ...Args, bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<!turbulenceModel, int> = 0> | ||
| 662 | ✗ | bool incorporateWallFunction_(Args&&... args) const | |
| 663 | ✗ | { return false; } | |
| 664 | |||
| 665 | //! if a turbulence model is used, ask the problem is a wall function shall be employed and get the flux accordingly | ||
| 666 | template<bool turbulenceModel = ModelTraits::usesTurbulenceModel(), std::enable_if_t<turbulenceModel, int> = 0> | ||
| 667 | 3516771 | bool incorporateWallFunction_(FacePrimaryVariables& lateralFlux, | |
| 668 | const Problem& problem, | ||
| 669 | const Element& element, | ||
| 670 | const FVElementGeometry& fvGeometry, | ||
| 671 | const SubControlVolumeFace& scvf, | ||
| 672 | const ElementVolumeVariables& elemVolVars, | ||
| 673 | const ElementFaceVariables& elemFaceVars, | ||
| 674 | const std::size_t localSubFaceIdx) const | ||
| 675 | { | ||
| 676 | 3516771 | const auto eIdx = scvf.insideScvIdx(); | |
| 677 | 7033542 | const auto& lateralFace = fvGeometry.scvf(eIdx, scvf.pairData(localSubFaceIdx).localLateralFaceIdx); | |
| 678 | |||
| 679 |
2/2✓ Branch 1 taken 508848 times.
✓ Branch 2 taken 3007923 times.
|
3516771 | if (problem.useWallFunction(element, lateralFace, Indices::velocity(scvf.directionIndex()))) |
| 680 | { | ||
| 681 | 508848 | FaceLateralSubControlVolumeFace lateralScvf(lateralStaggeredSCVFCenter_(lateralFace, scvf, localSubFaceIdx), 0.5*lateralFace.area()); | |
| 682 |
0/2✗ Branch 2 not taken.
✗ Branch 3 not taken.
|
1526544 | const auto lateralBoundaryFace = makeStaggeredBoundaryFace(lateralFace, lateralStaggeredFaceCenter_(scvf, localSubFaceIdx)); |
| 683 |
1/2✓ Branch 1 taken 508848 times.
✗ Branch 2 not taken.
|
508848 | lateralFlux += problem.wallFunction(element, fvGeometry, elemVolVars, elemFaceVars, scvf, lateralBoundaryFace)[Indices::velocity(scvf.directionIndex())] |
| 684 |
2/4✓ Branch 1 taken 508848 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 508848 times.
✗ Branch 4 not taken.
|
508848 | * extrusionFactor_(elemVolVars, lateralFace) * Extrusion::area(fvGeometry, lateralScvf); |
| 685 |
1/2✓ Branch 0 taken 508848 times.
✗ Branch 1 not taken.
|
508848 | return true; |
| 686 | } | ||
| 687 | else | ||
| 688 | return false; | ||
| 689 | } | ||
| 690 | }; | ||
| 691 | |||
| 692 | } // end namespace Dumux | ||
| 693 | |||
| 694 | #endif | ||
| 695 |