24#ifndef DUMUX_DISCRETIZATION_CC_TPFA_FACET_COUPLING_FICKS_LAW_HH
25#define DUMUX_DISCRETIZATION_CC_TPFA_FACET_COUPLING_FICKS_LAW_HH
30#include <dune/common/float_cmp.hh>
31#include <dune/common/fvector.hh>
47template<
class TypeTag,
58template<
class TypeTag, ReferenceSystemFormulation referenceSystem = ReferenceSystemFormulation::massAveraged>
68template<
class TypeTag, ReferenceSystemFormulation referenceSystem>
76 using FVElementGeometry =
typename GridGeometry::LocalView;
77 using SubControlVolume =
typename GridGeometry::SubControlVolume;
78 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
81 using GridView =
typename GridGeometry::GridView;
82 using Element =
typename GridView::template Codim<0>::Entity;
88 static const int numPhases = ModelTraits::numFluidPhases();
89 static const int numComponents = ModelTraits::numFluidComponents();
92 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
97 class FacetCouplingFicksLawCache
101 using Filler =
typename ParentType::Cache::Filler;
106 static constexpr int insideTijIdx = 0;
107 static constexpr int outsideTijIdx = 1;
108 static constexpr int facetTijIdx = 2;
111 using DiffusionTransmissibilityContainer = std::array<Scalar, 3>;
114 template<
class Problem,
class ElementVolumeVariables >
115 void updateDiffusion(
const Problem& problem,
116 const Element& element,
117 const FVElementGeometry& fvGeometry,
118 const ElementVolumeVariables& elemVolVars,
119 const SubControlVolumeFace &scvf,
120 unsigned int phaseIdx,
121 unsigned int compIdx)
123 tij_[phaseIdx][compIdx] = Implementation::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
130 Scalar diffusionTij(
unsigned int phaseIdx,
unsigned int compIdx)
const
131 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
134 Scalar diffusionTijInside(
unsigned int phaseIdx,
unsigned int compIdx)
const
135 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
138 Scalar diffusionTijOutside(
unsigned int phaseIdx,
unsigned int compIdx)
const
139 {
return tij_[phaseIdx][compIdx][outsideTijIdx]; }
142 Scalar diffusionTijFacet(
unsigned int phaseIdx,
unsigned int compIdx)
const
143 {
return tij_[phaseIdx][compIdx][facetTijIdx]; }
146 std::array< std::array<DiffusionTransmissibilityContainer, numComponents>, numPhases> tij_;
151 using Cache = FacetCouplingFicksLawCache;
159 {
return referenceSystem; }
162 template<
class Problem,
class ElementVolumeVariables,
class ElementFluxVarsCache >
163 static ComponentFluxVector
flux(
const Problem& problem,
164 const Element& element,
165 const FVElementGeometry& fvGeometry,
166 const ElementVolumeVariables& elemVolVars,
167 const SubControlVolumeFace& scvf,
169 const ElementFluxVarsCache& elemFluxVarsCache)
171 if (!problem.couplingManager().isOnInteriorBoundary(element, scvf))
172 return ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
174 ComponentFluxVector componentFlux(0.0);
175 for (
int compIdx = 0; compIdx < numComponents; compIdx++)
177 if(compIdx == FluidSystem::getMainComponent(phaseIdx))
181 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
182 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
183 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
184 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
187 const Scalar xInside =
massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx);
188 const Scalar xFacet =
massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx);
189 const Scalar xOutside =
massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx);
191 const Scalar rhoInside =
massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx);
193 const Scalar rho = 0.5*(rhoInside + rhoFacet);
195 componentFlux[compIdx] = fluxVarsCache.diffusionTijInside(phaseIdx, compIdx)*xInside
196 + fluxVarsCache.diffusionTijFacet(phaseIdx, compIdx)*xFacet;
198 if (!scvf.boundary())
199 componentFlux[compIdx] += fluxVarsCache.diffusionTijOutside(phaseIdx, compIdx)*xOutside;
200 componentFlux[compIdx] *= rho;
202 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
203 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
206 return componentFlux;
211 template<
class Problem,
class ElementVolumeVariables >
212 static typename Cache::DiffusionTransmissibilityContainer
214 const Element& element,
215 const FVElementGeometry& fvGeometry,
216 const ElementVolumeVariables& elemVolVars,
217 const SubControlVolumeFace& scvf,
218 unsigned int phaseIdx,
unsigned int compIdx)
220 typename Cache::DiffusionTransmissibilityContainer tij;
221 if (!problem.couplingManager().isCoupled(element, scvf))
224 tij[Cache::insideTijIdx] = ParentType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
229 static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(),
"FacetCoupling.Xi", 1.0);
231 const auto insideScvIdx = scvf.insideScvIdx();
232 const auto& insideScv = fvGeometry.scv(insideScvIdx);
233 const auto& insideVolVars = elemVolVars[insideScvIdx];
234 const auto wIn = Extrusion::area(scvf)
236 insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
237 insideVolVars.extrusionFactor());
240 const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf);
243 if (iBcTypes.hasOnlyNeumann())
245 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
246 const auto wFacet = 2.0*Extrusion::area(scvf)*insideVolVars.extrusionFactor()
247 /facetVolVars.extrusionFactor()
248 *
vtmv(scvf.unitOuterNormal(),
249 facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
250 scvf.unitOuterNormal());
259 if (!scvf.boundary())
261 const auto outsideScvIdx = scvf.outsideScvIdx();
262 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
263 const auto wOut = -1.0*Extrusion::area(scvf)
265 outsideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
266 outsideVolVars.extrusionFactor());
268 if ( !Dune::FloatCmp::eq(xi, 1.0, 1e-6) )
272 const Scalar factor = wIn * wFacet / ( wIn * wOut * ( 2.0 * xi - 1.0 ) + wFacet * ( xi * ( wIn + wOut ) + wFacet ) );
273 tij[Cache::insideTijIdx] = factor * ( wOut * xi + wFacet );
274 tij[Cache::outsideTijIdx] = factor * ( wOut * ( 1.0 - xi ) );
275 tij[Cache::facetTijIdx] = factor * ( - wOut - wFacet );
279 tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
280 tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
281 tij[Cache::outsideTijIdx] = 0.0;
286 tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
287 tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
288 tij[Cache::outsideTijIdx] = 0.0;
291 else if (iBcTypes.hasOnlyDirichlet())
293 tij[Cache::insideTijIdx] = wIn;
294 tij[Cache::outsideTijIdx] = 0.0;
295 tij[Cache::facetTijIdx] = -wIn;
298 DUNE_THROW(Dune::NotImplemented,
"Interior boundary types other than pure Dirichlet or Neumann");
308template<
class TypeTag, ReferenceSystemFormulation referenceSystem>
316 using FVElementGeometry =
typename GridGeometry::LocalView;
317 using SubControlVolume =
typename GridGeometry::SubControlVolume;
318 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
321 using GridView =
typename GridGeometry::GridView;
322 using Element =
typename GridView::template Codim<0>::Entity;
328 static const int numPhases = ModelTraits::numFluidPhases();
329 static const int numComponents = ModelTraits::numFluidComponents();
332 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
337 class FacetCouplingFicksLawCache
341 using Filler =
typename ParentType::Cache::Filler;
346 static constexpr int insideTijIdx = 0;
347 static constexpr int facetTijIdx = 1;
350 using DiffusionTransmissibilityContainer = std::array<Scalar, 2>;
353 template<
class Problem,
class ElementVolumeVariables >
354 void updateDiffusion(
const Problem& problem,
355 const Element& element,
356 const FVElementGeometry& fvGeometry,
357 const ElementVolumeVariables& elemVolVars,
358 const SubControlVolumeFace &scvf,
359 unsigned int phaseIdx,
360 unsigned int compIdx)
362 tij_[phaseIdx][compIdx] = Implementation::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
369 Scalar diffusionTij(
unsigned int phaseIdx,
unsigned int compIdx)
const
370 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
373 Scalar diffusionTijInside(
unsigned int phaseIdx,
unsigned int compIdx)
const
374 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
377 Scalar diffusionTijFacet(
unsigned int phaseIdx,
unsigned int compIdx)
const
378 {
return tij_[phaseIdx][compIdx][facetTijIdx]; }
381 std::array< std::array<DiffusionTransmissibilityContainer, numComponents>, numPhases> tij_;
386 using Cache = FacetCouplingFicksLawCache;
394 {
return referenceSystem; }
397 template<
class Problem,
class ElementVolumeVariables,
class ElementFluxVarsCache >
398 static ComponentFluxVector
flux(
const Problem& problem,
399 const Element& element,
400 const FVElementGeometry& fvGeometry,
401 const ElementVolumeVariables& elemVolVars,
402 const SubControlVolumeFace& scvf,
404 const ElementFluxVarsCache& elemFluxVarsCache)
406 if (!problem.couplingManager().isOnInteriorBoundary(element, scvf))
407 return ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
409 ComponentFluxVector componentFlux(0.0);
410 for (
int compIdx = 0; compIdx < numComponents; compIdx++)
412 if(compIdx == FluidSystem::getMainComponent(phaseIdx))
416 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
417 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
418 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
421 const Scalar xInside =
massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx);
422 const Scalar xFacet =
massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx);
424 const Scalar rhoInside =
massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx);
426 const Scalar rho = 0.5*(rhoInside + rhoFacet);
428 componentFlux[compIdx] = rho*(fluxVarsCache.diffusionTijInside(phaseIdx, compIdx)*xInside
429 + fluxVarsCache.diffusionTijFacet(phaseIdx, compIdx)*xFacet);
431 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
432 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
435 return componentFlux;
440 template<
class Problem,
class ElementVolumeVariables >
441 static typename Cache::DiffusionTransmissibilityContainer
443 const Element& element,
444 const FVElementGeometry& fvGeometry,
445 const ElementVolumeVariables& elemVolVars,
446 const SubControlVolumeFace& scvf,
447 unsigned int phaseIdx,
unsigned int compIdx)
449 typename Cache::DiffusionTransmissibilityContainer tij;
450 if (!problem.couplingManager().isCoupled(element, scvf))
453 tij[Cache::insideTijIdx] = ParentType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
458 static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(),
"FacetCoupling.Xi", 1.0);
463 if (Dune::FloatCmp::ne(xi, 1.0, 1e-6))
464 DUNE_THROW(Dune::InvalidStateException,
"Xi != 1.0 cannot be used on surface grids");
466 const auto insideScvIdx = scvf.insideScvIdx();
467 const auto& insideScv = fvGeometry.scv(insideScvIdx);
468 const auto& insideVolVars = elemVolVars[insideScvIdx];
469 const auto wIn = Extrusion::area(scvf)
471 insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
472 insideVolVars.extrusionFactor());
475 const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf);
478 if (iBcTypes.hasOnlyNeumann())
483 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
484 const auto wFacet = 2.0*Extrusion::area(scvf)*insideVolVars.extrusionFactor()
485 /sqrt(facetVolVars.extrusionFactor())
486 *
vtmv(scvf.unitOuterNormal(),
487 facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
488 scvf.unitOuterNormal());
490 tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
491 tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
493 else if (iBcTypes.hasOnlyDirichlet())
495 tij[Cache::insideTijIdx] = wIn;
496 tij[Cache::facetTijIdx] = -wIn;
499 DUNE_THROW(Dune::NotImplemented,
"Interior boundary types other than pure Dirichlet or Neumann");
Define some often used mathematical functions.
The infrastructure to retrieve run-time parameters from Dune::ParameterTrees.
Helper classes to compute the integration elements.
The available discretization methods in Dumux.
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
Tensor::field_type computeTpfaTransmissibility(const SubControlVolumeFace &scvf, const SubControlVolume &scv, const Tensor &T, typename SubControlVolume::Traits::Scalar extrusionFactor)
Free function to evaluate the Tpfa transmissibility associated with the flux (in the form of flux = T...
Definition: tpfa/computetransmissibility.hh:47
VolumeVariables::PrimaryVariables::value_type massOrMoleFraction(const VolumeVariables &volVars, ReferenceSystemFormulation referenceSys, const int phaseIdx, const int compIdx)
returns the mass or mole fraction to be used in Fick's law based on the reference system
Definition: referencesystemformulation.hh:66
VolumeVariables::PrimaryVariables::value_type massOrMolarDensity(const VolumeVariables &volVars, ReferenceSystemFormulation referenceSys, const int phaseIdx)
evaluates the density to be used in Fick's law based on the reference system
Definition: referencesystemformulation.hh:55
ReferenceSystemFormulation
The formulations available for Fick's law related to the reference system.
Definition: referencesystemformulation.hh:45
Dune::DenseMatrix< MAT >::value_type vtmv(const Dune::DenseVector< V1 > &v1, const Dune::DenseMatrix< MAT > &M, const Dune::DenseVector< V2 > &v2)
Evaluates the scalar product of a vector v2, projected by a matrix M, with a vector v1.
Definition: math.hh:863
typename Extrusion< T >::type Extrusion_t
Convenience alias for obtaining the extrusion type.
Definition: extrusion.hh:177
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property
Definition: propertysystem.hh:150
forward declaration of the method-specific implemetation
Definition: flux/box/fickslaw.hh:44
Fick's law for cell-centered finite volume schemes with two-point flux approximation.
Definition: flux/cctpfa/fickslaw.hh:51
TpfaFicksLawCache Cache
state the type for the corresponding cache and its filler
Definition: flux/cctpfa/fickslaw.hh:129
Forward declaration of the implementation.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:50
Specialization of CCTpfaFacetCouplingFicksLawImpl for dim=dimWorld.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:71
static ComponentFluxVector flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, int phaseIdx, const ElementFluxVarsCache &elemFluxVarsCache)
Compute the diffusive fluxes.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:163
static Cache::DiffusionTransmissibilityContainer calculateTransmissibility(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, unsigned int phaseIdx, unsigned int compIdx)
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:213
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
Return the reference system.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:158
Specialization of CCTpfaFacetCouplingFicksLawImpl for dim<dimWorld.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:311
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
Return the reference system.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:393
static Cache::DiffusionTransmissibilityContainer calculateTransmissibility(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, unsigned int phaseIdx, unsigned int compIdx)
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:442
static ComponentFluxVector flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, int phaseIdx, const ElementFluxVarsCache &elemFluxVarsCache)
Compute the diffusive fluxes.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:398
Declares all properties used in Dumux.
Free functions to evaluate the transmissibilities associated with flux evaluations across sub-control...
Fick's law for cell-centered finite volume schemes with two-point flux approximation.