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
33#include <opm/material/common/MathToolbox.hpp>
34#include <opm/material/common/Valgrind.hpp>
35
39
40#if HAVE_ECL_INPUT
41#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
42#include <opm/common/OpmLog/OpmLog.hpp>
43#endif
44
45#include <array>
46#include <cmath>
47#include <stdexcept>
48
49namespace Opm {
50
57template <class TypeTag, bool enableDispersion>
59
60template <class TypeTag, bool enableDispersion>
62
66template <class TypeTag>
67class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/false>
68{
74
75public:
77
78#if HAVE_ECL_INPUT
79 static void initFromState(const EclipseState&)
80 {}
81#endif
82
87 template <class Context>
88 static void addDispersiveFlux(RateVector&,
89 const Context&,
90 unsigned,
91 unsigned)
92 {}
93
94 template<class IntensiveQuantities, class Scalar>
95 static void addDispersiveFlux(RateVector&,
96 const IntensiveQuantities&,
97 const IntensiveQuantities&,
98 const Evaluation&,
99 const Scalar&)
100 {}
101};
102
106template <class TypeTag>
107class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/true>
108{
119
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 };
126
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;
131
132 using Toolbox = MathToolbox<Evaluation>;
133
134public:
136
137#if HAVE_ECL_INPUT
138 static void initFromState(const EclipseState& eclState)
139 {
140 if (!eclState.getSimulationConfig().rock_config().dispersion()) {
141 return;
142 }
143
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.");
148 }
149 }
150#endif
151
156 template <class Context>
157 static void addDispersiveFlux(RateVector& flux, const Context& context,
158 unsigned spaceIdx, unsigned timeIdx)
159 {
160 // Only work if dispersion is enabled by DISPERC in the deck
161 if (!context.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
162 return;
163 }
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);
170 }
171
192 template<class IntensiveQuantities, class Scalar>
193 static void addDispersiveFlux(RateVector& flux,
194 const IntensiveQuantities& inIq,
195 const IntensiveQuantities& exIq,
196 const Evaluation& dispersivity,
197 const Scalar& normVelocityAvg)
198 {
199 const auto& inFs = inIq.fluidState();
200 const auto& exFs = exIq.fluidState();
201 Evaluation diffR = 0.0;
202 if constexpr(enableBioeffects) {
203 // The dispersion coefficients are given for mass concentrations
204 const Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2;
205 diffR = inIq.microbialConcentration() - Toolbox::value(exIq.microbialConcentration());
206 flux[contiMicrobialEqIdx] +=
207 bAvg *
208 normVelocityAvg[waterPhaseIdx] *
209 dispersivity *
210 diffR;
211 if constexpr(enableMICP) {
212 diffR = inIq.oxygenConcentration() - Toolbox::value(exIq.oxygenConcentration());
213 flux[contiOxygenEqIdx] +=
214 bAvg *
215 normVelocityAvg[waterPhaseIdx] *
216 dispersivity *
217 diffR;
218 diffR = inIq.ureaConcentration() - Toolbox::value(exIq.ureaConcentration());
219 flux[contiUreaEqIdx] +=
220 bAvg *
221 normVelocityAvg[waterPhaseIdx] *
222 dispersivity *
223 diffR;
224 return;
225 }
226 }
227
228 unsigned pvtRegionIndex = inFs.pvtRegionIndex();
229 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
230 if (!FluidSystem::phaseIsActive(phaseIdx)) {
231 continue;
232 }
233
234 // no dispersion in water for blackoil models unless water can contain dissolved gas
235 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
236 continue;
237 }
238
239 // adding dispersion in the gas phase leads to
240 // convergence issues and unphysical results.
241 // we disable dispersion in the gas phase for now
242 if (FluidSystem::gasPhaseIdx == phaseIdx) {
243 continue;
244 }
245
246 // arithmetic mean of the phase's b factor
247 Evaluation bAvg = inFs.invB(phaseIdx);
248 bAvg += Toolbox::value(exFs.invB(phaseIdx));
249 bAvg /= 2;
250
251 Evaluation convFactor = 1.0;
252 if (FluidSystem::enableDissolvedGas() &&
253 FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) &&
254 phaseIdx == FluidSystem::oilPhaseIdx)
255 {
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());
259 }
260 if (FluidSystem::enableVaporizedOil() &&
261 FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) &&
262 phaseIdx == FluidSystem::gasPhaseIdx)
263 {
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());
268 }
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());
273 }
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());
279 }
280
281 // mass flux of solvent component
282 const unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
283 const unsigned activeSolventCompIdx = FluidSystem::canonicalToActiveCompIdx(solventCompIdx);
284 flux[conti0EqIdx + activeSolventCompIdx] +=
285 -bAvg *
286 normVelocityAvg[phaseIdx] *
287 convFactor *
288 dispersivity *
289 diffR;
290
291 // mass flux of solute component
292 const unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
293 const unsigned activeSoluteCompIdx = FluidSystem::canonicalToActiveCompIdx(soluteCompIdx);
294 flux[conti0EqIdx + activeSoluteCompIdx] +=
295 bAvg *
296 normVelocityAvg[phaseIdx] *
297 convFactor *
298 dispersivity *
299 diffR;
300 }
301 }
302
303private:
304 static Scalar toMassFractionGasOil (unsigned regionIdx)
305 {
306 const Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
307 const Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
308 return rhoO / rhoG;
309 }
310
311 static Scalar toMassFractionGasWater (unsigned regionIdx)
312 {
313 const Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
314 const Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
315 return rhoW / rhoG;
316 }
317};
318
326template <class TypeTag, bool enableDispersion>
328
332template <class TypeTag>
333class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/false>
334{
337
338public:
342 Scalar normVelocityCell(unsigned, unsigned) const
343 {
344 throw std::logic_error("Method normVelocityCell() "
345 "does not make sense if dispersion is disabled");
346 }
347
348protected:
353 template<class ElementContext>
354 void update_(ElementContext&,
355 unsigned,
356 unsigned)
357 {}
358};
359
363template <class TypeTag>
364class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/true>
365{
370
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>() };
381
382public:
386 Scalar normVelocityCell(unsigned phaseIdx) const
387 { return normVelocityCell_[phaseIdx]; }
388
389protected:
399 template<class ElementContext>
400 void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
401 {
402 // Only work if dispersion is enabled by DISPERC in the deck
403 if (!elemCtx.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
404 return;
405 }
406 const auto& problem = elemCtx.simulator().problem();
407 if (problem.model().linearizer().getVelocityInfo().empty()) {
408 return;
409 }
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;
417 }
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])]));
424 }
425 }
426 }
427 }
428
429private:
430 std::array<Scalar, numPhases> normVelocityCell_{};
431};
432
439template <class TypeTag, bool enableDispersion>
440class BlackOilDispersionExtensiveQuantities;
441
445template <class TypeTag>
446class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/false>
447{
452
453 enum { numPhases = FluidSystem::numPhases };
454
455protected:
460 void update_(const ElementContext&,
461 unsigned,
462 unsigned)
463 {}
464
465 template <class Context, class FluidState>
466 void updateBoundary_(const Context&,
467 unsigned,
468 unsigned,
469 const FluidState&)
470 {}
471
472public:
473 using ScalarArray = Scalar[numPhases];
474
475 static void update(ScalarArray&,
476 const IntensiveQuantities&,
477 const IntensiveQuantities&)
478 {}
479
484 Scalar dispersivity() const
485 {
486 throw std::logic_error("The method dispersivity() does not "
487 "make sense if dispersion is disabled.");
488 }
489
497 Scalar normVelocityAvg(unsigned) const
498 {
499 throw std::logic_error("The method normVelocityAvg() "
500 "does not make sense if dispersion is disabled.");
501 }
502};
503
507template <class TypeTag>
508class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/true>
509{
513
514 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
515
516public:
517 using ScalarArray = std::array<Scalar, numPhases>;
518 static void update(ScalarArray& normVelocityAvg,
519 const IntensiveQuantities& intQuantsInside,
520 const IntensiveQuantities& intQuantsOutside)
521 {
522 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
523 if (!FluidSystem::phaseIsActive(phaseIdx)) {
524 continue;
525 }
526 // no dispersion in water for blackoil models unless water can contain dissolved gas
527 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
528 continue;
529 }
530 // adding dispersion in the gas phase leads to
531 // convergence issues and unphysical results.
532 // we disable dispersion in the gas phase for now
533 if (FluidSystem::gasPhaseIdx == phaseIdx) {
534 continue;
535 }
536 // use the arithmetic average for the effective
537 // velocity coefficients at the face's integration point.
538 normVelocityAvg[phaseIdx] =
539 0.5 * (intQuantsInside.normVelocityCell(phaseIdx) +
540 intQuantsOutside.normVelocityCell(phaseIdx));
541 Valgrind::CheckDefined(normVelocityAvg[phaseIdx]);
542 }
543 }
544
545protected:
546 template <class Context, class FluidState>
547 void updateBoundary_(const Context&,
548 unsigned,
549 unsigned,
550 const FluidState&)
551 { throw std::runtime_error("Not implemented: Dispersion across boundary not implemented for blackoil"); }
552
553public:
560 Scalar dispersivity() const
561 { return dispersivity_; }
562
570 Scalar normVelocityAvg(unsigned phaseIdx) const
571 { return normVelocityAvg_[phaseIdx]; }
572
573 const auto& normVelocityAvg() const
574 { return normVelocityAvg_; }
575
576private:
577 Scalar dispersivity_;
578 ScalarArray normVelocityAvg_;
579};
580
581} // namespace Opm
582
583#endif
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