28 #ifndef EWOMS_DIFFUSION_MODULE_HH 29 #define EWOMS_DIFFUSION_MODULE_HH 34 #include <opm/common/Valgrind.hpp> 35 #include <opm/common/Unused.hpp> 36 #include <opm/common/ErrorMacros.hpp> 37 #include <opm/common/Exceptions.hpp> 39 #include <dune/common/fvector.hh> 42 namespace Properties {
52 template <
class TypeTag,
bool enableDiffusion>
58 template <
class TypeTag>
62 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
63 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
76 template <
class Context>
78 const Context& context OPM_UNUSED,
79 unsigned spaceIdx OPM_UNUSED,
80 unsigned timeIdx OPM_UNUSED)
87 template <
class TypeTag>
91 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
92 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
93 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
96 enum { numPhases = FluidSystem::numPhases };
97 enum { numComponents = FluidSystem::numComponents };
98 enum { conti0EqIdx = Indices::conti0EqIdx };
100 typedef Opm::MathToolbox<Evaluation> Toolbox;
113 template <
class Context>
115 unsigned spaceIdx,
unsigned timeIdx)
117 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
119 const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
120 const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
122 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
124 Evaluation rhoMolar = fluidStateI.molarDensity(phaseIdx);
125 rhoMolar += Toolbox::value(fluidStateJ.molarDensity(phaseIdx));
128 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx)
130 flux[conti0EqIdx + compIdx] +=
132 * extQuants.moleFractionGradientNormal(phaseIdx, compIdx)
133 * extQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx);
145 template <
class TypeTag,
bool enableDiffusion>
151 template <
class TypeTag>
155 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
156 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
165 OPM_THROW(std::logic_error,
"Method tortuosity() does not make sense " 166 "if diffusion is disabled");
175 OPM_THROW(std::logic_error,
"Method diffusionCoefficient() does not " 176 "make sense if diffusion is disabled");
185 OPM_THROW(std::logic_error,
"Method effectiveDiffusionCoefficient() " 186 "does not make sense if diffusion is " 195 template <
class Flu
idState>
197 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache OPM_UNUSED,
198 const ElementContext& elemCtx OPM_UNUSED,
199 unsigned dofIdx OPM_UNUSED,
200 unsigned timeIdx OPM_UNUSED)
207 template <
class TypeTag>
211 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
212 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
213 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
215 enum { numPhases = FluidSystem::numPhases };
216 enum { numComponents = FluidSystem::numComponents };
224 {
return diffusionCoefficient_[phaseIdx][compIdx]; }
231 {
return tortuosity_[phaseIdx]; }
238 {
return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; }
245 template <
class Flu
idState>
247 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
248 const ElementContext& elemCtx,
252 typedef Opm::MathToolbox<Evaluation> Toolbox;
254 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
255 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
256 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
262 const Evaluation& base =
265 * intQuants.fluidState().saturation(phaseIdx));
266 tortuosity_[phaseIdx] =
267 1.0 / (intQuants.porosity() * intQuants.porosity())
268 * Toolbox::pow(base, 7.0/3.0);
270 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
271 diffusionCoefficient_[phaseIdx][compIdx] =
272 FluidSystem::diffusionCoefficient(fluidState,
281 Evaluation tortuosity_[numPhases];
282 Evaluation diffusionCoefficient_[numPhases][numComponents];
291 template <
class TypeTag,
bool enableDiffusion>
297 template <
class TypeTag>
301 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
302 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
309 void update_(
const ElementContext& elemCtx OPM_UNUSED,
310 unsigned faceIdx OPM_UNUSED,
311 unsigned timeIdx OPM_UNUSED)
314 template <
class Context,
class Flu
idState>
315 void updateBoundary_(
const Context& context OPM_UNUSED,
316 unsigned bfIdx OPM_UNUSED,
317 unsigned timeIdx OPM_UNUSED,
318 const FluidState& fluidState OPM_UNUSED)
329 unsigned compIdx OPM_UNUSED)
const 331 OPM_THROW(std::logic_error,
332 "The method moleFractionGradient() does not " 333 "make sense if diffusion is disabled.");
344 unsigned compIdx OPM_UNUSED)
const 346 OPM_THROW(std::logic_error,
347 "The method effectiveDiffusionCoefficient() " 348 "does not make sense if diffusion is disabled.");
355 template <
class TypeTag>
359 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
360 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
363 enum { dimWorld = GridView::dimensionworld };
367 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
368 typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
375 void update_(
const ElementContext& elemCtx,
unsigned faceIdx,
unsigned timeIdx)
377 const auto& gradCalc = elemCtx.gradientCalculator();
380 DimEvalVector temperatureGrad;
382 const auto& face = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
383 const auto& normal = face.normal();
384 const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
386 const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
387 const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
389 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
390 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
394 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
397 DimEvalVector moleFractionGradient(0.0);
398 gradCalc.calculateGradient(moleFractionGradient,
401 moleFractionCallback);
403 moleFractionGradientNormal_[phaseIdx][compIdx] = 0.0;
404 for (
unsigned i = 0; i < normal.size(); ++i)
405 moleFractionGradientNormal_[phaseIdx][compIdx] +=
406 normal[i]*moleFractionGradient[i];
407 Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
411 effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
412 (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
414 intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
417 Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
422 template <
class Context,
class Flu
idState>
423 void updateBoundary_(
const Context& context,
426 const FluidState& fluidState)
428 const auto& stencil = context.stencil(timeIdx);
429 const auto& face = stencil.boundaryFace(bfIdx);
431 const auto& elemCtx = context.elementContext();
432 unsigned insideScvIdx = face.interiorIndex();
433 const auto& insideScv = stencil.subControlVolume(insideScvIdx);
435 const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
436 const auto& fluidStateInside = intQuantsInside.fluidState();
439 DimVector distVec = face.integrationPos();
440 distVec -= context.element().geometry().global(insideScv.localGeometry().center());
442 Scalar dist = distVec * face.normal();
448 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
449 if (!elemCtx.model().phaseIsConsidered(phaseIdx))
452 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
455 moleFractionGradientNormal_[phaseIdx][compIdx] =
456 (fluidState.moleFraction(phaseIdx, compIdx)
458 fluidStateInside.moleFraction(phaseIdx, compIdx))
460 Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
464 effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
465 intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx);
467 Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
480 {
return moleFractionGradientNormal_[phaseIdx][compIdx]; }
490 {
return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
493 Evaluation moleFractionGradientNormal_[numPhases][numComponents];
494 Evaluation effectiveDiffusionCoefficient_[numPhases][numComponents];
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes...
Definition: diffusionmodule.hh:146
Definition: baseauxiliarymodule.hh:37
Scalar effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: diffusionmodule.hh:183
static void addDiffusiveFlux(RateVector &flux, const Context &context, unsigned spaceIdx, unsigned timeIdx)
Adds the mass flux due to molecular diffusion to the flux vector over the flux integration point...
Definition: diffusionmodule.hh:114
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:230
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:479
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
const Evaluation & effectiveDiffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
The effective diffusion coeffcient of a component in a fluid phase at the face's integration point...
Definition: diffusionmodule.hh:343
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:375
const Evaluation & effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
The effective diffusion coeffcient of a component in a fluid phase at the face's integration point...
Definition: diffusionmodule.hh:489
Declare the properties used by the infrastructure code of the finite volume discretizations.
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:53
Scalar diffusionCoefficient(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:173
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:106
This method contains all callback classes for quantities that are required by some extensive quantiti...
void setComponentIndex(unsigned compIdx)
Set the index of the component for which the mole fraction should be returned.
Definition: quantitycallbacks.hh:460
void update_(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:246
static void addDiffusiveFlux(RateVector &flux OPM_UNUSED, const Context &context OPM_UNUSED, unsigned spaceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Adds the diffusive mass flux flux to the flux vector over a flux integration point.
Definition: diffusionmodule.hh:77
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which a mole fraction should be returned.
Definition: quantitycallbacks.hh:453
Callback class for a mole fraction of a component in a phase.
Definition: quantitycallbacks.hh:424
Scalar tortuosity(unsigned phaseIdx OPM_UNUSED) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:163
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx OPM_UNUSED, unsigned compIdx OPM_UNUSED) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:328
void update_(const ElementContext &elemCtx OPM_UNUSED, unsigned faceIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:309
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:69
Provides the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:292
#define NEW_PROP_TAG(PTagName)
Define a property tag.
Definition: propertysystem.hh:247
Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: diffusionmodule.hh:237
void update_(FluidState &fs OPM_UNUSED, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > ¶mCache OPM_UNUSED, const ElementContext &elemCtx OPM_UNUSED, unsigned dofIdx OPM_UNUSED, unsigned timeIdx OPM_UNUSED)
Update the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:196
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:223