blackoilextbomodules.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*/
32#ifndef EWOMS_BLACK_OIL_EXTBO_MODULE_HH
33#define EWOMS_BLACK_OIL_EXTBO_MODULE_HH
34
35#include "blackoilproperties.hh"
36
38
39//#include <opm/models/io/vtkBlackOilExtboModule.hh> //TODO: Missing ...
40
41#if HAVE_ECL_INPUT
42#include <opm/input/eclipse/EclipseState/EclipseState.hpp>
43#include <opm/input/eclipse/EclipseState/Tables/SsfnTable.hpp>
44#include <opm/input/eclipse/EclipseState/Tables/Sof2Table.hpp>
45#include <opm/input/eclipse/EclipseState/Tables/MsfnTable.hpp>
46#include <opm/input/eclipse/EclipseState/Tables/PmiscTable.hpp>
47#include <opm/input/eclipse/EclipseState/Tables/MiscTable.hpp>
48#include <opm/input/eclipse/EclipseState/Tables/SorwmisTable.hpp>
49#include <opm/input/eclipse/EclipseState/Tables/SgcwmisTable.hpp>
50#include <opm/input/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
51#endif
52
53#include <dune/common/fvector.hh>
54
55#include <cmath>
56#include <stdexcept>
57#include <string>
58#include <vector>
59
60namespace Opm {
61
67template <class TypeTag, bool enableExtboV = getPropValue<TypeTag, Properties::EnableExtbo>()>
69{
82
83 using Toolbox = MathToolbox<Evaluation>;
84
85 using TabulatedFunction = typename BlackOilExtboParams<Scalar>::TabulatedFunction;
86 using Tabulated2DFunction = typename BlackOilExtboParams<Scalar>::Tabulated2DFunction;
87
88 static constexpr unsigned zFractionIdx = Indices::zFractionIdx;
89 static constexpr unsigned contiZfracEqIdx = Indices::contiZfracEqIdx;
90 static constexpr unsigned enableExtbo = enableExtboV;
91 static constexpr unsigned numEq = getPropValue<TypeTag, Properties::NumEq>();
92 static constexpr unsigned numPhases = FluidSystem::numPhases;
93 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
94 static constexpr unsigned oilPhaseIdx = FluidSystem::oilPhaseIdx;
95 static constexpr unsigned waterPhaseIdx = FluidSystem::waterPhaseIdx;
96 static constexpr bool blackoilConserveSurfaceVolume = getPropValue<TypeTag, Properties::BlackoilConserveSurfaceVolume>();
97
98public:
99#if HAVE_ECL_INPUT
103 static void initFromState(const EclipseState& eclState)
104 {
105 // some sanity checks: if extended BO is enabled, the PVTSOL keyword must be
106 // present, if extended BO is disabled the keyword must not be present.
107 if (enableExtbo && !eclState.runspec().phases().active(Phase::ZFRACTION))
108 throw std::runtime_error("Extended black oil treatment requested at compile "
109 "time, but the deck does not contain the PVTSOL keyword");
110 else if (!enableExtbo && eclState.runspec().phases().active(Phase::ZFRACTION))
111 throw std::runtime_error("Extended black oil treatment disabled at compile time, but the deck "
112 "contains the PVTSOL keyword");
113
114 if (!eclState.runspec().phases().active(Phase::ZFRACTION))
115 return; // solvent treatment is supposed to be disabled
116
117 // pvt properties from kw PVTSOL:
118
119 const auto& tableManager = eclState.getTableManager();
120 const auto& pvtsolTables = tableManager.getPvtsolTables();
121
122 size_t numPvtRegions = pvtsolTables.size();
123
124 params_.BO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
125 params_.BG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
126 params_.RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
127 params_.RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
128 params_.X_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
129 params_.Y_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
130 params_.VISCO_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
131 params_.VISCG_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
132
133 params_.PBUB_RS_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
134 params_.PBUB_RV_.resize(numPvtRegions, Tabulated2DFunction{Tabulated2DFunction::InterpolationPolicy::LeftExtreme});
135
136 params_.zLim_.resize(numPvtRegions);
137
138 const bool extractCmpFromPvt = true; //<false>: Default values used in [*]
139 params_.oilCmp_.resize(numPvtRegions);
140 params_.gasCmp_.resize(numPvtRegions);
141
142 for (unsigned regionIdx = 0; regionIdx < numPvtRegions; ++ regionIdx) {
143 const auto& pvtsolTable = pvtsolTables[regionIdx];
144
145 const auto& saturatedTable = pvtsolTable.getSaturatedTable();
146 assert(saturatedTable.numRows() > 1);
147
148 std::vector<Scalar> oilCmp(saturatedTable.numRows(), -4.0e-9); //Default values used in [*]
149 std::vector<Scalar> gasCmp(saturatedTable.numRows(), -0.08); //-------------"-------------
150 params_.zLim_[regionIdx] = 0.7; //-------------"-------------
151 std::vector<Scalar> zArg(saturatedTable.numRows(), 0.0);
152
153 for (unsigned outerIdx = 0; outerIdx < saturatedTable.numRows(); ++ outerIdx) {
154 Scalar ZCO2 = saturatedTable.get("ZCO2", outerIdx);
155
156 zArg[outerIdx] = ZCO2;
157
158 params_.BO_[regionIdx].appendXPos(ZCO2);
159 params_.BG_[regionIdx].appendXPos(ZCO2);
160
161 params_.RS_[regionIdx].appendXPos(ZCO2);
162 params_.RV_[regionIdx].appendXPos(ZCO2);
163
164 params_.X_[regionIdx].appendXPos(ZCO2);
165 params_.Y_[regionIdx].appendXPos(ZCO2);
166
167 params_.VISCO_[regionIdx].appendXPos(ZCO2);
168 params_.VISCG_[regionIdx].appendXPos(ZCO2);
169
170 params_.PBUB_RS_[regionIdx].appendXPos(ZCO2);
171 params_.PBUB_RV_[regionIdx].appendXPos(ZCO2);
172
173 const auto& underSaturatedTable = pvtsolTable.getUnderSaturatedTable(outerIdx);
174 size_t numRows = underSaturatedTable.numRows();
175
176 Scalar bo0=0.0;
177 Scalar po0=0.0;
178 for (unsigned innerIdx = 0; innerIdx < numRows; ++ innerIdx) {
179 Scalar po = underSaturatedTable.get("P", innerIdx);
180 Scalar bo = underSaturatedTable.get("B_O", innerIdx);
181 Scalar bg = underSaturatedTable.get("B_G", innerIdx);
182 Scalar rs = underSaturatedTable.get("RS", innerIdx)+innerIdx*1.0e-10;
183 Scalar rv = underSaturatedTable.get("RV", innerIdx)+innerIdx*1.0e-10;
184 Scalar xv = underSaturatedTable.get("XVOL", innerIdx);
185 Scalar yv = underSaturatedTable.get("YVOL", innerIdx);
186 Scalar mo = underSaturatedTable.get("MU_O", innerIdx);
187 Scalar mg = underSaturatedTable.get("MU_G", innerIdx);
188
189 if (bo0 > bo) { // This is undersaturated oil-phase for ZCO2 <= zLim ...
190 // Here we assume tabulated bo to decay beyond boiling point
191 if (extractCmpFromPvt) {
192 Scalar cmpFactor = (bo-bo0)/(po-po0);
193 oilCmp[outerIdx] = cmpFactor;
194 params_.zLim_[regionIdx] = ZCO2;
195 //std::cout << "### cmpFactorOil: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl;
196 }
197 break;
198 } else if (bo0 == bo) { // This is undersaturated gas-phase for ZCO2 > zLim ...
199 // Here we assume tabulated bo to be constant extrapolated beyond dew point
200 if (innerIdx+1 < numRows && ZCO2<1.0 && extractCmpFromPvt) {
201 Scalar rvNxt = underSaturatedTable.get("RV", innerIdx+1)+innerIdx*1.0e-10;
202 Scalar bgNxt = underSaturatedTable.get("B_G", innerIdx+1);
203 Scalar cmpFactor = (bgNxt-bg)/(rvNxt-rv);
204 gasCmp[outerIdx] = cmpFactor;
205 //std::cout << "### cmpFactorGas: " << cmpFactor << " zLim: " << zLim_[regionIdx] << std::endl;
206 }
207
208 params_.BO_[regionIdx].appendSamplePoint(outerIdx,po,bo);
209 params_.BG_[regionIdx].appendSamplePoint(outerIdx,po,bg);
210 params_.RS_[regionIdx].appendSamplePoint(outerIdx,po,rs);
211 params_.RV_[regionIdx].appendSamplePoint(outerIdx,po,rv);
212 params_.X_[regionIdx].appendSamplePoint(outerIdx,po,xv);
213 params_.Y_[regionIdx].appendSamplePoint(outerIdx,po,yv);
214 params_.VISCO_[regionIdx].appendSamplePoint(outerIdx,po,mo);
215 params_.VISCG_[regionIdx].appendSamplePoint(outerIdx,po,mg);
216 break;
217 }
218
219 bo0=bo;
220 po0=po;
221
222 params_.BO_[regionIdx].appendSamplePoint(outerIdx,po,bo);
223 params_.BG_[regionIdx].appendSamplePoint(outerIdx,po,bg);
224
225 params_.RS_[regionIdx].appendSamplePoint(outerIdx,po,rs);
226 params_.RV_[regionIdx].appendSamplePoint(outerIdx,po,rv);
227
228 params_.X_[regionIdx].appendSamplePoint(outerIdx,po,xv);
229 params_.Y_[regionIdx].appendSamplePoint(outerIdx,po,yv);
230
231 params_.VISCO_[regionIdx].appendSamplePoint(outerIdx,po,mo);
232 params_.VISCG_[regionIdx].appendSamplePoint(outerIdx,po,mg);
233
234 // rs,rv -> pressure
235 params_.PBUB_RS_[regionIdx].appendSamplePoint(outerIdx, rs, po);
236 params_.PBUB_RV_[regionIdx].appendSamplePoint(outerIdx, rv, po);
237
238 }
239 }
240 params_.oilCmp_[regionIdx].setXYContainers(zArg, oilCmp, /*sortInput=*/false);
241 params_.gasCmp_[regionIdx].setXYContainers(zArg, gasCmp, /*sortInput=*/false);
242 }
243
244 // Reference density for pure z-component taken from kw SDENSITY
245 const auto& sdensityTables = eclState.getTableManager().getSolventDensityTables();
246 if (sdensityTables.size() == numPvtRegions) {
247 params_.zReferenceDensity_.resize(numPvtRegions);
248 for (unsigned regionIdx = 0; regionIdx < numPvtRegions; ++ regionIdx) {
249 Scalar rhoRefS = sdensityTables[regionIdx].getSolventDensityColumn().front();
250 params_.zReferenceDensity_[regionIdx]=rhoRefS;
251 }
252 }
253 else
254 throw std::runtime_error("Extbo: kw SDENSITY is missing or not aligned with NTPVT\n");
255 }
256#endif
257
261 static void registerParameters()
262 {
263 }
264
268 static void registerOutputModules(Model&,
269 Simulator&)
270 {
271 }
272
273 static bool primaryVarApplies(unsigned pvIdx)
274 {
275 if constexpr (enableExtbo)
276 return pvIdx == zFractionIdx;
277 else
278 return false;
279 }
280
281 static std::string primaryVarName([[maybe_unused]] unsigned pvIdx)
282 {
283 assert(primaryVarApplies(pvIdx));
284
285 return "z_fraction";
286 }
287
288 static Scalar primaryVarWeight([[maybe_unused]] unsigned pvIdx)
289 {
290 assert(primaryVarApplies(pvIdx));
291
292 // TODO: it may be beneficial to chose this differently.
293 return static_cast<Scalar>(1.0);
294 }
295
296 static bool eqApplies(unsigned eqIdx)
297 {
298 if constexpr (enableExtbo)
299 return eqIdx == contiZfracEqIdx;
300 else
301 return false;
302 }
303
304 static std::string eqName([[maybe_unused]] unsigned eqIdx)
305 {
306 assert(eqApplies(eqIdx));
307
308 return "conti^solvent";
309 }
310
311 static Scalar eqWeight([[maybe_unused]] unsigned eqIdx)
312 {
313 assert(eqApplies(eqIdx));
314
315 // TODO: it may be beneficial to chose this differently.
316 return static_cast<Scalar>(1.0);
317 }
318
319 template <class LhsEval>
320 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
321 const IntensiveQuantities& intQuants)
322 {
323 if constexpr (enableExtbo) {
324 if constexpr (blackoilConserveSurfaceVolume) {
325 storage[contiZfracEqIdx] =
326 Toolbox::template decay<LhsEval>(intQuants.porosity())
327 * Toolbox::template decay<LhsEval>(intQuants.yVolume())
328 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(gasPhaseIdx))
329 * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(gasPhaseIdx));
330 if (FluidSystem::enableDissolvedGas()) { // account for dissolved z in oil phase
331 storage[contiZfracEqIdx] +=
332 Toolbox::template decay<LhsEval>(intQuants.porosity())
333 * Toolbox::template decay<LhsEval>(intQuants.xVolume())
334 * Toolbox::template decay<LhsEval>(intQuants.fluidState().Rs())
335 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(oilPhaseIdx))
336 * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(oilPhaseIdx));
337 }
338 // Reg. terms: Preliminary attempt to avoid singular behaviour when solvent is invading a pure water
339 // region. Results seems insensitive to the weighting factor.
340 // TODO: Further investigations ...
341 const Scalar regWghtFactor = 1.0e-6;
342 storage[contiZfracEqIdx] += regWghtFactor*(1.0-Toolbox::template decay<LhsEval>(intQuants.zFraction()))
343 + regWghtFactor*Toolbox::template decay<LhsEval>(intQuants.porosity())
344 * Toolbox::template decay<LhsEval>(intQuants.fluidState().saturation(gasPhaseIdx))
345 * Toolbox::template decay<LhsEval>(intQuants.fluidState().invB(gasPhaseIdx));
346 storage[contiZfracEqIdx-1] += regWghtFactor*Toolbox::template decay<LhsEval>(intQuants.zFraction());
347 }
348 else {
349 throw std::runtime_error("Only component conservation in terms of surface volumes is implemented. ");
350 }
351 }
352 }
353
354 static void computeFlux([[maybe_unused]] RateVector& flux,
355 [[maybe_unused]] const ElementContext& elemCtx,
356 [[maybe_unused]] unsigned scvfIdx,
357 [[maybe_unused]] unsigned timeIdx)
358
359 {
360 if constexpr (enableExtbo) {
361 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
362
363 if constexpr (blackoilConserveSurfaceVolume) {
364 unsigned inIdx = extQuants.interiorIndex();
365
366 unsigned upIdxGas = static_cast<unsigned>(extQuants.upstreamIndex(gasPhaseIdx));
367 const auto& upGas = elemCtx.intensiveQuantities(upIdxGas, timeIdx);
368 const auto& fsGas = upGas.fluidState();
369 if (upIdxGas == inIdx) {
370 flux[contiZfracEqIdx] =
371 extQuants.volumeFlux(gasPhaseIdx)
372 * (upGas.yVolume())
373 * fsGas.invB(gasPhaseIdx);
374 }
375 else {
376 flux[contiZfracEqIdx] =
377 extQuants.volumeFlux(gasPhaseIdx)
378 * (decay<Scalar>(upGas.yVolume()))
379 * decay<Scalar>(fsGas.invB(gasPhaseIdx));
380 }
381 if (FluidSystem::enableDissolvedGas()) { // account for dissolved z in oil phase
382 unsigned upIdxOil = static_cast<unsigned>(extQuants.upstreamIndex(oilPhaseIdx));
383 const auto& upOil = elemCtx.intensiveQuantities(upIdxOil, timeIdx);
384 const auto& fsOil = upOil.fluidState();
385 if (upIdxOil == inIdx) {
386 flux[contiZfracEqIdx] +=
387 extQuants.volumeFlux(oilPhaseIdx)
388 * upOil.xVolume()
389 * fsOil.Rs()
390 * fsOil.invB(oilPhaseIdx);
391 }
392 else {
393 flux[contiZfracEqIdx] +=
394 extQuants.volumeFlux(oilPhaseIdx)
395 * decay<Scalar>(upOil.xVolume())
396 * decay<Scalar>(fsOil.Rs())
397 * decay<Scalar>(fsOil.invB(oilPhaseIdx));
398 }
399 }
400 }
401 else {
402 throw std::runtime_error("Only component conservation in terms of surface volumes is implemented. ");
403 }
404 }
405 }
406
410 static void assignPrimaryVars(PrimaryVariables& priVars,
411 Scalar zFraction)
412 {
413 if constexpr (enableExtbo)
414 priVars[zFractionIdx] = zFraction;
415 }
416
420 static void updatePrimaryVars(PrimaryVariables& newPv,
421 const PrimaryVariables& oldPv,
422 const EqVector& delta)
423 {
424 if constexpr (enableExtbo)
425 // do a plain unchopped Newton update
426 newPv[zFractionIdx] = oldPv[zFractionIdx] - delta[zFractionIdx];
427 }
428
432 static Scalar computeUpdateError(const PrimaryVariables&,
433 const EqVector&)
434 {
435 // do not consider consider the cange of solvent primary variables for
436 // convergence
437 // TODO: maybe this should be changed
438 return static_cast<Scalar>(0.0);
439 }
440
444 static Scalar computeResidualError(const EqVector& resid)
445 {
446 // do not weight the residual of solvents when it comes to convergence
447 return std::abs(Toolbox::scalarValue(resid[contiZfracEqIdx]));
448 }
449
450 template <class DofEntity>
451 static void serializeEntity(const Model& model, std::ostream& outstream, const DofEntity& dof)
452 {
453 if constexpr (enableExtbo) {
454 unsigned dofIdx = model.dofMapper().index(dof);
455
456 const PrimaryVariables& priVars = model.solution(/*timeIdx=*/0)[dofIdx];
457 outstream << priVars[zFractionIdx];
458 }
459 }
460
461 template <class DofEntity>
462 static void deserializeEntity(Model& model, std::istream& instream, const DofEntity& dof)
463 {
464 if constexpr (enableExtbo) {
465 unsigned dofIdx = model.dofMapper().index(dof);
466
467 PrimaryVariables& priVars0 = model.solution(/*timeIdx=*/0)[dofIdx];
468 PrimaryVariables& priVars1 = model.solution(/*timeIdx=*/1)[dofIdx];
469
470 instream >> priVars0[zFractionIdx];
471
472 // set the primary variables for the beginning of the current time step.
473 priVars1 = priVars0[zFractionIdx];
474 }
475 }
476
477 template <typename Value>
478 static Value xVolume(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
479 return params_.X_[pvtRegionIdx].eval(z, pressure, true);
480 }
481
482 template <typename Value>
483 static Value yVolume(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
484 return params_.Y_[pvtRegionIdx].eval(z, pressure, true);
485 }
486
487 template <typename Value>
488 static Value pbubRs(unsigned pvtRegionIdx, const Value& z, const Value& rs) {
489 return params_.PBUB_RS_[pvtRegionIdx].eval(z, rs, true);
490 }
491
492 template <typename Value>
493 static Value pbubRv(unsigned pvtRegionIdx, const Value& z, const Value& rv) {
494 return params_.PBUB_RV_[pvtRegionIdx].eval(z, rv, true);
495 }
496
497 template <typename Value>
498 static Value oilViscosity(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
499 return params_.VISCO_[pvtRegionIdx].eval(z, pressure, true);
500 }
501
502 template <typename Value>
503 static Value gasViscosity(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
504 return params_.VISCG_[pvtRegionIdx].eval(z, pressure, true);
505 }
506
507 template <typename Value>
508 static Value bo(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
509 return params_.BO_[pvtRegionIdx].eval(z, pressure, true);
510 }
511
512 template <typename Value>
513 static Value bg(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
514 return params_.BG_[pvtRegionIdx].eval(z, pressure, true);
515 }
516
517 template <typename Value>
518 static Value rs(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
519 return params_.RS_[pvtRegionIdx].eval(z, pressure, true);
520 }
521
522 template <typename Value>
523 static Value rv(unsigned pvtRegionIdx, const Value& pressure, const Value& z) {
524 return params_.RV_[pvtRegionIdx].eval(z, pressure, true);
525 }
526
527 static Scalar referenceDensity(unsigned regionIdx) {
528 return params_.zReferenceDensity_[regionIdx];
529 }
530
531 static Scalar zLim(unsigned regionIdx) {
532 return params_.zLim_[regionIdx];
533 }
534
535 template <typename Value>
536 static Value oilCmp(unsigned pvtRegionIdx, const Value& z) {
537 return params_.oilCmp_[pvtRegionIdx].eval(z, true);
538 }
539
540 template <typename Value>
541 static Value gasCmp(unsigned pvtRegionIdx, const Value& z) {
542 return params_.gasCmp_[pvtRegionIdx].eval(z, true);
543 }
544
545private:
546 static BlackOilExtboParams<Scalar> params_;
547};
548
549template <class TypeTag, bool enableExtboV>
550BlackOilExtboParams<typename BlackOilExtboModule<TypeTag, enableExtboV>::Scalar>
551BlackOilExtboModule<TypeTag, enableExtboV>::params_;
552
560template <class TypeTag, bool enableExtboV = getPropValue<TypeTag, Properties::EnableExtbo>()>
562{
564
572
574
575 enum { numPhases = getPropValue<TypeTag, Properties::NumPhases>() };
576 static constexpr int zFractionIdx = Indices::zFractionIdx;
577 static constexpr int oilPhaseIdx = FluidSystem::oilPhaseIdx;
578 static constexpr int gasPhaseIdx = FluidSystem::gasPhaseIdx;
579 static constexpr int waterPhaseIdx = FluidSystem::waterPhaseIdx;
580 static constexpr double cutOff = 1e-12;
581
582
583public:
589 void zFractionUpdate_(const ElementContext& elemCtx,
590 unsigned dofIdx,
591 unsigned timeIdx)
592 {
593 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
594 unsigned pvtRegionIdx = priVars.pvtRegionIndex();
595 auto& fs = asImp_().fluidState_;
596
597 zFraction_ = priVars.makeEvaluation(zFractionIdx, timeIdx);
598
599 oilViscosity_ = ExtboModule::oilViscosity(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
600 gasViscosity_ = ExtboModule::gasViscosity(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
601
602 bo_ = ExtboModule::bo(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
603 bg_ = ExtboModule::bg(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
604
605 bz_ = ExtboModule::bg(pvtRegionIdx, fs.pressure(oilPhaseIdx), Evaluation{0.99});
606
607 if (FluidSystem::enableDissolvedGas())
608 rs_ = ExtboModule::rs(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
609 else
610 rs_ = 0.0;
611
612 if (FluidSystem::enableVaporizedOil())
613 rv_ = ExtboModule::rv(pvtRegionIdx, fs.pressure(gasPhaseIdx), zFraction_);
614 else
615 rv_ = 0.0;
616
617 xVolume_ = ExtboModule::xVolume(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
618 yVolume_ = ExtboModule::yVolume(pvtRegionIdx, fs.pressure(oilPhaseIdx), zFraction_);
619
620 Evaluation pbub = fs.pressure(oilPhaseIdx);
621
622 if (priVars.primaryVarsMeaningWater() == PrimaryVariables::WaterMeaning::Sw) {
623 static const Scalar thresholdWaterFilledCell = 1.0 - 1e-6;
624 Scalar sw = priVars.makeEvaluation(Indices::waterSwitchIdx, timeIdx).value();
625 if (sw >= thresholdWaterFilledCell)
626 rs_ = 0.0; // water only, zero rs_ ...
627 }
628
629 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rs) {
630 rs_ = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
631 const Evaluation zLim = ExtboModule::zLim(pvtRegionIdx);
632 if (zFraction_ > zLim) {
633 pbub = ExtboModule::pbubRs(pvtRegionIdx, zLim, rs_);
634 } else {
635 pbub = ExtboModule::pbubRs(pvtRegionIdx, zFraction_, rs_);
636 }
637 bo_ = ExtboModule::bo(pvtRegionIdx, pbub, zFraction_) + ExtboModule::oilCmp(pvtRegionIdx, zFraction_)*(fs.pressure(oilPhaseIdx)-pbub);
638
639 xVolume_ = ExtboModule::xVolume(pvtRegionIdx, pbub, zFraction_);
640 }
641
642 if (priVars.primaryVarsMeaningGas() == PrimaryVariables::GasMeaning::Rv) {
643 rv_ = priVars.makeEvaluation(Indices::compositionSwitchIdx, timeIdx);
644 Evaluation rvsat = ExtboModule::rv(pvtRegionIdx, pbub, zFraction_);
645 bg_ = ExtboModule::bg(pvtRegionIdx, pbub, zFraction_) + ExtboModule::gasCmp(pvtRegionIdx, zFraction_)*(rv_-rvsat);
646
647 yVolume_ = ExtboModule::yVolume(pvtRegionIdx, pbub, zFraction_);
648 }
649 }
650
657 {
658 const auto& iq = asImp_();
659 auto& fs = asImp_().fluidState_;
660
661 unsigned pvtRegionIdx = iq.pvtRegionIndex();
663
664 fs.setInvB(oilPhaseIdx, 1.0/bo_);
665 fs.setInvB(gasPhaseIdx, 1.0/bg_);
666
667 fs.setDensity(oilPhaseIdx,
668 fs.invB(oilPhaseIdx)
669 *(FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)
670 + (1.0-xVolume_)*fs.Rs()*FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx)
671 + xVolume_*fs.Rs()*zRefDensity_ ));
672 fs.setDensity(gasPhaseIdx,
673 fs.invB(gasPhaseIdx)
674 *(FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx)*(1.0-yVolume_)+yVolume_*zRefDensity_
675 + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx)*fs.Rv()));
676 }
677
678 const Evaluation& zFraction() const
679 { return zFraction_; }
680
681 const Evaluation& xVolume() const
682 { return xVolume_; }
683
684 const Evaluation& yVolume() const
685 { return yVolume_; }
686
687 const Evaluation& oilViscosity() const
688 { return oilViscosity_; }
689
690 const Evaluation& gasViscosity() const
691 { return gasViscosity_; }
692
693 const Evaluation& bo() const
694 { return bo_; }
695
696 const Evaluation& bg() const
697 { return bg_; }
698
699 const Evaluation& rs() const
700 { return rs_; }
701
702 const Evaluation& rv() const
703 { return rv_; }
704
705 const Evaluation zPureInvFormationVolumeFactor() const
706 { return 1.0/bz_; }
707
708 const Scalar& zRefDensity() const
709 { return zRefDensity_; }
710
711private:
712
713
714protected:
715 Implementation& asImp_()
716 { return *static_cast<Implementation*>(this); }
717
718 // Abstract "mass fraction" accounting for the solvent component. The relation between this
719 // quantity and the actual mass fraction of solvent, is implicitly defined from the specific
720 // pvt measurements as provided by kw PVTSOL.
721 Evaluation zFraction_;
722
723 // The solvent component is assumed gas at surface conditions
724 Evaluation xVolume_; // Solvent volume fraction of Rs
725 Evaluation yVolume_; // Solvent volume fraction of Sg/Bg
726
727 // Standard black oil parameters modified for presence of solvent
728 Evaluation oilViscosity_;
729 Evaluation gasViscosity_;
730 Evaluation bo_;
731 Evaluation bg_;
732 Evaluation rs_;
733 Evaluation rv_;
734
735 // Properties of pure solvent
736 Evaluation bz_;
738};
739
740template <class TypeTag>
742{
746
747public:
748
750 { }
751
752 void zFractionUpdate_(const ElementContext&,
753 unsigned,
754 unsigned)
755 { }
756
757 const Evaluation& xVolume() const
758 { throw std::runtime_error("xVolume() called but extbo is disabled"); }
759
760 const Evaluation& yVolume() const
761 { throw std::runtime_error("yVolume() called but extbo is disabled"); }
762
763 const Evaluation& oilViscosity() const
764 { throw std::runtime_error("oilViscosity() called but extbo is disabled"); }
765
766 const Evaluation& gasViscosity() const
767 { throw std::runtime_error("gasViscosity() called but extbo is disabled"); }
768
769 const Evaluation& rs() const
770 { throw std::runtime_error("rs() called but extbo is disabled"); }
771
772 const Evaluation& rv() const
773 { throw std::runtime_error("rv() called but extbo is disabled"); }
774
775 const Evaluation& zPureInvFormationVolumeFactor() const
776 { throw std::runtime_error("zPureInvFormationVolumeFactor() called but extbo is disabled"); }
777
778 const Evaluation& zFraction() const
779 { throw std::runtime_error("zFraction() called but extbo is disabled"); }
780
781 const Evaluation& zInverseFormationVolumeFactor() const
782 { throw std::runtime_error("zInverseFormationVolumeFactor() called but extbo is disabled"); }
783
784 const Scalar& zRefDensity() const
785 { throw std::runtime_error("zRefDensity() called but extbo is disabled"); }
786};
787
795template <class TypeTag, bool enableExtboV = getPropValue<TypeTag, Properties::EnableExtbo>()>
797{
799
807
808 using Toolbox = MathToolbox<Evaluation>;
809
810 static constexpr unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
811 static constexpr int dimWorld = GridView::dimensionworld;
812
813 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
814 typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
815
816public:
817
818private:
819 Implementation& asImp_()
820 { return *static_cast<Implementation*>(this); }
821
822};
823
824template <class TypeTag>
826{
829
830public:
831
832};
833
834} // namespace Opm
835
836#endif
Contains the parameters required to extend the black-oil model by solvent component....
Declares the properties required by the black oil model.
Provides the solvent specific extensive quantities to the generic black-oil module's extensive quanti...
Definition: blackoilextbomodules.hh:797
const Scalar & zRefDensity() const
Definition: blackoilextbomodules.hh:784
const Evaluation & xVolume() const
Definition: blackoilextbomodules.hh:757
const Evaluation & rs() const
Definition: blackoilextbomodules.hh:769
const Evaluation & zInverseFormationVolumeFactor() const
Definition: blackoilextbomodules.hh:781
const Evaluation & zPureInvFormationVolumeFactor() const
Definition: blackoilextbomodules.hh:775
const Evaluation & oilViscosity() const
Definition: blackoilextbomodules.hh:763
void zFractionUpdate_(const ElementContext &, unsigned, unsigned)
Definition: blackoilextbomodules.hh:752
const Evaluation & zFraction() const
Definition: blackoilextbomodules.hh:778
void zPvtUpdate_()
Definition: blackoilextbomodules.hh:749
const Evaluation & rv() const
Definition: blackoilextbomodules.hh:772
const Evaluation & gasViscosity() const
Definition: blackoilextbomodules.hh:766
const Evaluation & yVolume() const
Definition: blackoilextbomodules.hh:760
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilextbomodules.hh:562
const Evaluation & rv() const
Definition: blackoilextbomodules.hh:702
Evaluation zFraction_
Definition: blackoilextbomodules.hh:721
const Evaluation & oilViscosity() const
Definition: blackoilextbomodules.hh:687
const Evaluation & zFraction() const
Definition: blackoilextbomodules.hh:678
const Evaluation & bo() const
Definition: blackoilextbomodules.hh:693
void zFractionUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Compute extended pvt properties from table lookups.
Definition: blackoilextbomodules.hh:589
Scalar zRefDensity_
Definition: blackoilextbomodules.hh:737
void zPvtUpdate_()
Re-compute face densities to account for zFraction dependency.
Definition: blackoilextbomodules.hh:656
const Evaluation & yVolume() const
Definition: blackoilextbomodules.hh:684
Evaluation gasViscosity_
Definition: blackoilextbomodules.hh:729
Evaluation bg_
Definition: blackoilextbomodules.hh:731
Evaluation rs_
Definition: blackoilextbomodules.hh:732
const Evaluation & gasViscosity() const
Definition: blackoilextbomodules.hh:690
const Evaluation & rs() const
Definition: blackoilextbomodules.hh:699
Implementation & asImp_()
Definition: blackoilextbomodules.hh:715
Evaluation rv_
Definition: blackoilextbomodules.hh:733
const Evaluation & bg() const
Definition: blackoilextbomodules.hh:696
const Evaluation zPureInvFormationVolumeFactor() const
Definition: blackoilextbomodules.hh:705
const Scalar & zRefDensity() const
Definition: blackoilextbomodules.hh:708
Evaluation bo_
Definition: blackoilextbomodules.hh:730
Evaluation bz_
Definition: blackoilextbomodules.hh:736
Evaluation yVolume_
Definition: blackoilextbomodules.hh:725
Evaluation xVolume_
Definition: blackoilextbomodules.hh:724
Evaluation oilViscosity_
Definition: blackoilextbomodules.hh:728
const Evaluation & xVolume() const
Definition: blackoilextbomodules.hh:681
Contains the high level supplements required to extend the black oil model.
Definition: blackoilextbomodules.hh:69
static void addStorage(Dune::FieldVector< LhsEval, numEq > &storage, const IntensiveQuantities &intQuants)
Definition: blackoilextbomodules.hh:320
static bool primaryVarApplies(unsigned pvIdx)
Definition: blackoilextbomodules.hh:273
static Scalar zLim(unsigned regionIdx)
Definition: blackoilextbomodules.hh:531
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilextbomodules.hh:444
static Value yVolume(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:483
static void updatePrimaryVars(PrimaryVariables &newPv, const PrimaryVariables &oldPv, const EqVector &delta)
Do a Newton-Raphson update the primary variables of the solvents.
Definition: blackoilextbomodules.hh:420
static Value gasViscosity(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:503
static std::string eqName(unsigned eqIdx)
Definition: blackoilextbomodules.hh:304
static bool eqApplies(unsigned eqIdx)
Definition: blackoilextbomodules.hh:296
static Value pbubRv(unsigned pvtRegionIdx, const Value &z, const Value &rv)
Definition: blackoilextbomodules.hh:493
static Scalar primaryVarWeight(unsigned pvIdx)
Definition: blackoilextbomodules.hh:288
static Value rs(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:518
static void registerOutputModules(Model &, Simulator &)
Register all solvent specific VTK and ECL output modules.
Definition: blackoilextbomodules.hh:268
static Scalar computeUpdateError(const PrimaryVariables &, const EqVector &)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilextbomodules.hh:432
static Value xVolume(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:478
static void serializeEntity(const Model &model, std::ostream &outstream, const DofEntity &dof)
Definition: blackoilextbomodules.hh:451
static void computeFlux(RateVector &flux, const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Definition: blackoilextbomodules.hh:354
static void registerParameters()
Register all run-time parameters for the black-oil solvent module.
Definition: blackoilextbomodules.hh:261
static void deserializeEntity(Model &model, std::istream &instream, const DofEntity &dof)
Definition: blackoilextbomodules.hh:462
static std::string primaryVarName(unsigned pvIdx)
Definition: blackoilextbomodules.hh:281
static Value rv(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:523
static Value pbubRs(unsigned pvtRegionIdx, const Value &z, const Value &rs)
Definition: blackoilextbomodules.hh:488
static Value oilCmp(unsigned pvtRegionIdx, const Value &z)
Definition: blackoilextbomodules.hh:536
static Scalar referenceDensity(unsigned regionIdx)
Definition: blackoilextbomodules.hh:527
static Value bg(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:513
static Value oilViscosity(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:498
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar zFraction)
Assign the solvent specific primary variables to a PrimaryVariables object.
Definition: blackoilextbomodules.hh:410
static Scalar eqWeight(unsigned eqIdx)
Definition: blackoilextbomodules.hh:311
static Value bo(unsigned pvtRegionIdx, const Value &pressure, const Value &z)
Definition: blackoilextbomodules.hh:508
static Value gasCmp(unsigned pvtRegionIdx, const Value &z)
Definition: blackoilextbomodules.hh:541
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
Struct holding the parameters for the BlackoilExtboModule class.
Definition: blackoilextboparams.hh:44
UniformXTabulated2DFunction< Scalar > Tabulated2DFunction
Definition: blackoilextboparams.hh:46
Tabulated1DFunction< Scalar > TabulatedFunction
Definition: blackoilextboparams.hh:45