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];
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
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes...
Definition: diffusionmodule.hh:146
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 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
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
#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, unsigned compIdx) const
The effective diffusion coeffcient of a component in a fluid phase at the face's integration point...
Definition: diffusionmodule.hh:489
#define GET_PROP_TYPE(TypeTag, PropTagName)
Access the type attribute of a property for a type tag.
Definition: propertysystem.hh:486
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:479
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:375
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
Declare the properties used by the infrastructure code of the finite volume discretizations.
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:230
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:53
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
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...
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
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:223
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
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
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