| 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 CCTpfaFlux | ||
| 10 | * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH | ||
| 13 | #define DUMUX_DISCRETIZATION_CC_TPFA_DARCYS_LAW_HH | ||
| 14 | |||
| 15 | #include <dumux/common/math.hh> | ||
| 16 | #include <dumux/common/parameters.hh> | ||
| 17 | #include <dumux/common/properties.hh> | ||
| 18 | |||
| 19 | #include <dumux/discretization/method.hh> | ||
| 20 | #include <dumux/discretization/extrusion.hh> | ||
| 21 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 22 | |||
| 23 | namespace Dumux { | ||
| 24 | |||
| 25 | // forward declarations | ||
| 26 | template<class TypeTag, class DiscretizationMethod> | ||
| 27 | class DarcysLawImplementation; | ||
| 28 | |||
| 29 | /*! | ||
| 30 | * \ingroup CCTpfaFlux | ||
| 31 | * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation | ||
| 32 | * \note Darcy's law is specialized for network and surface grids (i.e. if grid dim < dimWorld) | ||
| 33 | * \tparam Scalar the scalar type for scalar physical quantities | ||
| 34 | * \tparam GridGeometry the grid geometry | ||
| 35 | * \tparam isNetwork whether we are computing on a network grid embedded in a higher world dimension | ||
| 36 | */ | ||
| 37 | template<class Scalar, class GridGeometry, bool isNetwork> | ||
| 38 | class CCTpfaDarcysLaw; | ||
| 39 | |||
| 40 | /*! | ||
| 41 | * \ingroup CCTpfaFlux | ||
| 42 | * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation | ||
| 43 | * \note Darcy's law is specialized for network and surface grids (i.e. if grid dim < dimWorld) | ||
| 44 | */ | ||
| 45 | template <class TypeTag> | ||
| 46 | class DarcysLawImplementation<TypeTag, DiscretizationMethods::CCTpfa> | ||
| 47 | : public CCTpfaDarcysLaw<GetPropType<TypeTag, Properties::Scalar>, | ||
| 48 | GetPropType<TypeTag, Properties::GridGeometry>, | ||
| 49 | (GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension < GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld)> | ||
| 50 | {}; | ||
| 51 | |||
| 52 | /*! | ||
| 53 | * \ingroup CCTpfaFlux | ||
| 54 | * \brief Class that fills the cache corresponding to tpfa Darcy's Law | ||
| 55 | */ | ||
| 56 | template<class GridGeometry> | ||
| 57 | class TpfaDarcysLawCacheFiller | ||
| 58 | { | ||
| 59 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 60 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 61 | using Element = typename GridGeometry::GridView::template Codim<0>::Entity; | ||
| 62 | |||
| 63 | public: | ||
| 64 | //! Function to fill a TpfaDarcysLawCache of a given scvf | ||
| 65 | //! This interface has to be met by any advection-related cache filler class | ||
| 66 | //! TODO: Probably get cache type out of the filler | ||
| 67 | template<class FluxVariablesCache, class Problem, class ElementVolumeVariables, class FluxVariablesCacheFiller> | ||
| 68 | ✗ | static void fill(FluxVariablesCache& scvfFluxVarsCache, | |
| 69 | const Problem& problem, | ||
| 70 | const Element& element, | ||
| 71 | const FVElementGeometry& fvGeometry, | ||
| 72 | const ElementVolumeVariables& elemVolVars, | ||
| 73 | const SubControlVolumeFace& scvf, | ||
| 74 | const FluxVariablesCacheFiller& fluxVarsCacheFiller) | ||
| 75 | { | ||
| 76 | 406166438 | scvfFluxVarsCache.updateAdvection(problem, element, fvGeometry, elemVolVars, scvf); | |
| 77 | ✗ | } | |
| 78 | }; | ||
| 79 | |||
| 80 | /*! | ||
| 81 | * \ingroup CCTpfaFlux | ||
| 82 | * \brief The cache corresponding to tpfa Darcy's Law | ||
| 83 | */ | ||
| 84 | template<class AdvectionType, class GridGeometry> | ||
| 85 | 4560 | class TpfaDarcysLawCache | |
| 86 | { | ||
| 87 | using Scalar = typename AdvectionType::Scalar; | ||
| 88 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 89 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 90 | using Element = typename GridGeometry::GridView::template Codim<0>::Entity; | ||
| 91 | |||
| 92 | public: | ||
| 93 | using Filler = TpfaDarcysLawCacheFiller<GridGeometry>; | ||
| 94 | |||
| 95 | template<class Problem, class ElementVolumeVariables> | ||
| 96 | ✗ | void updateAdvection(const Problem& problem, | |
| 97 | const Element& element, | ||
| 98 | const FVElementGeometry& fvGeometry, | ||
| 99 | const ElementVolumeVariables& elemVolVars, | ||
| 100 | const SubControlVolumeFace &scvf) | ||
| 101 | { | ||
| 102 | 196781737 | tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf); | |
| 103 | ✗ | } | |
| 104 | |||
| 105 | const Scalar& advectionTij() const | ||
| 106 | 399249161 | { return tij_; } | |
| 107 | |||
| 108 | private: | ||
| 109 | Scalar tij_; | ||
| 110 | }; | ||
| 111 | |||
| 112 | /*! | ||
| 113 | * \ingroup CCTpfaFlux | ||
| 114 | * \brief Specialization of the CCTpfaDarcysLaw grids where dim=dimWorld | ||
| 115 | */ | ||
| 116 | template<class ScalarType, class GridGeometry> | ||
| 117 | class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false> | ||
| 118 | { | ||
| 119 | using ThisType = CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ false>; | ||
| 120 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 121 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 122 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 123 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 124 | using GridView = typename GridGeometry::GridView; | ||
| 125 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 126 | |||
| 127 | static constexpr int dim = GridView::dimension; | ||
| 128 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 129 | |||
| 130 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 131 | |||
| 132 | public: | ||
| 133 | //! state the scalar type of the law | ||
| 134 | using Scalar = ScalarType; | ||
| 135 | |||
| 136 | using DiscretizationMethod = DiscretizationMethods::CCTpfa; | ||
| 137 | //! state the discretization method this implementation belongs to | ||
| 138 | static constexpr DiscretizationMethod discMethod{}; | ||
| 139 | |||
| 140 | //! state the type for the corresponding cache | ||
| 141 | using Cache = TpfaDarcysLawCache<ThisType, GridGeometry>; | ||
| 142 | |||
| 143 | /*! | ||
| 144 | * \brief Returns the advective flux of a fluid phase | ||
| 145 | * across the given sub-control volume face. | ||
| 146 | * \note This assembles the term | ||
| 147 | * \f$-|\sigma| \mathbf{n}^T \mathbf{K} \left( \nabla p - \rho \mathbf{g} \right)\f$, | ||
| 148 | * where \f$|\sigma|\f$ is the area of the face and \f$\mathbf{n}\f$ is the outer | ||
| 149 | * normal vector. Thus, the flux is given in N*m, and can be converted | ||
| 150 | * into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme | ||
| 151 | * for the mobility or the product of density and mobility, respectively. | ||
| 152 | */ | ||
| 153 | template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache> | ||
| 154 | 458278406 | static Scalar flux(const Problem& problem, | |
| 155 | const Element& element, | ||
| 156 | const FVElementGeometry& fvGeometry, | ||
| 157 | const ElementVolumeVariables& elemVolVars, | ||
| 158 | const SubControlVolumeFace& scvf, | ||
| 159 | int phaseIdx, | ||
| 160 | const ElementFluxVarsCache& elemFluxVarsCache) | ||
| 161 | { | ||
| 162 |
6/8✓ Branch 0 taken 303 times.
✓ Branch 1 taken 454429975 times.
✓ Branch 3 taken 188 times.
✓ Branch 4 taken 115 times.
✓ Branch 6 taken 188 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 188 times.
✗ Branch 10 not taken.
|
458278406 | static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity"); |
| 163 | |||
| 164 | 458278406 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 165 | |||
| 166 | // Get the inside and outside volume variables | ||
| 167 |
4/4✓ Branch 0 taken 131736831 times.
✓ Branch 1 taken 88490198 times.
✓ Branch 2 taken 131736831 times.
✓ Branch 3 taken 88490198 times.
|
916556812 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 168 | 458278406 | const auto& insideVolVars = elemVolVars[insideScv]; | |
| 169 | 916556812 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 170 | |||
| 171 |
2/2✓ Branch 0 taken 367072980 times.
✓ Branch 1 taken 87357298 times.
|
458278406 | if (enableGravity) |
| 172 | { | ||
| 173 | // do averaging for the density over all neighboring elements | ||
| 174 |
3/4✓ Branch 0 taken 935421 times.
✓ Branch 1 taken 366137559 times.
✓ Branch 2 taken 21760 times.
✗ Branch 3 not taken.
|
367082196 | const auto rho = scvf.boundary() ? outsideVolVars.density(phaseIdx) |
| 175 |
2/4✓ Branch 0 taken 8078264 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 8078264 times.
✗ Branch 3 not taken.
|
1098441765 | : (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5; |
| 176 | |||
| 177 | // Obtain inside and outside pressures | ||
| 178 |
2/2✓ Branch 0 taken 26704012 times.
✓ Branch 1 taken 12 times.
|
367082196 | const auto pInside = insideVolVars.pressure(phaseIdx); |
| 179 |
2/2✓ Branch 0 taken 26704012 times.
✓ Branch 1 taken 12 times.
|
367082196 | const auto pOutside = outsideVolVars.pressure(phaseIdx); |
| 180 | |||
| 181 | 734005272 | const auto& tij = fluxVarsCache.advectionTij(); | |
| 182 |
6/6✓ Branch 0 taken 18603988 times.
✓ Branch 1 taken 12 times.
✓ Branch 2 taken 18603988 times.
✓ Branch 3 taken 12 times.
✓ Branch 4 taken 18603988 times.
✓ Branch 5 taken 12 times.
|
1101246588 | const auto& g = problem.spatialParams().gravity(scvf.ipGlobal()); |
| 183 | |||
| 184 | //! compute alpha := n^T*K*g | ||
| 185 |
4/4✓ Branch 0 taken 18605068 times.
✓ Branch 1 taken 72 times.
✓ Branch 2 taken 18603988 times.
✓ Branch 3 taken 12 times.
|
738158772 | const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor(); |
| 186 | |||
| 187 |
12/12✓ Branch 0 taken 366137559 times.
✓ Branch 1 taken 935421 times.
✓ Branch 2 taken 1080 times.
✓ Branch 3 taken 60 times.
✓ Branch 4 taken 1080 times.
✓ Branch 5 taken 60 times.
✓ Branch 6 taken 1080 times.
✓ Branch 7 taken 60 times.
✓ Branch 8 taken 1080 times.
✓ Branch 9 taken 60 times.
✓ Branch 10 taken 1080 times.
✓ Branch 11 taken 60 times.
|
367087896 | Scalar flux = tij*(pInside - pOutside) + rho*Extrusion::area(fvGeometry, scvf)*alpha_inside; |
| 188 | |||
| 189 | //! On interior faces we have to add K-weighted gravitational contributions | ||
| 190 |
2/2✓ Branch 0 taken 366137559 times.
✓ Branch 1 taken 935421 times.
|
367082196 | if (!scvf.boundary()) |
| 191 | { | ||
| 192 |
4/4✓ Branch 0 taken 66719886 times.
✓ Branch 1 taken 299417673 times.
✓ Branch 2 taken 66719886 times.
✓ Branch 3 taken 299417673 times.
|
732293070 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 193 | 366146535 | const auto outsideK = outsideVolVars.permeability(); | |
| 194 | 732293070 | const auto outsideTi = fvGeometry.gridGeometry().isPeriodic() | |
| 195 |
2/2✓ Branch 0 taken 120592 times.
✓ Branch 1 taken 366016967 times.
|
713567410 | ? computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, outsideK, outsideVolVars.extrusionFactor()) |
| 196 | 384631011 | : -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, outsideK, outsideVolVars.extrusionFactor()); | |
| 197 | 732294150 | const auto alpha_outside = vtmv(scvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor(); | |
| 198 | |||
| 199 | 366151935 | flux -= rho*tij/outsideTi*(alpha_inside - alpha_outside); | |
| 200 | } | ||
| 201 | |||
| 202 | 367082196 | return flux; | |
| 203 | } | ||
| 204 | else | ||
| 205 | { | ||
| 206 | // Obtain inside and outside pressures | ||
| 207 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
91196210 | const auto pInside = insideVolVars.pressure(phaseIdx); |
| 208 |
0/2✗ Branch 0 not taken.
✗ Branch 1 not taken.
|
91196210 | const auto pOutside = outsideVolVars.pressure(phaseIdx); |
| 209 | |||
| 210 | // return flux | ||
| 211 | 94882471 | return fluxVarsCache.advectionTij()*(pInside - pOutside); | |
| 212 | } | ||
| 213 | } | ||
| 214 | |||
| 215 | // The flux variables cache has to be bound to an element prior to flux calculations | ||
| 216 | // During the binding, the transmissibility will be computed and stored using the method below. | ||
| 217 | template<class Problem, class ElementVolumeVariables> | ||
| 218 | 206267031 | static Scalar calculateTransmissibility(const Problem& problem, | |
| 219 | const Element& element, | ||
| 220 | const FVElementGeometry& fvGeometry, | ||
| 221 | const ElementVolumeVariables& elemVolVars, | ||
| 222 | const SubControlVolumeFace& scvf) | ||
| 223 | { | ||
| 224 | 2280 | Scalar tij; | |
| 225 | |||
| 226 |
2/2✓ Branch 0 taken 92509915 times.
✓ Branch 1 taken 82924260 times.
|
206267031 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 227 |
2/2✓ Branch 0 taken 92509915 times.
✓ Branch 1 taken 82924260 times.
|
206267031 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 228 | 206267031 | const auto& insideVolVars = elemVolVars[insideScvIdx]; | |
| 229 | |||
| 230 | 618796533 | const Scalar ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, | |
| 231 |
2/2✓ Branch 0 taken 38930 times.
✓ Branch 1 taken 32703500 times.
|
206267031 | getPermeability_(problem, insideVolVars, scvf.ipGlobal()), |
| 232 | insideVolVars.extrusionFactor()); | ||
| 233 | |||
| 234 | // on the boundary (dirichlet) we only need ti | ||
| 235 |
2/2✓ Branch 0 taken 10029847 times.
✓ Branch 1 taken 192122148 times.
|
206267031 | if (scvf.boundary()) |
| 236 | 20271086 | tij = Extrusion::area(fvGeometry, scvf)*ti; | |
| 237 | |||
| 238 | // otherwise we compute a tpfa harmonic mean | ||
| 239 | else | ||
| 240 | { | ||
| 241 |
2/2✓ Branch 0 taken 82924260 times.
✓ Branch 1 taken 82924260 times.
|
196131548 | const auto outsideScvIdx = scvf.outsideScvIdx(); |
| 242 | // as we assemble fluxes from the neighbor to our element the outside index | ||
| 243 | // refers to the scv of our element, so we use the scv method | ||
| 244 |
2/2✓ Branch 0 taken 82924260 times.
✓ Branch 1 taken 82924260 times.
|
196131548 | const auto& outsideScv = fvGeometry.scv(outsideScvIdx); |
| 245 | 196131548 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | |
| 246 | 392263096 | const Scalar tj = fvGeometry.gridGeometry().isPeriodic() | |
| 247 |
2/2✓ Branch 0 taken 241184 times.
✓ Branch 1 taken 191880964 times.
|
359557436 | ? computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor()) |
| 248 | 359075068 | : -1.0*computeTpfaTransmissibility(fvGeometry, scvf, outsideScv, getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), outsideVolVars.extrusionFactor()); | |
| 249 | |||
| 250 | // harmonic mean (check for division by zero!) | ||
| 251 | // TODO: This could lead to problems!? Is there a better way to do this? | ||
| 252 |
3/6✓ Branch 0 taken 192122148 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2160 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2160 times.
✗ Branch 5 not taken.
|
196135868 | if (ti*tj <= 0.0) |
| 253 | tij = 0; | ||
| 254 | else | ||
| 255 | 392271736 | tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); | |
| 256 | } | ||
| 257 | |||
| 258 | 206267031 | return tij; | |
| 259 | } | ||
| 260 | |||
| 261 | private: | ||
| 262 | template<class Problem, class VolumeVariables, | ||
| 263 | std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0> | ||
| 264 | ✗ | static decltype(auto) getPermeability_(const Problem& problem, | |
| 265 | const VolumeVariables& volVars, | ||
| 266 | const GlobalPosition& scvfIpGlobal) | ||
| 267 |
2/2✓ Branch 0 taken 38930 times.
✓ Branch 1 taken 32703500 times.
|
446148219 | { return volVars.permeability(); } |
| 268 | |||
| 269 | template<class Problem, class VolumeVariables, | ||
| 270 | std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0> | ||
| 271 | static decltype(auto) getPermeability_(const Problem& problem, | ||
| 272 | const VolumeVariables& volVars, | ||
| 273 | const GlobalPosition& scvfIpGlobal) | ||
| 274 | { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); } | ||
| 275 | }; | ||
| 276 | |||
| 277 | /*! | ||
| 278 | * \ingroup CCTpfaFlux | ||
| 279 | * \brief Specialization of the CCTpfaDarcysLaw grids where dim < dimWorld (network/surface grids) | ||
| 280 | */ | ||
| 281 | template<class ScalarType, class GridGeometry> | ||
| 282 | class CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ true> | ||
| 283 | { | ||
| 284 | using ThisType = CCTpfaDarcysLaw<ScalarType, GridGeometry, /*isNetwork*/ true>; | ||
| 285 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 286 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 287 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 288 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 289 | using GridView = typename GridGeometry::GridView; | ||
| 290 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 291 | |||
| 292 | static constexpr int dim = GridView::dimension; | ||
| 293 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 294 | |||
| 295 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 296 | |||
| 297 | public: | ||
| 298 | //! state the scalar type of the law | ||
| 299 | using Scalar = ScalarType; | ||
| 300 | |||
| 301 | using DiscretizationMethod = DiscretizationMethods::CCTpfa; | ||
| 302 | //! state the discretization method this implementation belongs to | ||
| 303 | static constexpr DiscretizationMethod discMethod{}; | ||
| 304 | |||
| 305 | //! state the type for the corresponding cache | ||
| 306 | using Cache = TpfaDarcysLawCache<ThisType, GridGeometry>; | ||
| 307 | |||
| 308 | /*! | ||
| 309 | * \brief Returns the advective flux of a fluid phase | ||
| 310 | * across the given sub-control volume face. | ||
| 311 | * \note This assembles the term | ||
| 312 | * \f$-|\sigma| \mathbf{n}^T \mathbf{K} \left( \nabla p - \rho \mathbf{g} \right)\f$, | ||
| 313 | * where \f$|\sigma|\f$ is the area of the face and \f$\mathbf{n}\f$ is the outer | ||
| 314 | * normal vector. Thus, the flux is given in N*m, and can be converted | ||
| 315 | * into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme | ||
| 316 | * for the mobility or the product of density and mobility, respectively. | ||
| 317 | */ | ||
| 318 | template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache> | ||
| 319 | 56118600 | static Scalar flux(const Problem& problem, | |
| 320 | const Element& element, | ||
| 321 | const FVElementGeometry& fvGeometry, | ||
| 322 | const ElementVolumeVariables& elemVolVars, | ||
| 323 | const SubControlVolumeFace& scvf, | ||
| 324 | int phaseIdx, | ||
| 325 | const ElementFluxVarsCache& elemFluxVarsCache) | ||
| 326 | { | ||
| 327 |
6/8✓ Branch 0 taken 127 times.
✓ Branch 1 taken 53967409 times.
✓ Branch 3 taken 74 times.
✓ Branch 4 taken 53 times.
✓ Branch 6 taken 74 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 74 times.
✗ Branch 10 not taken.
|
56118600 | static const bool gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity"); |
| 328 | |||
| 329 | 56118600 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 330 | |||
| 331 | // Get the inside and outside volume variables | ||
| 332 |
4/4✓ Branch 0 taken 1146990 times.
✓ Branch 1 taken 1129947 times.
✓ Branch 2 taken 1146990 times.
✓ Branch 3 taken 1129947 times.
|
112237200 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 333 | 56118600 | const auto& insideVolVars = elemVolVars[insideScv]; | |
| 334 | 112237200 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 335 | |||
| 336 |
2/2✓ Branch 0 taken 32415621 times.
✓ Branch 1 taken 21551915 times.
|
56118600 | if (gravity) |
| 337 | { | ||
| 338 | // do averaging for the density over all neighboring elements | ||
| 339 | 64911250 | const auto rho = [&]() | |
| 340 | { | ||
| 341 | // boundaries | ||
| 342 |
2/2✓ Branch 0 taken 71820 times.
✓ Branch 1 taken 32263793 times.
|
32415621 | if (scvf.boundary()) |
| 343 | 72052 | return outsideVolVars.density(phaseIdx); | |
| 344 | |||
| 345 | // inner faces with two neighboring elements | ||
| 346 |
4/4✓ Branch 0 taken 30446278 times.
✓ Branch 1 taken 1817515 times.
✓ Branch 2 taken 30446278 times.
✓ Branch 3 taken 1817515 times.
|
64687138 | else if (scvf.numOutsideScvs() == 1) |
| 347 | 61052108 | return (insideVolVars.density(phaseIdx) + outsideVolVars.density(phaseIdx))*0.5; | |
| 348 | |||
| 349 | // inner faces in networks (general case) | ||
| 350 | else | ||
| 351 | { | ||
| 352 | 3635030 | Scalar rho(insideVolVars.density(phaseIdx)); | |
| 353 |
4/4✓ Branch 0 taken 5982596 times.
✓ Branch 1 taken 1817515 times.
✓ Branch 2 taken 5982596 times.
✓ Branch 3 taken 1817515 times.
|
13782707 | for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) |
| 354 | { | ||
| 355 | 5982596 | const auto outsideScvIdx = scvf.outsideScvIdx(i); | |
| 356 | 5982596 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | |
| 357 | 5982596 | rho += outsideVolVars.density(phaseIdx); | |
| 358 | } | ||
| 359 | 3635030 | return rho/(scvf.numOutsideScvs()+1); | |
| 360 | } | ||
| 361 | 70893846 | }(); | |
| 362 | |||
| 363 | 64832138 | const auto& tij = fluxVarsCache.advectionTij(); | |
| 364 | 97486887 | const auto& g = problem.spatialParams().gravity(scvf.ipGlobal()); | |
| 365 | |||
| 366 | // Obtain inside and outside pressures | ||
| 367 | 32495629 | const auto pInside = insideVolVars.pressure(phaseIdx); | |
| 368 | 32415621 | const auto pOutside = [&]() | |
| 369 | { | ||
| 370 | // Dirichlet boundaries and inner faces with one neighbor | ||
| 371 |
4/4✓ Branch 0 taken 30518098 times.
✓ Branch 1 taken 1817515 times.
✓ Branch 2 taken 30518098 times.
✓ Branch 3 taken 1817515 times.
|
64831242 | if (scvf.numOutsideScvs() == 1) |
| 372 | 30598106 | return outsideVolVars.pressure(phaseIdx); | |
| 373 | |||
| 374 | // inner faces in networks (general case) | ||
| 375 | else | ||
| 376 | { | ||
| 377 | 1817515 | Scalar sumTi(tij); | |
| 378 | 1817515 | Scalar sumPTi(tij*pInside); | |
| 379 | |||
| 380 | // add inside gravitational contribution | ||
| 381 | 1817515 | sumPTi += rho*Extrusion::area(fvGeometry, scvf) | |
| 382 | 3635030 | *insideVolVars.extrusionFactor() | |
| 383 | 3635030 | *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g); | |
| 384 | |||
| 385 |
4/4✓ Branch 0 taken 5982596 times.
✓ Branch 1 taken 1817515 times.
✓ Branch 2 taken 5982596 times.
✓ Branch 3 taken 1817515 times.
|
13782707 | for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) |
| 386 | { | ||
| 387 | 5982596 | const auto outsideScvIdx = scvf.outsideScvIdx(i); | |
| 388 | 5982596 | const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i); | |
| 389 | 5982596 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | |
| 390 | 5982596 | const auto& outsideFluxVarsCache = elemFluxVarsCache[flippedScvf]; | |
| 391 | 5982596 | sumTi += outsideFluxVarsCache.advectionTij(); | |
| 392 | 17947788 | sumPTi += outsideFluxVarsCache.advectionTij()*outsideVolVars.pressure(phaseIdx); | |
| 393 | |||
| 394 | // add outside gravitational contribution | ||
| 395 | 5982596 | sumPTi += rho*Extrusion::area(fvGeometry, scvf) | |
| 396 | 5982596 | *outsideVolVars.extrusionFactor() | |
| 397 | 17947788 | *vtmv(flippedScvf.unitOuterNormal(), outsideVolVars.permeability(), g); | |
| 398 | } | ||
| 399 | 1817515 | return sumPTi/sumTi; | |
| 400 | } | ||
| 401 | 69076331 | }(); | |
| 402 | |||
| 403 | //! precompute alpha := n^T*K*g | ||
| 404 | 64991258 | const auto alpha_inside = vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g)*insideVolVars.extrusionFactor(); | |
| 405 | |||
| 406 |
2/2✓ Branch 0 taken 32343569 times.
✓ Branch 1 taken 72052 times.
|
32495629 | Scalar flux = tij*(pInside - pOutside) + Extrusion::area(fvGeometry, scvf)*rho*alpha_inside; |
| 407 | |||
| 408 | //! On interior faces with one neighbor we have to add K-weighted gravitational contributions | ||
| 409 |
6/6✓ Branch 0 taken 32343569 times.
✓ Branch 1 taken 72052 times.
✓ Branch 2 taken 30526054 times.
✓ Branch 3 taken 1817515 times.
✓ Branch 4 taken 30526054 times.
✓ Branch 5 taken 1817515 times.
|
32495629 | if (!scvf.boundary() && scvf.numOutsideScvs() == 1) |
| 410 | { | ||
| 411 |
4/4✓ Branch 0 taken 41008 times.
✓ Branch 1 taken 41008 times.
✓ Branch 2 taken 41008 times.
✓ Branch 3 taken 41008 times.
|
61211660 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 412 | 30605830 | const auto& outsideScvf = fvGeometry.flipScvf(scvf.index()); | |
| 413 | 30605830 | const auto outsideK = outsideVolVars.permeability(); | |
| 414 | 30605830 | const auto outsideTi = computeTpfaTransmissibility(fvGeometry, outsideScvf, outsideScv, outsideK, outsideVolVars.extrusionFactor()); | |
| 415 | 61211660 | const auto alpha_outside = vtmv(outsideScvf.unitOuterNormal(), outsideK, g)*outsideVolVars.extrusionFactor(); | |
| 416 | |||
| 417 | 30605830 | flux -= rho*tij/outsideTi*(alpha_inside + alpha_outside); | |
| 418 | } | ||
| 419 | |||
| 420 | return flux; | ||
| 421 | } | ||
| 422 | else | ||
| 423 | { | ||
| 424 | // Obtain inside and outside pressures | ||
| 425 | 23622971 | const auto pInside = insideVolVars.pressure(phaseIdx); | |
| 426 | 21551915 | const auto pOutside = [&]() | |
| 427 | { | ||
| 428 | // Dirichlet boundaries and inner faces with two neighboring elements | ||
| 429 |
4/4✓ Branch 0 taken 18461414 times.
✓ Branch 1 taken 1019445 times.
✓ Branch 2 taken 18461414 times.
✓ Branch 3 taken 1019445 times.
|
43103830 | if (scvf.numOutsideScvs() <= 1) |
| 430 | 20532470 | return outsideVolVars.pressure(phaseIdx); | |
| 431 | |||
| 432 | // inner faces in networks (general case) | ||
| 433 | else | ||
| 434 | { | ||
| 435 | 1019445 | const auto& insideFluxVarsCache = elemFluxVarsCache[scvf]; | |
| 436 | 1019445 | Scalar sumTi(insideFluxVarsCache.advectionTij()); | |
| 437 | 1019445 | Scalar sumPTi(insideFluxVarsCache.advectionTij()*pInside); | |
| 438 | |||
| 439 |
4/4✓ Branch 0 taken 3040045 times.
✓ Branch 1 taken 1019445 times.
✓ Branch 2 taken 3040045 times.
✓ Branch 3 taken 1019445 times.
|
7099535 | for (unsigned int i = 0; i < scvf.numOutsideScvs(); ++i) |
| 440 | { | ||
| 441 | 3040045 | const auto outsideScvIdx = scvf.outsideScvIdx(i); | |
| 442 | 3040045 | const auto& flippedScvf = fvGeometry.flipScvf(scvf.index(), i); | |
| 443 | 3040045 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | |
| 444 | 3040045 | const auto& outsideFluxVarsCache = elemFluxVarsCache[flippedScvf]; | |
| 445 | 3040045 | sumTi += outsideFluxVarsCache.advectionTij(); | |
| 446 | 3040045 | sumPTi += outsideFluxVarsCache.advectionTij()*outsideVolVars.pressure(phaseIdx); | |
| 447 | } | ||
| 448 | 1019445 | return sumPTi/sumTi; | |
| 449 | } | ||
| 450 | 47195486 | }(); | |
| 451 | |||
| 452 | // return flux | ||
| 453 | 27729107 | return fluxVarsCache.advectionTij()*(pInside - pOutside); | |
| 454 | } | ||
| 455 | } | ||
| 456 | |||
| 457 | // The flux variables cache has to be bound to an element prior to flux calculations | ||
| 458 | // During the binding, the transmissibility will be computed and stored using the method below. | ||
| 459 | template<class Problem, class ElementVolumeVariables> | ||
| 460 | 6441630 | static Scalar calculateTransmissibility(const Problem& problem, | |
| 461 | const Element& element, | ||
| 462 | const FVElementGeometry& fvGeometry, | ||
| 463 | const ElementVolumeVariables& elemVolVars, | ||
| 464 | const SubControlVolumeFace& scvf) | ||
| 465 | { | ||
| 466 | Scalar tij; | ||
| 467 | |||
| 468 |
2/2✓ Branch 0 taken 1623835 times.
✓ Branch 1 taken 1570503 times.
|
6441630 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 469 |
2/2✓ Branch 0 taken 1623835 times.
✓ Branch 1 taken 1570503 times.
|
6441630 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 470 | 6441630 | const auto& insideVolVars = elemVolVars[insideScvIdx]; | |
| 471 | |||
| 472 | 19325540 | const Scalar ti = computeTpfaTransmissibility(fvGeometry, scvf, insideScv, | |
| 473 | 6441630 | getPermeability_(problem, insideVolVars, scvf.ipGlobal()), | |
| 474 | insideVolVars.extrusionFactor()); | ||
| 475 | |||
| 476 | // for the boundary (dirichlet) or at branching points we only need ti | ||
| 477 |
6/6✓ Branch 0 taken 3340075 times.
✓ Branch 1 taken 58343 times.
✓ Branch 2 taken 8791 times.
✓ Branch 3 taken 3331284 times.
✓ Branch 4 taken 8791 times.
✓ Branch 5 taken 3331284 times.
|
6441630 | if (scvf.boundary() || scvf.numOutsideScvs() > 1) |
| 478 | 236552 | tij = Extrusion::area(fvGeometry, scvf)*ti; | |
| 479 | |||
| 480 | // otherwise we compute a tpfa harmonic mean | ||
| 481 | else | ||
| 482 | { | ||
| 483 |
2/2✓ Branch 0 taken 1569063 times.
✓ Branch 1 taken 1569063 times.
|
6323354 | const auto outsideScvIdx = scvf.outsideScvIdx(); |
| 484 | // as we assemble fluxes from the neighbor to our element the outside index | ||
| 485 | // refers to the scv of our element, so we use the scv method | ||
| 486 |
2/2✓ Branch 0 taken 1569063 times.
✓ Branch 1 taken 1569063 times.
|
6323354 | const auto& outsideScv = fvGeometry.scv(outsideScvIdx); |
| 487 | 6323354 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | |
| 488 | 6323354 | const Scalar tj = computeTpfaTransmissibility(fvGeometry, fvGeometry.flipScvf(scvf.index()), outsideScv, | |
| 489 | 12646092 | getPermeability_(problem, outsideVolVars, scvf.ipGlobal()), | |
| 490 | outsideVolVars.extrusionFactor()); | ||
| 491 | |||
| 492 | // harmonic mean (check for division by zero!) | ||
| 493 | // TODO: This could lead to problems!? Is there a better way to do this? | ||
| 494 |
1/2✓ Branch 0 taken 3331284 times.
✗ Branch 1 not taken.
|
6323354 | if (ti*tj <= 0.0) |
| 495 | tij = 0; | ||
| 496 | else | ||
| 497 | 12646708 | tij = Extrusion::area(fvGeometry, scvf)*(ti * tj)/(ti + tj); | |
| 498 | } | ||
| 499 | |||
| 500 | 6441630 | return tij; | |
| 501 | } | ||
| 502 | |||
| 503 | private: | ||
| 504 | template<class Problem, class VolumeVariables, | ||
| 505 | std::enable_if_t<!Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0> | ||
| 506 | ✗ | static decltype(auto) getPermeability_(const Problem& problem, | |
| 507 | const VolumeVariables& volVars, | ||
| 508 | const GlobalPosition& scvfIpGlobal) | ||
| 509 | 6728436 | { return volVars.permeability(); } | |
| 510 | |||
| 511 | template<class Problem, class VolumeVariables, | ||
| 512 | std::enable_if_t<Problem::SpatialParams::evaluatePermeabilityAtScvfIP(), int> = 0> | ||
| 513 | ✗ | static decltype(auto) getPermeability_(const Problem& problem, | |
| 514 | const VolumeVariables& volVars, | ||
| 515 | const GlobalPosition& scvfIpGlobal) | ||
| 516 | 3798 | { return problem.spatialParams().permeabilityAtPos(scvfIpGlobal); } | |
| 517 | }; | ||
| 518 | |||
| 519 | } // end namespace Dumux | ||
| 520 | |||
| 521 | #endif | ||
| 522 |