28#ifndef OPM_BLACKOIL_DIFFUSION_MODULE_HH
29#define OPM_BLACKOIL_DIFFUSION_MODULE_HH
31#include <dune/common/fvector.hh>
33#include <opm/material/common/Valgrind.hpp>
49template <
class TypeTag,
bool enableDiffusion>
52template <
class TypeTag,
bool enableDiffusion>
58template <
class TypeTag>
68 static void initFromState(
const EclipseState&)
82 template <
class Context>
93template <
class TypeTag>
105 enum { numPhases = FluidSystem::numPhases };
106 enum { conti0EqIdx = Indices::conti0EqIdx };
108 enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
109 enum { enableMICP = Indices::enableMICP };
111 static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
112 static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
113 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
114 static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
116 using Toolbox = MathToolbox<Evaluation>;
125 static void initFromState(
const EclipseState& eclState)
127 use_mole_fraction_ = eclState.getTableManager().diffMoleFraction();
159 template <
class Context>
161 unsigned spaceIdx,
unsigned timeIdx)
164 if (!FluidSystem::enableDiffusion()) {
168 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
169 const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
170 const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
171 const auto& diffusivity = extQuants.diffusivity();
172 const auto& effectiveDiffusionCoefficient = extQuants.effectiveDiffusionCoefficient();
173 addDiffusiveFlux(flux, inIq, exIq, diffusivity, effectiveDiffusionCoefficient);
174 if constexpr(enableBioeffects) {
175 const auto& effectiveBioDiffCoefficient = extQuants.effectiveBioDiffCoefficient();
176 addBioDiffFlux(flux, inIq, exIq, diffusivity, effectiveBioDiffCoefficient);
180 template<
class IntensiveQuantities,
class EvaluationArray>
182 const IntensiveQuantities& inIq,
183 const IntensiveQuantities& exIq,
184 const Evaluation& diffusivity,
185 const EvaluationArray& effectiveDiffusionCoefficient)
187 const auto& inFs = inIq.fluidState();
188 const auto& exFs = exIq.fluidState();
189 unsigned pvtRegionIndex = inFs.pvtRegionIndex();
190 Evaluation diffR = 0.0;
192 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
193 if (!FluidSystem::phaseIsActive(phaseIdx)) {
198 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
203 if (FluidSystem::gasPhaseIdx == phaseIdx && !FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
208 Evaluation bSAvg = inFs.saturation(phaseIdx) * inFs.invB(phaseIdx);
209 bSAvg += Toolbox::value(exFs.saturation(phaseIdx)) * Toolbox::value(exFs.invB(phaseIdx));
213 if (bSAvg < 1.0e-6) {
216 Evaluation convFactor = 1.0;
217 if (FluidSystem::enableDissolvedGas() &&
218 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
219 phaseIdx == FluidSystem::oilPhaseIdx)
221 const Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2;
222 convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex) + rsAvg);
223 diffR = inFs.Rs() - Toolbox::value(exFs.Rs());
225 if (FluidSystem::enableVaporizedOil() &&
226 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
227 phaseIdx == FluidSystem::gasPhaseIdx)
229 const Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2;
230 convFactor = toFractionGasOil(pvtRegionIndex) /
231 (1.0 + rvAvg * toFractionGasOil(pvtRegionIndex));
232 diffR = inFs.Rv() - Toolbox::value(exFs.Rv());
234 if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
235 const Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2;
236 convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex) + rsAvg);
237 diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw());
239 if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
240 const Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2;
241 convFactor = toFractionGasWater(pvtRegionIndex) /
242 (1.0 + rvAvg * toFractionGasWater(pvtRegionIndex));
243 diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw());
247 const unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
248 const unsigned activeSolventCompIdx = FluidSystem::canonicalToActiveCompIdx(solventCompIdx);
249 flux[conti0EqIdx + activeSolventCompIdx] +=
254 effectiveDiffusionCoefficient[phaseIdx][solventCompIdx];
257 const unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
258 const unsigned activeSoluteCompIdx = FluidSystem::canonicalToActiveCompIdx(soluteCompIdx);
259 flux[conti0EqIdx + activeSoluteCompIdx] +=
264 effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx];
268 template<
class IntensiveQuantities,
class EvaluationArray>
270 const IntensiveQuantities& inIq,
271 const IntensiveQuantities& exIq,
272 const Evaluation& diffusivity,
273 const EvaluationArray& effectiveBioDiffCoefficient)
275 const auto& inFs = inIq.fluidState();
276 const auto& exFs = exIq.fluidState();
277 Evaluation diffR = 0.0;
280 Evaluation bAvg = (inFs.saturation(waterPhaseIdx) * inFs.invB(waterPhaseIdx) +
281 Toolbox::value(exFs.saturation(waterPhaseIdx)) * Toolbox::value(exFs.invB(waterPhaseIdx))) / 2;
282 diffR = inIq.microbialConcentration() - Toolbox::value(exIq.microbialConcentration());
283 flux[contiMicrobialEqIdx] +=
287 effectiveBioDiffCoefficient[BioeffectsParams::micrDiffIdx];
288 if constexpr(enableMICP) {
289 diffR = inIq.oxygenConcentration() - Toolbox::value(exIq.oxygenConcentration());
290 flux[contiOxygenEqIdx] +=
294 effectiveBioDiffCoefficient[BioeffectsParams::oxygDiffIdx];
295 diffR = inIq.ureaConcentration() - Toolbox::value(exIq.ureaConcentration());
296 flux[contiUreaEqIdx] +=
300 effectiveBioDiffCoefficient[BioeffectsParams::ureaDiffIdx];
305 static Scalar toFractionGasOil (
unsigned regionIdx)
307 const Scalar mMOil = use_mole_fraction_ ? FluidSystem::molarMass(FluidSystem::oilCompIdx, regionIdx) : 1;
308 const Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
309 const Scalar mMGas = use_mole_fraction_ ? FluidSystem::molarMass(FluidSystem::gasCompIdx, regionIdx) : 1;
310 const Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
311 return rhoO * mMGas / (rhoG * mMOil);
314 static Scalar toFractionGasWater (
unsigned regionIdx)
316 const Scalar mMWater = use_mole_fraction_ ? FluidSystem::molarMass(FluidSystem::waterCompIdx, regionIdx) : 1;
317 const Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
318 const Scalar mMGas = use_mole_fraction_ ? FluidSystem::molarMass(FluidSystem::gasCompIdx, regionIdx) : 1;
319 const Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
320 return rhoW * mMGas / (rhoG * mMWater);
323 static bool use_mole_fraction_;
326template <
typename TypeTag>
328BlackOilDiffusionModule<TypeTag, true>::use_mole_fraction_;
337template <
class TypeTag,
bool enableDiffusion>
343template <
class TypeTag>
356 throw std::logic_error(
"Method tortuosity() does not make sense "
357 "if diffusion is disabled");
366 throw std::logic_error(
"Method diffusionCoefficient() does not "
367 "make sense if diffusion is disabled");
376 throw std::logic_error(
"Method effectiveDiffusionCoefficient() "
377 "does not make sense if diffusion is disabled");
385 template <
class Flu
idState>
388 const ElementContext&,
397template <
class TypeTag>
408 enum { numPhases = FluidSystem::numPhases };
409 enum { numComponents = FluidSystem::numComponents };
410 enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
411 enum { enableMICP = Indices::enableMICP };
412 enum { numBioInWat = Indices::numBioInWat };
414 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
430 if (FluidSystem::enableDiffusion()) {
431 tortuosity_ = rhs.tortuosity_;
432 diffusionCoefficient_ = rhs.diffusionCoefficient_;
442 {
return diffusionCoefficient_[phaseIdx][compIdx]; }
449 {
return tortuosity_[phaseIdx]; }
459 static constexpr bool enableTortuosity =
false;
460 if constexpr (enableTortuosity) {
461 return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx];
463 return diffusionCoefficient_[phaseIdx][compIdx];
472 {
return bioDiffCoefficient_[compIdx]; }
481 static bool enableTortuosity =
false;
482 if (enableTortuosity)
483 return tortuosity_[waterPhaseIdx] * bioDiffCoefficient_[compIdx];
485 return bioDiffCoefficient_[compIdx];
493 template <
class Flu
idState>
495 const unsigned regionIdx,
496 const ElementContext& elemCtx,
501 if (!FluidSystem::enableDiffusion()) {
505 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
506 update_(fluidState, regionIdx, intQuants);
509 template<
class Flu
idState>
511 const unsigned regionIdx,
512 const IntensiveQuantities& intQuants)
514 using Toolbox = MathToolbox<Evaluation>;
516 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
517 if (!FluidSystem::phaseIsActive(phaseIdx)) {
522 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
529 constexpr double myeps = 0.0001;
530 const Evaluation& base =
532 intQuants.porosity() *
533 intQuants.fluidState().saturation(phaseIdx));
534 tortuosity_[phaseIdx] =
535 1.0 / (intQuants.porosity() * intQuants.porosity()) *
536 Toolbox::pow(base, 10.0 / 3.0);
538 using PCache =
typename FluidSystem::template ParameterCache<Scalar>;
539 PCache pcache(regionIdx);
540 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
541 diffusionCoefficient_[phaseIdx][compIdx] =
542 FluidSystem::diffusionCoefficient(fluidState,
549 if constexpr(enableBioeffects) {
550 unsigned pvtRegionIndex = intQuants.pvtRegionIndex();
551 for (
unsigned compIdx = 0; compIdx < numBioInWat; ++compIdx) {
552 bioDiffCoefficient_[compIdx] =
553 BioeffectsModule::bioDiffCoefficient(pvtRegionIndex, compIdx);
559 std::array<Evaluation, numPhases> tortuosity_{};
560 std::array<std::array<Evaluation, numComponents>, numPhases> diffusionCoefficient_{};
561 std::array<Evaluation, numBioInWat> bioDiffCoefficient_{};
570template <
class TypeTag,
bool enableDiffusion>
571class BlackOilDiffusionExtensiveQuantities;
576template <
class TypeTag>
593 template <
class Context,
class Flu
idState>
607 throw std::logic_error(
"The method diffusivity() does not "
608 "make sense if diffusion is disabled.");
621 throw std::logic_error(
"The method effectiveDiffusionCoefficient() "
622 "does not make sense if diffusion is disabled.");
629template <
class TypeTag>
639 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
640 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
641 enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
642 enum { enableMICP = Indices::enableMICP };
643 enum { numBioInWat = Indices::numBioInWat };
653 void update_(
const ElementContext& elemCtx,
unsigned faceIdx,
unsigned timeIdx)
656 if (!FluidSystem::enableDiffusion()) {
660 const auto& stencil = elemCtx.stencil(timeIdx);
661 const auto& face = stencil.interiorFace(faceIdx);
662 const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
663 const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
664 const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
666 const Scalar diff = elemCtx.problem().diffusivity(elemCtx, face.interiorIndex(), face.exteriorIndex());
667 const Scalar faceArea = face.area();
668 diffusivity_ = diff / faceArea;
669 update(effectiveDiffusionCoefficient_, intQuantsInside, intQuantsOutside);
670 Valgrind::CheckDefined(diffusivity_);
671 if constexpr(enableBioeffects) {
672 updateBio(effectiveBioDiffCoefficient_, intQuantsInside, intQuantsOutside);
678 const IntensiveQuantities& intQuantsInside,
679 const IntensiveQuantities& intQuantsOutside)
684 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
685 if (!FluidSystem::phaseIsActive(phaseIdx)) {
689 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
692 for (
unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
693 effectiveDiffusionCoefficient[phaseIdx][compIdx] =
694 0.5 * (intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) +
695 intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx));
696 Valgrind::CheckDefined(effectiveDiffusionCoefficient[phaseIdx][compIdx]);
701 static void updateBio(std::array<Evaluation, numBioInWat>& effectiveBioDiffCoefficient,
702 const IntensiveQuantities& intQuantsInside,
703 const IntensiveQuantities& intQuantsOutside) {
705 for (
unsigned compIdx = 0; compIdx < numBioInWat; ++compIdx) {
706 effectiveBioDiffCoefficient[compIdx] =
707 0.5 * (intQuantsInside.effectiveBioDiffCoefficient(compIdx) +
708 intQuantsOutside.effectiveBioDiffCoefficient(compIdx));
709 Valgrind::CheckDefined(effectiveBioDiffCoefficient[compIdx]);
714 template <
class Context,
class Flu
idState>
720 throw std::runtime_error(
"Not implemented: Diffusion across boundary not implemented for blackoil");
731 {
return diffusivity_; }
741 {
return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
744 {
return effectiveDiffusionCoefficient_; }
747 {
return effectiveBioDiffCoefficient_[compIdx]; }
750 return effectiveBioDiffCoefficient_;
755 EvaluationArray effectiveDiffusionCoefficient_;
756 std::array<Evaluation, numBioInWat> effectiveBioDiffCoefficient_;
Contains the classes required to extend the black-oil model by bioeffects.
Contains the high level supplements required to extend the black oil model by bioeffects.
Definition: blackoilbioeffectsmodules.hh:93
Scalar diffusivity() const
The diffusivity the face.
Definition: blackoildiffusionmodule.hh:605
const Evaluation & effectiveDiffusionCoefficient(unsigned, unsigned) const
The effective diffusion coeffcient of a component in a fluid phase at the face's integration point.
Definition: blackoildiffusionmodule.hh:618
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:588
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildiffusionmodule.hh:594
Provides the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:631
Scalar diffusivity() const
The diffusivity of the face.
Definition: blackoildiffusionmodule.hh:730
static void update(EvaluationArray &effectiveDiffusionCoefficient, const IntensiveQuantities &intQuantsInside, const IntensiveQuantities &intQuantsOutside)
Definition: blackoildiffusionmodule.hh:677
const Evaluation & effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
The effective diffusion coefficient of a component in a fluid phase at the face's integration point.
Definition: blackoildiffusionmodule.hh:740
const Evaluation & effectiveBioDiffCoefficient(unsigned compIdx) const
Definition: blackoildiffusionmodule.hh:746
static void updateBio(std::array< Evaluation, numBioInWat > &effectiveBioDiffCoefficient, const IntensiveQuantities &intQuantsInside, const IntensiveQuantities &intQuantsOutside)
Definition: blackoildiffusionmodule.hh:701
std::array< std::array< Evaluation, numComponents >, numPhases > EvaluationArray
Definition: blackoildiffusionmodule.hh:646
const auto & effectiveBioDiffCoefficient() const
Definition: blackoildiffusionmodule.hh:749
const auto & effectiveDiffusionCoefficient() const
Definition: blackoildiffusionmodule.hh:743
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildiffusionmodule.hh:715
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:653
Provides the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:53
void update_(FluidState &, const unsigned, const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:386
Scalar effectiveDiffusionCoefficient(unsigned, unsigned) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: blackoildiffusionmodule.hh:374
Scalar diffusionCoefficient(unsigned, unsigned) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: blackoildiffusionmodule.hh:364
Scalar tortuosity(unsigned) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: blackoildiffusionmodule.hh:354
Evaluation bioDiffCoefficient(unsigned compIdx) const
Returns the molecular diffusion coefficient for a biocomponent in the water phase.
Definition: blackoildiffusionmodule.hh:471
Evaluation effectiveBioDiffCoefficient(unsigned compIdx) const
Returns the effective molecular diffusion coefficient of the porous medium for a biocomponent in the ...
Definition: blackoildiffusionmodule.hh:478
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: blackoildiffusionmodule.hh:441
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: blackoildiffusionmodule.hh:448
BlackOilDiffusionIntensiveQuantities(BlackOilDiffusionIntensiveQuantities &&) noexcept=default
Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: blackoildiffusionmodule.hh:455
BlackOilDiffusionIntensiveQuantities()=default
void update_(FluidState &fluidState, const unsigned regionIdx, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:494
void update_(FluidState &fluidState, const unsigned regionIdx, const IntensiveQuantities &intQuants)
Definition: blackoildiffusionmodule.hh:510
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes.
Definition: blackoildiffusionmodule.hh:338
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: blackoildiffusionmodule.hh:75
static void addDiffusiveFlux(RateVector &, const Context &, unsigned, unsigned)
Adds the diffusive mass flux flux to the flux vector over a flux integration point.
Definition: blackoildiffusionmodule.hh:83
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 integration point....
Definition: blackoildiffusionmodule.hh:160
static void addDiffusiveFlux(RateVector &flux, const IntensiveQuantities &inIq, const IntensiveQuantities &exIq, const Evaluation &diffusivity, const EvaluationArray &effectiveDiffusionCoefficient)
Definition: blackoildiffusionmodule.hh:181
static void addBioDiffFlux(RateVector &flux, const IntensiveQuantities &inIq, const IntensiveQuantities &exIq, const Evaluation &diffusivity, const EvaluationArray &effectiveBioDiffCoefficient)
Definition: blackoildiffusionmodule.hh:269
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: blackoildiffusionmodule.hh:134
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: blackoildiffusionmodule.hh:50
Declare the properties used by the infrastructure code of the finite volume discretizations.
Definition: blackoilbioeffectsmodules.hh:43
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
Struct holding the parameters for the BlackOilBioeffectsModule class.
Definition: blackoilbioeffectsparams.hpp:44