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{
66
67 enum { numPhases = FluidSystem::numPhases };
68
69public:
71
72#if HAVE_ECL_INPUT
73 static void initFromState(const EclipseState&)
74 {
75 }
76#endif
77
82 template <class Context>
83 static void addDispersiveFlux(RateVector&,
84 const Context&,
85 unsigned,
86 unsigned)
87 {}
88
89 template<class IntensiveQuantities, class Scalar>
90 static void addDispersiveFlux(RateVector&,
91 const IntensiveQuantities&,
92 const IntensiveQuantities&,
93 const Evaluation&,
94 const Scalar&)
95 {}
96};
97
101template <class TypeTag>
102class BlackOilDispersionModule<TypeTag, /*enableDispersion=*/true>
103{
115
116 enum { numPhases = FluidSystem::numPhases };
117 enum { numComponents = FluidSystem::numComponents };
118 enum { conti0EqIdx = Indices::conti0EqIdx };
119 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
120 enum { enableMICP = getPropValue<TypeTag, Properties::EnableMICP>() };
121
122 static constexpr unsigned contiMicrobialEqIdx = Indices::contiMicrobialEqIdx;
123 static constexpr unsigned contiOxygenEqIdx = Indices::contiOxygenEqIdx;
124 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
125 static constexpr unsigned contiUreaEqIdx = Indices::contiUreaEqIdx;
126
127 using Toolbox = MathToolbox<Evaluation>;
128
129public:
131#if HAVE_ECL_INPUT
132 static void initFromState(const EclipseState& eclState)
133 {
134 if (!eclState.getSimulationConfig().rock_config().dispersion()) {
135 return;
136 }
137
138 if (eclState.getSimulationConfig().hasVAPWAT() || eclState.getSimulationConfig().hasVAPOIL()) {
139 OpmLog::warning("Dispersion is activated in combination with VAPWAT/VAPOIL. \n"
140 "Water/oil is still allowed to vaporize, but dispersion in the "
141 "gas phase is ignored.");
142 }
143 }
144#endif
149 template <class Context>
150 static void addDispersiveFlux(RateVector& flux, const Context& context,
151 unsigned spaceIdx, unsigned timeIdx)
152 {
153 // Only work if dispersion is enabled by DISPERC in the deck
154 if (!context.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
155 return;
156 }
157 const auto& extQuants = context.extensiveQuantities(spaceIdx, timeIdx);
158 const auto& inIq = context.intensiveQuantities(extQuants.interiorIndex(), timeIdx);
159 const auto& exIq = context.intensiveQuantities(extQuants.exteriorIndex(), timeIdx);
160 const auto& dispersivity = extQuants.dispersivity();
161 const auto& normVelocityAvg = extQuants.normVelocityAvg();
162 addDispersiveFlux(flux, inIq, exIq, dispersivity, normVelocityAvg);
163 }
164
185 template<class IntensiveQuantities, class Scalar>
186 static void addDispersiveFlux(RateVector& flux,
187 const IntensiveQuantities& inIq,
188 const IntensiveQuantities& exIq,
189 const Evaluation& dispersivity,
190 const Scalar& normVelocityAvg)
191 {
192 const auto& inFs = inIq.fluidState();
193 const auto& exFs = exIq.fluidState();
194 Evaluation diffR = 0.0;
195 if constexpr(enableMICP) {
196 // The dispersion coefficients are given for mass concentrations
197 Evaluation bAvg = (inFs.invB(waterPhaseIdx) + Toolbox::value(exFs.invB(waterPhaseIdx))) / 2;
198 diffR = inIq.microbialConcentration()- Toolbox::value(exIq.microbialConcentration());
199 flux[contiMicrobialEqIdx] +=
200 bAvg
201 * normVelocityAvg[waterPhaseIdx]
202 * dispersivity
203 * diffR;
204 diffR = inIq.oxygenConcentration()- Toolbox::value(exIq.oxygenConcentration());
205 flux[contiOxygenEqIdx] +=
206 bAvg
207 * normVelocityAvg[waterPhaseIdx]
208 * dispersivity
209 * diffR;
210 diffR = inIq.ureaConcentration()- Toolbox::value(exIq.ureaConcentration());
211 flux[contiUreaEqIdx] +=
212 bAvg
213 * normVelocityAvg[waterPhaseIdx]
214 * dispersivity
215 * diffR;
216 return;
217 }
218
219 unsigned pvtRegionIndex = inFs.pvtRegionIndex();
220 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
221 if (!FluidSystem::phaseIsActive(phaseIdx)) {
222 continue;
223 }
224
225 // no dispersion in water for blackoil models unless water can contain dissolved gas
226 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
227 continue;
228 }
229
230 // adding dispersion in the gas phase leads to
231 // convergence issues and unphysical results.
232 // we disable dispersion in the gas phase for now
233 if (FluidSystem::gasPhaseIdx == phaseIdx) {
234 continue;
235 }
236
237 // arithmetic mean of the phase's b factor
238 Evaluation bAvg = inFs.invB(phaseIdx);
239 bAvg += Toolbox::value(exFs.invB(phaseIdx));
240 bAvg /= 2;
241
242 Evaluation convFactor = 1.0;
243 if (FluidSystem::enableDissolvedGas() && FluidSystem::phaseIsActive(FluidSystem::gasPhaseIdx) && phaseIdx == FluidSystem::oilPhaseIdx) {
244 Evaluation rsAvg = (inFs.Rs() + Toolbox::value(exFs.Rs())) / 2;
245 convFactor = 1.0 / (toMassFractionGasOil(pvtRegionIndex) + rsAvg);
246 diffR = inFs.Rs() - Toolbox::value(exFs.Rs());
247 }
248 if (FluidSystem::enableVaporizedOil() && FluidSystem::phaseIsActive(FluidSystem::oilPhaseIdx) && phaseIdx == FluidSystem::gasPhaseIdx) {
249 Evaluation rvAvg = (inFs.Rv() + Toolbox::value(exFs.Rv())) / 2;
250 convFactor = toMassFractionGasOil(pvtRegionIndex) / (1.0 + rvAvg*toMassFractionGasOil(pvtRegionIndex));
251 diffR = inFs.Rv() - Toolbox::value(exFs.Rv());
252 }
253 if (FluidSystem::enableDissolvedGasInWater() && phaseIdx == FluidSystem::waterPhaseIdx) {
254 Evaluation rsAvg = (inFs.Rsw() + Toolbox::value(exFs.Rsw())) / 2;
255 convFactor = 1.0 / (toMassFractionGasWater(pvtRegionIndex) + rsAvg);
256 diffR = inFs.Rsw() - Toolbox::value(exFs.Rsw());
257 }
258 if (FluidSystem::enableVaporizedWater() && phaseIdx == FluidSystem::gasPhaseIdx) {
259 Evaluation rvAvg = (inFs.Rvw() + Toolbox::value(exFs.Rvw())) / 2;
260 convFactor = toMassFractionGasWater(pvtRegionIndex)/ (1.0 + rvAvg*toMassFractionGasWater(pvtRegionIndex));
261 diffR = inFs.Rvw() - Toolbox::value(exFs.Rvw());
262 }
263
264 // mass flux of solvent component
265 unsigned solventCompIdx = FluidSystem::solventComponentIndex(phaseIdx);
266 unsigned activeSolventCompIdx = Indices::canonicalToActiveComponentIndex(solventCompIdx);
267 flux[conti0EqIdx + activeSolventCompIdx] +=
268 - bAvg
269 * normVelocityAvg[phaseIdx]
270 * convFactor
271 * dispersivity
272 * diffR;
273
274 // mass flux of solute component
275 unsigned soluteCompIdx = FluidSystem::soluteComponentIndex(phaseIdx);
276 unsigned activeSoluteCompIdx = Indices::canonicalToActiveComponentIndex(soluteCompIdx);
277 flux[conti0EqIdx + activeSoluteCompIdx] +=
278 bAvg
279 * normVelocityAvg[phaseIdx]
280 * convFactor
281 * dispersivity
282 * diffR;
283 }
284 }
285
286private:
287
288 static Scalar toMassFractionGasOil (unsigned regionIdx) {
289 Scalar rhoO = FluidSystem::referenceDensity(FluidSystem::oilPhaseIdx, regionIdx);
290 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
291 return rhoO / rhoG;
292 }
293 static Scalar toMassFractionGasWater (unsigned regionIdx) {
294 Scalar rhoW = FluidSystem::referenceDensity(FluidSystem::waterPhaseIdx, regionIdx);
295 Scalar rhoG = FluidSystem::referenceDensity(FluidSystem::gasPhaseIdx, regionIdx);
296 return rhoW / rhoG;
297 }
298};
299
307template <class TypeTag, bool enableDispersion>
309
313template <class TypeTag>
314class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/false>
315{
319
320public:
324 Scalar normVelocityCell(unsigned, unsigned) const
325 {
326 throw std::logic_error("Method normVelocityCell() "
327 "does not make sense if dispersion is disabled");
328 }
329
330protected:
335 template<class ElementContext>
336 void update_(ElementContext&,
337 unsigned,
338 unsigned)
339 { }
340};
341
345template <class TypeTag>
346class BlackOilDispersionIntensiveQuantities<TypeTag, /*enableDispersion=*/true>
347{
354 enum { numPhases = FluidSystem::numPhases };
355 enum { numComponents = FluidSystem::numComponents };
356 enum { oilPhaseIdx = FluidSystem::oilPhaseIdx };
357 enum { gasPhaseIdx = FluidSystem::gasPhaseIdx };
358 enum { waterPhaseIdx = FluidSystem::waterPhaseIdx };
359 enum { gasCompIdx = FluidSystem::gasCompIdx };
360 enum { oilCompIdx = FluidSystem::oilCompIdx };
361 enum { waterCompIdx = FluidSystem::waterCompIdx };
362 enum { conti0EqIdx = Indices::conti0EqIdx };
363 enum { enableDispersion = getPropValue<TypeTag, Properties::EnableDispersion>() };
364
365public:
369 Scalar normVelocityCell(unsigned phaseIdx) const
370 {
371
372 return normVelocityCell_[phaseIdx];
373 }
374
375protected:
385 template<class ElementContext>
386 void update_(const ElementContext& elemCtx, unsigned dofIdx, unsigned timeIdx)
387 {
388 // Only work if dispersion is enabled by DISPERC in the deck
389 if (!elemCtx.simulator().vanguard().eclState().getSimulationConfig().rock_config().dispersion()) {
390 return;
391 }
392 const auto& problem = elemCtx.simulator().problem();
393 if (problem.model().linearizer().getVelocityInfo().empty()) {
394 return;
395 }
396 const std::array<int, 3> phaseIdxs = { gasPhaseIdx, oilPhaseIdx, waterPhaseIdx };
397 const std::array<int, 3> compIdxs = { gasCompIdx, oilCompIdx, waterCompIdx };
398 const auto& velocityInf = problem.model().linearizer().getVelocityInfo();
399 unsigned globalDofIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
400 auto velocityInfos = velocityInf[globalDofIdx];
401 for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
402 normVelocityCell_[i] = 0;
403 }
404 for (const auto& velocityInfo : velocityInfos) {
405 for (unsigned i = 0; i < phaseIdxs.size(); ++i) {
406 if (FluidSystem::phaseIsActive(phaseIdxs[i])) {
407 normVelocityCell_[phaseIdxs[i]] = max( normVelocityCell_[phaseIdxs[i]],
408 std::abs( velocityInfo.velocity[conti0EqIdx
409 + Indices::canonicalToActiveComponentIndex(compIdxs[i])] ));
410 }
411 }
412 }
413 }
414
415private:
416 Scalar normVelocityCell_[numPhases];
417};
418
425template <class TypeTag, bool enableDispersion>
426class BlackOilDispersionExtensiveQuantities;
427
431template <class TypeTag>
432class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/false>
433{
439
440 enum { numPhases = FluidSystem::numPhases };
441
442protected:
447 void update_(const ElementContext&,
448 unsigned,
449 unsigned)
450 {}
451
452 template <class Context, class FluidState>
453 void updateBoundary_(const Context&,
454 unsigned,
455 unsigned,
456 const FluidState&)
457 {}
458
459public:
460 using ScalarArray = Scalar[numPhases];
461
462 static void update(ScalarArray&,
463 const IntensiveQuantities&,
464 const IntensiveQuantities&)
465 {}
466
471 const Scalar& dispersivity() const
472 {
473 throw std::logic_error("The method dispersivity() does not "
474 "make sense if dispersion is disabled.");
475 }
476
484 const Scalar& normVelocityAvg(unsigned) const
485 {
486 throw std::logic_error("The method normVelocityAvg() "
487 "does not make sense if dispersion is disabled.");
488 }
489
490};
491
495template <class TypeTag>
496class BlackOilDispersionExtensiveQuantities<TypeTag, /*enableDispersion=*/true>
497{
503 using Toolbox = MathToolbox<Evaluation>;
505
506 enum { dimWorld = GridView::dimensionworld };
508 enum { numComponents = getPropValue<TypeTag, Properties::NumComponents>() };
509
510 using DimVector = Dune::FieldVector<Scalar, dimWorld>;
511 using DimEvalVector = Dune::FieldVector<Evaluation, dimWorld>;
512public:
513 using ScalarArray = Scalar[numPhases];
514 static void update(ScalarArray& normVelocityAvg,
515 const IntensiveQuantities& intQuantsInside,
516 const IntensiveQuantities& intQuantsOutside)
517 {
518 for (unsigned phaseIdx = 0; phaseIdx < numPhases; ++phaseIdx) {
519 if (!FluidSystem::phaseIsActive(phaseIdx)) {
520 continue;
521 }
522 // no dispersion in water for blackoil models unless water can contain dissolved gas
523 if (!FluidSystem::enableDissolvedGasInWater() && FluidSystem::waterPhaseIdx == phaseIdx) {
524 continue;
525 }
526 // adding dispersion in the gas phase leads to
527 // convergence issues and unphysical results.
528 // we disable dispersion in the gas phase for now
529 if (FluidSystem::gasPhaseIdx == phaseIdx) {
530 continue;
531 }
532 // use the arithmetic average for the effective
533 // velocity coefficients at the face's integration point.
534 normVelocityAvg[phaseIdx] = 0.5 *
535 ( intQuantsInside.normVelocityCell(phaseIdx) +
536 intQuantsOutside.normVelocityCell(phaseIdx) );
537 Valgrind::CheckDefined(normVelocityAvg[phaseIdx]);
538 }
539 }
540protected:
541 template <class Context, class FluidState>
542 void updateBoundary_(const Context&,
543 unsigned,
544 unsigned,
545 const FluidState&)
546 {
547 throw std::runtime_error("Not implemented: Dispersion across boundary not implemented for blackoil");
548 }
549
550public:
557 const Scalar& dispersivity() const
558 { return dispersivity_; }
559
567 const Scalar& normVelocityAvg(unsigned phaseIdx) const
568 { return normVelocityAvg_[phaseIdx]; }
569
570 const auto& normVelocityAvg() const{
571 return normVelocityAvg_;
572 }
573
574private:
575 Scalar dispersivity_;
576 ScalarArray normVelocityAvg_;
577};
578
579} // namespace Opm
580
581#endif
Definition blackoildispersionmodule.hh:433
void update_(const ElementContext &, unsigned, unsigned)
Update the quantities required to calculate the dispersive fluxes.
Definition blackoildispersionmodule.hh:447
const Scalar & dispersivity() const
The dispersivity the face.
Definition blackoildispersionmodule.hh:471
const Scalar & normVelocityAvg(unsigned) const
The effective filter velocity coefficient in a fluid phase at the face's integration point.
Definition blackoildispersionmodule.hh:484
Provides the quantities required to calculate dispersive mass fluxes.
Definition blackoildispersionmodule.hh:497
const Scalar & dispersivity() const
The dispersivity of the face.
Definition blackoildispersionmodule.hh:557
const Scalar & normVelocityAvg(unsigned phaseIdx) const
The effective velocity coefficient in a fluid phase at the face's integration point.
Definition blackoildispersionmodule.hh:567
Provides the quantities required to calculate dispersive mass fluxes.
Definition blackoildispersionmodule.hh:53
Scalar normVelocityCell(unsigned, unsigned) const
Returns the max.
Definition blackoildispersionmodule.hh:324
void update_(ElementContext &, unsigned, unsigned)
Update the quantities required to calculate dispersive fluxes.
Definition blackoildispersionmodule.hh:336
void update_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Update the quantities required to calculate dispersive mass fluxes.
Definition blackoildispersionmodule.hh:386
Scalar normVelocityCell(unsigned phaseIdx) const
Returns the max.
Definition blackoildispersionmodule.hh:369
Provides the volumetric quantities required for the calculation of dispersive fluxes.
Definition blackoildispersionmodule.hh:308
static void addDispersiveFlux(RateVector &, const Context &, unsigned, unsigned)
Adds the dispersive flux to the flux vector over a flux integration point.
Definition blackoildispersionmodule.hh:83
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:186
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:150
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