| 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-FileCopyrightText: 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::CCTpfaFacetCouplingFicksLaw | ||
| 11 | */ | ||
| 12 | #ifndef DUMUX_DISCRETIZATION_CC_TPFA_FACET_COUPLING_FICKS_LAW_HH | ||
| 13 | #define DUMUX_DISCRETIZATION_CC_TPFA_FACET_COUPLING_FICKS_LAW_HH | ||
| 14 | |||
| 15 | #include <array> | ||
| 16 | #include <cmath> | ||
| 17 | |||
| 18 | #include <dune/common/float_cmp.hh> | ||
| 19 | #include <dune/common/fvector.hh> | ||
| 20 | |||
| 21 | #include <dumux/common/math.hh> | ||
| 22 | #include <dumux/common/parameters.hh> | ||
| 23 | #include <dumux/common/properties.hh> | ||
| 24 | |||
| 25 | #include <dumux/discretization/method.hh> | ||
| 26 | #include <dumux/discretization/extrusion.hh> | ||
| 27 | #include <dumux/discretization/cellcentered/tpfa/computetransmissibility.hh> | ||
| 28 | |||
| 29 | #include <dumux/flux/referencesystemformulation.hh> | ||
| 30 | #include <dumux/flux/cctpfa/fickslaw.hh> | ||
| 31 | |||
| 32 | namespace Dumux { | ||
| 33 | |||
| 34 | //! Forward declaration of the implementation | ||
| 35 | template<class TypeTag, | ||
| 36 | ReferenceSystemFormulation referenceSystem, | ||
| 37 | bool isNetwork> | ||
| 38 | class CCTpfaFacetCouplingFicksLawImpl; | ||
| 39 | |||
| 40 | /*! | ||
| 41 | * \ingroup FacetCoupling | ||
| 42 | * \brief Fick's law for cell-centered finite volume schemes with two-point flux approximation | ||
| 43 | * in the context of coupled models where the coupling occurs across the facets of the bulk | ||
| 44 | * domain elements with a lower-dimensional domain defined on these facets. | ||
| 45 | */ | ||
| 46 | template<class TypeTag, ReferenceSystemFormulation referenceSystem = ReferenceSystemFormulation::massAveraged> | ||
| 47 | using CCTpfaFacetCouplingFicksLaw = | ||
| 48 | CCTpfaFacetCouplingFicksLawImpl< TypeTag, referenceSystem, | ||
| 49 | ( int(GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimension) | ||
| 50 | < int(GetPropType<TypeTag, Properties::GridGeometry>::GridView::dimensionworld) ) >; | ||
| 51 | |||
| 52 | /*! | ||
| 53 | * \ingroup FacetCoupling | ||
| 54 | * \brief Specialization of CCTpfaFacetCouplingFicksLawImpl for dim=dimWorld | ||
| 55 | */ | ||
| 56 | template<class TypeTag, ReferenceSystemFormulation referenceSystem> | ||
| 57 | class CCTpfaFacetCouplingFicksLawImpl<TypeTag, referenceSystem, /*isNetwork*/false> | ||
| 58 | : public FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem> | ||
| 59 | { | ||
| 60 | using Implementation = CCTpfaFacetCouplingFicksLawImpl<TypeTag, referenceSystem, false>; | ||
| 61 | using ParentType = FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem>; | ||
| 62 | |||
| 63 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 64 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 65 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 66 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 67 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 68 | |||
| 69 | using GridView = typename GridGeometry::GridView; | ||
| 70 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 71 | |||
| 72 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 73 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 74 | using BalanceEqOpts = GetPropType<TypeTag, Properties::BalanceEqOpts>; | ||
| 75 | |||
| 76 | static const int numPhases = ModelTraits::numFluidPhases(); | ||
| 77 | static const int numComponents = ModelTraits::numFluidComponents(); | ||
| 78 | |||
| 79 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 80 | using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>; | ||
| 81 | |||
| 82 | /*! | ||
| 83 | * \brief The cache class used with this specialization of Fick's law. | ||
| 84 | */ | ||
| 85 | class FacetCouplingFicksLawCache | ||
| 86 | { | ||
| 87 | public: | ||
| 88 | //! export the corresponding filler class | ||
| 89 | using Filler = typename ParentType::Cache::Filler; | ||
| 90 | |||
| 91 | //! we store the transmissibilities associated with the interior | ||
| 92 | //! cell, outside cell, and the fracture facet in an array. Access | ||
| 93 | //! to this array should be done using the following indices: | ||
| 94 | static constexpr int insideTijIdx = 0; | ||
| 95 | static constexpr int outsideTijIdx = 1; | ||
| 96 | static constexpr int facetTijIdx = 2; | ||
| 97 | |||
| 98 | //! Export transmissibility storage type | ||
| 99 | using DiffusionTransmissibilityContainer = std::array<Scalar, 3>; | ||
| 100 | |||
| 101 | //! update subject to a given problem | ||
| 102 | template< class Problem, class ElementVolumeVariables > | ||
| 103 | 2656312 | void updateDiffusion(const Problem& problem, | |
| 104 | const Element& element, | ||
| 105 | const FVElementGeometry& fvGeometry, | ||
| 106 | const ElementVolumeVariables& elemVolVars, | ||
| 107 | const SubControlVolumeFace &scvf, | ||
| 108 | unsigned int phaseIdx, | ||
| 109 | unsigned int compIdx) | ||
| 110 | { | ||
| 111 | 2656312 | tij_[phaseIdx][compIdx] = Implementation::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx); | |
| 112 | } | ||
| 113 | |||
| 114 | //! We use the same name as in the TpfaFicksLawCache so | ||
| 115 | //! that this cache and the law implementation for non-coupled | ||
| 116 | //! models can be reused here on facets that do not lie on an | ||
| 117 | //! interior boundary, i.e. do not coincide with a facet element | ||
| 118 | 1850640 | Scalar diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const | |
| 119 | 1850640 | { return tij_[phaseIdx][compIdx][insideTijIdx]; } | |
| 120 | |||
| 121 | //! returns the transmissibility associated with the inside cell | ||
| 122 | 55440 | Scalar diffusionTijInside(unsigned int phaseIdx, unsigned int compIdx) const | |
| 123 | 55440 | { return tij_[phaseIdx][compIdx][insideTijIdx]; } | |
| 124 | |||
| 125 | //! returns the transmissibility associated with the outside cell | ||
| 126 | 55440 | Scalar diffusionTijOutside(unsigned int phaseIdx, unsigned int compIdx) const | |
| 127 | 55440 | { return tij_[phaseIdx][compIdx][outsideTijIdx]; } | |
| 128 | |||
| 129 | //! returns the transmissibility associated with the outside cell | ||
| 130 | 55440 | Scalar diffusionTijFacet(unsigned int phaseIdx, unsigned int compIdx) const | |
| 131 | 55440 | { return tij_[phaseIdx][compIdx][facetTijIdx]; } | |
| 132 | |||
| 133 | private: | ||
| 134 | std::array< std::array<DiffusionTransmissibilityContainer, numComponents>, numPhases> tij_; | ||
| 135 | }; | ||
| 136 | |||
| 137 | public: | ||
| 138 | //! export the type for the corresponding cache | ||
| 139 | using Cache = FacetCouplingFicksLawCache; | ||
| 140 | |||
| 141 | using DiscretizationMethod = DiscretizationMethods::CCTpfa; | ||
| 142 | //! state the discretization method this implementation belongs to | ||
| 143 | static constexpr DiscretizationMethod discMethod{}; | ||
| 144 | |||
| 145 | //! Return the reference system | ||
| 146 | static constexpr ReferenceSystemFormulation referenceSystemFormulation() | ||
| 147 | { return referenceSystem; } | ||
| 148 | |||
| 149 | //! Compute the diffusive fluxes | ||
| 150 | template< class Problem, class ElementVolumeVariables, class ElementFluxVarsCache > | ||
| 151 |
2/2✓ Branch 0 taken 1850640 times.
✓ Branch 1 taken 55440 times.
|
1906080 | static ComponentFluxVector flux(const Problem& problem, |
| 152 | const Element& element, | ||
| 153 | const FVElementGeometry& fvGeometry, | ||
| 154 | const ElementVolumeVariables& elemVolVars, | ||
| 155 | const SubControlVolumeFace& scvf, | ||
| 156 | int phaseIdx, | ||
| 157 | const ElementFluxVarsCache& elemFluxVarsCache) | ||
| 158 | { | ||
| 159 |
2/2✓ Branch 0 taken 1850640 times.
✓ Branch 1 taken 55440 times.
|
1906080 | if (!problem.couplingManager().isOnInteriorBoundary(element, scvf)) |
| 160 | 1850640 | return ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache); | |
| 161 | |||
| 162 | 55440 | ComponentFluxVector componentFlux(0.0); | |
| 163 |
2/2✓ Branch 0 taken 110880 times.
✓ Branch 1 taken 55440 times.
|
166320 | for (int compIdx = 0; compIdx < numComponents; compIdx++) |
| 164 | { | ||
| 165 |
2/2✓ Branch 0 taken 55440 times.
✓ Branch 1 taken 55440 times.
|
110880 | if(compIdx == FluidSystem::getMainComponent(phaseIdx)) |
| 166 | 55440 | continue; | |
| 167 | |||
| 168 | // get inside/outside volume variables | ||
| 169 | 55440 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 170 | 55440 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 171 | 55440 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 172 | 55440 | const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()]; | |
| 173 | |||
| 174 | // the inside and outside mass/mole fractions fractions | ||
| 175 | 55440 | const Scalar xInside = massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx); | |
| 176 | 55440 | const Scalar xFacet = massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx); | |
| 177 |
1/2✓ Branch 0 taken 55440 times.
✗ Branch 1 not taken.
|
55440 | const Scalar xOutside = massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx); |
| 178 | |||
| 179 | 55440 | const Scalar rhoInside = massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx); | |
| 180 | 55440 | const Scalar rhoFacet = massOrMolarDensity(facetVolVars, referenceSystem, phaseIdx); | |
| 181 | 55440 | const Scalar rho = 0.5*(rhoInside + rhoFacet); | |
| 182 | |||
| 183 |
1/2✓ Branch 0 taken 55440 times.
✗ Branch 1 not taken.
|
55440 | componentFlux[compIdx] = fluxVarsCache.diffusionTijInside(phaseIdx, compIdx)*xInside |
| 184 |
1/2✓ Branch 0 taken 55440 times.
✗ Branch 1 not taken.
|
55440 | + fluxVarsCache.diffusionTijFacet(phaseIdx, compIdx)*xFacet; |
| 185 | |||
| 186 |
1/2✓ Branch 0 taken 55440 times.
✗ Branch 1 not taken.
|
55440 | if (!scvf.boundary()) |
| 187 | 55440 | componentFlux[compIdx] += fluxVarsCache.diffusionTijOutside(phaseIdx, compIdx)*xOutside; | |
| 188 | 55440 | componentFlux[compIdx] *= rho; | |
| 189 | |||
| 190 | if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem()) | ||
| 191 | 55440 | componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx]; | |
| 192 | } | ||
| 193 | |||
| 194 | 55440 | return componentFlux; | |
| 195 | } | ||
| 196 | |||
| 197 | // The flux variables cache has to be bound to an element prior to flux calculations | ||
| 198 | // During the binding, the transmissibility will be computed and stored using the method below. | ||
| 199 | template< class Problem, class ElementVolumeVariables > | ||
| 200 | static typename Cache::DiffusionTransmissibilityContainer | ||
| 201 |
2/2✓ Branch 0 taken 2571640 times.
✓ Branch 1 taken 84672 times.
|
2656312 | calculateTransmissibility(const Problem& problem, |
| 202 | const Element& element, | ||
| 203 | const FVElementGeometry& fvGeometry, | ||
| 204 | const ElementVolumeVariables& elemVolVars, | ||
| 205 | const SubControlVolumeFace& scvf, | ||
| 206 | unsigned int phaseIdx, unsigned int compIdx) | ||
| 207 | { | ||
| 208 | typename Cache::DiffusionTransmissibilityContainer tij; | ||
| 209 |
2/2✓ Branch 0 taken 2571640 times.
✓ Branch 1 taken 84672 times.
|
2656312 | if (!problem.couplingManager().isCoupled(element, scvf)) |
| 210 | { | ||
| 211 | //! use the standard Fick's law and only compute one transmissibility | ||
| 212 | 2571640 | tij[Cache::insideTijIdx] = ParentType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx); | |
| 213 | 2571640 | return tij; | |
| 214 | } | ||
| 215 | |||
| 216 | //! xi factor for the coupling conditions | ||
| 217 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 84670 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
84672 | static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0); |
| 218 | |||
| 219 |
2/2✓ Branch 0 taken 42336 times.
✓ Branch 1 taken 42336 times.
|
84672 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 220 | 84672 | const auto& insideScv = fvGeometry.scv(insideScvIdx); | |
| 221 | 84672 | const auto& insideVolVars = elemVolVars[insideScvIdx]; | |
| 222 | 169344 | const auto wIn = Extrusion::area(fvGeometry, scvf) | |
| 223 | 84672 | *computeTpfaTransmissibility(fvGeometry, scvf, insideScv, | |
| 224 | insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx), | ||
| 225 | insideVolVars.extrusionFactor()); | ||
| 226 | |||
| 227 | // proceed depending on the interior BC types used | ||
| 228 |
1/2✓ Branch 0 taken 84672 times.
✗ Branch 1 not taken.
|
84672 | const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf); |
| 229 | |||
| 230 | // neumann-coupling | ||
| 231 |
1/2✓ Branch 0 taken 84672 times.
✗ Branch 1 not taken.
|
84672 | if (iBcTypes.hasOnlyNeumann()) |
| 232 | { | ||
| 233 | 84672 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 234 | 84672 | const auto wFacet = 2.0*Extrusion::area(fvGeometry, scvf)*insideVolVars.extrusionFactor() | |
| 235 | 84672 | /facetVolVars.extrusionFactor() | |
| 236 |
1/2✓ Branch 0 taken 84672 times.
✗ Branch 1 not taken.
|
84672 | *vtmv(scvf.unitOuterNormal(), |
| 237 | facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx), | ||
| 238 | 84672 | scvf.unitOuterNormal()); | |
| 239 | |||
| 240 | // The fluxes across this face and the outside face can be expressed in matrix form: | ||
| 241 | // \f$\mathbf{C} \bar{\mathbf{u}} + \mathbf{D} \mathbf{u} + \mathbf{E} \mathbf{u}_\gamma\f$, | ||
| 242 | // where \f$\gamma$\f denotes the domain living on the facets and \f$\bar{\mathbf{u}}$\f are | ||
| 243 | // intermediate face unknowns in the matrix domain. Equivalently, flux continuity reads: | ||
| 244 | // \f$\mathbf{A} \bar{\mathbf{u}} = \mathbf{B} \mathbf{u} + \mathbf{M} \mathbf{u}_\gamma\f$. | ||
| 245 | // Combining the two, we can eliminate the intermediate unknowns and compute the transmissibilities | ||
| 246 | // that allow the description of the fluxes as functions of the cell and Dirichlet mass/mole fractions only. | ||
| 247 |
1/2✓ Branch 0 taken 84672 times.
✗ Branch 1 not taken.
|
84672 | if (!scvf.boundary()) |
| 248 | { | ||
| 249 | 84672 | const auto outsideScvIdx = scvf.outsideScvIdx(); | |
| 250 | 84672 | const auto& outsideVolVars = elemVolVars[outsideScvIdx]; | |
| 251 | 169344 | const auto wOut = -1.0*Extrusion::area(fvGeometry, scvf) | |
| 252 |
2/2✓ Branch 1 taken 42336 times.
✓ Branch 2 taken 42336 times.
|
84672 | *computeTpfaTransmissibility(fvGeometry, scvf, fvGeometry.scv(outsideScvIdx), |
| 253 | outsideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx), | ||
| 254 | outsideVolVars.extrusionFactor()); | ||
| 255 | |||
| 256 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 84672 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 84672 times.
|
84672 | if ( !Dune::FloatCmp::eq(xi, 1.0, 1e-6) ) |
| 257 | { | ||
| 258 | // optimized implementation: factorization obtained using sympy | ||
| 259 | // see CCTpfaFacetCouplingDarcysLaw for more details | ||
| 260 | ✗ | const Scalar factor = wIn * wFacet / ( wIn * wOut * ( 2.0 * xi - 1.0 ) + wFacet * ( xi * ( wIn + wOut ) + wFacet ) ); | |
| 261 | ✗ | tij[Cache::insideTijIdx] = factor * ( wOut * xi + wFacet ); | |
| 262 | ✗ | tij[Cache::outsideTijIdx] = factor * ( wOut * ( 1.0 - xi ) ); | |
| 263 | ✗ | tij[Cache::facetTijIdx] = factor * ( - wOut - wFacet ); | |
| 264 | } | ||
| 265 | else | ||
| 266 | { | ||
| 267 | 84672 | tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet); | |
| 268 | 84672 | tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx]; | |
| 269 | 84672 | tij[Cache::outsideTijIdx] = 0.0; | |
| 270 | } | ||
| 271 | } | ||
| 272 | else | ||
| 273 | { | ||
| 274 | ✗ | tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet); | |
| 275 | ✗ | tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx]; | |
| 276 | ✗ | tij[Cache::outsideTijIdx] = 0.0; | |
| 277 | } | ||
| 278 | } | ||
| 279 | ✗ | else if (iBcTypes.hasOnlyDirichlet()) | |
| 280 | { | ||
| 281 | ✗ | tij[Cache::insideTijIdx] = wIn; | |
| 282 | ✗ | tij[Cache::outsideTijIdx] = 0.0; | |
| 283 | ✗ | tij[Cache::facetTijIdx] = -wIn; | |
| 284 | } | ||
| 285 | else | ||
| 286 | ✗ | DUNE_THROW(Dune::NotImplemented, "Interior boundary types other than pure Dirichlet or Neumann"); | |
| 287 | |||
| 288 | return tij; | ||
| 289 | } | ||
| 290 | }; | ||
| 291 | |||
| 292 | /*! | ||
| 293 | * \ingroup FacetCoupling | ||
| 294 | * \brief Specialization of CCTpfaFacetCouplingFicksLawImpl for dim<dimWorld | ||
| 295 | */ | ||
| 296 | template<class TypeTag, ReferenceSystemFormulation referenceSystem> | ||
| 297 | class CCTpfaFacetCouplingFicksLawImpl<TypeTag, referenceSystem, /*isNetwork*/true> | ||
| 298 | : public FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem> | ||
| 299 | { | ||
| 300 | using Implementation = CCTpfaFacetCouplingFicksLawImpl<TypeTag, referenceSystem, true>; | ||
| 301 | using ParentType = FicksLawImplementation<TypeTag, DiscretizationMethods::CCTpfa, referenceSystem>; | ||
| 302 | |||
| 303 | using GridGeometry = GetPropType<TypeTag, Properties::GridGeometry>; | ||
| 304 | using FVElementGeometry = typename GridGeometry::LocalView; | ||
| 305 | using SubControlVolume = typename GridGeometry::SubControlVolume; | ||
| 306 | using SubControlVolumeFace = typename GridGeometry::SubControlVolumeFace; | ||
| 307 | using Extrusion = Extrusion_t<GridGeometry>; | ||
| 308 | |||
| 309 | using GridView = typename GridGeometry::GridView; | ||
| 310 | using Element = typename GridView::template Codim<0>::Entity; | ||
| 311 | |||
| 312 | using FluidSystem = GetPropType<TypeTag, Properties::FluidSystem>; | ||
| 313 | using ModelTraits = GetPropType<TypeTag, Properties::ModelTraits>; | ||
| 314 | using BalanceEqOpts = GetPropType<TypeTag, Properties::BalanceEqOpts>; | ||
| 315 | |||
| 316 | static const int numPhases = ModelTraits::numFluidPhases(); | ||
| 317 | static const int numComponents = ModelTraits::numFluidComponents(); | ||
| 318 | |||
| 319 | using Scalar = GetPropType<TypeTag, Properties::Scalar>; | ||
| 320 | using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>; | ||
| 321 | |||
| 322 | /*! | ||
| 323 | * \brief The cache class used with this specialization of Fick's law. | ||
| 324 | */ | ||
| 325 | class FacetCouplingFicksLawCache | ||
| 326 | { | ||
| 327 | public: | ||
| 328 | //! export the corresponding filler class | ||
| 329 | using Filler = typename ParentType::Cache::Filler; | ||
| 330 | |||
| 331 | //! we store the transmissibilities associated with the interior | ||
| 332 | //! cell and the fracture facet in an array. Access to this array | ||
| 333 | //! should be done using the following indices: | ||
| 334 | static constexpr int insideTijIdx = 0; | ||
| 335 | static constexpr int facetTijIdx = 1; | ||
| 336 | |||
| 337 | //! Export transmissibility storage type | ||
| 338 | using DiffusionTransmissibilityContainer = std::array<Scalar, 2>; | ||
| 339 | |||
| 340 | //! update subject to a given problem | ||
| 341 | template< class Problem, class ElementVolumeVariables > | ||
| 342 | 4495956 | void updateDiffusion(const Problem& problem, | |
| 343 | const Element& element, | ||
| 344 | const FVElementGeometry& fvGeometry, | ||
| 345 | const ElementVolumeVariables& elemVolVars, | ||
| 346 | const SubControlVolumeFace &scvf, | ||
| 347 | unsigned int phaseIdx, | ||
| 348 | unsigned int compIdx) | ||
| 349 | { | ||
| 350 | 4495956 | tij_[phaseIdx][compIdx] = Implementation::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx); | |
| 351 | } | ||
| 352 | |||
| 353 | //! We use the same name as in the TpfaFicksLawCache so | ||
| 354 | //! that this cache and the law implementation for non-coupled | ||
| 355 | //! models can be reused here on facets that do not lie on an | ||
| 356 | //! interior boundary, i.e. do not coincide with a facet element | ||
| 357 | 3157304 | Scalar diffusionTij(unsigned int phaseIdx, unsigned int compIdx) const | |
| 358 | 3157304 | { return tij_[phaseIdx][compIdx][insideTijIdx]; } | |
| 359 | |||
| 360 | //! returns the transmissibility associated with the inside cell | ||
| 361 | 94584 | Scalar diffusionTijInside(unsigned int phaseIdx, unsigned int compIdx) const | |
| 362 | 94584 | { return tij_[phaseIdx][compIdx][insideTijIdx]; } | |
| 363 | |||
| 364 | //! returns the transmissibility associated with the outside cell | ||
| 365 | 94584 | Scalar diffusionTijFacet(unsigned int phaseIdx, unsigned int compIdx) const | |
| 366 | 94584 | { return tij_[phaseIdx][compIdx][facetTijIdx]; } | |
| 367 | |||
| 368 | private: | ||
| 369 | std::array< std::array<DiffusionTransmissibilityContainer, numComponents>, numPhases> tij_; | ||
| 370 | }; | ||
| 371 | |||
| 372 | public: | ||
| 373 | //! export the type for the corresponding cache | ||
| 374 | using Cache = FacetCouplingFicksLawCache; | ||
| 375 | |||
| 376 | using DiscretizationMethod = DiscretizationMethods::CCTpfa; | ||
| 377 | //! state the discretization method this implementation belongs to | ||
| 378 | static constexpr DiscretizationMethod discMethod{}; | ||
| 379 | |||
| 380 | //! Return the reference system | ||
| 381 | static constexpr ReferenceSystemFormulation referenceSystemFormulation() | ||
| 382 | { return referenceSystem; } | ||
| 383 | |||
| 384 | //! Compute the diffusive fluxes | ||
| 385 | template< class Problem, class ElementVolumeVariables, class ElementFluxVarsCache > | ||
| 386 |
2/2✓ Branch 0 taken 3157304 times.
✓ Branch 1 taken 94584 times.
|
3251888 | static ComponentFluxVector flux(const Problem& problem, |
| 387 | const Element& element, | ||
| 388 | const FVElementGeometry& fvGeometry, | ||
| 389 | const ElementVolumeVariables& elemVolVars, | ||
| 390 | const SubControlVolumeFace& scvf, | ||
| 391 | int phaseIdx, | ||
| 392 | const ElementFluxVarsCache& elemFluxVarsCache) | ||
| 393 | { | ||
| 394 |
2/2✓ Branch 0 taken 3157304 times.
✓ Branch 1 taken 94584 times.
|
3251888 | if (!problem.couplingManager().isOnInteriorBoundary(element, scvf)) |
| 395 | 3157304 | return ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache); | |
| 396 | |||
| 397 | 94584 | ComponentFluxVector componentFlux(0.0); | |
| 398 |
2/2✓ Branch 0 taken 189168 times.
✓ Branch 1 taken 94584 times.
|
283752 | for (int compIdx = 0; compIdx < numComponents; compIdx++) |
| 399 | { | ||
| 400 |
2/2✓ Branch 0 taken 94584 times.
✓ Branch 1 taken 94584 times.
|
189168 | if(compIdx == FluidSystem::getMainComponent(phaseIdx)) |
| 401 | 94584 | continue; | |
| 402 | |||
| 403 | // get inside/outside volume variables | ||
| 404 | 94584 | const auto& fluxVarsCache = elemFluxVarsCache[scvf]; | |
| 405 | 94584 | const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()]; | |
| 406 | 94584 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 407 | |||
| 408 | // the inside and outside mass/mole fractions fractions | ||
| 409 | 94584 | const Scalar xInside = massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx); | |
| 410 | 94584 | const Scalar xFacet = massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx); | |
| 411 | |||
| 412 | 94584 | const Scalar rhoInside = massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx); | |
| 413 | 94584 | const Scalar rhoFacet = massOrMolarDensity(facetVolVars, referenceSystem, phaseIdx); | |
| 414 | 94584 | const Scalar rho = 0.5*(rhoInside + rhoFacet); | |
| 415 | |||
| 416 | 94584 | componentFlux[compIdx] = rho*(fluxVarsCache.diffusionTijInside(phaseIdx, compIdx)*xInside | |
| 417 | 94584 | + fluxVarsCache.diffusionTijFacet(phaseIdx, compIdx)*xFacet); | |
| 418 | |||
| 419 | if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem()) | ||
| 420 | 94584 | componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx]; | |
| 421 | } | ||
| 422 | |||
| 423 | 94584 | return componentFlux; | |
| 424 | } | ||
| 425 | |||
| 426 | // The flux variables cache has to be bound to an element prior to flux calculations | ||
| 427 | // During the binding, the transmissibility will be computed and stored using the method below. | ||
| 428 | template< class Problem, class ElementVolumeVariables > | ||
| 429 | static typename Cache::DiffusionTransmissibilityContainer | ||
| 430 |
2/2✓ Branch 0 taken 4350916 times.
✓ Branch 1 taken 145040 times.
|
4495956 | calculateTransmissibility(const Problem& problem, |
| 431 | const Element& element, | ||
| 432 | const FVElementGeometry& fvGeometry, | ||
| 433 | const ElementVolumeVariables& elemVolVars, | ||
| 434 | const SubControlVolumeFace& scvf, | ||
| 435 | unsigned int phaseIdx, unsigned int compIdx) | ||
| 436 | { | ||
| 437 | typename Cache::DiffusionTransmissibilityContainer tij; | ||
| 438 |
2/2✓ Branch 0 taken 4350916 times.
✓ Branch 1 taken 145040 times.
|
4495956 | if (!problem.couplingManager().isCoupled(element, scvf)) |
| 439 | { | ||
| 440 | //! use the standard Fick's law and only compute one transmissibility | ||
| 441 | 4350916 | tij[Cache::insideTijIdx] = ParentType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx); | |
| 442 | 4350916 | return tij; | |
| 443 | } | ||
| 444 | |||
| 445 | //! xi factor for the coupling conditions | ||
| 446 |
4/6✓ Branch 0 taken 2 times.
✓ Branch 1 taken 145038 times.
✓ Branch 3 taken 2 times.
✗ Branch 4 not taken.
✓ Branch 6 taken 2 times.
✗ Branch 7 not taken.
|
145040 | static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(), "FacetCoupling.Xi", 1.0); |
| 447 | |||
| 448 | // On surface grids only xi = 1.0 can be used, as the coupling condition | ||
| 449 | // for xi != 1.0 does not generalize for surface grids where the normal | ||
| 450 | // vectors of the inside/outside elements have different orientations. | ||
| 451 |
2/4✗ Branch 0 not taken.
✓ Branch 1 taken 145040 times.
✗ Branch 2 not taken.
✓ Branch 3 taken 145040 times.
|
145040 | if (Dune::FloatCmp::ne(xi, 1.0, 1e-6)) |
| 452 |
2/22✗ Branch 1 not taken.
✗ Branch 2 not taken.
✗ Branch 4 not taken.
✗ Branch 5 not taken.
✗ Branch 7 not taken.
✗ Branch 8 not taken.
✗ Branch 10 not taken.
✗ Branch 11 not taken.
✗ Branch 13 not taken.
✗ Branch 14 not taken.
✗ Branch 16 not taken.
✗ Branch 17 not taken.
✗ Branch 19 not taken.
✗ Branch 20 not taken.
✗ Branch 22 not taken.
✗ Branch 23 not taken.
✗ Branch 25 not taken.
✗ Branch 26 not taken.
✗ Branch 28 not taken.
✗ Branch 29 not taken.
✓ Branch 33 taken 72520 times.
✓ Branch 34 taken 72520 times.
|
145040 | DUNE_THROW(Dune::InvalidStateException, "Xi != 1.0 cannot be used on surface grids"); |
| 453 | |||
| 454 |
2/2✓ Branch 0 taken 72520 times.
✓ Branch 1 taken 72520 times.
|
145040 | const auto insideScvIdx = scvf.insideScvIdx(); |
| 455 | 145040 | const auto& insideScv = fvGeometry.scv(insideScvIdx); | |
| 456 | 145040 | const auto& insideVolVars = elemVolVars[insideScvIdx]; | |
| 457 | 290080 | const auto wIn = Extrusion::area(fvGeometry, scvf) | |
| 458 | 145040 | *computeTpfaTransmissibility(fvGeometry, scvf, insideScv, | |
| 459 | insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx), | ||
| 460 | insideVolVars.extrusionFactor()); | ||
| 461 | |||
| 462 | // proceed depending on the interior BC types used | ||
| 463 |
1/2✓ Branch 0 taken 145040 times.
✗ Branch 1 not taken.
|
145040 | const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf); |
| 464 | |||
| 465 | // neumann-coupling | ||
| 466 |
1/2✓ Branch 0 taken 145040 times.
✗ Branch 1 not taken.
|
145040 | if (iBcTypes.hasOnlyNeumann()) |
| 467 | { | ||
| 468 | // Here we use the square root of the facet extrusion factor | ||
| 469 | // as an approximate average distance from scvf ip to facet center | ||
| 470 | using std::sqrt; | ||
| 471 | 145040 | const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf); | |
| 472 | 145040 | const auto wFacet = 2.0*Extrusion::area(fvGeometry, scvf)*insideVolVars.extrusionFactor() | |
| 473 | 145040 | /sqrt(facetVolVars.extrusionFactor()) | |
| 474 | 290080 | *vtmv(scvf.unitOuterNormal(), | |
| 475 | facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx), | ||
| 476 | 145040 | scvf.unitOuterNormal()); | |
| 477 | |||
| 478 | 145040 | tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet); | |
| 479 | 145040 | tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx]; | |
| 480 | } | ||
| 481 | ✗ | else if (iBcTypes.hasOnlyDirichlet()) | |
| 482 | { | ||
| 483 | ✗ | tij[Cache::insideTijIdx] = wIn; | |
| 484 | ✗ | tij[Cache::facetTijIdx] = -wIn; | |
| 485 | } | ||
| 486 | else | ||
| 487 | ✗ | DUNE_THROW(Dune::NotImplemented, "Interior boundary types other than pure Dirichlet or Neumann"); | |
| 488 | |||
| 489 | 145040 | return tij; | |
| 490 | } | ||
| 491 | }; | ||
| 492 | |||
| 493 | } // end namespace Dumux | ||
| 494 | |||
| 495 | #endif | ||
| 496 |