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
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 enableDispersion>
49
50template <class TypeTag, bool enableDispersion>
52
56template <class TypeTag>
57class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/false>
58{
63
64 enum { numPhases = FluidSystem::numPhases };
65
66public:
68
69#if HAVE_ECL_INPUT
70 static void initFromState(const EclipseState&)
71 {
72 }
73#endif
74
79 template <class Context>
80 static void addDispersiveFlux(RateVector&,
81 const Context&,
82 unsigned,
83 unsigned)
84 {}
85
86 template<class FluidState, class Scalar>
87 static void addDispersiveFlux(RateVector&,
88 const FluidState&,
89 const FluidState&,
90 const Evaluation&,
91 const Scalar&)
92 {}
93};
94
98template <class TypeTag>
99class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/true>
100{
112
113 enum { numPhases = FluidSystem::numPhases };
114 enum { numComponents = FluidSystem::numComponents };
115 enum { conti0EqIdx = Indices::conti0EqIdx };
116 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
117
118 using Toolbox = MathToolbox<Evaluation>;
119
120public:
122#if HAVE_ECL_INPUT
123 static void initFromState(const EclipseState& eclState)
124 {
125 if (!eclState.getSimulationConfig().rock_config().dispersion()) {
126 return;
127 }
128
129 if (eclState.getSimulationConfig().hasVAPWAT() || eclState.getSimulationConfig().hasVAPOIL()) {
130 OpmLog::warning("Dispersion is activated in combination with VAPWAT/VAPOIL. \n"
131 "Water/oil is still allowed to vaporize, but dispersion in the "
132 "gas phase is ignored.");
133 }
134 }
135#endif
140 template <class Context>
141 static void addDispersiveFlux(RateVector& flux, const Context& context,
142 unsigned spaceIdx, unsigned timeIdx)
143 {
144 // Only work if dispersion is enabled by DISPERC in the deck
145 if (!context.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
146 return;
147 }
148 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
149 const auto& fluidStateI = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx).fluidState();
150 const auto& fluidStateJ = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx).fluidState();
151 const auto& dispersivity = extQuants.dispersivity();
152 const auto& normVelocityAvg = extQuants.normVelocityAvg();
153 addDispersiveFlux(flux, fluidStateI, fluidStateJ, dispersivity, normVelocityAvg);
154 }
155
176 template<class FluidState, class Scalar>
177 static void addDispersiveFlux(RateVector& flux,
178 const FluidState& fluidStateI,
179 const FluidState& fluidStateJ,
180 const Evaluation& dispersivity,
181 const Scalar& normVelocityAvg)
182 {
183 unsigned pvtRegionIndex = fluidStateI.pvtRegionIndex();
184 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
185 if (!FluidSystem::phaseIsActive(phaseIdx)) {
186 continue;
187 }
188
189 // no dispersion in water for blackoil models unless water can contain dissolved gas
190 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
191 continue;
192 }
193
194 // Adding dispersion in the gas phase leads to
195 // convergence issues and unphysical results.
196 // We disable dispersion in the gas phase for now
197 if (FluidSystem::gasPhaseIdx == phaseIdx) {
198 continue;
199 }
200
201 // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
202 if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) {
203 continue;
204 }
205
206 // arithmetic mean of the phase's b factor
207 Evaluation bAvg = fluidStateI.invB(phaseIdx);
208 bAvg += Toolbox::value(fluidStateJ.invB(phaseIdx));
209 bAvg /= 2;
210
211 Evaluation convFactor = 1.0;
212 Evaluation diffR = 0.0;
213 if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) {
214 Evaluation rsAvg = (fluidStateI.Rs() + Toolbox::value(fluidStateJ.Rs())) / 2;
215 convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg);
216 diffR = fluidStateI.Rs() - Toolbox::value(fluidStateJ.Rs());
217 }
218 if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) {
219 Evaluation rvAvg = (fluidStateI.Rv() + Toolbox::value(fluidStateJ.Rv())) / 2;
220 convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex));
221 diffR = fluidStateI.Rv() - Toolbox::value(fluidStateJ.Rv());
222 }
223 if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
224 Evaluation rsAvg = (fluidStateI.Rsw() + Toolbox::value(fluidStateJ.Rsw())) / 2;
225 convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg);
226 diffR = fluidStateI.Rsw() - Toolbox::value(fluidStateJ.Rsw());
227 }
228 if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
229 Evaluation rvAvg = (fluidStateI.Rvw() + Toolbox::value(fluidStateJ.Rvw())) / 2;
230 convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex));
231 diffR = fluidStateI.Rvw() - Toolbox::value(fluidStateJ.Rvw());
232 }
233
234 // mass flux of solvent component
235 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
236 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
237 flux[conti0EqIdx + activeSolventCompIdx] +=
238 - bAvg
239 * normVelocityAvg[phaseIdx]
240 * convFactor
241 * dispersivity
242 * diffR;
243
244 // mass flux of solute component
245 unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
246 unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx);
247 flux[conti0EqIdx + activeSoluteCompIdx] +=
248 bAvg
249 * normVelocityAvg[phaseIdx]
250 * convFactor
251 * dispersivity
252 * diffR;
253 }
254 }
255
256private:
257
258 static Scalar toMassFractionGasOil (unsigned regionIdx) {
259 Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
260 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
261 return rhoO / rhoG;
262 }
263 static Scalar toMassFractionGasWater (unsigned regionIdx) {
264 Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
265 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
266 return rhoW / rhoG;
267 }
268};
269
277template <class TypeTag, bool enableDispersion>
279
283template <class TypeTag>
284class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/false>
285{
289
290public:
294 Scalar normVelocityCell(unsigned, unsigned) const
295 {
296 throw std::logic_error("Method normVelocityCell() "
297 "does not make sense if dispersion is disabled");
298 }
299
300protected:
305 template<class ElementContext>
306 void update_(ElementContext&,
307 unsigned,
308 unsigned)
309 { }
310};
311
315template <class TypeTag>
316class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/true>
317{
324 enum { numPhases = FluidSystem::numPhases };
325 enum { numComponents = FluidSystem::numComponents };
326 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
327 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
328 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
329 enum { gasCompIdx = FluidSystem::gasCompIdx };
330 enum { oilCompIdx = FluidSystem::oilCompIdx };
331 enum { waterCompIdx = FluidSystem::waterCompIdx };
332 enum { conti0EqIdx = Indices::conti0EqIdx };
333 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
334
335public:
339 Scalar normVelocityCell(unsigned phaseIdx) const
340 {
341
342 return normVelocityCell_[phaseIdx];
343 }
344
345protected:
355 template<class ElementContext>
356 void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
357 {
358 // Only work if dispersion is enabled by DISPERC in the deck
359 if (!elemCtx.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
360 return;
361 }
362 const auto& problem = elemCtx.simulator().problem();
363 if (problem.model().linearizer().getVelocityInfo().empty()) {
364 return;
365 }
366 const std::array<int, 3> phaseIdxs = { gasPhaseIdx, oilPhaseIdx, waterPhaseIdx };
367 const std::array<int, 3> compIdxs = { gasCompIdx, oilCompIdx, waterCompIdx };
368 const auto& velocityInf = problem.model().linearizer().getVelocityInfo();
369 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
370 auto velocityInfos = velocityInf[globalDofIdx];
371 for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
372 normVelocityCell_[i] = 0;
373 }
374 for (auto& velocityInfo : velocityInfos) {
375 for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
376 if (FluidSystem::phaseIsActive(phaseIdxs[i])) {
377 normVelocityCell_[phaseIdxs[i]] = max( normVelocityCell_[phaseIdxs[i]],
378 std::abs( velocityInfo.velocity[conti0EqIdx
379 + Indices::canonicalToActiveComponentIndex(compIdxs[i])] ));
380 }
381 }
382 }
383 }
384
385private:
386 Scalar normVelocityCell_[numPhases];
387};
388
395template <class TypeTag, bool enableDispersion>
396class BlackOilDispersionExtensiveQuantities;
397
401template <class TypeTag>
402class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/false>
403{
409
410 enum { numPhases = FluidSystem::numPhases };
411
412protected:
417 void update_(const ElementContext&,
418 unsigned,
419 unsigned)
420 {}
421
422 template <class Context, class FluidState>
423 void updateBoundary_(const Context&,
424 unsigned,
425 unsigned,
426 const FluidState&)
427 {}
428
429public:
430 using ScalarArray = Scalar[numPhases];
431
432 static void update(ScalarArray&,
433 const IntensiveQuantities&,
434 const IntensiveQuantities&)
435 {}
436
441 const Scalar& dispersivity() const
442 {
443 throw std::logic_error("The method dispersivity() does not "
444 "make sense if dispersion is disabled.");
445 }
446
454 const Scalar& normVelocityAvg(unsigned) const
455 {
456 throw std::logic_error("The method normVelocityAvg() "
457 "does not make sense if dispersion is disabled.");
458 }
459
460};
461
465template <class TypeTag>
466class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/true>
467{
473 using Toolbox = MathToolbox<Evaluation>;
475
476 enum { dimWorld = GridView::dimensionworld };
477 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
478 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
479
480 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
481 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
482public:
483 using ScalarArray = Scalar[numPhases];
484 static void update(ScalarArray& normVelocityAvg,
485 const IntensiveQuantities& intQuantsInside,
486 const IntensiveQuantities& intQuantsOutside)
487 {
488 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
489 if (!FluidSystem::phaseIsActive(phaseIdx)) {
490 continue;
491 }
492 // no dispersion in water for blackoil models unless water can contain dissolved gas
493 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
494 continue;
495 }
496 // no dispersion in gas for blackoil models unless gas can contain evaporated water or oil
497 if ((!FluidSystem::enableVaporizedWater() && !FluidSystem::enableVaporizedOil()) && FluidSystem::gasPhaseIdx == phaseIdx) {
498 continue;
499 }
500 // use the arithmetic average for the effective
501 // velocity coefficients at the face's integration point.
502 normVelocityAvg[phaseIdx] = 0.5 *
503 ( intQuantsInside.normVelocityCell(phaseIdx) +
504 intQuantsOutside.normVelocityCell(phaseIdx) );
505 Valgrind::CheckDefined(normVelocityAvg[phaseIdx]);
506 }
507 }
508protected:
509 template <class Context, class FluidState>
510 void updateBoundary_(const Context&,
511 unsigned,
512 unsigned,
513 const FluidState&)
514 {
515 throw std::runtime_error("Not implemented: Dispersion across boundary not implemented for blackoil");
516 }
517
518public:
525 const Scalar& dispersivity() const
526 { return dispersivity_; }
527
535 const Scalar& normVelocityAvg(unsigned phaseIdx) const
536 { return normVelocityAvg_[phaseIdx]; }
537
538 const auto& normVelocityAvg() const{
539 return normVelocityAvg_;
540 }
541
542private:
543 Scalar dispersivity_;
544 ScalarArray normVelocityAvg_;
545};
546
547} // namespace Opm
548
549#endif
Definition: blackoildispersionmodule.hh:403
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildispersionmodule.hh:423
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the dispersive fluxes.
Definition: blackoildispersionmodule.hh:417
const Scalar & dispersivity() const
The dispersivity the face.
Definition: blackoildispersionmodule.hh:441
Scalar[numPhases] ScalarArray
Definition: blackoildispersionmodule.hh:430
static void update(ScalarArray &, const IntensiveQuantities &, const IntensiveQuantities &)
Definition: blackoildispersionmodule.hh:432
const Scalar & normVelocityAvg(unsigned) const
The effective filter velocity coefficient in a fluid phase at the face's integration point.
Definition: blackoildispersionmodule.hh:454
Provides the quantities required to calculate dispersive mass fluxes.
Definition: blackoildispersionmodule.hh:467
static void update(ScalarArray &normVelocityAvg, const IntensiveQuantities &intQuantsInside, const IntensiveQuantities &intQuantsOutside)
Definition: blackoildispersionmodule.hh:484
const Scalar & dispersivity() const
The dispersivity of the face.
Definition: blackoildispersionmodule.hh:525
Scalar[numPhases] ScalarArray
Definition: blackoildispersionmodule.hh:483
void updateBoundary_(const Context &, unsigned, unsigned, const FluidState &)
Definition: blackoildispersionmodule.hh:510
const auto & normVelocityAvg() const
Definition: blackoildispersionmodule.hh:538
const Scalar & normVelocityAvg(unsigned phaseIdx) const
The effective velocity coefficient in a fluid phase at the face's integration point.
Definition: blackoildispersionmodule.hh:535
Provides the quantities required to calculate dispersive mass fluxes.
Definition: blackoildispersionmodule.hh:51
Scalar normVelocityCell(unsigned, unsigned) const
Returns the max. norm of the filter velocity of the cell.
Definition: blackoildispersionmodule.hh:294
void update_(ElementContext &, unsigned, unsigned)
Update the quantities required to calculate dispersive fluxes.
Definition: blackoildispersionmodule.hh:306
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:356
Scalar normVelocityCell(unsigned phaseIdx) const
Returns the max. norm of the filter velocity of the cell.
Definition: blackoildispersionmodule.hh:339
Provides the volumetric quantities required for the calculation of dispersive fluxes.
Definition: blackoildispersionmodule.hh:278
static void addDispersiveFlux(RateVector &, const Context &, unsigned, unsigned)
Adds the dispersive flux to the flux vector over a flux integration point.
Definition: blackoildispersionmodule.hh:80
static void addDispersiveFlux(RateVector &, const FluidState &, const FluidState &, const Evaluation &, const Scalar &)
Definition: blackoildispersionmodule.hh:87
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:177
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:141
Provides the auxiliary methods required for consideration of the dispersion equation.
Definition: blackoildispersionmodule.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