blackoildiffusionmodule.hh
Go to the documentation of this file.
1// -*- mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*-
2// vi: set et ts=4 sw=4 sts=4:
3/*
4 This file is part of the Open Porous Media project (OPM).
5
6 OPM is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 2 of the License, or
9 (at your option) any later version.
10
11 OPM is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with OPM. If not, see <http://www.gnu.org/licenses/>.
18
19 Consult the COPYING file in the top-level source directory of this
20 module for the precise wording of the license and the list of
21 copyright holders.
22*/
28#ifndef OPM_BLACKOIL_DIFFUSION_MODULE_HH
29#define OPM_BLACKOIL_DIFFUSION_MODULE_HH
30
32
33#include <opm/material/common/Valgrind.hpp>
34
35#include <dune/common/fvector.hh>
36
37#include <stdexcept>
38
39namespace Opm {
40
47template <class TypeTag, bool enableDiffusion>
49
50template <class TypeTag, bool enableDiffusion>
52
56template <class TypeTag>
57class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/false>
58{
62
63public:
64
65 #if HAVE_ECL_INPUT
69 static void initFromState(const EclipseState&)
70 {
71 }
72 #endif
73
77 static void registerParameters()
78 {}
79
84 template <class Context>
85 static void addDiffusiveFlux(RateVector&,
86 const Context&,
87 unsigned,
88 unsigned)
89 {}
90};
91
95template <class TypeTag>
96class BlackOilDiffusionModule<TypeTag, /*enableDiffusion=*/true>
97{
103
104 enum { numPhases = FluidSystem::numPhases };
105 enum { numComponents = FluidSystem::numComponents };
106 enum { conti0EqIdx = Indices::conti0EqIdx };
107
108 using Toolbox = MathToolbox<Evaluation>;
109
110public:
112
113
114 #if HAVE_ECL_INPUT
118 static void initFromState(const EclipseState& eclState)
119 {
120 use_mole_fraction_ = eclState.getTableManager().diffMoleFraction();
121 }
122 #endif
123
127 static void registerParameters()
128 {}
129
152 template <class Context>
153 static void addDiffusiveFlux(RateVector& flux, const Context& context,
154 unsigned spaceIdx, unsigned timeIdx)
155 {
156 // Only work if diffusion is enabled run-time by DIFFUSE in the deck
157 if(!FluidSystem::enableDiffusion())
158 return;
159 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
160 const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
161 const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
162 const auto& diffusivity = extQuants.diffusivity();
163 const auto& effectiveDiffusionCoefficient = extQuants.effectiveDiffusionCoefficient();
164 addDiffusiveFlux(flux, fluidStateI, fluidStateJ, diffusivity, effectiveDiffusionCoefficient);
165 }
166
167 template<class FluidState,class EvaluationArray>
168 static void addDiffusiveFlux(RateVector& flux,
169 const FluidState& fluidStateI,
170 const FluidState& fluidStateJ,
171 const Evaluation& diffusivity,
172 const EvaluationArray& effectiveDiffusionCoefficient)
173 {
174 unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex();
175 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
176 if (!FluidSystem::phaseIsActive(phaseIdx)) {
177 continue;
178 }
179
180 // no diffusion in water for blackoil models
181 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
182 continue;
183 }
184
185 // no diffusion in gas phase in water + gas system.
186 if (FluidSystem::gasPhaseIdx == phaseIdx && !FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx)) {
187 continue;
188 }
189
190 // arithmetic mean of the phase's b factor weighed by saturation
191 Evaluation bSAvg = fluidStateI.saturation(phaseIdx) * fluidStateI.invB(phaseIdx);
192 bSAvg += Toolbox::value(fluidStateJ.saturation(phaseIdx)) * Toolbox::value(fluidStateJ.invB(phaseIdx));
193 bSAvg /= 2;
194
195 // phase not present, skip
196 if(bSAvg < 1.0e-6)
197 continue;
198 Evaluation convFactor = 1.0;
199 Evaluation diffR = 0.0;
200 if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) {
201 Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2;
202 convFactor = 1.0 / (toFractionGasOil(pvtRegionIndex) + rsAvg);
203 diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs());
204 }
205 if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) {
206 Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2;
207 convFactor = toFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toFractionGasOil(pvtRegionIndex));
208 diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv());
209 }
210 if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
211 Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2;
212 convFactor = 1.0 / (toFractionGasWater(pvtRegionIndex) + rsAvg);
213 diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw());
214 }
215 if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
216 Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2;
217 convFactor = toFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toFractionGasWater(pvtRegionIndex));
218 diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw());
219 }
220
221 // mass flux of solvent component (oil in oil or gas in gas)
222 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
223 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
224 flux[conti0EqIdx + activeSolventCompIdx] +=
225 - bSAvg
226 * convFactor
227 * diffR
228 * diffusivity
229 * effectiveDiffusionCoefficient[phaseIdx][solventCompIdx];
230
231 // mass flux of solute component (gas in oil or oil in gas)
232 unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
233 unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx);
234 flux[conti0EqIdx + activeSoluteCompIdx] +=
235 bSAvg
236 * diffR
237 * convFactor
238 * diffusivity
239 * effectiveDiffusionCoefficient[phaseIdx][soluteCompIdx];
240 }
241 }
242
243private:
244 static Scalar toFractionGasOil (unsigned regionIdx) {
245 Scalar mMOil = use_mole_fraction_? FluidSystem::molarMass(FluidSystem::oilCompIdx, regionIdx) : 1;
246 Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
247 Scalar mMGas = use_mole_fraction_? FluidSystem::molarMass(FluidSystem::gasCompIdx, regionIdx) : 1;
248 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
249 return rhoO * mMGas / (rhoG * mMOil);
250 }
251 static Scalar toFractionGasWater (unsigned regionIdx) {
252 Scalar mMWater = use_mole_fraction_? FluidSystem::molarMass(FluidSystem::waterCompIdx, regionIdx) : 1;
253 Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
254 Scalar mMGas = use_mole_fraction_? FluidSystem::molarMass(FluidSystem::gasCompIdx, regionIdx) : 1;
255 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
256 return rhoW * mMGas / (rhoG * mMWater);
257 }
258
259 static bool use_mole_fraction_;
260};
261
262template <typename TypeTag>
263bool
264BlackOilDiffusionModule<TypeTag, true>::use_mole_fraction_;
265
273template <class TypeTag, bool enableDiffusion>
275
279template <class TypeTag>
280class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/false>
281{
285
286public:
291 Scalar tortuosity(unsigned) const
292 {
293 throw std::logic_error("Method tortuosity() does not make sense "
294 "if diffusion is disabled");
295 }
296
301 Scalar diffusionCoefficient(unsigned, unsigned) const
302 {
303 throw std::logic_error("Method diffusionCoefficient() does not "
304 "make sense if diffusion is disabled");
305 }
306
311 Scalar effectiveDiffusionCoefficient(unsigned, unsigned) const
312 {
313 throw std::logic_error("Method effectiveDiffusionCoefficient() "
314 "does not make sense if diffusion is disabled");
315 }
316
317protected:
322 template <class FluidState>
323 void update_(FluidState&,
324 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>&,
325 const ElementContext&,
326 unsigned,
327 unsigned)
328 { }
329};
330
334template <class TypeTag>
335class BlackOilDiffusionIntensiveQuantities<TypeTag, /*enableDiffusion=*/true>
336{
342 enum { numPhases = FluidSystem::numPhases };
343 enum { numComponents = FluidSystem::numComponents };
344
345public:
349
351
354 {
355 if (this == &rhs) return *this;
356
357 if (FluidSystem::enableDiffusion()) {
358 std::copy(rhs.tortuosity_, rhs.tortuosity_ + numPhases, tortuosity_);
359 for (size_t i = 0; i < numPhases; ++i) {
360 std::copy(rhs.diffusionCoefficient_[i],
361 rhs.diffusionCoefficient_[i]+numComponents,
362 diffusionCoefficient_[i]);
363 }
364 }
365 return *this;
366 }
371 Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
372 { return diffusionCoefficient_[phaseIdx][compIdx]; }
373
378 Evaluation tortuosity(unsigned phaseIdx) const
379 { return tortuosity_[phaseIdx]; }
380
385 Evaluation effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
386 {
387 // For the blackoil model tortuosity is disabled.
388 // TODO add a run-time parameter to enable tortuosity
389 static bool enableTortuosity = false;
390 if (enableTortuosity)
391 return tortuosity_[phaseIdx] * diffusionCoefficient_[phaseIdx][compIdx];
392
393 return diffusionCoefficient_[phaseIdx][compIdx];
394 }
395
396protected:
401 template <class FluidState>
402 void update_(FluidState& fluidState,
403 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
404 const ElementContext& elemCtx,
405 unsigned dofIdx,
406 unsigned timeIdx)
407 {
408 // Only work if diffusion is enabled run-time by DIFFUSE in the deck
409 if(!FluidSystem::enableDiffusion())
410 return;
411 const auto& intQuants = elemCtx.intensiveQuantities(dofIdx, timeIdx);
412 update_(fluidState, paramCache, intQuants);
413 }
414
415 template<class FluidState>
416 void update_(FluidState& fluidState,
417 typename FluidSystem::template ParameterCache<typename FluidState::Scalar>& paramCache,
418 const IntensiveQuantities& intQuants) {
419 using Toolbox = MathToolbox<Evaluation>;
420
421 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
422 if (!FluidSystem::phaseIsActive(phaseIdx)) {
423 continue;
424 }
425
426 // no diffusion in water for blackoil models
427 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
428 continue;
429 }
430
431 // Based on Millington, R. J., & Quirk, J. P. (1961).
432 // \Note: it is possible to use NumericalConstants later
433 // constexpr auto& numconst = GetPropValue<TypeTag, Properties::NumericalConstants>;
434 constexpr double myeps = 0.0001; //numconst.blackoildiffusionmoduleeps;
435 const Evaluation& base =
436 Toolbox::max(myeps, //0.0001,
437 intQuants.porosity()
438 * intQuants.fluidState().saturation(phaseIdx));
439 tortuosity_[phaseIdx] =
440 1.0 / (intQuants.porosity() * intQuants.porosity())
441 * Toolbox::pow(base, 10.0/3.0);
442
443 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
444 diffusionCoefficient_[phaseIdx][compIdx] =
445 FluidSystem::diffusionCoefficient(fluidState,
446 paramCache,
447 phaseIdx,
448 compIdx);
449 }
450 }
451 }
452
453private:
454 Evaluation tortuosity_[numPhases];
455 Evaluation diffusionCoefficient_[numPhases][numComponents];
456};
457
464template <class TypeTag, bool enableDiffusion>
465class BlackOilDiffusionExtensiveQuantities;
466
470template <class TypeTag>
471class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/false>
472{
476
477protected:
482 void update_(const ElementContext&,
483 unsigned,
484 unsigned)
485 {}
486
487 template <class Context, class FluidState>
488 void updateBoundary_(const Context&,
489 unsigned,
490 unsigned,
491 const FluidState&)
492 {}
493
494public:
499 const Scalar& diffusivity() const
500 {
501 throw std::logic_error("The method diffusivity() does not "
502 "make sense if diffusion is disabled.");
503 }
504
512 const Evaluation& effectiveDiffusionCoefficient(unsigned,
513 unsigned) const
514 {
515 throw std::logic_error("The method effectiveDiffusionCoefficient() "
516 "does not make sense if diffusion is disabled.");
517 }
518
519};
520
524template <class TypeTag>
525class BlackOilDiffusionExtensiveQuantities<TypeTag, /*enableDiffusion=*/true>
526{
532 using Toolbox = MathToolbox<Evaluation>;
534
535 enum { dimWorld = GridView::dimensionworld };
536 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
537 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
538
539 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
540 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
541public:
542 using EvaluationArray = Evaluation[numPhases][numComponents];
543protected:
548 void update_(const ElementContext& elemCtx, unsigned faceIdx, unsigned timeIdx)
549 {
550 // Only work if diffusion is enabled run-time by DIFFUSE in the deck
551 if(!FluidSystem::enableDiffusion())
552 return;
553
554 const auto& stencil = elemCtx.stencil(timeIdx);
555 const auto& face = stencil.interiorFace(faceIdx);
556 const auto& extQuants = elemCtx.extensiveQuantities(faceIdx, timeIdx);
557 const auto& intQuantsInside = elemCtx.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
558 const auto& intQuantsOutside = elemCtx.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
559
560 const Scalar diffusivity = elemCtx.problem().diffusivity(elemCtx, face.interiorIndex(), face.exteriorIndex());
561 const Scalar faceArea = face.area();
562 diffusivity_ = diffusivity / faceArea;
563 update(effectiveDiffusionCoefficient_, intQuantsInside, intQuantsOutside);
564 Valgrind::CheckDefined(diffusivity_);
565 }
566
567public:
568 static void update(EvaluationArray& effectiveDiffusionCoefficient,
569 const IntensiveQuantities& intQuantsInside,
570 const IntensiveQuantities& intQuantsOutside) {
571 // opm-models expects per area flux
572 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
573 if (!FluidSystem::phaseIsActive(phaseIdx)) {
574 continue;
575 }
576 // no diffusion in water for blackoil models
577 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
578 continue;
579 }
580 for (unsigned compIdx = 0; compIdx < numComponents; ++compIdx) {
581 // use the arithmetic average for the effective
582 // diffusion coefficients.
583 effectiveDiffusionCoefficient[phaseIdx][compIdx] = 0.5 *
584 ( intQuantsInside.effectiveDiffusionCoefficient(phaseIdx, compIdx) +
585 intQuantsOutside.effectiveDiffusionCoefficient(phaseIdx, compIdx) );
586 Valgrind::CheckDefined(effectiveDiffusionCoefficient[phaseIdx][compIdx]);
587 }
588 }
589 }
590protected:
591 template <class Context, class FluidState>
592 void updateBoundary_(const Context&,
593 unsigned,
594 unsigned,
595 const FluidState&)
596 {
597 throw std::runtime_error("Not implemented: Diffusion across boundary not implemented for blackoil");
598 }
599
600public:
607 const Scalar& diffusivity() const
608 { return diffusivity_; }
609
617 const Evaluation& effectiveDiffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
618 { return effectiveDiffusionCoefficient_[phaseIdx][compIdx]; }
619
620 const auto& effectiveDiffusionCoefficient() const{
621 return effectiveDiffusionCoefficient_;
622 }
623
624private:
625 Scalar diffusivity_;
626 EvaluationArray effectiveDiffusionCoefficient_;
627};
628
629} // namespace Opm
630
631#endif
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:512
const Scalar & diffusivity() const
The diffusivity the face.
Definition: blackoildiffusionmodule.hh:499
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:482
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildiffusionmodule.hh:488
Provides the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:526
static void update(EvaluationArray &effectiveDiffusionCoefficient, const IntensiveQuantities &intQuantsInside, const IntensiveQuantities &intQuantsOutside)
Definition: blackoildiffusionmodule.hh:568
Evaluation[numPhases][numComponents] EvaluationArray
Definition: blackoildiffusionmodule.hh:542
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: blackoildiffusionmodule.hh:617
const auto & effectiveDiffusionCoefficient() const
Definition: blackoildiffusionmodule.hh:620
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildiffusionmodule.hh:592
void update_(const ElementContext &elemCtx, unsigned faceIdx, unsigned timeIdx)
Update the quantities required to calculate the diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:548
const Scalar & diffusivity() const
The diffusivity of the face.
Definition: blackoildiffusionmodule.hh:607
Provides the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:51
Scalar effectiveDiffusionCoefficient(unsigned, unsigned) const
Returns the effective molecular diffusion coefficient of the porous medium for a component in a phase...
Definition: blackoildiffusionmodule.hh:311
Scalar diffusionCoefficient(unsigned, unsigned) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: blackoildiffusionmodule.hh:301
void update_(FluidState &, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &, const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:323
Scalar tortuosity(unsigned) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: blackoildiffusionmodule.hh:291
void update_(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache, const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the quantities required to calculate diffusive mass fluxes.
Definition: blackoildiffusionmodule.hh:402
Evaluation diffusionCoefficient(unsigned phaseIdx, unsigned compIdx) const
Returns the molecular diffusion coefficient for a component in a phase.
Definition: blackoildiffusionmodule.hh:371
Evaluation tortuosity(unsigned phaseIdx) const
Returns the tortuousity of the sub-domain of a fluid phase in the porous medium.
Definition: blackoildiffusionmodule.hh:378
BlackOilDiffusionIntensiveQuantities(BlackOilDiffusionIntensiveQuantities &&) noexcept=default
void update_(FluidState &fluidState, typename FluidSystem::template ParameterCache< typename FluidState::Scalar > &paramCache, const IntensiveQuantities &intQuants)
Definition: blackoildiffusionmodule.hh:416
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:385
Provides the volumetric quantities required for the calculation of molecular diffusive fluxes.
Definition: blackoildiffusionmodule.hh:274
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: blackoildiffusionmodule.hh:77
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:85
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:153
static void addDiffusiveFlux(RateVector &flux, const FluidState &fluidStateI, const FluidState &fluidStateJ, const Evaluation &diffusivity, const EvaluationArray &effectiveDiffusionCoefficient)
Definition: blackoildiffusionmodule.hh:168
static void registerParameters()
Register all run-time parameters for the diffusion module.
Definition: blackoildiffusionmodule.hh:127
Provides the auxiliary methods required for consideration of the diffusion equation.
Definition: blackoildiffusionmodule.hh:48
Declare the properties used by the infrastructure code of the finite volume discretizations.
Definition: blackoilboundaryratevector.hh:37
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:235