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>
46template<
class TypeTag,
57template<
class TypeTag, ReferenceSystemFormulation referenceSystem = ReferenceSystemFormulation::massAveraged>
67template<
class TypeTag, ReferenceSystemFormulation referenceSystem>
75 using FVElementGeometry =
typename GridGeometry::LocalView;
76 using SubControlVolume =
typename GridGeometry::SubControlVolume;
77 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
79 using GridView =
typename GridGeometry::GridView;
80 using Element =
typename GridView::template Codim<0>::Entity;
86 static const int numPhases = ModelTraits::numFluidPhases();
87 static const int numComponents = ModelTraits::numFluidComponents();
90 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
95 class FacetCouplingFicksLawCache
99 using Filler =
typename ParentType::Cache::Filler;
104 static constexpr int insideTijIdx = 0;
105 static constexpr int outsideTijIdx = 1;
106 static constexpr int facetTijIdx = 2;
109 using DiffusionTransmissibilityContainer = std::array<Scalar, 3>;
112 template<
class Problem,
class ElementVolumeVariables >
113 void updateDiffusion(
const Problem& problem,
114 const Element& element,
115 const FVElementGeometry& fvGeometry,
116 const ElementVolumeVariables& elemVolVars,
117 const SubControlVolumeFace &scvf,
118 unsigned int phaseIdx,
119 unsigned int compIdx)
121 tij_[phaseIdx][compIdx] = Implementation::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
128 Scalar diffusionTij(
unsigned int phaseIdx,
unsigned int compIdx)
const
129 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
132 Scalar diffusionTijInside(
unsigned int phaseIdx,
unsigned int compIdx)
const
133 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
136 Scalar diffusionTijOutside(
unsigned int phaseIdx,
unsigned int compIdx)
const
137 {
return tij_[phaseIdx][compIdx][outsideTijIdx]; }
140 Scalar diffusionTijFacet(
unsigned int phaseIdx,
unsigned int compIdx)
const
141 {
return tij_[phaseIdx][compIdx][facetTijIdx]; }
144 std::array< std::array<DiffusionTransmissibilityContainer, numComponents>, numPhases> tij_;
149 using Cache = FacetCouplingFicksLawCache;
156 {
return referenceSystem; }
159 template<
class Problem,
class ElementVolumeVariables,
class ElementFluxVarsCache >
160 static ComponentFluxVector
flux(
const Problem& problem,
161 const Element& element,
162 const FVElementGeometry& fvGeometry,
163 const ElementVolumeVariables& elemVolVars,
164 const SubControlVolumeFace& scvf,
166 const ElementFluxVarsCache& elemFluxVarsCache)
168 if (!problem.couplingManager().isOnInteriorBoundary(element, scvf))
169 return ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
171 ComponentFluxVector componentFlux(0.0);
172 for (
int compIdx = 0; compIdx < numComponents; compIdx++)
174 if(compIdx == FluidSystem::getMainComponent(phaseIdx))
178 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
179 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
180 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
181 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
184 const Scalar xInside =
massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx);
185 const Scalar xFacet =
massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx);
186 const Scalar xOutside =
massOrMoleFraction(outsideVolVars, referenceSystem, phaseIdx, compIdx);
188 const Scalar rhoInside =
massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx);
190 const Scalar rho = 0.5*(rhoInside + rhoFacet);
192 componentFlux[compIdx] = fluxVarsCache.diffusionTijInside(phaseIdx, compIdx)*xInside
193 + fluxVarsCache.diffusionTijFacet(phaseIdx, compIdx)*xFacet;
195 if (!scvf.boundary())
196 componentFlux[compIdx] += fluxVarsCache.diffusionTijOutside(phaseIdx, compIdx)*xOutside;
197 componentFlux[compIdx] *= rho;
199 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
200 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
203 return componentFlux;
208 template<
class Problem,
class ElementVolumeVariables >
209 static typename Cache::DiffusionTransmissibilityContainer
211 const Element& element,
212 const FVElementGeometry& fvGeometry,
213 const ElementVolumeVariables& elemVolVars,
214 const SubControlVolumeFace& scvf,
215 unsigned int phaseIdx,
unsigned int compIdx)
217 typename Cache::DiffusionTransmissibilityContainer tij;
218 if (!problem.couplingManager().isCoupled(element, scvf))
221 tij[Cache::insideTijIdx] = ParentType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
226 static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(),
"FacetCoupling.Xi", 1.0);
228 const auto insideScvIdx = scvf.insideScvIdx();
229 const auto& insideScv = fvGeometry.scv(insideScvIdx);
230 const auto& insideVolVars = elemVolVars[insideScvIdx];
233 insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
234 insideVolVars.extrusionFactor());
237 const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf);
240 if (iBcTypes.hasOnlyNeumann())
242 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
243 const auto wFacet = 2.0*scvf.area()*insideVolVars.extrusionFactor()
244 /facetVolVars.extrusionFactor()
245 *
vtmv(scvf.unitOuterNormal(),
246 facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
247 scvf.unitOuterNormal());
256 if (!scvf.boundary())
258 const auto outsideScvIdx = scvf.outsideScvIdx();
259 const auto& outsideVolVars = elemVolVars[outsideScvIdx];
261 fvGeometry.scv(outsideScvIdx),
262 outsideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
263 outsideVolVars.extrusionFactor());
265 if ( !Dune::FloatCmp::eq(xi, 1.0, 1e-6) )
269 const Scalar factor = wIn * wFacet / ( wIn * wOut * ( 2.0 * xi - 1.0 ) + wFacet * ( xi * ( wIn + wOut ) + wFacet ) );
270 tij[Cache::insideTijIdx] = factor * ( wOut * xi + wFacet );
271 tij[Cache::outsideTijIdx] = factor * ( wOut * ( 1.0 - xi ) );
272 tij[Cache::facetTijIdx] = factor * ( - wOut - wFacet );
276 tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
277 tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
278 tij[Cache::outsideTijIdx] = 0.0;
283 tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
284 tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
285 tij[Cache::outsideTijIdx] = 0.0;
288 else if (iBcTypes.hasOnlyDirichlet())
290 tij[Cache::insideTijIdx] = wIn;
291 tij[Cache::outsideTijIdx] = 0.0;
292 tij[Cache::facetTijIdx] = -wIn;
295 DUNE_THROW(Dune::NotImplemented,
"Interior boundary types other than pure Dirichlet or Neumann");
305template<
class TypeTag, ReferenceSystemFormulation referenceSystem>
313 using FVElementGeometry =
typename GridGeometry::LocalView;
314 using SubControlVolume =
typename GridGeometry::SubControlVolume;
315 using SubControlVolumeFace =
typename GridGeometry::SubControlVolumeFace;
317 using GridView =
typename GridGeometry::GridView;
318 using Element =
typename GridView::template Codim<0>::Entity;
324 static const int numPhases = ModelTraits::numFluidPhases();
325 static const int numComponents = ModelTraits::numFluidComponents();
328 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
333 class FacetCouplingFicksLawCache
337 using Filler =
typename ParentType::Cache::Filler;
342 static constexpr int insideTijIdx = 0;
343 static constexpr int facetTijIdx = 1;
346 using DiffusionTransmissibilityContainer = std::array<Scalar, 2>;
349 template<
class Problem,
class ElementVolumeVariables >
350 void updateDiffusion(
const Problem& problem,
351 const Element& element,
352 const FVElementGeometry& fvGeometry,
353 const ElementVolumeVariables& elemVolVars,
354 const SubControlVolumeFace &scvf,
355 unsigned int phaseIdx,
356 unsigned int compIdx)
358 tij_[phaseIdx][compIdx] = Implementation::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
365 Scalar diffusionTij(
unsigned int phaseIdx,
unsigned int compIdx)
const
366 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
369 Scalar diffusionTijInside(
unsigned int phaseIdx,
unsigned int compIdx)
const
370 {
return tij_[phaseIdx][compIdx][insideTijIdx]; }
373 Scalar diffusionTijFacet(
unsigned int phaseIdx,
unsigned int compIdx)
const
374 {
return tij_[phaseIdx][compIdx][facetTijIdx]; }
377 std::array< std::array<DiffusionTransmissibilityContainer, numComponents>, numPhases> tij_;
382 using Cache = FacetCouplingFicksLawCache;
389 {
return referenceSystem; }
392 template<
class Problem,
class ElementVolumeVariables,
class ElementFluxVarsCache >
393 static ComponentFluxVector
flux(
const Problem& problem,
394 const Element& element,
395 const FVElementGeometry& fvGeometry,
396 const ElementVolumeVariables& elemVolVars,
397 const SubControlVolumeFace& scvf,
399 const ElementFluxVarsCache& elemFluxVarsCache)
401 if (!problem.couplingManager().isOnInteriorBoundary(element, scvf))
402 return ParentType::flux(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, elemFluxVarsCache);
404 ComponentFluxVector componentFlux(0.0);
405 for (
int compIdx = 0; compIdx < numComponents; compIdx++)
407 if(compIdx == FluidSystem::getMainComponent(phaseIdx))
411 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
412 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
413 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
416 const Scalar xInside =
massOrMoleFraction(insideVolVars, referenceSystem, phaseIdx, compIdx);
417 const Scalar xFacet =
massOrMoleFraction(facetVolVars, referenceSystem, phaseIdx, compIdx);
419 const Scalar rhoInside =
massOrMolarDensity(insideVolVars, referenceSystem, phaseIdx);
421 const Scalar rho = 0.5*(rhoInside + rhoFacet);
423 componentFlux[compIdx] = rho*(fluxVarsCache.diffusionTijInside(phaseIdx, compIdx)*xInside
424 + fluxVarsCache.diffusionTijFacet(phaseIdx, compIdx)*xFacet);
426 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx) && !FluidSystem::isTracerFluidSystem())
427 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
430 return componentFlux;
435 template<
class Problem,
class ElementVolumeVariables >
436 static typename Cache::DiffusionTransmissibilityContainer
438 const Element& element,
439 const FVElementGeometry& fvGeometry,
440 const ElementVolumeVariables& elemVolVars,
441 const SubControlVolumeFace& scvf,
442 unsigned int phaseIdx,
unsigned int compIdx)
444 typename Cache::DiffusionTransmissibilityContainer tij;
445 if (!problem.couplingManager().isCoupled(element, scvf))
448 tij[Cache::insideTijIdx] = ParentType::calculateTransmissibility(problem, element, fvGeometry, elemVolVars, scvf, phaseIdx, compIdx);
453 static const Scalar xi = getParamFromGroup<Scalar>(problem.paramGroup(),
"FacetCoupling.Xi", 1.0);
458 if (Dune::FloatCmp::ne(xi, 1.0, 1e-6))
459 DUNE_THROW(Dune::InvalidStateException,
"Xi != 1.0 cannot be used on surface grids");
461 const auto insideScvIdx = scvf.insideScvIdx();
462 const auto& insideScv = fvGeometry.scv(insideScvIdx);
463 const auto& insideVolVars = elemVolVars[insideScvIdx];
466 insideVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
467 insideVolVars.extrusionFactor());
470 const auto iBcTypes = problem.interiorBoundaryTypes(element, scvf);
473 if (iBcTypes.hasOnlyNeumann())
478 const auto& facetVolVars = problem.couplingManager().getLowDimVolVars(element, scvf);
479 const auto wFacet = 2.0*scvf.area()*insideVolVars.extrusionFactor()
480 /sqrt(facetVolVars.extrusionFactor())
481 *
vtmv(scvf.unitOuterNormal(),
482 facetVolVars.effectiveDiffusionCoefficient(phaseIdx, phaseIdx, compIdx),
483 scvf.unitOuterNormal());
485 tij[Cache::insideTijIdx] = wFacet*wIn/(wIn+wFacet);
486 tij[Cache::facetTijIdx] = -tij[Cache::insideTijIdx];
488 else if (iBcTypes.hasOnlyDirichlet())
490 tij[Cache::insideTijIdx] = wIn;
491 tij[Cache::facetTijIdx] = -wIn;
494 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.
The available discretization methods in Dumux.
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
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:840
typename Properties::Detail::GetPropImpl< TypeTag, Property >::type::type GetPropType
get the type alias defined in the property (equivalent to old macro GET_PROP_TYPE(....
Definition: propertysystem.hh:149
forward declaration of the method-specific implemetation
Definition: flux/box/fickslaw.hh:43
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:125
Forward declaration of the implementation.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:49
Specialization of CCTpfaFacetCouplingFicksLawImpl for dim=dimWorld.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:70
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:160
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:210
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
Return the reference system.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:155
Specialization of CCTpfaFacetCouplingFicksLawImpl for dim<dimWorld.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:308
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
Return the reference system.
Definition: multidomain/facet/cellcentered/tpfa/fickslaw.hh:388
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:437
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:393
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.