| 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 FacetCoupling | ||
| 10 | * \copydoc Dumux::CCTpfaFacetCouplingDarcysLaw | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_DISCRETIZATION_CC_TPFA_FACET_COUPLING_DARCYS_LAW_HH | ||
| 13 | #define DUMUX_DISCRETIZATION_CC_TPFA_FACET_COUPLING_DARCYS_LAW_HH | ||
| 14 | |||
| 15 | #include <array> | ||
| 16 | #include <cmath> | ||
| 17 | |||
| 18 | #include <dune/common/float_cmp.hh> | ||
| 19 | |||
| 20 | #include <dumux/common/math.hh> | ||
| 21 | #include <dumux/common/parameters.hh> | ||
| 22 | #include <dumux/common/properties.hh> | ||
| 23 | |||
| 24 | #include <dumux/discretization/method.hh> | ||
| 25 | #include <dumux/discretization/extrusion.hh> | ||
| 26 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 27 | #include <dumux/flux/cctpfa/darcyslaw.hh> | ||
| 28 | |||
| 29 | namespace Dumux { | ||
| 30 | |||
| 31 | //! Forward declaration of the implementation | ||
| 32 | template<class ScalarType, class GridGeometry, bool isNetwork> | ||
| 33 | class CCTpfaFacetCouplingDarcysLawImpl; | ||
| 34 | |||
| 35 | /*! | ||
| 36 | * \ingroup FacetCoupling | ||
| 37 | * \brief The cache corresponding to tpfa Darcy's Law with facet coupling | ||
| 38 | * \note We distinguish between network and non-network grids here. Specializations | ||
| 39 | * for the two cases can be found below. | ||
| 40 | */ | ||
| 41 | template<class AdvectionType, class GridGeometry, bool isNetwork> | ||
| 42 | class CCTpfaFacetCouplingDarcysLawCache; | ||
| 43 | |||
| 44 | /*! | ||
| 45 | * \ingroup FacetCoupling | ||
| 46 | * \brief Darcy's law for cell-centered finite volume schemes with two-point flux approximation | ||
| 47 | * in the context of coupled models where the coupling occurs across the facets of the bulk | ||
| 48 | * domain elements with a lower-dimensional domain living on these facets. | ||
| 49 | */ | ||
| 50 | template<class ScalarType, class GridGeometry> | ||
| 51 | using CCTpfaFacetCouplingDarcysLaw = | ||
| 52 | CCTpfaFacetCouplingDarcysLawImpl< ScalarType, GridGeometry, ( int(GridGeometry::GridView::dimension) < | ||
| 53 | int(GridGeometry::GridView::dimensionworld) ) >; | ||
| 54 | |||
| 55 | /*! | ||
| 56 | * \ingroup FacetCoupling | ||
| 57 | * \brief Specialization of FacetCouplingTpfaDarcysLawCache for non-network grids. | ||
| 58 | */ | ||
| 59 | template<class AdvectionType, class GridGeometry> | ||
| 60 | class CCTpfaFacetCouplingDarcysLawCache<AdvectionType, GridGeometry, /*isNetwork*/false> | ||
| 61 | { | ||
| 62 | using Scalar = typename AdvectionType::Scalar; | ||
| 63 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 64 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 65 | using Element = typename GridGeometry::GridView::template Codim<0>::Entity; | ||
| 66 | |||
| 67 | public: | ||
| 68 | //! export the corresponding filler class | ||
| 69 | using Filler = TpfaDarcysLawCacheFiller<GridGeometry>; | ||
| 70 | |||
| 71 | //! we store the transmissibilities associated with the interior | ||
| 72 | //! cell, outside cell, and the fracture facet in an array. Access | ||
| 73 | //! to this array should be done using the following indices: | ||
| 74 | static constexpr int insideTijIdx = 0; | ||
| 75 | static constexpr int outsideTijIdx = 1; | ||
| 76 | static constexpr int facetTijIdx = 2; | ||
| 77 | |||
| 78 | //! Export transmissibility storage type | ||
| 79 | using AdvectionTransmissibilityContainer = std::array<Scalar, 3>; | ||
| 80 | |||
| 81 | //! Export the type used for the gravity coefficients | ||
| 82 | using GravityCoefficients = std::array<Scalar, 2>; | ||
| 83 | |||
| 84 | //! update subject to a given problem | ||
| 85 | template< class Problem, class ElementVolumeVariables > | ||
| 86 | void updateAdvection(const Problem& problem, | ||
| 87 | const Element& element, | ||
| 88 | const FVElementGeometry& fvGeometry, | ||
| 89 | const ElementVolumeVariables& elemVolVars, | ||
| 90 | const SubControlVolumeFace &scvf) | ||
| 91 | { | ||
| 92 | 5904262 | tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g_); | |
| 93 | } | ||
| 94 | |||
| 95 | //! We use the same name as in the TpfaDarcysLawCache so | ||
| 96 | //! that this cache and the law implementation for non-coupled | ||
| 97 | //! models can be reused here on facets that do not lie on an | ||
| 98 | //! interior boundary, i.e. do not coincide with a facet element | ||
| 99 | Scalar advectionTij() const | ||
| 100 | 7690762 | { return tij_[insideTijIdx]; } | |
| 101 | |||
| 102 | //! returns the transmissibility associated with the inside cell | ||
| 103 | Scalar advectionTijInside() const | ||
| 104 | 208912 | { return tij_[insideTijIdx]; } | |
| 105 | |||
| 106 | //! returns the transmissibility associated with the outside cell | ||
| 107 | Scalar advectionTijOutside() const | ||
| 108 | 207952 | { return tij_[outsideTijIdx]; } | |
| 109 | |||
| 110 | //! returns the transmissibility associated with the outside cell | ||
| 111 | Scalar advectionTijFacet() const | ||
| 112 | 208912 | { return tij_[facetTijIdx]; } | |
| 113 | |||
| 114 | //! return the coefficients for the computation of gravity at the scvf | ||
| 115 | const GravityCoefficients& gravityCoefficients() const | ||
| 116 | 5568 | { return g_; } | |
| 117 | |||
| 118 | private: | ||
| 119 | std::array<Scalar, 3> tij_; | ||
| 120 | GravityCoefficients g_; | ||
| 121 | }; | ||
| 122 | |||
| 123 | /*! | ||
| 124 | * \ingroup FacetCoupling | ||
| 125 | * \brief Specialization of CCTpfaFacetCouplingDarcysLawImpl for dim=dimWorld | ||
| 126 | */ | ||
| 127 | template<class ScalarType, class GridGeometry> | ||
| 128 | class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/false> | ||
| 129 | { | ||
| 130 | using ThisType = CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, false>; | ||
| 131 | using TpfaDarcysLaw = CCTpfaDarcysLaw<ScalarType, GridGeometry, false>; | ||
| 132 | |||
| 133 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 134 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 135 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 136 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 137 | |||
| 138 | using GridView = typename GridGeometry::GridView; | ||
| 139 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 140 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 141 | |||
| 142 | public: | ||
| 143 | //! state the scalar type of the law | ||
| 144 | using Scalar = ScalarType; | ||
| 145 | //! export the discretization method this implementation belongs to | ||
| 146 | using DiscretizationMethod = DiscretizationMethods::CCTpfa; | ||
| 147 | static constexpr DiscretizationMethod discMethod{}; | ||
| 148 | //! export the type for the corresponding cache | ||
| 149 | using Cache = CCTpfaFacetCouplingDarcysLawCache<ThisType, GridGeometry, /*isNetwork*/false>; | ||
| 150 | //! export the type used to store transmissibilities | ||
| 151 | using TijContainer = typename Cache::AdvectionTransmissibilityContainer; | ||
| 152 | |||
| 153 | |||
| 154 | //! Compute the advective flux | ||
| 155 | template< class Problem, class ElementVolumeVariables, class ElementFluxVarsCache > | ||
| 156 | 3949837 | static Scalar flux(const Problem& problem, | |
| 157 | const Element& element, | ||
| 158 | const FVElementGeometry& fvGeometry, | ||
| 159 | const ElementVolumeVariables& elemVolVars, | ||
| 160 | const SubControlVolumeFace& scvf, | ||
| 161 | int phaseIdx, | ||
| 162 | const ElementFluxVarsCache& elemFluxVarsCache) | ||
| 163 | { | ||
| 164 |
6/6✓ Branch 0 taken 3845381 times.
✓ Branch 1 taken 104456 times.
✓ Branch 2 taken 3845381 times.
✓ Branch 3 taken 104456 times.
✓ Branch 4 taken 3845381 times.
✓ Branch 5 taken 104456 times.
|
11849511 | if (!problem.couplingManager().isOnInteriorBoundary(element, scvf)) |
| 165 | 3845381 | return TpfaDarcysLaw::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache); | |
| 166 | |||
| 167 | // Obtain inside and fracture pressures | ||
| 168 | 208912 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 169 | 208912 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 170 | 104456 | const auto pInside = insideVolVars.pressure(phaseIdx); | |
| 171 | 104456 | const auto pFacet = facetVolVars.pressure(phaseIdx); | |
| 172 | |||
| 173 | // compute and return flux | ||
| 174 | 104456 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 175 |
4/4✓ Branch 0 taken 19 times.
✓ Branch 1 taken 104437 times.
✓ Branch 2 taken 19 times.
✓ Branch 3 taken 104437 times.
|
208912 | Scalar flux = fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet; |
| 176 | |||
| 177 | // maybe add gravitational acceleration | ||
| 178 |
5/8✓ Branch 0 taken 19 times.
✓ Branch 1 taken 104437 times.
✓ Branch 3 taken 19 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 19 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 19 times.
✗ Branch 10 not taken.
|
104456 | static const Scalar gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity"); |
| 179 |
2/2✓ Branch 0 taken 2784 times.
✓ Branch 1 taken 101672 times.
|
104456 | if (gravity) |
| 180 | { | ||
| 181 | // compute alpha := n^T*K*g and add to flux (use arithmetic mean for density) | ||
| 182 | 8352 | const auto& g = problem.spatialParams().gravity(scvf.ipGlobal()); | |
| 183 | 5568 | const auto rho = 0.5*(insideVolVars.density(phaseIdx) + facetVolVars.density(phaseIdx)); | |
| 184 | 2784 | const auto rhoTimesArea = rho*Extrusion::area(fvGeometry, scvf); | |
| 185 | 5568 | const auto alpha_inside = rhoTimesArea*insideVolVars.extrusionFactor() | |
| 186 | 5568 | *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g); | |
| 187 | |||
| 188 | 2784 | flux += alpha_inside; | |
| 189 |
1/2✓ Branch 0 taken 2784 times.
✗ Branch 1 not taken.
|
2784 | if (!scvf.boundary()) |
| 190 | { | ||
| 191 | 5568 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 192 | |||
| 193 | // add further gravitational contributions | ||
| 194 |
3/6✓ Branch 0 taken 2784 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 2784 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 2784 times.
✗ Branch 5 not taken.
|
8352 | if ( problem.interiorBoundaryTypes(element, scvf).hasOnlyNeumann() ) |
| 195 | { | ||
| 196 |
5/8✓ Branch 0 taken 2 times.
✓ Branch 1 taken 2782 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 2 times.
✗ Branch 10 not taken.
|
2784 | static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0); |
| 197 | 5568 | const auto alpha_facet = rhoTimesArea*insideVolVars.extrusionFactor() | |
| 198 | 5568 | *vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), g); | |
| 199 | 5568 | const auto alpha_outside = rhoTimesArea*outsideVolVars.extrusionFactor() | |
| 200 | 5568 | *vtmv(scvf.unitOuterNormal(), outsideVolVars.permeability(), g); | |
| 201 | |||
| 202 | 2784 | flux -= fluxVarsCache.gravityCoefficients()[0]*(xi*alpha_inside - alpha_facet + (1.0 - xi)*alpha_outside); | |
| 203 | 5568 | flux += fluxVarsCache.gravityCoefficients()[1]*(xi*alpha_outside - alpha_facet + (1.0 - xi)*alpha_inside); | |
| 204 | } | ||
| 205 | |||
| 206 | // add outside contribution | ||
| 207 | 8352 | flux += fluxVarsCache.advectionTijOutside()*outsideVolVars.pressure(phaseIdx); | |
| 208 | } | ||
| 209 | |||
| 210 | 2784 | return flux; | |
| 211 | } | ||
| 212 | else | ||
| 213 |
2/2✓ Branch 0 taken 101192 times.
✓ Branch 1 taken 480 times.
|
101672 | return scvf.boundary() ? flux |
| 214 | 303576 | : flux + fluxVarsCache.advectionTijOutside()*elemVolVars[scvf.outsideScvIdx()].pressure(phaseIdx); | |
| 215 | } | ||
| 216 | |||
| 217 | // The flux variables cache has to be bound to an element prior to flux calculations | ||
| 218 | // During the binding, the transmissibility will be computed and stored using the method below. | ||
| 219 | template< class Problem, class ElementVolumeVariables > | ||
| 220 | static TijContainer calculateTransmissibility(const Problem& problem, | ||
| 221 | const Element& element, | ||
| 222 | const FVElementGeometry& fvGeometry, | ||
| 223 | const ElementVolumeVariables& elemVolVars, | ||
| 224 | const SubControlVolumeFace& scvf) | ||
| 225 | { | ||
| 226 | typename Cache::GravityCoefficients g; | ||
| 227 | return calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g); | ||
| 228 | } | ||
| 229 | |||
| 230 | // This overload additionally receives a container in which the coefficients required | ||
| 231 | // for the computation of the gravitational acceleration ar the scvf are stored | ||
| 232 | template< class Problem, class ElementVolumeVariables > | ||
| 233 | 5904262 | static TijContainer calculateTransmissibility(const Problem& problem, | |
| 234 | const Element& element, | ||
| 235 | const FVElementGeometry& fvGeometry, | ||
| 236 | const ElementVolumeVariables& elemVolVars, | ||
| 237 | const SubControlVolumeFace& scvf, | ||
| 238 | typename Cache::GravityCoefficients& g) | ||
| 239 | { | ||
| 240 | TijContainer tij; | ||
| 241 |
6/6✓ Branch 0 taken 5747554 times.
✓ Branch 1 taken 156708 times.
✓ Branch 2 taken 5747554 times.
✓ Branch 3 taken 156708 times.
✓ Branch 4 taken 5747554 times.
✓ Branch 5 taken 156708 times.
|
17712786 | if (!problem.couplingManager().isCoupled(element, scvf)) |
| 242 | { | ||
| 243 | //! use the standard darcy's law and only compute one transmissibility | ||
| 244 | 5747554 | tij[Cache::insideTijIdx] = TpfaDarcysLaw::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf); | |
| 245 | 5747554 | return tij; | |
| 246 | } | ||
| 247 | |||
| 248 | //! xi factor for the coupling conditions | ||
| 249 |
5/8✓ Branch 0 taken 19 times.
✓ Branch 1 taken 156689 times.
✓ Branch 3 taken 19 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 19 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 19 times.
✗ Branch 10 not taken.
|
156708 | static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0); |
| 250 | |||
| 251 |
2/2✓ Branch 0 taken 78654 times.
✓ Branch 1 taken 78054 times.
|
156708 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 252 |
2/2✓ Branch 0 taken 78654 times.
✓ Branch 1 taken 78054 times.
|
156708 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 253 | 156708 | const auto& insideVolVars = elemVolVars[insideScvIdx]; | |
| 254 | 156708 | const auto wIn = Extrusion::area(fvGeometry, scvf) | |
| 255 | 156708 | *computeTpfaTransmissibility(fvGeometry, scvf, insideScv, | |
| 256 | insideVolVars.permeability(), | ||
| 257 | insideVolVars.extrusionFactor()); | ||
| 258 | |||
| 259 | // proceed depending on the interior BC types used | ||
| 260 |
2/2✓ Branch 0 taken 71436 times.
✓ Branch 1 taken 600 times.
|
156708 | const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf); |
| 261 | |||
| 262 | // neumann-coupling | ||
| 263 |
4/4✓ Branch 0 taken 156108 times.
✓ Branch 1 taken 600 times.
✓ Branch 2 taken 156108 times.
✓ Branch 3 taken 600 times.
|
313416 | if (iBcTypes.hasOnlyNeumann()) |
| 264 | { | ||
| 265 | 312216 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 266 | 156108 | const auto wFacet = 2.0*Extrusion::area(fvGeometry, scvf)*insideVolVars.extrusionFactor() | |
| 267 | 156108 | /facetVolVars.extrusionFactor() | |
| 268 | 468324 | *vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), scvf.unitOuterNormal()); | |
| 269 | |||
| 270 | // The fluxes across this face and the outside face can be expressed in matrix form: | ||
| 271 | // \f$\mathbf{C} \bar{\mathbf{u}} + \mathbf{D} \mathbf{u} + \mathbf{E} \mathbf{u}_\gamma\f$, | ||
| 272 | // where \f$\gamma$\f denotes the domain living on the facets and \f$\bar{\mathbf{u}}$\f are | ||
| 273 | // intermediate face unknowns in the matrix domain. Equivalently, flux continuity reads: | ||
| 274 | // \f$\mathbf{A} \bar{\mathbf{u}} = \mathbf{B} \mathbf{u} + \mathbf{M} \mathbf{u}_\gamma\f$. | ||
| 275 | // Combining the two, we can eliminate the intermediate unknowns and compute the transmissibilities | ||
| 276 | // that allow the description of the fluxes as functions of the cell and Dirichlet pressures only. | ||
| 277 |
1/2✓ Branch 0 taken 156108 times.
✗ Branch 1 not taken.
|
156108 | if (!scvf.boundary()) |
| 278 | { | ||
| 279 | 156108 | const auto outsideScvIdx = scvf.outsideScvIdx(); | |
| 280 | 156108 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | |
| 281 |
2/2✓ Branch 0 taken 78054 times.
✓ Branch 1 taken 78054 times.
|
156108 | const auto wOut = -1.0*Extrusion::area(fvGeometry, scvf) |
| 282 |
2/2✓ Branch 0 taken 78054 times.
✓ Branch 1 taken 78054 times.
|
312216 | *computeTpfaTransmissibility(fvGeometry, scvf, fvGeometry.scv(outsideScvIdx), |
| 283 | outsideVolVars.permeability(), | ||
| 284 | outsideVolVars.extrusionFactor()); | ||
| 285 | |||
| 286 |
4/4✓ Branch 0 taken 5950 times.
✓ Branch 1 taken 150158 times.
✓ Branch 2 taken 5950 times.
✓ Branch 3 taken 150158 times.
|
162058 | if ( !Dune::FloatCmp::eq(xi, 1.0, 1e-6) ) |
| 287 | { | ||
| 288 | // The gravity coefficients are the first row of the inverse of the A matrix in the local eq system | ||
| 289 | // multiplied with wIn. Note that we never compute the inverse but use an optimized implementation below. | ||
| 290 | // The A matrix has the following coefficients: | ||
| 291 | // A = | xi*wIn + wFacet, (xi - 1.0)*wOut | -> AInv = 1/detA | xi*wOut + wFacet, -(xi - 1.0)*wOut | | ||
| 292 | // | wIn*(xi - 1.0) , xi*wOut + wFacet | | -wIn*(xi - 1.0) , xi*wIn + wFacet | | ||
| 293 | 5950 | const Scalar xiMinusOne = (xi - 1.0); | |
| 294 | 5950 | const Scalar a01 = xiMinusOne*wOut; | |
| 295 | 5950 | const Scalar a11 = xi*wOut + wFacet; | |
| 296 | 5950 | const Scalar detA = (xi*wIn + wFacet)*a11 - xiMinusOne*wIn*a01; | |
| 297 | 11900 | g[0] = wIn*a11/detA; g[1] = -wIn*a01/detA; | |
| 298 | |||
| 299 | // optimized implementation: factorization obtained using sympy | ||
| 300 | 5950 | const Scalar factor = wIn * wFacet / ( wIn * wOut * ( 2.0 * xi - 1.0 ) + wFacet * ( xi * ( wIn + wOut ) + wFacet ) ); | |
| 301 | 5950 | tij[Cache::insideTijIdx] = factor * ( wOut * xi + wFacet ); | |
| 302 | 5950 | tij[Cache::outsideTijIdx] = factor * ( wOut * ( 1.0 - xi ) ); | |
| 303 | 11900 | tij[Cache::facetTijIdx] = factor * ( - wOut - wFacet ); | |
| 304 | } | ||
| 305 | else | ||
| 306 | { | ||
| 307 | 300316 | g[0] = wIn/(wIn+wFacet); g[1] = 0.0; | |
| 308 | 300316 | tij[Cache::insideTijIdx] = wFacet*g[0]; | |
| 309 | 300316 | tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx]; | |
| 310 | 300316 | tij[Cache::outsideTijIdx] = 0.0; | |
| 311 | } | ||
| 312 | } | ||
| 313 | else | ||
| 314 | { | ||
| 315 | // TODO: check for division by zero?? | ||
| 316 | ✗ | tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet); | |
| 317 | ✗ | tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx]; | |
| 318 | ✗ | tij[Cache::outsideTijIdx] = 0.0; | |
| 319 | } | ||
| 320 | } | ||
| 321 |
2/4✓ Branch 0 taken 600 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 600 times.
✗ Branch 3 not taken.
|
1200 | else if (iBcTypes.hasOnlyDirichlet()) |
| 322 | { | ||
| 323 | 600 | tij[Cache::insideTijIdx] = wIn; | |
| 324 | 600 | tij[Cache::outsideTijIdx] = 0.0; | |
| 325 | 1200 | tij[Cache::facetTijIdx] = -wIn; | |
| 326 | } | ||
| 327 | else | ||
| 328 | ✗ | DUNE_THROW(Dune::NotImplemented, "Interior boundary types other than pure Dirichlet or Neumann"); | |
| 329 | |||
| 330 | return tij; | ||
| 331 | } | ||
| 332 | }; | ||
| 333 | |||
| 334 | /*! | ||
| 335 | * \ingroup FacetCoupling | ||
| 336 | * \brief Specialization of FacetCouplingTpfaDarcysLawCache for network grids | ||
| 337 | */ | ||
| 338 | template<class AdvectionType, class GridGeometry> | ||
| 339 | class CCTpfaFacetCouplingDarcysLawCache<AdvectionType, GridGeometry, /*isNetwork*/true> | ||
| 340 | { | ||
| 341 | using Scalar = typename AdvectionType::Scalar; | ||
| 342 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 343 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 344 | using Element = typename GridGeometry::GridView::template Codim<0>::Entity; | ||
| 345 | |||
| 346 | public: | ||
| 347 | //! export the corresponding filler class | ||
| 348 | using Filler = TpfaDarcysLawCacheFiller<GridGeometry>; | ||
| 349 | |||
| 350 | //! we store the transmissibilities associated with the interior | ||
| 351 | //! cell and the fracture facet in an array. Access to this array | ||
| 352 | //! should be done using the following indices: | ||
| 353 | static constexpr int insideTijIdx = 0; | ||
| 354 | static constexpr int facetTijIdx = 1; | ||
| 355 | |||
| 356 | //! Export transmissibility storage type | ||
| 357 | using AdvectionTransmissibilityContainer = std::array<Scalar, 2>; | ||
| 358 | |||
| 359 | //! Export the type used for the gravity coefficients | ||
| 360 | using GravityCoefficients = std::array<Scalar, 1>; | ||
| 361 | |||
| 362 | //! update subject to a given problem | ||
| 363 | template< class Problem, class ElementVolumeVariables > | ||
| 364 | void updateAdvection(const Problem& problem, | ||
| 365 | const Element& element, | ||
| 366 | const FVElementGeometry& fvGeometry, | ||
| 367 | const ElementVolumeVariables& elemVolVars, | ||
| 368 | const SubControlVolumeFace &scvf) | ||
| 369 | { | ||
| 370 | 3110716 | tij_ = AdvectionType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g_); | |
| 371 | } | ||
| 372 | |||
| 373 | //! We use the same name as in the TpfaDarcysLawCache so | ||
| 374 | //! that this cache and the law implementation for non-coupled | ||
| 375 | //! models can be reused here on facets that do not lie on an | ||
| 376 | //! interior boundary, i.e. do not coincide with a facet element | ||
| 377 | Scalar advectionTij() const | ||
| 378 | 4265256 | { return tij_[insideTijIdx]; } | |
| 379 | |||
| 380 | //! returns the transmissibility associated with the inside cell | ||
| 381 | Scalar advectionTijInside() const | ||
| 382 | 120192 | { return tij_[insideTijIdx]; } | |
| 383 | |||
| 384 | //! returns the transmissibility associated with the outside cell | ||
| 385 | Scalar advectionTijFacet() const | ||
| 386 | 120192 | { return tij_[facetTijIdx]; } | |
| 387 | |||
| 388 | //! return the coefficients for the computation of gravity at the scvf | ||
| 389 | const GravityCoefficients& gravityCoefficients() const | ||
| 390 | 1392 | { return g_; } | |
| 391 | |||
| 392 | private: | ||
| 393 | AdvectionTransmissibilityContainer tij_; | ||
| 394 | GravityCoefficients g_; | ||
| 395 | }; | ||
| 396 | |||
| 397 | /*! | ||
| 398 | * \ingroup FacetCoupling | ||
| 399 | * \brief Specialization of CCTpfaFacetCouplingDarcysLawImpl for dim<dimWorld | ||
| 400 | */ | ||
| 401 | template<class ScalarType, class GridGeometry> | ||
| 402 | class CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, /*isNetwork*/true> | ||
| 403 | { | ||
| 404 | using ThisType = CCTpfaFacetCouplingDarcysLawImpl<ScalarType, GridGeometry, true>; | ||
| 405 | using TpfaDarcysLaw = CCTpfaDarcysLaw<ScalarType, GridGeometry, true>; | ||
| 406 | |||
| 407 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 408 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 409 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 410 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 411 | |||
| 412 | using GridView = typename GridGeometry::GridView; | ||
| 413 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 414 | using GlobalPosition = typename Element::Geometry::GlobalCoordinate; | ||
| 415 | |||
| 416 | public: | ||
| 417 | //! state the scalar type of the law | ||
| 418 | using Scalar = ScalarType; | ||
| 419 | //! state the discretization method this implementation belongs to | ||
| 420 | using DiscretizationMethod = DiscretizationMethods::CCTpfa; | ||
| 421 | static constexpr DiscretizationMethod discMethod{}; | ||
| 422 | //! state the type for the corresponding cache | ||
| 423 | using Cache = CCTpfaFacetCouplingDarcysLawCache<ThisType, GridGeometry, /*isNetwork*/true>; | ||
| 424 | //! export the type used to store transmissibilities | ||
| 425 | using TijContainer = typename Cache::AdvectionTransmissibilityContainer; | ||
| 426 | |||
| 427 | //! Compute the advective flux | ||
| 428 | template< class Problem, class ElementVolumeVariables, class ElementFluxVarsCache > | ||
| 429 | 2192724 | static Scalar flux(const Problem& problem, | |
| 430 | const Element& element, | ||
| 431 | const FVElementGeometry& fvGeometry, | ||
| 432 | const ElementVolumeVariables& elemVolVars, | ||
| 433 | const SubControlVolumeFace& scvf, | ||
| 434 | int phaseIdx, | ||
| 435 | const ElementFluxVarsCache& elemFluxVarsCache) | ||
| 436 | { | ||
| 437 |
6/6✓ Branch 0 taken 2132628 times.
✓ Branch 1 taken 60096 times.
✓ Branch 2 taken 2132628 times.
✓ Branch 3 taken 60096 times.
✓ Branch 4 taken 2132628 times.
✓ Branch 5 taken 60096 times.
|
6578172 | if (!problem.couplingManager().isOnInteriorBoundary(element, scvf)) |
| 438 | 2132628 | return TpfaDarcysLaw::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache); | |
| 439 | |||
| 440 | // On surface grids only xi = 1.0 can be used, as the coupling condition | ||
| 441 | // for xi != 1.0 does not generalize for surface grids where there can be | ||
| 442 | // seveal neighbor meeting at a branching point. | ||
| 443 |
5/8✓ Branch 0 taken 5 times.
✓ Branch 1 taken 60091 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
|
60096 | static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0); |
| 444 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 60096 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 60096 times.
|
60096 | if (Dune::FloatCmp::ne(xi, 1.0, 1e-6)) |
| 445 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Xi != 1.0 cannot be used on surface grids"); | |
| 446 | |||
| 447 | // Obtain inside and fracture pressures | ||
| 448 | 120192 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 449 | 120192 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 450 | 60096 | const auto pInside = insideVolVars.pressure(phaseIdx); | |
| 451 | 60096 | const auto pFacet = facetVolVars.pressure(phaseIdx); | |
| 452 | |||
| 453 | // compute and return flux | ||
| 454 | 60096 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 455 |
4/4✓ Branch 0 taken 5 times.
✓ Branch 1 taken 60091 times.
✓ Branch 2 taken 5 times.
✓ Branch 3 taken 60091 times.
|
120192 | Scalar flux = fluxVarsCache.advectionTijInside()*pInside + fluxVarsCache.advectionTijFacet()*pFacet; |
| 456 | |||
| 457 |
5/8✓ Branch 0 taken 5 times.
✓ Branch 1 taken 60091 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
|
60096 | static const Scalar gravity = getParamFromGroup<bool>(problem.paramGroup(), "Problem.EnableGravity"); |
| 458 |
2/2✓ Branch 0 taken 1392 times.
✓ Branch 1 taken 58704 times.
|
60096 | if (gravity) |
| 459 | { | ||
| 460 | // compute alpha := n^T*K*g and add to flux (use arithmetic mean for density) | ||
| 461 | 4176 | const auto& g = problem.spatialParams().gravity(scvf.ipGlobal()); | |
| 462 | 2784 | const auto rho = 0.5*(insideVolVars.density(phaseIdx) + facetVolVars.density(phaseIdx)); | |
| 463 | 1392 | const auto rhoTimesArea = rho*Extrusion::area(fvGeometry, scvf); | |
| 464 | 2784 | const auto alpha_inside = rhoTimesArea*insideVolVars.extrusionFactor() | |
| 465 | 2784 | *vtmv(scvf.unitOuterNormal(), insideVolVars.permeability(), g); | |
| 466 | |||
| 467 | 1392 | flux += alpha_inside; | |
| 468 | |||
| 469 | // maybe add further gravitational contributions | ||
| 470 |
4/8✓ Branch 0 taken 1392 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 1392 times.
✗ Branch 3 not taken.
✓ Branch 4 taken 1392 times.
✗ Branch 5 not taken.
✓ Branch 6 taken 1392 times.
✗ Branch 7 not taken.
|
1392 | if ( !scvf.boundary() && problem.interiorBoundaryTypes(element, scvf).hasOnlyNeumann() ) |
| 471 | { | ||
| 472 | 2784 | const auto alpha_facet = rhoTimesArea*insideVolVars.extrusionFactor() | |
| 473 | 2784 | *vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), g); | |
| 474 | |||
| 475 | 2784 | flux -= fluxVarsCache.gravityCoefficients()[0]*(alpha_inside - alpha_facet); | |
| 476 | } | ||
| 477 | } | ||
| 478 | |||
| 479 | return flux; | ||
| 480 | } | ||
| 481 | |||
| 482 | // The flux variables cache has to be bound to an element prior to flux calculations | ||
| 483 | // During the binding, the transmissibility will be computed and stored using the method below. | ||
| 484 | template< class Problem, class ElementVolumeVariables > | ||
| 485 | static TijContainer calculateTransmissibility(const Problem& problem, | ||
| 486 | const Element& element, | ||
| 487 | const FVElementGeometry& fvGeometry, | ||
| 488 | const ElementVolumeVariables& elemVolVars, | ||
| 489 | const SubControlVolumeFace& scvf) | ||
| 490 | { | ||
| 491 | typename Cache::GravityCoefficients g; | ||
| 492 | return calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, g); | ||
| 493 | } | ||
| 494 | |||
| 495 | // This overload additionally receives a container in which the coefficients required | ||
| 496 | // for the computation of the gravitational acceleration ar the scvf are stored | ||
| 497 | template< class Problem, class ElementVolumeVariables > | ||
| 498 | 3110716 | static TijContainer calculateTransmissibility(const Problem& problem, | |
| 499 | const Element& element, | ||
| 500 | const FVElementGeometry& fvGeometry, | ||
| 501 | const ElementVolumeVariables& elemVolVars, | ||
| 502 | const SubControlVolumeFace& scvf, | ||
| 503 | typename Cache::GravityCoefficients& g) | ||
| 504 | { | ||
| 505 | TijContainer tij; | ||
| 506 |
6/6✓ Branch 0 taken 3018402 times.
✓ Branch 1 taken 92314 times.
✓ Branch 2 taken 3018402 times.
✓ Branch 3 taken 92314 times.
✓ Branch 4 taken 3018402 times.
✓ Branch 5 taken 92314 times.
|
9332148 | if (!problem.couplingManager().isCoupled(element, scvf)) |
| 507 | { | ||
| 508 | //! use the standard darcy's law and only compute one transmissibility | ||
| 509 | 3018402 | tij[Cache::insideTijIdx] = TpfaDarcysLaw::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf); | |
| 510 | 3018402 | return tij; | |
| 511 | } | ||
| 512 | |||
| 513 | //! xi factor for the coupling conditions | ||
| 514 |
5/8✓ Branch 0 taken 5 times.
✓ Branch 1 taken 92309 times.
✓ Branch 3 taken 5 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 5 times.
✗ Branch 7 not taken.
✓ Branch 9 taken 5 times.
✗ Branch 10 not taken.
|
92314 | static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0); |
| 515 | |||
| 516 | // On surface grids only xi = 1.0 can be used, as the coupling condition | ||
| 517 | // for xi != 1.0 does not generalize for surface grids where the normal | ||
| 518 | // vectors of the inside/outside elements have different orientations. | ||
| 519 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 92314 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 92314 times.
|
92314 | if (Dune::FloatCmp::ne(xi, 1.0, 1e-6)) |
| 520 | ✗ | DUNE_THROW(Dune::InvalidStateException, "Xi != 1.0 cannot be used on surface grids"); | |
| 521 | |||
| 522 |
2/2✓ Branch 0 taken 45545 times.
✓ Branch 1 taken 46769 times.
|
92314 | const auto area = Extrusion::area(fvGeometry, scvf); |
| 523 |
2/2✓ Branch 0 taken 45545 times.
✓ Branch 1 taken 46769 times.
|
92314 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 524 |
2/2✓ Branch 0 taken 45545 times.
✓ Branch 1 taken 46769 times.
|
92314 | const auto& insideScv = fvGeometry.scv(insideScvIdx); |
| 525 | 92314 | const auto& insideVolVars = elemVolVars[insideScvIdx]; | |
| 526 | 92314 | const auto wIn = area*computeTpfaTransmissibility(fvGeometry, scvf, insideScv, insideVolVars.permeability(), insideVolVars.extrusionFactor()); | |
| 527 | |||
| 528 | // proceed depending on the interior BC types used | ||
| 529 |
1/2✓ Branch 0 taken 8398 times.
✗ Branch 1 not taken.
|
92314 | const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf); |
| 530 | |||
| 531 | // neumann-coupling | ||
| 532 |
2/4✓ Branch 0 taken 92314 times.
✗ Branch 1 not taken.
✓ Branch 2 taken 92314 times.
✗ Branch 3 not taken.
|
184628 | if (iBcTypes.hasOnlyNeumann()) |
| 533 | { | ||
| 534 | 184628 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 535 | |||
| 536 | // Here we use the square root of the facet extrusion factor | ||
| 537 | // as an approximate average distance from scvf ip to facet center | ||
| 538 | using std::sqrt; | ||
| 539 | 184628 | const auto wFacet = 2.0*area*insideVolVars.extrusionFactor() | |
| 540 | 92314 | /sqrt(facetVolVars.extrusionFactor()) | |
| 541 | 276942 | *vtmv(scvf.unitOuterNormal(), facetVolVars.permeability(), scvf.unitOuterNormal()); | |
| 542 | |||
| 543 | // TODO: check for division by zero?? | ||
| 544 | 92314 | g[0] = wIn/(wIn+wFacet); | |
| 545 | 184628 | tij[Cache::insideTijIdx] = wFacet*g[0]; | |
| 546 | 276942 | tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx]; | |
| 547 | } | ||
| 548 | ✗ | else if (iBcTypes.hasOnlyDirichlet()) | |
| 549 | { | ||
| 550 | ✗ | tij[Cache::insideTijIdx] = wIn; | |
| 551 | ✗ | tij[Cache::facetTijIdx] = -wIn; | |
| 552 | } | ||
| 553 | else | ||
| 554 | ✗ | DUNE_THROW(Dune::NotImplemented, "Interior boundary types other than pure Dirichlet or Neumann"); | |
| 555 | |||
| 556 | 92314 | return tij; | |
| 557 | } | ||
| 558 | }; | ||
| 559 | |||
| 560 | } // end namespace Dumux | ||
| 561 | |||
| 562 | #endif | ||
| 563 |