28#ifndef EWOMS_DISPERSION_MODULE_HH
29#define EWOMS_DISPERSION_MODULE_HH
31#include <dune/common/fvector.hh>
33#include <opm/material/common/MathToolbox.hpp>
34#include <opm/material/common/Valgrind.hpp>
41#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
42#include <opm/common/OpmLog/OpmLog.hpp>
57template <
class TypeTag,
bool enableDispersion>
60template <
class TypeTag,
bool enableDispersion>
66template <
class TypeTag>
79 static void initFromState(
const EclipseState&)
87 template <
class Context>
94 template<
class IntensiveQuantities,
class Scalar>
96 const IntensiveQuantities&,
97 const IntensiveQuantities&,
106template <
class TypeTag>
120 enum { numPhases = FluidSystem::numPhases };
121 enum { numComponents = FluidSystem::numComponents };
122 enum { conti0EqIdx = Indices::conti0EqIdx };
123 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
124 enum { enableBioeffects = getPropValue<TypeTag, Properties::EnableBioeffects>() };
125 enum { enableMICP = Indices::enableMICP };
127 static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
128 static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
129 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
130 static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
132 using Toolbox = MathToolbox<Evaluation>;
138 static void initFromState(
const EclipseState& eclState)
140 if (!eclState.getSimulationConfig().rock_config().dispersion()) {
144 if (eclState.getSimulationConfig().hasVAPWAT() || eclState.getSimulationConfig().hasVAPOIL()) {
145 OpmLog::warning(
"Dispersion is activated in combination with VAPWAT/VAPOIL. \n"
146 "Water/oil is still allowed to vaporize, but dispersion in the "
147 "gas phase is ignored.");
156 template <
class Context>
158 unsigned spaceIdx,
unsigned timeIdx)
161 if (!context.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
164 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
165 const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
166 const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
167 const auto& dispersivity = extQuants.dispersivity();
168 const auto& normVelocityAvg = extQuants.normVelocityAvg();
169 addDispersiveFlux(flux, inIq, exIq, dispersivity, normVelocityAvg);
192 template<
class IntensiveQuantities,
class Scalar>
194 const IntensiveQuantities& inIq,
195 const IntensiveQuantities& exIq,
196 const Evaluation& dispersivity,
197 const Scalar& normVelocityAvg)
199 const auto& inFs = inIq.fluidState();
200 const auto& exFs = exIq.fluidState();
201 Evaluation diffR = 0.0;
202 if constexpr(enableBioeffects) {
204 const Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2;
205 diffR = inIq.microbialConcentration() - Toolbox::value(exIq.microbialConcentration());
206 flux[contiMicrobialEqIdx] +=
208 normVelocityAvg[waterPhaseIdx] *
211 if constexpr(enableMICP) {
212 diffR = inIq.oxygenConcentration() - Toolbox::value(exIq.oxygenConcentration());
213 flux[contiOxygenEqIdx] +=
215 normVelocityAvg[waterPhaseIdx] *
218 diffR = inIq.ureaConcentration() - Toolbox::value(exIq.ureaConcentration());
219 flux[contiUreaEqIdx] +=
221 normVelocityAvg[waterPhaseIdx] *
228 unsigned pvtRegionIndex = inFs.pvtRegionIndex();
229 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
230 if (!FluidSystem::phaseIsActive(phaseIdx)) {
235 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
242 if (FluidSystem::gasPhaseIdx == phaseIdx) {
247 Evaluation bAvg = inFs.invB(phaseIdx);
248 bAvg += Toolbox::value(exFs.invB(phaseIdx));
251 Evaluation convFactor = 1.0;
252 if (FluidSystem::enableDissolvedGas() &&
253 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
254 phaseIdx == FluidSystem::oilPhaseIdx)
256 const Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2;
257 convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg);
258 diffR = inFs.Rs() - Toolbox::value(exFs.Rs());
260 if (FluidSystem::enableVaporizedOil() &&
261 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
262 phaseIdx == FluidSystem::gasPhaseIdx)
264 const Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2;
265 convFactor = toMassFractionGasOil(pvtRegionIndex) /
266 (1.0 + rvAvg * toMassFractionGasOil(pvtRegionIndex));
267 diffR = inFs.Rv() - Toolbox::value(exFs.Rv());
269 if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
270 const Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2;
271 convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg);
272 diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw());
274 if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
275 const Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2;
276 convFactor = toMassFractionGasWater(pvtRegionIndex) /
277 (1.0 + rvAvg * toMassFractionGasWater(pvtRegionIndex));
278 diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw());
282 const unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
283 const unsigned activeSolventCompIdx = FluidSystem::canonicalToActiveCompIdx(solventCompIdx);
284 flux[conti0EqIdx + activeSolventCompIdx] +=
286 normVelocityAvg[phaseIdx] *
292 const unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
293 const unsigned activeSoluteCompIdx = FluidSystem::canonicalToActiveCompIdx(soluteCompIdx);
294 flux[conti0EqIdx + activeSoluteCompIdx] +=
296 normVelocityAvg[phaseIdx] *
304 static Scalar toMassFractionGasOil (
unsigned regionIdx)
306 const Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
307 const Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
311 static Scalar toMassFractionGasWater (
unsigned regionIdx)
313 const Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
314 const Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
326template <
class TypeTag,
bool enableDispersion>
332template <
class TypeTag>
344 throw std::logic_error(
"Method normVelocityCell() "
345 "does not make sense if dispersion is disabled");
353 template<
class ElementContext>
363template <
class TypeTag>
371 enum { numPhases = FluidSystem::numPhases };
372 enum { numComponents = FluidSystem::numComponents };
373 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
374 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
375 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
376 enum { gasCompIdx = FluidSystem::gasCompIdx };
377 enum { oilCompIdx = FluidSystem::oilCompIdx };
378 enum { waterCompIdx = FluidSystem::waterCompIdx };
379 enum { conti0EqIdx = Indices::conti0EqIdx };
380 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
387 {
return normVelocityCell_[phaseIdx]; }
399 template<
class ElementContext>
400 void update_(
const ElementContext& elemCtx,
unsigned dofIdx,
unsigned timeIdx)
403 if (!elemCtx.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
406 const auto& problem = elemCtx.simulator().problem();
407 if (problem.model().linearizer().getVelocityInfo().empty()) {
410 const std::array<int, 3> phaseIdxs = {gasPhaseIdx, oilPhaseIdx, waterPhaseIdx};
411 const std::array<int, 3> compIdxs = {gasCompIdx, oilCompIdx, waterCompIdx};
412 const auto& velocityInf = problem.model().linearizer().getVelocityInfo();
413 const unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
414 const auto& velocityInfos = velocityInf[globalDofIdx];
415 for (
unsigned i = 0; i < phaseIdxs.size(); ++i) {
416 normVelocityCell_[i] = 0;
418 for (
const auto& velocityInfo : velocityInfos) {
419 for (
unsigned i = 0; i < phaseIdxs.size(); ++i) {
420 if (FluidSystem::phaseIsActive(phaseIdxs[i])) {
421 normVelocityCell_[phaseIdxs[i]] = std::max(normVelocityCell_[phaseIdxs[i]],
422 std::abs(velocityInfo.velocity[conti0EqIdx +
423 FluidSystem::canonicalToActiveCompIdx(compIdxs[i])]));
430 std::array<Scalar, numPhases> normVelocityCell_{};
439template <
class TypeTag,
bool enableDispersion>
440class BlackOilDispersionExtensiveQuantities;
445template <
class TypeTag>
453 enum { numPhases = FluidSystem::numPhases };
465 template <
class Context,
class Flu
idState>
476 const IntensiveQuantities&,
477 const IntensiveQuantities&)
486 throw std::logic_error(
"The method dispersivity() does not "
487 "make sense if dispersion is disabled.");
499 throw std::logic_error(
"The method normVelocityAvg() "
500 "does not make sense if dispersion is disabled.");
507template <
class TypeTag>
514 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
519 const IntensiveQuantities& intQuantsInside,
520 const IntensiveQuantities& intQuantsOutside)
522 for (
unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
523 if (!FluidSystem::phaseIsActive(phaseIdx)) {
527 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
533 if (FluidSystem::gasPhaseIdx == phaseIdx) {
538 normVelocityAvg[phaseIdx] =
539 0.5 * (intQuantsInside.normVelocityCell(phaseIdx) +
540 intQuantsOutside.normVelocityCell(phaseIdx));
541 Valgrind::CheckDefined(normVelocityAvg[phaseIdx]);
546 template <
class Context,
class Flu
idState>
551 {
throw std::runtime_error(
"Not implemented: Dispersion across boundary not implemented for blackoil"); }
561 {
return dispersivity_; }
571 {
return normVelocityAvg_[phaseIdx]; }
574 {
return normVelocityAvg_; }
577 Scalar dispersivity_;
578 ScalarArray normVelocityAvg_;
Declares the properties required by the black oil model.
Definition: blackoildispersionmodule.hh:447
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildispersionmodule.hh:466
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the dispersive fluxes.
Definition: blackoildispersionmodule.hh:460
Scalar[numPhases] ScalarArray
Definition: blackoildispersionmodule.hh:473
static void update(ScalarArray &, const IntensiveQuantities &, const IntensiveQuantities &)
Definition: blackoildispersionmodule.hh:475
Scalar normVelocityAvg(unsigned) const
The effective filter velocity coefficient in a fluid phase at the face's integration point.
Definition: blackoildispersionmodule.hh:497
Scalar dispersivity() const
The dispersivity the face.
Definition: blackoildispersionmodule.hh:484
Provides the quantities required to calculate dispersive mass fluxes.
Definition: blackoildispersionmodule.hh:509
Scalar dispersivity() const
The dispersivity of the face.
Definition: blackoildispersionmodule.hh:560
static void update(ScalarArray &normVelocityAvg, const IntensiveQuantities &intQuantsInside, const IntensiveQuantities &intQuantsOutside)
Definition: blackoildispersionmodule.hh:518
std::array< Scalar, numPhases > ScalarArray
Definition: blackoildispersionmodule.hh:517
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildispersionmodule.hh:547
const auto & normVelocityAvg() const
Definition: blackoildispersionmodule.hh:573
Scalar normVelocityAvg(unsigned phaseIdx) const
The effective velocity coefficient in a fluid phase at the face's integration point.
Definition: blackoildispersionmodule.hh:570
Provides the quantities required to calculate dispersive mass fluxes.
Definition: blackoildispersionmodule.hh:61
Scalar normVelocityCell(unsigned, unsigned) const
Returns the max. norm of the filter velocity of the cell.
Definition: blackoildispersionmodule.hh:342
void update_(ElementContext &, unsigned, unsigned)
Update the quantities required to calculate dispersive fluxes.
Definition: blackoildispersionmodule.hh:354
void update_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the quantities required to calculate dispersive mass fluxes. This considers the linear disperi...
Definition: blackoildispersionmodule.hh:400
Scalar normVelocityCell(unsigned phaseIdx) const
Returns the max. norm of the filter velocity of the cell.
Definition: blackoildispersionmodule.hh:386
Provides the volumetric quantities required for the calculation of dispersive fluxes.
Definition: blackoildispersionmodule.hh:327
static void addDispersiveFlux(RateVector &, const Context &, unsigned, unsigned)
Adds the dispersive flux to the flux vector over a flux integration point.
Definition: blackoildispersionmodule.hh:88
static void addDispersiveFlux(RateVector &, const IntensiveQuantities &, const IntensiveQuantities &, const Evaluation &, const Scalar &)
Definition: blackoildispersionmodule.hh:95
static void addDispersiveFlux(RateVector &flux, const IntensiveQuantities &inIq, const IntensiveQuantities &exIq, const Evaluation &dispersivity, const Scalar &normVelocityAvg)
Adds the mass flux due to dispersion to the flux vector over the integration point....
Definition: blackoildispersionmodule.hh:193
static void addDispersiveFlux(RateVector &flux, const Context &context, unsigned spaceIdx, unsigned timeIdx)
Adds the mass flux due to dispersion to the flux vector over the flux integration point.
Definition: blackoildispersionmodule.hh:157
Provides the auxiliary methods required for consideration of the dispersion equation.
Definition: blackoildispersionmodule.hh:58
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: 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