| 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 CVFEFlux | ||
| 10 | * \brief Specialization of Darcy's Law for control-volume finite element schemes | ||
| 11 | * | ||
| 12 | * This file contains the data which is required to calculate | ||
| 13 | * volume and mass fluxes of fluid phases over a face of a finite volume by means | ||
| 14 | * of the Darcy approximation. | ||
| 15 | */ | ||
| 16 | #ifndef DUMUX_DISCRETIZATION_CVFE_DARCYS_LAW_HH | ||
| 17 | #define DUMUX_DISCRETIZATION_CVFE_DARCYS_LAW_HH | ||
| 18 | |||
| 19 | #include <dumux/common/math.hh> | ||
| 20 | #include <dumux/common/parameters.hh> | ||
| 21 | #include <dumux/common/properties.hh> | ||
| 22 | #include <dumux/discretization/method.hh> | ||
| 23 | #include <dumux/discretization/extrusion.hh> | ||
| 24 | #include <dumux/flux/facetensoraverage.hh> | ||
| 25 | |||
| 26 | namespace Dumux { | ||
| 27 | |||
| 28 | /*! | ||
| 29 | * \ingroup CVFEFlux | ||
| 30 | * \brief Darcy's law for control-volume finite element schemes | ||
| 31 | * \tparam Scalar the scalar type for scalar physical quantities | ||
| 32 | * \tparam GridGeometry the grid geometry | ||
| 33 | */ | ||
| 34 | template<class Scalar, class GridGeometry> | ||
| 35 | class CVFEDarcysLaw | ||
| 36 | { | ||
| 37 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 38 | using SubControlVolumeFace = typename FVElementGeometry::SubControlVolumeFace; | ||
| 39 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 40 | using GridView = typename GridGeometry::GridView; | ||
| 41 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 42 | |||
| 43 | static constexpr int dim = GridView::dimension; | ||
| 44 | static constexpr int dimWorld = GridView::dimensionworld; | ||
| 45 | |||
| 46 | public: | ||
| 47 | |||
| 48 | /*! | ||
| 49 | * \brief Returns the advective flux of a fluid phase | ||
| 50 | * across the given sub-control volume face. | ||
| 51 | * \note This assembles the term | ||
| 52 | * \f$-|\sigma| \mathbf{n}^T \mathbf{K} \left( \nabla p - \rho \mathbf{g} \right)\f$, | ||
| 53 | * where \f$|\sigma|\f$ is the area of the face and \f$\mathbf{n}\f$ is the outer | ||
| 54 | * normal vector. Thus, the flux is given in N*m, and can be converted | ||
| 55 | * into a volume flux (m^3/s) or mass flux (kg/s) by applying an upwind scheme | ||
| 56 | * for the mobility or the product of density and mobility, respectively. | ||
| 57 | */ | ||
| 58 | template<class Problem, class ElementVolumeVariables, class ElementFluxVarsCache> | ||
| 59 | 645855798 | static Scalar flux(const Problem& problem, | |
| 60 | const Element& element, | ||
| 61 | const FVElementGeometry& fvGeometry, | ||
| 62 | const ElementVolumeVariables& elemVolVars, | ||
| 63 | const SubControlVolumeFace& scvf, | ||
| 64 | const int phaseIdx, | ||
| 65 | const ElementFluxVarsCache& elemFluxVarCache) | ||
| 66 | { | ||
| 67 |
2/2✓ Branch 0 taken 20400 times.
✓ Branch 1 taken 635763664 times.
|
645855798 | const auto& fluxVarCache = elemFluxVarCache[scvf]; |
| 68 |
4/4✓ Branch 0 taken 20400 times.
✓ Branch 1 taken 635763664 times.
✓ Branch 2 taken 20400 times.
✓ Branch 3 taken 635763664 times.
|
1291711596 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 69 |
3/4✓ Branch 0 taken 20400 times.
✓ Branch 1 taken 635763664 times.
✓ Branch 2 taken 624293152 times.
✗ Branch 3 not taken.
|
645876198 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 70 |
1/2✓ Branch 0 taken 624293152 times.
✗ Branch 1 not taken.
|
645855798 | const auto& insideVolVars = elemVolVars[insideScv]; |
| 71 |
1/2✓ Branch 0 taken 624293152 times.
✗ Branch 1 not taken.
|
645855798 | const auto& outsideVolVars = elemVolVars[outsideScv]; |
| 72 | |||
| 73 | 645855798 | auto insideK = insideVolVars.permeability(); | |
| 74 | 645855798 | auto outsideK = outsideVolVars.permeability(); | |
| 75 | |||
| 76 | // scale with correct extrusion factor | ||
| 77 | 645855798 | insideK *= insideVolVars.extrusionFactor(); | |
| 78 | 645855798 | outsideK *= outsideVolVars.extrusionFactor(); | |
| 79 | |||
| 80 |
2/4✓ Branch 0 taken 624293152 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 624272752 times.
✗ Branch 3 not taken.
|
1291691196 | const auto K = faceTensorAverage(insideK, outsideK, scvf.unitOuterNormal()); |
| 81 |
6/8✓ Branch 0 taken 167 times.
✓ Branch 1 taken 635783897 times.
✓ Branch 3 taken 123 times.
✓ Branch 4 taken 44 times.
✓ Branch 6 taken 123 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 123 times.
✗ Branch 10 not taken.
|
645855798 | static const bool enableGravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity"); |
| 82 | |||
| 83 | 645855798 | const auto& shapeValues = fluxVarCache.shapeValues(); | |
| 84 | |||
| 85 | // evaluate gradP - rho*g at integration point | ||
| 86 | 645855798 | Dune::FieldVector<Scalar, dimWorld> gradP(0.0); | |
| 87 | 645855798 | Scalar rho(0.0); | |
| 88 |
4/4✓ Branch 0 taken 2442773396 times.
✓ Branch 1 taken 635784064 times.
✓ Branch 2 taken 2442773396 times.
✓ Branch 3 taken 635784064 times.
|
6262993260 | for (auto&& scv : scvs(fvGeometry)) |
| 89 | { | ||
| 90 |
2/2✓ Branch 0 taken 1735222304 times.
✓ Branch 1 taken 707551092 times.
|
2485640832 | const auto& volVars = elemVolVars[scv]; |
| 91 | |||
| 92 |
2/2✓ Branch 0 taken 1735222304 times.
✓ Branch 1 taken 707551092 times.
|
2485640832 | if (enableGravity) |
| 93 |
1/2✓ Branch 0 taken 28800000 times.
✗ Branch 1 not taken.
|
5205666912 | rho += volVars.density(phaseIdx)*shapeValues[scv.indexInElement()][0]; |
| 94 | |||
| 95 | // the global shape function gradient | ||
| 96 |
2/4✓ Branch 0 taken 28800000 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 28800000 times.
✗ Branch 3 not taken.
|
9942563328 | gradP.axpy(volVars.pressure(phaseIdx), fluxVarCache.gradN(scv.indexInElement())); |
| 97 | } | ||
| 98 | |||
| 99 |
2/2✓ Branch 0 taken 449942144 times.
✓ Branch 1 taken 185841920 times.
|
645855798 | if (enableGravity) |
| 100 | 1799768576 | gradP.axpy(-rho, problem.spatialParams().gravity(scvf.center())); | |
| 101 | |||
| 102 | // apply the permeability and return the flux | ||
| 103 | 2583402792 | return -1.0*vtmv(scvf.unitOuterNormal(), K, gradP)*Extrusion::area(fvGeometry, scvf); | |
| 104 | } | ||
| 105 | |||
| 106 | // compute transmissibilities ti for analytical Jacobians | ||
| 107 | template<class Problem, class ElementVolumeVariables, class FluxVarCache> | ||
| 108 | 589556 | static std::vector<Scalar> calculateTransmissibilities(const Problem& problem, | |
| 109 | const Element& element, | ||
| 110 | const FVElementGeometry& fvGeometry, | ||
| 111 | const ElementVolumeVariables& elemVolVars, | ||
| 112 | const SubControlVolumeFace& scvf, | ||
| 113 | const FluxVarCache& fluxVarCache) | ||
| 114 | { | ||
| 115 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 589556 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 589556 times.
|
1179112 | const auto& insideScv = fvGeometry.scv(scvf.insideScvIdx()); |
| 116 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 589556 times.
✓ Branch 2 taken 589556 times.
✗ Branch 3 not taken.
|
589556 | const auto& outsideScv = fvGeometry.scv(scvf.outsideScvIdx()); |
| 117 |
1/2✓ Branch 0 taken 589556 times.
✗ Branch 1 not taken.
|
589556 | const auto& insideVolVars = elemVolVars[insideScv]; |
| 118 |
1/2✓ Branch 0 taken 589556 times.
✗ Branch 1 not taken.
|
589556 | const auto& outsideVolVars = elemVolVars[outsideScv]; |
| 119 | |||
| 120 | 589556 | auto insideK = insideVolVars.permeability(); | |
| 121 | 589556 | auto outsideK = outsideVolVars.permeability(); | |
| 122 | |||
| 123 | // scale with correct extrusion factor | ||
| 124 | 589556 | insideK *= insideVolVars.extrusionFactor(); | |
| 125 | 589556 | outsideK *= outsideVolVars.extrusionFactor(); | |
| 126 | |||
| 127 |
2/4✓ Branch 0 taken 589556 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 589556 times.
✗ Branch 3 not taken.
|
1179112 | const auto K = faceTensorAverage(insideK, outsideK, scvf.unitOuterNormal()); |
| 128 | |||
| 129 | 1768668 | std::vector<Scalar> ti(fvGeometry.numScv()); | |
| 130 |
4/4✓ Branch 0 taken 2089224 times.
✓ Branch 1 taken 589556 times.
✓ Branch 2 taken 2089224 times.
✓ Branch 3 taken 589556 times.
|
5357560 | for (const auto& scv : scvs(fvGeometry)) |
| 131 | 2089224 | ti[scv.indexInElement()] = | |
| 132 | 10446120 | -1.0*Extrusion::area(fvGeometry, scvf)*vtmv(scvf.unitOuterNormal(), K, fluxVarCache.gradN(scv.indexInElement())); | |
| 133 | |||
| 134 | 589556 | return ti; | |
| 135 | } | ||
| 136 | }; | ||
| 137 | |||
| 138 | // forward declaration | ||
| 139 | template<class TypeTag, class DiscretizationMethod> | ||
| 140 | class DarcysLawImplementation; | ||
| 141 | |||
| 142 | /*! | ||
| 143 | * \ingroup CVFEFlux | ||
| 144 | * \brief Specialization of Darcy's Law for control-volume finite element schemes | ||
| 145 | */ | ||
| 146 | template<class TypeTag, class DM> | ||
| 147 | class DarcysLawImplementation<TypeTag, DiscretizationMethods::CVFE<DM>> | ||
| 148 | : public CVFEDarcysLaw<GetPropType<TypeTag, Properties::Scalar>, GetPropType<TypeTag, Properties::GridGeometry>> | ||
| 149 | {}; | ||
| 150 | |||
| 151 | } // end namespace Dumux | ||
| 152 | |||
| 153 | #endif | ||
| 154 |