28#ifndef OPM_DIFFUSION_MODULE_HH
29#define OPM_DIFFUSION_MODULE_HH
31#include <dune/common/fvector.hh>
33#include <opm/material/common/Valgrind.hpp>
50template <
class TypeTag,
bool enableDiffusion>
56template <
class TypeTag>
74 template <
class Context>
85template <
class TypeTag>
94 enum { numPhases = FluidSystem::numPhases };
95 enum { numComponents = FluidSystem::numComponents };
96 enum { conti0EqIdx = Indices::conti0EqIdx };
98 using Toolbox = Opm::MathToolbox<Evaluation>;
111 template <
class Context>
113 unsigned spaceIdx,
unsigned timeIdx)
115 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
117 const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
118 const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
120 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
122 Evaluation rhoMolar = fluidStateI.molarDensity(phaseIdx);
123 rhoMolar += Toolbox::value(fluidStateJ.molarDensity(phaseIdx));
126 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
128 flux[conti0EqIdx + compIdx] +=
130 * extQuants.moleFractionGradientNormal(phaseIdx, compIdx)
131 * extQuants.effectiveDiffusionCoefficient(phaseIdx, compIdx);
144template <
class TypeTag,
bool enableDiffusion>
150template <
class TypeTag>
164 throw std::logic_error(
"Method tortuosity() does not make sense "
165 "if diffusion is disabled");
174 throw std::logic_error(
"Method diffusionCoefficient() does not "
175 "make sense if diffusion is disabled");
184 throw std::logic_error(
"Method effectiveDiffusionCoefficient() "
185 "does not make sense if diffusion is disabled");
193 template <
class Flu
idState>
195 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>&,
196 const ElementContext&,
205template <
class TypeTag>
213 enum { numPhases = FluidSystem::numPhases };
214 enum { numComponents = FluidSystem::numComponents };
222 {
return diffusionCoefficient_[phaseIdx][compIdx]; }
229 {
return tortuosity_[phaseIdx]; }
236 {
return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx]; }
243 template <
class Flu
idState>
245 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
246 const ElementContext& elemCtx,
250 using Toolbox = Opm::MathToolbox<Evaluation>;
252 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
253 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
254 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, 10.0/3.0);
270 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
271 diffusionCoefficient_[phaseIdx][compIdx] =
272 FluidSystem::diffusionCoefficient(fluidState,
281 std::array<Evaluation, numPhases> tortuosity_;
282 std::array<std::array<Evaluation, numComponents>, numPhases> diffusionCoefficient_;
291template <
class TypeTag,
bool enableDiffusion>
297template <
class TypeTag>
314 template <
class Context,
class Flu
idState>
331 throw std::logic_error(
"The method moleFractionGradient() does not "
332 "make sense if diffusion is disabled.");
345 throw std::logic_error(
"The method effectiveDiffusionCoefficient() "
346 "does not make sense if diffusion is disabled.");
353template <
class TypeTag>
361 enum { dimWorld = GridView::dimensionworld };
362 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
363 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
365 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
366 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
373 void update_(
const ElementContext& elemCtx,
unsigned faceIdx,
unsigned timeIdx)
375 const auto& gradCalc = elemCtx.gradientCalculator();
378 const auto& face = elemCtx.stencil(timeIdx).interiorFace(faceIdx);
379 const auto& normal = face.normal();
380 const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
382 const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
383 const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
385 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
386 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
391 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
394 DimEvalVector moleFractionGradient(0.0);
395 gradCalc.calculateGradient(moleFractionGradient,
398 moleFractionCallback);
400 moleFractionGradientNormal_[phaseIdx][compIdx] = 0.0;
401 for (
unsigned i = 0; i < normal.size(); ++i) {
402 moleFractionGradientNormal_[phaseIdx][compIdx] +=
403 normal[i]*moleFractionGradient[i];
405 Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
409 effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
410 (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx)
412 intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx))
415 Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
420 template <
class Context,
class Flu
idState>
424 const FluidState& fluidState)
426 const auto& stencil = context.stencil(timeIdx);
427 const auto& face = stencil.boundaryFace(bfIdx);
429 const auto& elemCtx = context.elementContext();
430 const unsigned insideScvIdx = face.interiorIndex();
431 const auto& insideScv = stencil.subControlVolume(insideScvIdx);
433 const auto& intQuantsInside = elemCtx.intensiveQuantities(insideScvIdx, timeIdx);
434 const auto& fluidStateInside = intQuantsInside.fluidState();
437 DimVector distVec = face.integrationPos();
438 distVec -= context.element().geometry().global(insideScv.localGeometry().center());
440 const Scalar dist = distVec * face.normal();
446 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
447 if (!elemCtx.model().phaseIsConsidered(phaseIdx)) {
451 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
454 moleFractionGradientNormal_[phaseIdx][compIdx] =
455 (fluidState.moleFraction(phaseIdx, compIdx)
457 fluidStateInside.moleFraction(phaseIdx, compIdx))
459 Opm::Valgrind::CheckDefined(moleFractionGradientNormal_[phaseIdx][compIdx]);
463 effectiveDiffusionCoefficient_[phaseIdx][compIdx] =
464 intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx);
466 Opm::Valgrind::CheckDefined(effectiveDiffusionCoefficient_[phaseIdx][compIdx]);
479 {
return moleFractionGradientNormal_[phaseIdx][compIdx]; }
489 {
return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
492 std::array<std::array<Evaluation, numComponents>, numPhases> moleFractionGradientNormal_;
493 std::array<std::array<Evaluation, numComponents>, numPhases> effectiveDiffusionCoefficient_;
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:309
const Evaluation & effectiveDiffusionCoefficient(unsigned, unsigned) const
The effective diffusion coeffcient of a component in a fluid phase at the face's integration point.
Definition: diffusionmodule.hh:342
const Evaluation & moleFractionGradientNormal(unsigned, unsigned) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:328
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: diffusionmodule.hh:315
void updateBoundary_(const Context &context, unsigned bfIdx, unsigned timeIdx, const FluidState &fluidState)
Definition: diffusionmodule.hh:421
const Evaluation & moleFractionGradientNormal(unsigned phaseIdx, unsigned compIdx) const
The the gradient of the mole fraction times the face normal.
Definition: diffusionmodule.hh:478
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:488
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: diffusionmodule.hh:373
Provides the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:292
Scalar diffusionCoefficient(unsigned, unsigned) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:172
Scalar effectiveDiffusionCoefficient(unsigned, unsigned) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: diffusionmodule.hh:182
void update_(FluidState &, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &, const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate diffusive mass fluxes.
Definition: diffusionmodule.hh:194
Scalar tortuosity(unsigned) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:162
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: diffusionmodule.hh:228
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:244
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:235
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: diffusionmodule.hh:221
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes.
Definition: diffusionmodule.hh:145
static void addDiffusiveFlux(RateVector &, const Context &, unsigned, unsigned)
Adds the diffusive mass flux flux to the flux vector over a flux integration point.
Definition: diffusionmodule.hh:75
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:67
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:112
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: diffusionmodule.hh:104
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: diffusionmodule.hh:51
Callback class for a mole fraction of a component in a phase.
Definition: quantitycallbacks.hh:426
void setComponentIndex(unsigned compIdx)
Set the index of the component for which the mole fraction should be returned.
Definition: quantitycallbacks.hh:461
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which a mole fraction should be returned.
Definition: quantitycallbacks.hh:454
Declare the properties used by the infrastructure code of the finite volume discretizations.
Defines the common properties required by the porous medium multi-phase models.
Definition: blackoilboundaryratevector.hh:39
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:233
This method contains all callback classes for quantities that are required by some extensive quantiti...