My Project
Loading...
Searching...
No Matches
blackoildispersionmodule.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 EWOMS_DISPERSION_MODULE_HH
29#define EWOMS_DISPERSION_MODULE_HH
30
31#include <dune/common/fvector.hh>
32
35
36#include <opm/material/common/MathToolbox.hpp>
37#include <opm/material/common/Valgrind.hpp>
38
39#include <stdexcept>
40
41namespace Opm {
42
49template <class TypeTag, bool enableDispersion>
51
52template <class TypeTag, bool enableDispersion>
54
58template <class TypeTag>
59class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/false>
60{
65
66 enum { numPhases = FluidSystem::numPhases };
67
68public:
70
71#if HAVE_ECL_INPUT
72 static void initFromState(const EclipseState&)
73 {
74 }
75#endif
76
81 template <class Context>
82 static void addDispersiveFlux(RateVector&,
83 const Context&,
84 unsigned,
85 unsigned)
86 {}
87
88 template<class FluidState, class Scalar>
89 static void addDispersiveFlux(RateVector&,
90 const FluidState&,
91 const FluidState&,
92 const Evaluation&,
93 const Scalar&)
94 {}
95};
96
100template <class TypeTag>
101class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/true>
102{
114
115 enum { numPhases = FluidSystem::numPhases };
116 enum { numComponents = FluidSystem::numComponents };
117 enum { conti0EqIdx = Indices::conti0EqIdx };
118 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
119
120 using Toolbox = MathToolbox<Evaluation>;
121
122public:
124#if HAVE_ECL_INPUT
125 static void initFromState(const EclipseState& eclState)
126 {
127 if (!eclState.getSimulationConfig().rock_config().dispersion()) {
128 return;
129 }
130
131 if (eclState.getSimulationConfig().hasVAPWAT() || eclState.getSimulationConfig().hasVAPOIL()) {
132 OpmLog::warning("Dispersion is activated in combination with VAPWAT/VAPOIL. \n"
133 "Water/oil is still allowed to vaporize, but dispersion in the "
134 "gas phase is ignored.");
135 }
136 }
137#endif
142 template <class Context>
143 static void addDispersiveFlux(RateVector& flux, const Context& context,
144 unsigned spaceIdx, unsigned timeIdx)
145 {
146 // Only work if dispersion is enabled by DISPERC in the deck
147 if (!context.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
148 return;
149 }
150 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
151 const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
152 const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
153 const auto& dispersivity = extQuants.dispersivity();
154 const auto& normVelocityAvg = extQuants.normVelocityAvg();
155 addDispersiveFlux(flux, fluidStateI, fluidStateJ, dispersivity, normVelocityAvg);
156 }
157
178 template<class FluidState, class Scalar>
179 static void addDispersiveFlux(RateVector& flux,
180 const FluidState& fluidStateI,
181 const FluidState& fluidStateJ,
182 const Evaluation& dispersivity,
183 const Scalar& normVelocityAvg)
184 {
185 unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex();
186 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
187 if (!FluidSystem::phaseIsActive(phaseIdx)) {
188 continue;
189 }
190
191 // no dispersion in water for blackoil models unless water can contain dissolved gas
192 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
193 continue;
194 }
195
196 // Adding dispersion in the gas phase leads to
197 // convergence issues and unphysical results.
198 // We disable dispersion in the gas phase for now
199 // See comment below
200 if (FluidSystem::gasPhaseIdx == phaseIdx) {
201 continue;
202 }
203
204 // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
205 // phase check disabled due to if above, reenable when removing unconditional gas phase disablement
206 if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil())
207 /*&& FluidSystem::gasPhaseIdx == phaseIdx*/) {
208 continue;
209 }
210
211 // arithmetic mean of the phase's b factor
212 Evaluation bAvg = fluidStateI.invB(phaseIdx);
213 bAvg += Toolbox::value(fluidStateJ.invB(phaseIdx));
214 bAvg /= 2;
215
216 Evaluation convFactor = 1.0;
217 Evaluation diffR = 0.0;
218 if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) {
219 Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2;
220 convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg);
221 diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs());
222 }
223 if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) {
224 Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2;
225 convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex));
226 diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv());
227 }
228 if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
229 Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2;
230 convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg);
231 diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw());
232 }
233 if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
234 Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2;
235 convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex));
236 diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw());
237 }
238
239 // mass flux of solvent component
240 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
241 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
242 flux[conti0EqIdx + activeSolventCompIdx] +=
243 - bAvg
244 * normVelocityAvg[phaseIdx]
245 * convFactor
246 * dispersivity
247 * diffR;
248
249 // mass flux of solute component
250 unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
251 unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx);
252 flux[conti0EqIdx + activeSoluteCompIdx] +=
253 bAvg
254 * normVelocityAvg[phaseIdx]
255 * convFactor
256 * dispersivity
257 * diffR;
258 }
259 }
260
261private:
262
263 static Scalar toMassFractionGasOil (unsigned regionIdx) {
264 Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
265 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
266 return rhoO / rhoG;
267 }
268 static Scalar toMassFractionGasWater (unsigned regionIdx) {
269 Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
270 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
271 return rhoW / rhoG;
272 }
273};
274
282template <class TypeTag, bool enableDispersion>
284
288template <class TypeTag>
289class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/false>
290{
294
295public:
299 Scalar normVelocityCell(unsigned, unsigned) const
300 {
301 throw std::logic_error("Method normVelocityCell() "
302 "does not make sense if dispersion is disabled");
303 }
304
305protected:
310 template<class ElementContext>
311 void update_(ElementContext&,
312 unsigned,
313 unsigned)
314 { }
315};
316
320template <class TypeTag>
321class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/true>
322{
329 enum { numPhases = FluidSystem::numPhases };
330 enum { numComponents = FluidSystem::numComponents };
331 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
332 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
333 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
334 enum { gasCompIdx = FluidSystem::gasCompIdx };
335 enum { oilCompIdx = FluidSystem::oilCompIdx };
336 enum { waterCompIdx = FluidSystem::waterCompIdx };
337 enum { conti0EqIdx = Indices::conti0EqIdx };
338 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
339
340public:
344 Scalar normVelocityCell(unsigned phaseIdx) const
345 {
346
347 return normVelocityCell_[phaseIdx];
348 }
349
350protected:
360 template<class ElementContext>
361 void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
362 {
363 // Only work if dispersion is enabled by DISPERC in the deck
364 if (!elemCtx.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
365 return;
366 }
367 const auto& problem = elemCtx.simulator().problem();
368 if (problem.model().linearizer().getVelocityInfo().empty()) {
369 return;
370 }
371 const std::array<int, 3> phaseIdxs = { gasPhaseIdx, oilPhaseIdx, waterPhaseIdx };
372 const std::array<int, 3> compIdxs = { gasCompIdx, oilCompIdx, waterCompIdx };
373 const auto& velocityInf = problem.model().linearizer().getVelocityInfo();
374 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
375 auto velocityInfos = velocityInf[globalDofIdx];
376 for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
377 normVelocityCell_[i] = 0;
378 }
379 for (const auto& velocityInfo : velocityInfos) {
380 for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
381 if (FluidSystem::phaseIsActive(phaseIdxs[i])) {
382 normVelocityCell_[phaseIdxs[i]] = max( normVelocityCell_[phaseIdxs[i]],
383 std::abs( velocityInfo.velocity[conti0EqIdx
384 + Indices::canonicalToActiveComponentIndex(compIdxs[i])] ));
385 }
386 }
387 }
388 }
389
390private:
391 Scalar normVelocityCell_[numPhases];
392};
393
400template <class TypeTag, bool enableDispersion>
401class BlackOilDispersionExtensiveQuantities;
402
406template <class TypeTag>
407class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/false>
408{
414
415 enum { numPhases = FluidSystem::numPhases };
416
417protected:
422 void update_(const ElementContext&,
423 unsigned,
424 unsigned)
425 {}
426
427 template <class Context, class FluidState>
428 void updateBoundary_(const Context&,
429 unsigned,
430 unsigned,
431 const FluidState&)
432 {}
433
434public:
435 using ScalarArray = Scalar[numPhases];
436
437 static void update(ScalarArray&,
438 const IntensiveQuantities&,
439 const IntensiveQuantities&)
440 {}
441
446 const Scalar& dispersivity() const
447 {
448 throw std::logic_error("The method dispersivity() does not "
449 "make sense if dispersion is disabled.");
450 }
451
459 const Scalar& normVelocityAvg(unsigned) const
460 {
461 throw std::logic_error("The method normVelocityAvg() "
462 "does not make sense if dispersion is disabled.");
463 }
464
465};
466
470template <class TypeTag>
471class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/true>
472{
478 using Toolbox = MathToolbox<Evaluation>;
480
481 enum { dimWorld = GridView::dimensionworld };
483 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
484
485 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
486 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
487public:
488 using ScalarArray = Scalar[numPhases];
489 static void update(ScalarArray& normVelocityAvg,
490 const IntensiveQuantities& intQuantsInside,
491 const IntensiveQuantities& intQuantsOutside)
492 {
493 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
494 if (!FluidSystem::phaseIsActive(phaseIdx)) {
495 continue;
496 }
497 // no dispersion in water for blackoil models unless water can contain dissolved gas
498 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
499 continue;
500 }
501 // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
502 if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) {
503 continue;
504 }
505 // use the arithmetic average for the effective
506 // velocity coefficients at the face's integration point.
507 normVelocityAvg[phaseIdx] = 0.5 *
508 ( intQuantsInside.normVelocityCell(phaseIdx) +
509 intQuantsOutside.normVelocityCell(phaseIdx) );
510 Valgrind::CheckDefined(normVelocityAvg[phaseIdx]);
511 }
512 }
513protected:
514 template <class Context, class FluidState>
515 void updateBoundary_(const Context&,
516 unsigned,
517 unsigned,
518 const FluidState&)
519 {
520 throw std::runtime_error("Not implemented: Dispersion across boundary not implemented for blackoil");
521 }
522
523public:
530 const Scalar& dispersivity() const
531 { return dispersivity_; }
532
540 const Scalar& normVelocityAvg(unsigned phaseIdx) const
541 { return normVelocityAvg_[phaseIdx]; }
542
543 const auto& normVelocityAvg() const{
544 return normVelocityAvg_;
545 }
546
547private:
548 Scalar dispersivity_;
549 ScalarArray normVelocityAvg_;
550};
551
552} // namespace Opm
553
554#endif
Definition blackoildispersionmodule.hh:408
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the dispersive fluxes.
Definition blackoildispersionmodule.hh:422
const Scalar & dispersivity() const
The dispersivity the face.
Definition blackoildispersionmodule.hh:446
const Scalar & normVelocityAvg(unsigned) const
The effective filter velocity coefficient in a fluid phase at the face's integration point.
Definition blackoildispersionmodule.hh:459
Provides the quantities required to calculate dispersive mass fluxes.
Definition blackoildispersionmodule.hh:472
const Scalar & dispersivity() const
The dispersivity of the face.
Definition blackoildispersionmodule.hh:530
const Scalar & normVelocityAvg(unsigned phaseIdx) const
The effective velocity coefficient in a fluid phase at the face's integration point.
Definition blackoildispersionmodule.hh:540
Provides the quantities required to calculate dispersive mass fluxes.
Definition blackoildispersionmodule.hh:53
Scalar normVelocityCell(unsigned, unsigned) const
Returns the max.
Definition blackoildispersionmodule.hh:299
void update_(ElementContext &, unsigned, unsigned)
Update the quantities required to calculate dispersive fluxes.
Definition blackoildispersionmodule.hh:311
void update_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the quantities required to calculate dispersive mass fluxes.
Definition blackoildispersionmodule.hh:361
Scalar normVelocityCell(unsigned phaseIdx) const
Returns the max.
Definition blackoildispersionmodule.hh:344
Provides the volumetric quantities required for the calculation of dispersive fluxes.
Definition blackoildispersionmodule.hh:283
static void addDispersiveFlux(RateVector &, const Context &, unsigned, unsigned)
Adds the dispersive flux to the flux vector over a flux integration point.
Definition blackoildispersionmodule.hh:82
static void addDispersiveFlux(RateVector &flux, const FluidState &fluidStateI, const FluidState &fluidStateJ, const Evaluation &dispersivity, const Scalar &normVelocityAvg)
Adds the mass flux due to dispersion to the flux vector over the integration point.
Definition blackoildispersionmodule.hh:179
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:143
Provides the auxiliary methods required for consideration of the dispersion equation.
Definition blackoildispersionmodule.hh:50
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.
This file contains a set of helper functions used by VFPProd / VFPInj.
Definition blackoilboundaryratevector.hh:37
constexpr auto getPropValue()
get the value data member of a property
Definition propertysystem.hh:242
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