25#ifndef DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
26#define DUMUX_DISCRETIZATION_BOX_FICKS_LAW_HH
28#include <dune/common/fvector.hh>
29#include <dune/common/fmatrix.hh>
42template<
class TypeTag, DiscretizationMethod discMethod, ReferenceSystemFormulation referenceSystem>
49template <
class TypeTag, ReferenceSystemFormulation referenceSystem>
57 using SubControlVolume =
typename FVElementGeometry::SubControlVolume;
58 using SubControlVolumeFace =
typename FVElementGeometry::SubControlVolumeFace;
64 using Element =
typename GridView::template Codim<0>::Entity;
67 enum { dim = GridView::dimension} ;
68 enum { dimWorld = GridView::dimensionworld} ;
71 numPhases = ModelTraits::numFluidPhases(),
72 numComponents = ModelTraits::numFluidComponents()
74 using DimWorldMatrix = Dune::FieldMatrix<Scalar, dimWorld, dimWorld>;
75 using ComponentFluxVector = Dune::FieldVector<Scalar, numComponents>;
82 {
return referenceSystem; }
84 static ComponentFluxVector
flux(
const Problem& problem,
85 const Element& element,
86 const FVElementGeometry& fvGeometry,
87 const ElementVolumeVariables& elemVolVars,
88 const SubControlVolumeFace& scvf,
90 const ElementFluxVariablesCache& elemFluxVarsCache)
93 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
94 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
97 const auto& fluxVarsCache = elemFluxVarsCache[scvf];
98 const auto& shapeValues = fluxVarsCache.shapeValues();
102 for (
auto&& scv : scvs(fvGeometry))
103 rho +=
massOrMolarDensity(elemVolVars[scv], referenceSystem, phaseIdx)*shapeValues[scv.indexInElement()][0];
105 ComponentFluxVector componentFlux(0.0);
106 for (
int compIdx = 0; compIdx < numComponents; compIdx++)
108 if constexpr (!FluidSystem::isTracerFluidSystem())
109 if (compIdx == FluidSystem::getMainComponent(phaseIdx))
112 const auto D = averageDiffusionCoefficient_(phaseIdx, compIdx, insideVolVars, outsideVolVars, problem, scvf);
115 const auto massOrMoleFrac = [&](
const SubControlVolume& scv){
return massOrMoleFraction(elemVolVars[scv], referenceSystem, phaseIdx, compIdx); };
116 componentFlux[compIdx] = discreteFlux_(fvGeometry, scvf, fluxVarsCache, massOrMoleFrac, D, rho);
119 if constexpr (!FluidSystem::isTracerFluidSystem())
120 if (BalanceEqOpts::mainComponentIsBalanced(phaseIdx))
121 componentFlux[FluidSystem::getMainComponent(phaseIdx)] -= componentFlux[compIdx];
123 return componentFlux;
127 static std::array<std::vector<Scalar>, numComponents>
129 const Element& element,
130 const FVElementGeometry& fvGeometry,
131 const ElementVolumeVariables& elemVolVars,
132 const SubControlVolumeFace& scvf,
133 const FluxVarCache& fluxVarCache,
137 const auto& shapeValues = fluxVarCache.shapeValues();
138 for (
auto&& scv : scvs(fvGeometry))
139 rho +=
massOrMolarDensity(elemVolVars[scv], referenceSystem, phaseIdx)*shapeValues[scv.indexInElement()][0];
141 const auto& insideVolVars = elemVolVars[scvf.insideScvIdx()];
142 const auto& outsideVolVars = elemVolVars[scvf.outsideScvIdx()];
144 std::array<std::vector<Scalar>, numComponents> ti;
145 for (
int compIdx = 0; compIdx < numComponents; compIdx++)
147 if constexpr (!FluidSystem::isTracerFluidSystem())
148 if (compIdx == FluidSystem::getMainComponent(phaseIdx))
151 const auto D = averageDiffusionCoefficient_(phaseIdx, compIdx, insideVolVars, outsideVolVars, problem, scvf);
153 ti[compIdx].resize(fvGeometry.numScv());
154 for (
auto&& scv : scvs(fvGeometry))
155 ti[compIdx][scv.indexInElement()] = -rho*
vtmv(scvf.unitOuterNormal(), D, fluxVarCache.gradN(scv.indexInElement()))*scvf.area();
162 static Scalar averageDiffusionCoefficient_(
const int phaseIdx,
const int compIdx,
163 const VolumeVariables& insideVV,
const VolumeVariables& outsideVV,
164 const Problem& problem,
165 const SubControlVolumeFace& scvf)
168 auto [insideD, outsideD] = diffusionCoefficientsAtInterface_(phaseIdx, compIdx, insideVV, outsideVV);
171 insideD *= insideVV.extrusionFactor();
172 outsideD *= outsideVV.extrusionFactor();
175 return problem.spatialParams().harmonicMean(insideD, outsideD, scvf.unitOuterNormal());
178 static std::pair<Scalar, Scalar>
179 diffusionCoefficientsAtInterface_([[maybe_unused]]
const int phaseIdx,
const int compIdx,
180 const VolumeVariables& insideVV,
const VolumeVariables& outsideVV)
182 if constexpr (!FluidSystem::isTracerFluidSystem())
184 using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
185 const auto mainCompIdx = FluidSystem::getMainComponent(phaseIdx);
186 const auto insideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(insideVV, phaseIdx, mainCompIdx, compIdx);
187 const auto outsideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(outsideVV, phaseIdx, mainCompIdx, compIdx);
188 return { std::move(insideD), std::move(outsideD) };
192 using EffDiffModel = GetPropType<TypeTag, Properties::EffectiveDiffusivityModel>;
193 const auto insideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(insideVV, 0, 0, compIdx);
194 const auto outsideD = Deprecated::template effectiveDiffusionCoefficient<EffDiffModel>(outsideVV, 0, 0, compIdx);
195 return { std::move(insideD), std::move(outsideD) };
199 template<
class EvaluateVariable,
class Tensor>
200 static Scalar discreteFlux_(
const FVElementGeometry& fvGeometry,
201 const SubControlVolumeFace& scvf,
202 const FluxVarCache& fluxVarsCache,
204 const Tensor& D,
const Scalar preFactor)
206 Dune::FieldVector<Scalar, dimWorld> gradX(0.0);
207 for (
auto&& scv : scvs(fvGeometry))
209 return -1.0*preFactor*
vtmv(scvf.unitOuterNormal(), D, gradX)*scvf.area();
Define some often used mathematical functions.
The available discretization methods in Dumux.
The reference frameworks and formulations available for splitting total fluxes into a advective and d...
Container storing the diffusion coefficients required by Fick's law. Uses the minimal possible contai...
DiscretizationMethod
The available discretization methods in Dumux.
Definition: method.hh:37
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
static constexpr ReferenceSystemFormulation referenceSystemFormulation()
Definition: flux/box/fickslaw.hh:81
static ComponentFluxVector flux(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const int phaseIdx, const ElementFluxVariablesCache &elemFluxVarsCache)
Definition: flux/box/fickslaw.hh:84
static std::array< std::vector< Scalar >, numComponents > calculateTransmissibilities(const Problem &problem, const Element &element, const FVElementGeometry &fvGeometry, const ElementVolumeVariables &elemVolVars, const SubControlVolumeFace &scvf, const FluxVarCache &fluxVarCache, const int phaseIdx)
Definition: flux/box/fickslaw.hh:128
Container storing the diffusion coefficients required by Fick's law. Uses the minimal possible contai...
Definition: fickiandiffusioncoefficients.hh:44
Declares all properties used in Dumux.