28 #ifndef EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
29 #define EWOMS_BLACK_OIL_SOLVENT_MODULE_HH
35 #include <opm/material/fluidsystems/blackoilpvt/SolventPvt.hpp>
36 #include <opm/material/common/Tabulated1DFunction.hpp>
39 #include <opm/parser/eclipse/Deck/Deck.hpp>
40 #include <opm/parser/eclipse/EclipseState/EclipseState.hpp>
41 #include <opm/parser/eclipse/EclipseState/Tables/SsfnTable.hpp>
42 #include <opm/parser/eclipse/EclipseState/Tables/Sof2Table.hpp>
43 #include <opm/parser/eclipse/EclipseState/Tables/MsfnTable.hpp>
44 #include <opm/parser/eclipse/EclipseState/Tables/PmiscTable.hpp>
45 #include <opm/parser/eclipse/EclipseState/Tables/MiscTable.hpp>
46 #include <opm/parser/eclipse/EclipseState/Tables/SorwmisTable.hpp>
47 #include <opm/parser/eclipse/EclipseState/Tables/SgcwmisTable.hpp>
48 #include <opm/parser/eclipse/EclipseState/Tables/TlpmixpaTable.hpp>
51 #include <opm/common/Valgrind.hpp>
52 #include <opm/common/Unused.hpp>
53 #include <opm/common/ErrorMacros.hpp>
54 #include <opm/common/Exceptions.hpp>
56 #include <dune/common/fvector.hh>
66 template <
class TypeTag,
bool enableSolventV = GET_PROP_VALUE(TypeTag, EnableSolvent)>
67 class BlackOilSolventModule
69 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
70 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
71 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
72 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
73 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
74 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
75 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
76 typedef typename GET_PROP_TYPE(TypeTag, Model) Model;
77 typedef typename GET_PROP_TYPE(TypeTag, Simulator) Simulator;
78 typedef typename GET_PROP_TYPE(TypeTag, EqVector) EqVector;
79 typedef typename GET_PROP_TYPE(TypeTag, RateVector) RateVector;
80 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
82 typedef Opm::MathToolbox<Evaluation> Toolbox;
83 typedef Opm::SolventPvt<Scalar> SolventPvt;
85 typedef typename Opm::Tabulated1DFunction<Scalar> TabulatedFunction;
87 static constexpr
unsigned solventSaturationIdx = Indices::solventSaturationIdx;
88 static constexpr
unsigned contiSolventEqIdx = Indices::contiSolventEqIdx;
89 static constexpr
unsigned enableSolvent = enableSolventV;
91 static constexpr
unsigned numPhases = FluidSystem::numPhases;
92 static constexpr
bool blackoilConserveSurfaceVolume =
GET_PROP_VALUE(TypeTag, BlackoilConserveSurfaceVolume);
100 static void initFromDeck(
const Opm::Deck& deck,
const Opm::EclipseState& eclState)
104 if (enableSolvent && !deck.hasKeyword(
"SOLVENT")) {
105 OPM_THROW(std::runtime_error,
106 "Non-trivial solvent treatment requested at compile time, but "
107 "the deck does not contain the SOLVENT keyword");
109 else if (!enableSolvent && deck.hasKeyword(
"SOLVENT")) {
110 OPM_THROW(std::runtime_error,
111 "Solvent treatment disabled at compile time, but the deck "
112 "contains the SOLVENT keyword");
115 if (!deck.hasKeyword(
"SOLVENT"))
118 solventPvt_.initFromDeck(deck, eclState);
120 const auto& tableManager = eclState.getTableManager();
122 const auto& ssfnTables = tableManager.getSsfnTables();
123 unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
125 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
126 const auto& ssfnTable = ssfnTables.template getTable<Opm::SsfnTable>(satRegionIdx);
127 ssfnKrg_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(),
128 ssfnTable.getGasRelPermMultiplierColumn(),
130 ssfnKrs_[satRegionIdx].setXYContainers(ssfnTable.getSolventFractionColumn(),
131 ssfnTable.getSolventRelPermMultiplierColumn(),
137 if (deck.hasKeyword(
"MISCIBLE")) {
140 unsigned numSatRegions = tableManager.getTabdims().getNumSatTables();
141 unsigned numMiscRegions = 1;
144 const auto& sof2Tables = tableManager.getSof2Tables();
145 if (!sof2Tables.empty()) {
148 sof2Krn_.resize(numSatRegions);
149 for (
unsigned satRegionIdx = 0; satRegionIdx < numSatRegions; ++ satRegionIdx) {
150 const auto& sof2Table = sof2Tables.template getTable<Opm::Sof2Table>(satRegionIdx);
151 sof2Krn_[satRegionIdx].setXYContainers(sof2Table.getSoColumn(),
152 sof2Table.getKroColumn(),
157 OPM_THROW(std::runtime_error,
"SOF2 must be specified in MISCIBLE (SOLVENT) runs\n");
160 const auto& miscTables = tableManager.getMiscTables();
161 if (!miscTables.empty()) {
163 assert(numMiscRegions == miscTables.size());
166 misc_.resize(numMiscRegions);
167 for (
unsigned miscRegionIdx = 0; miscRegionIdx < numMiscRegions; ++miscRegionIdx) {
168 const auto& miscTable = miscTables.template getTable<Opm::MiscTable>(miscRegionIdx);
171 const auto& solventFraction = miscTable.getSolventFractionColumn();
172 const auto& misc = miscTable.getMiscibilityColumn();
173 misc_[miscRegionIdx].setXYContainers(solventFraction, misc);
177 OPM_THROW(std::runtime_error,
"MISC must be specified in MISCIBLE (SOLVENT) runs\n");
181 pmisc_.resize(numMiscRegions);
182 const auto& pmiscTables = tableManager.getPmiscTables();
183 if (!pmiscTables.empty()) {
185 assert(numMiscRegions == pmiscTables.size());
187 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
188 const auto& pmiscTable = pmiscTables.template getTable<Opm::PmiscTable>(regionIdx);
191 const auto& po = pmiscTable.getOilPhasePressureColumn();
192 const auto& pmisc = pmiscTable.getMiscibilityColumn();
194 pmisc_[regionIdx].setXYContainers(po, pmisc);
198 std::vector<double> x = {0.0,1.0e20};
199 std::vector<double> y = {1.0,1.0};
200 TabulatedFunction constant = TabulatedFunction(2, x, y);
201 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
208 msfnKrsg_.resize(numSatRegions);
209 msfnKro_.resize(numSatRegions);
210 const auto& msfnTables = tableManager.getMsfnTables();
211 if (!msfnTables.empty()) {
213 assert(numSatRegions == msfnTables.size());
216 for (
unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) {
217 const Opm::MsfnTable& msfnTable = msfnTables.template getTable<Opm::MsfnTable>(regionIdx);
221 const auto& Ssg = msfnTable.getGasPhaseFractionColumn();
222 const auto& krsg = msfnTable.getGasSolventRelpermMultiplierColumn();
223 const auto& kro = msfnTable.getOilRelpermMultiplierColumn();
225 msfnKrsg_[regionIdx].setXYContainers(Ssg, krsg);
226 msfnKro_[regionIdx].setXYContainers(Ssg, kro);
230 std::vector<double> x = {0.0,1.0};
231 std::vector<double> y = {1.0,0.0};
232 TabulatedFunction unit = TabulatedFunction(2, x, x);
233 TabulatedFunction inv_unit = TabulatedFunction(2, x, y);
235 for (
unsigned regionIdx = 0; regionIdx < numSatRegions; ++regionIdx) {
236 setMsfn(regionIdx, unit, inv_unit);
240 sorwmis_.resize(numMiscRegions);
241 const auto& sorwmisTables = tableManager.getSorwmisTables();
242 if (!sorwmisTables.empty()) {
243 assert(numMiscRegions == sorwmisTables.size());
245 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
246 const auto& sorwmisTable = sorwmisTables.template getTable<Opm::SorwmisTable>(regionIdx);
249 const auto& sw = sorwmisTable.getWaterSaturationColumn();
250 const auto& sorwmis = sorwmisTable.getMiscibleResidualOilColumn();
252 sorwmis_[regionIdx].setXYContainers(sw, sorwmis);
256 std::vector<double> x = {0.0,1.0};
257 std::vector<double> y = {0.0,0.0};
258 TabulatedFunction zero = TabulatedFunction(2, x, y);
259 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
265 sgcwmis_.resize(numMiscRegions);
266 const auto& sgcwmisTables = tableManager.getSgcwmisTables();
267 if (!sgcwmisTables.empty()) {
269 assert(numMiscRegions ==sgcwmisTables.size());
271 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
272 const auto& sgcwmisTable = sgcwmisTables.template getTable<Opm::SgcwmisTable>(regionIdx);
275 const auto& sw = sgcwmisTable.getWaterSaturationColumn();
276 const auto& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn();
278 sgcwmis_[regionIdx].setXYContainers(sw, sgcwmis);
282 std::vector<double> x = {0.0,1.0};
283 std::vector<double> y = {0.0,0.0};
284 TabulatedFunction zero = TabulatedFunction(2, x, y);
285 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
291 if (deck.hasKeyword(
"TLMIXPAR")) {
293 tlMixParamViscosity_.resize(numMiscRegions);
294 tlMixParamDensity_.resize(numMiscRegions);
296 assert(numMiscRegions == deck.getKeyword(
"TLMIXPAR").size() );
298 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
299 const auto& tlmixparRecord = deck.getKeyword(
"TLMIXPAR").getRecord(regionIdx);
300 const auto& mix_params_viscosity = tlmixparRecord.getItem(
"TL_VISCOSITY_PARAMETER").getSIDoubleData();
301 tlMixParamViscosity_[regionIdx] = mix_params_viscosity[0];
302 const auto& mix_params_density = tlmixparRecord.getItem(
"TL_DENSITY_PARAMETER").getSIDoubleData();
303 const int numDensityItems = mix_params_density.size();
304 if (numDensityItems == 0) {
305 tlMixParamDensity_[regionIdx] = tlMixParamViscosity_[regionIdx];
306 }
else if (numDensityItems == 1) {
307 tlMixParamDensity_[regionIdx] = mix_params_density[0];
309 OPM_THROW(std::runtime_error,
"Only one value can be entered for the TL parameter pr MISC region.");
313 OPM_THROW(std::runtime_error,
"TLMIXPAR must be specified in MISCIBLE (SOLVENT) runs\n");
317 tlPMixTable_.resize(numMiscRegions);
318 if (deck.hasKeyword(
"TLPMIXPA")) {
319 const auto& tlpmixparTables = tableManager.getTlpmixpaTables();
320 if (!tlpmixparTables.empty()) {
322 assert(numMiscRegions == tlpmixparTables.size());
323 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
324 const auto& tlpmixparTable = tlpmixparTables.template getTable<Opm::TlpmixpaTable>(regionIdx);
327 const auto& po = tlpmixparTable.getOilPhasePressureColumn();
328 const auto& tlpmixpa = tlpmixparTable.getMiscibilityColumn();
330 tlPMixTable_[regionIdx].setXYContainers(po, tlpmixpa);
335 if (pmisc_.size() > 0) {
336 tlPMixTable_ = pmisc_;
338 OPM_THROW(std::invalid_argument,
"If the pressure dependent TL values in TLPMIXPA is defaulted (no entries), then the PMISC tables must be specified.");
343 std::vector<double> x = {0.0,1.0e20};
344 std::vector<double> y = {1.0,1.0};
345 TabulatedFunction ones = TabulatedFunction(2, x, y);
346 for (
unsigned regionIdx = 0; regionIdx < numMiscRegions; ++regionIdx) {
366 ssfnKrg_.resize(numRegions);
367 ssfnKrs_.resize(numRegions);
376 const TabulatedFunction& ssfnKrg,
377 const TabulatedFunction& ssfnKrs)
379 ssfnKrg_[satRegionIdx] = ssfnKrg;
380 ssfnKrs_[satRegionIdx] = ssfnKrs;
389 const TabulatedFunction& sof2Krn)
391 sof2Krn_[satRegionIdx] = sof2Krn;
400 const TabulatedFunction& misc)
402 misc_[miscRegionIdx] = misc;
411 const TabulatedFunction& pmisc)
413 pmisc_[miscRegionIdx] = pmisc;
422 const TabulatedFunction& msfnKrsg,
423 const TabulatedFunction& msfnKro)
425 msfnKrsg_[satRegionIdx] = msfnKrsg;
426 msfnKro_[satRegionIdx] = msfnKro;
435 const TabulatedFunction& sorwmis)
437 sorwmis_[miscRegionIdx] = sorwmis;
446 const TabulatedFunction& sgcwmis)
448 sgcwmis_[miscRegionIdx] = sgcwmis;
457 const Scalar& tlMixParamViscosity,
458 const Scalar& tlMixParamDensity)
460 tlMixParamViscosity_[miscRegionIdx] = tlMixParamViscosity;
461 tlMixParamDensity_[miscRegionIdx] = tlMixParamDensity;
470 const TabulatedFunction& tlPMixTable)
472 tlPMixTable_[miscRegionIdx] = tlPMixTable;
479 { solventPvt_ = value; }
482 static void setIsMiscible(
const bool isMiscible)
483 { isMiscible_ = isMiscible; }
501 Simulator& simulator)
510 static bool primaryVarApplies(
unsigned pvIdx)
516 return pvIdx == solventSaturationIdx;
519 static std::string primaryVarName(
unsigned pvIdx OPM_OPTIM_UNUSED)
521 assert(primaryVarApplies(pvIdx));
523 return "saturation_solvent";
526 static Scalar primaryVarWeight(
unsigned pvIdx OPM_OPTIM_UNUSED)
528 assert(primaryVarApplies(pvIdx));
531 return static_cast<Scalar
>(1.0);
534 static bool eqApplies(
unsigned eqIdx)
539 return eqIdx == contiSolventEqIdx;
542 static std::string eqName(
unsigned eqIdx OPM_OPTIM_UNUSED)
544 assert(eqApplies(eqIdx));
546 return "conti^solvent";
549 static Scalar eqWeight(
unsigned eqIdx OPM_OPTIM_UNUSED)
551 assert(eqApplies(eqIdx));
554 return static_cast<Scalar
>(1.0);
557 template <
class LhsEval>
558 static void addStorage(Dune::FieldVector<LhsEval, numEq>& storage,
559 const IntensiveQuantities& intQuants)
563 if (blackoilConserveSurfaceVolume) {
564 storage[contiSolventEqIdx] +=
565 Toolbox::template decay<LhsEval>(intQuants.porosity())
566 * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
567 * Toolbox::template decay<LhsEval>(intQuants.solventInverseFormationVolumeFactor());
569 storage[contiSolventEqIdx] +=
570 Toolbox::template decay<LhsEval>(intQuants.porosity())
571 * Toolbox::template decay<LhsEval>(intQuants.solventSaturation())
572 * Toolbox::template decay<LhsEval>(intQuants.solventDensity());
576 static void computeFlux(RateVector& flux,
577 const ElementContext& elemCtx,
585 const auto& extQuants = elemCtx.extensiveQuantities(scvfIdx, timeIdx);
587 unsigned upIdx = extQuants.solventUpstreamIndex();
588 unsigned inIdx = extQuants.interiorIndex();
589 const auto& up = elemCtx.intensiveQuantities(upIdx, timeIdx);
591 if (blackoilConserveSurfaceVolume) {
593 flux[contiSolventEqIdx] =
594 extQuants.solventVolumeFlux()
595 *up.solventInverseFormationVolumeFactor();
597 flux[contiSolventEqIdx] =
598 extQuants.solventVolumeFlux()
599 *Opm::decay<Scalar>(up.solventInverseFormationVolumeFactor());
602 flux[contiSolventEqIdx] =
603 extQuants.solventVolumeFlux()
604 *up.solventDensity();
606 flux[contiSolventEqIdx] =
607 extQuants.solventVolumeFlux()
608 *Opm::decay<Scalar>(up.solventDensity());
616 Scalar solventSaturation)
621 priVars[solventSaturationIdx] = solventSaturation;
628 const PrimaryVariables& oldPv,
629 const EqVector& delta)
635 newPv[solventSaturationIdx] = oldPv[solventSaturationIdx] - delta[solventSaturationIdx];
642 const EqVector& delta OPM_UNUSED)
647 return static_cast<Scalar
>(0.0);
656 return std::abs(Toolbox::scalarValue(resid[contiSolventEqIdx]));
659 template <
class DofEntity>
660 static void serializeEntity(
const Model& model, std::ostream& outstream,
const DofEntity& dof)
665 unsigned dofIdx = model.dofMapper().index(dof);
667 const PrimaryVariables& priVars = model.solution(0)[dofIdx];
668 outstream << priVars[solventSaturationIdx];
671 template <
class DofEntity>
672 static void deserializeEntity(Model& model, std::istream& instream,
const DofEntity& dof)
677 unsigned dofIdx = model.dofMapper().index(dof);
679 PrimaryVariables& priVars0 = model.solution(0)[dofIdx];
680 PrimaryVariables& priVars1 = model.solution(1)[dofIdx];
682 instream >> priVars0[solventSaturationIdx];
685 priVars1 = priVars0[solventSaturationIdx];
688 static const SolventPvt& solventPvt()
689 {
return solventPvt_; }
691 static const TabulatedFunction& ssfnKrg(
const ElementContext& elemCtx,
695 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
696 return ssfnKrg_[satnumRegionIdx];
699 static const TabulatedFunction& ssfnKrs(
const ElementContext& elemCtx,
703 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
704 return ssfnKrs_[satnumRegionIdx];
707 static const TabulatedFunction& sof2Krn(
const ElementContext& elemCtx,
711 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
712 return sof2Krn_[satnumRegionIdx];
715 static const TabulatedFunction& misc(
const ElementContext& elemCtx,
719 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
720 return misc_[miscnumRegionIdx];
723 static const TabulatedFunction& pmisc(
const ElementContext& elemCtx,
727 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
728 return pmisc_[miscnumRegionIdx];
731 static const TabulatedFunction& msfnKrsg(
const ElementContext& elemCtx,
735 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
736 return msfnKrsg_[satnumRegionIdx];
739 static const TabulatedFunction& msfnKro(
const ElementContext& elemCtx,
743 unsigned satnumRegionIdx = elemCtx.problem().satnumRegionIndex(elemCtx, scvIdx, timeIdx);
744 return msfnKro_[satnumRegionIdx];
747 static const TabulatedFunction& sorwmis(
const ElementContext& elemCtx,
751 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
752 return sorwmis_[miscnumRegionIdx];
755 static const TabulatedFunction& sgcwmis(
const ElementContext& elemCtx,
759 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
760 return sgcwmis_[miscnumRegionIdx];
763 static const TabulatedFunction& tlPMixTable(
const ElementContext& elemCtx,
767 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
768 return tlPMixTable_[miscnumRegionIdx];
771 static const Scalar& tlMixParamViscosity(
const ElementContext& elemCtx,
775 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
776 return tlMixParamViscosity_[miscnumRegionIdx];
779 static const Scalar& tlMixParamDensity(
const ElementContext& elemCtx,
783 unsigned miscnumRegionIdx = elemCtx.problem().miscnumRegionIndex(elemCtx, scvIdx, timeIdx);
784 return tlMixParamDensity_[miscnumRegionIdx];
787 static bool isMiscible()
794 static SolventPvt solventPvt_;
796 static std::vector<TabulatedFunction> ssfnKrg_;
797 static std::vector<TabulatedFunction> ssfnKrs_;
798 static std::vector<TabulatedFunction> sof2Krn_;
799 static std::vector<TabulatedFunction> misc_;
800 static std::vector<TabulatedFunction> pmisc_;
801 static std::vector<TabulatedFunction> msfnKrsg_;
802 static std::vector<TabulatedFunction> msfnKro_;
803 static std::vector<TabulatedFunction> sorwmis_;
804 static std::vector<TabulatedFunction> sgcwmis_;
806 static std::vector<Scalar> tlMixParamViscosity_;
807 static std::vector<Scalar> tlMixParamDensity_;
808 static std::vector<TabulatedFunction> tlPMixTable_;
810 static bool isMiscible_;
813 template <
class TypeTag,
bool enableSolventV>
814 typename BlackOilSolventModule<TypeTag, enableSolventV>::SolventPvt
815 BlackOilSolventModule<TypeTag, enableSolventV>::solventPvt_;
817 template <
class TypeTag,
bool enableSolventV>
818 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
819 BlackOilSolventModule<TypeTag, enableSolventV>::ssfnKrg_;
821 template <
class TypeTag,
bool enableSolventV>
822 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
823 BlackOilSolventModule<TypeTag, enableSolventV>::ssfnKrs_;
825 template <
class TypeTag,
bool enableSolventV>
826 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
827 BlackOilSolventModule<TypeTag, enableSolventV>::sof2Krn_;
829 template <
class TypeTag,
bool enableSolventV>
830 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
831 BlackOilSolventModule<TypeTag, enableSolventV>::misc_;
833 template <
class TypeTag,
bool enableSolventV>
834 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
835 BlackOilSolventModule<TypeTag, enableSolventV>::pmisc_;
837 template <
class TypeTag,
bool enableSolventV>
838 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
839 BlackOilSolventModule<TypeTag, enableSolventV>::msfnKrsg_;
841 template <
class TypeTag,
bool enableSolventV>
842 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
843 BlackOilSolventModule<TypeTag, enableSolventV>::msfnKro_;
845 template <
class TypeTag,
bool enableSolventV>
846 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
847 BlackOilSolventModule<TypeTag, enableSolventV>::sorwmis_;
849 template <
class TypeTag,
bool enableSolventV>
850 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
851 BlackOilSolventModule<TypeTag, enableSolventV>::sgcwmis_;
853 template <
class TypeTag,
bool enableSolventV>
854 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::Scalar>
855 BlackOilSolventModule<TypeTag, enableSolventV>::tlMixParamViscosity_;
857 template <
class TypeTag,
bool enableSolventV>
858 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::Scalar>
859 BlackOilSolventModule<TypeTag, enableSolventV>::tlMixParamDensity_;
862 template <
class TypeTag,
bool enableSolventV>
863 std::vector<typename BlackOilSolventModule<TypeTag, enableSolventV>::TabulatedFunction>
864 BlackOilSolventModule<TypeTag, enableSolventV>::tlPMixTable_;
866 template <
class TypeTag,
bool enableSolventV>
868 BlackOilSolventModule<TypeTag, enableSolventV>::isMiscible_;
878 template <
class TypeTag,
bool enableSolventV = GET_PROP_VALUE(TypeTag, EnableSolvent)>
881 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) Implementation;
883 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
884 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
885 typedef typename GET_PROP_TYPE(TypeTag, PrimaryVariables) PrimaryVariables;
886 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
887 typedef typename GET_PROP_TYPE(TypeTag, MaterialLaw) MaterialLaw;
888 typedef typename GET_PROP_TYPE(TypeTag, Indices) Indices;
889 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
894 static constexpr
int solventSaturationIdx = Indices::solventSaturationIdx;
895 static constexpr
int oilPhaseIdx = FluidSystem::oilPhaseIdx;
896 static constexpr
int gasPhaseIdx = FluidSystem::gasPhaseIdx;
897 static constexpr
int waterPhaseIdx = FluidSystem::waterPhaseIdx;
898 static constexpr
double cutOff = 1e-16;
912 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
913 auto& fs = asImp_().fluidState_;
914 solventSaturation_ = priVars.makeEvaluation(solventSaturationIdx, timeIdx);
915 hydrocarbonSaturation_ = fs.saturation(gasPhaseIdx);
918 if (std::abs(solventSaturation().value()) < cutOff) {
924 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_ + solventSaturation_);
940 auto& fs = asImp_().fluidState_;
941 fs.setSaturation(gasPhaseIdx, hydrocarbonSaturation_);
943 solventMobility_ = 0.0;
946 if (std::abs(solventSaturation().value()) < cutOff) {
951 if(SolventModule::isMiscible()) {
953 const Evaluation& p = fs.pressure(oilPhaseIdx);
954 const Evaluation pmisc = SolventModule::pmisc(elemCtx, dofIdx, timeIdx).eval(p,
true);
955 const Evaluation& pgImisc = fs.pressure(gasPhaseIdx);
958 const auto& problem = elemCtx.problem();
959 const PrimaryVariables& priVars = elemCtx.primaryVars(dofIdx, timeIdx);
960 Evaluation pgMisc = 0.0;
961 Evaluation pC[numPhases];
962 const auto& materialParams = problem.materialLawParams(elemCtx, dofIdx, timeIdx);
963 MaterialLaw::capillaryPressures(pC, materialParams, fs);
966 if (priVars.primaryVarsMeaning() == PrimaryVariables::Sw_pg_Rv) {
967 pgMisc = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
969 const Evaluation& po = priVars.makeEvaluation(Indices::pressureSwitchIdx, timeIdx);
970 pgMisc = po + (pC[gasPhaseIdx] - pC[oilPhaseIdx]);
973 fs.setPressure(gasPhaseIdx, pmisc * pgMisc + (1.0 - pmisc) * pgImisc);
977 Evaluation gasSolventSat = hydrocarbonSaturation_ + solventSaturation_;
979 if (std::abs(gasSolventSat.value()) < cutOff) {
983 Evaluation Fhydgas = hydrocarbonSaturation_/gasSolventSat;
984 Evaluation Fsolgas = solventSaturation_/gasSolventSat;
987 if(SolventModule::isMiscible()) {
988 const auto& misc = SolventModule::misc(elemCtx, dofIdx, timeIdx);
989 const auto& pmisc = SolventModule::pmisc(elemCtx, dofIdx, timeIdx);
990 const Evaluation& p = fs.pressure(oilPhaseIdx);
991 const Evaluation miscibility = misc.eval(Fsolgas,
true) * pmisc.eval(p,
true);
994 unsigned cellIdx = elemCtx.globalSpaceIndex(dofIdx, timeIdx);
995 const auto& materialLawManager = elemCtx.problem().materialLawManager();
996 const auto& scaledDrainageInfo =
997 materialLawManager->oilWaterScaledEpsInfoDrainage(cellIdx);
999 const Scalar& sgcr = scaledDrainageInfo.Sgcr;
1000 const Scalar& sogcr = scaledDrainageInfo.Sogcr;
1001 const Evaluation& sw = fs.saturation(waterPhaseIdx);
1002 const auto& sorwmis = SolventModule::sorwmis(elemCtx, dofIdx, timeIdx);
1003 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, dofIdx, timeIdx);
1005 Evaluation sor = miscibility * sorwmis.eval(sw,
true) + ( 1.0 - miscibility) * sogcr;
1006 Evaluation sgc = miscibility * sgcwmis.eval(sw,
true) + ( 1.0 - miscibility) * sgcr;
1008 const Evaluation oilGasSolventSat = gasSolventSat + fs.saturation(oilPhaseIdx);
1009 const Evaluation zero = 0.0;
1010 const Evaluation oilGasSolventEffSat = std::max(oilGasSolventSat - sor - sgc, zero);
1012 Evaluation F_totalGas = 0.0;
1013 if (std::abs(oilGasSolventEffSat.value()) > cutOff) {
1014 const Evaluation gasSolventEffSat = std::max(gasSolventSat - sgc, zero);
1015 F_totalGas = gasSolventEffSat / oilGasSolventEffSat;
1017 const auto& msfnKro = SolventModule::msfnKro(elemCtx, dofIdx, timeIdx);
1018 const auto& msfnKrsg = SolventModule::msfnKrsg(elemCtx, dofIdx, timeIdx);
1019 const auto& sof2Krn = SolventModule::sof2Krn(elemCtx, dofIdx, timeIdx);
1021 const Evaluation mkrgt = msfnKrsg.eval( F_totalGas,
true ) * sof2Krn.eval( oilGasSolventSat,
true);
1022 const Evaluation mkro = msfnKro.eval( F_totalGas,
true ) * sof2Krn.eval( oilGasSolventSat,
true );
1024 Evaluation& kro = asImp_().mobility_[oilPhaseIdx];
1025 Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
1028 krg *= (1.0 - miscibility);
1029 krg += miscibility * mkrgt;
1030 kro *= (1.0 - miscibility);
1031 kro += miscibility * mkro;
1037 const auto& ssfnKrg = SolventModule::ssfnKrg(elemCtx, dofIdx, timeIdx);
1038 const auto& ssfnKrs = SolventModule::ssfnKrs(elemCtx, dofIdx, timeIdx);
1040 Evaluation& krg = asImp_().mobility_[gasPhaseIdx];
1041 solventMobility_ = krg * ssfnKrs.eval(Fsolgas,
true);
1042 krg *= ssfnKrg.eval(Fhydgas,
true);
1056 const auto& iq = asImp_();
1057 const auto& fs = iq.fluidState();
1058 const auto& solventPvt = SolventModule::solventPvt();
1060 unsigned pvtRegionIdx = iq.pvtRegionIndex();
1061 solventRefDensity_ = solventPvt.referenceDensity(pvtRegionIdx);
1062 const Evaluation& T = fs.temperature(gasPhaseIdx);
1063 const Evaluation& p = fs.pressure(gasPhaseIdx);
1064 solventInvFormationVolumeFactor_ = solventPvt.inverseFormationVolumeFactor(pvtRegionIdx, T, p);
1066 solventDensity_ = solventInvFormationVolumeFactor_*solventRefDensity_;
1067 solventViscosity_ = solventPvt.viscosity(pvtRegionIdx, T, p);
1069 effectiveProperties(elemCtx, scvIdx, timeIdx);
1071 solventMobility_ /= solventViscosity_;
1076 const Evaluation& solventSaturation()
const
1077 {
return solventSaturation_; }
1079 const Evaluation& solventDensity()
const
1080 {
return solventDensity_; }
1082 const Evaluation& solventViscosity()
const
1083 {
return solventViscosity_; }
1085 const Evaluation& solventMobility()
const
1086 {
return solventMobility_; }
1088 const Evaluation& solventInverseFormationVolumeFactor()
const
1089 {
return solventInvFormationVolumeFactor_; }
1092 const Scalar& solventRefDensity()
const
1093 {
return solventRefDensity_; }
1098 void effectiveProperties(
const ElementContext& elemCtx,
1102 if (!SolventModule::isMiscible()) {
1108 if ( solventSaturation() < cutOff) {
1112 auto& fs = asImp_().fluidState_;
1115 const auto& sorwmis = SolventModule::sorwmis(elemCtx, scvIdx, timeIdx);
1116 const auto& sgcwmis = SolventModule::sgcwmis(elemCtx, scvIdx, timeIdx);
1117 const Evaluation& sw = fs.saturation(waterPhaseIdx);
1118 const Evaluation oilEffSat = fs.saturation(oilPhaseIdx) - sorwmis.eval(sw,
true);
1119 const Evaluation gasEffSat = fs.saturation(gasPhaseIdx) - sgcwmis.eval(sw,
true);
1120 const Evaluation solventEffSat = solventSaturation() - sgcwmis.eval(sw,
true);
1121 const Evaluation oilGasSolventEffSat = oilEffSat + gasEffSat + solventEffSat;
1122 const Evaluation oilSolventEffSat = oilEffSat + solventEffSat;
1123 const Evaluation solventGasEffSat = solventEffSat + gasEffSat;
1126 const Evaluation& muGas = fs.viscosity(gasPhaseIdx);
1127 const Evaluation& muOil = fs.viscosity(oilPhaseIdx);
1128 const Evaluation& muSolvent = solventViscosity_;
1130 const Evaluation muOilPow = pow(muOil, 0.25);
1131 const Evaluation muGasPow = pow(muGas, 0.25);
1132 const Evaluation muSolventPow = pow(muSolvent, 0.25);
1134 Evaluation muMixOilSolvent = muOil;
1135 if ( std::abs(oilSolventEffSat.value()) > cutOff)
1136 muMixOilSolvent *= muSolvent / pow( ( (oilEffSat / oilSolventEffSat) * muSolventPow) + ( (solventEffSat / oilSolventEffSat) * muOilPow) , 4.0);
1137 Evaluation muMixSolventGas = muGas;
1138 if ( std::abs(solventGasEffSat.value()) > cutOff)
1139 muMixSolventGas *= muSolvent / pow( ( (gasEffSat / solventGasEffSat) * muSolventPow) + ( (solventEffSat / solventGasEffSat) * muGasPow) , 4.0);
1141 Evaluation muMixSolventGasOil = muOil;
1142 if (std::abs(oilGasSolventEffSat.value()) > cutOff)
1143 muMixSolventGasOil *= muSolvent * muGas / pow( ( (oilEffSat / oilGasSolventEffSat) * muSolventPow * muGasPow)
1144 + ( (solventEffSat / oilGasSolventEffSat) * muOilPow * muGasPow) + ( (gasEffSat / oilGasSolventEffSat) * muSolventPow * muOilPow), 4.0);
1149 const Evaluation& po = fs.pressure(oilPhaseIdx);
1150 const auto& tlPMixTable = SolventModule::tlPMixTable(elemCtx, scvIdx, timeIdx);
1151 const Evaluation tlMixParamMu = SolventModule::tlMixParamViscosity(elemCtx, scvIdx, timeIdx) * tlPMixTable.eval(po,
true);
1153 Evaluation muOilEff = pow(muOil,1.0 - tlMixParamMu) * pow(muMixOilSolvent, tlMixParamMu);
1154 Evaluation muGasEff = pow(muGas,1.0 - tlMixParamMu) * pow(muMixSolventGas, tlMixParamMu);
1155 Evaluation muSolventEff = pow(muSolvent,1.0 - tlMixParamMu) * pow(muMixSolventGasOil, tlMixParamMu);
1158 const Evaluation& rhoGas = fs.density(gasPhaseIdx);
1159 const Evaluation& rhoOil = fs.density(oilPhaseIdx);
1160 const Evaluation& rhoSolvent = solventDensity_;
1165 const Evaluation tlMixParamRho = SolventModule::tlMixParamDensity(elemCtx, scvIdx, timeIdx) * tlPMixTable.eval(po,
true);
1169 const Evaluation muOilEffPow = pow(pow(muOil, 1.0 - tlMixParamRho) * pow(muMixOilSolvent, tlMixParamRho), 0.25);
1170 const Evaluation muGasEffPow = pow(pow(muGas, 1.0 - tlMixParamRho) * pow(muMixSolventGas, tlMixParamRho), 0.25);
1171 const Evaluation muSolventEffPow = pow(pow(muSolvent, 1.0 - tlMixParamRho) * pow(muMixSolventGasOil, tlMixParamRho), 0.25);
1173 const Evaluation oilGasEffSaturation = oilEffSat + gasEffSat;
1174 Evaluation sof = 0.0;
1175 Evaluation sgf = 0.0;
1176 if ( std::abs(oilGasEffSaturation.value()) > cutOff ) {
1177 sof = oilEffSat / oilGasEffSaturation;
1178 sgf = gasEffSat / oilGasEffSaturation;
1181 const Evaluation muSolventOilGasPow = muSolventPow * ( (sgf * muOilPow) + (sof * muGasPow) );
1183 Evaluation rhoMixSolventGasOil = 0.0;
1184 if (std::abs(oilGasSolventEffSat.value()) > cutOff )
1185 rhoMixSolventGasOil = (rhoOil * oilEffSat / oilGasSolventEffSat) + (rhoGas * gasEffSat / oilGasSolventEffSat) + (rhoSolvent * solventEffSat / oilGasSolventEffSat);
1187 Evaluation rhoGasEff = 0.0;
1188 if ( std::abs(muSolvent.value() - muGas.value()) < cutOff ) {
1189 rhoGasEff = ( (1.0 - tlMixParamRho) * rhoGas) + (tlMixParamRho * rhoMixSolventGasOil);
1191 const Evaluation solventGasEffFraction = (muGasPow * (muSolventPow - muGasEffPow)) / (muGasEffPow * (muSolventPow - muGasPow));
1192 rhoGasEff = (rhoGas * solventGasEffFraction) + (rhoSolvent * (1.0 - solventGasEffFraction));
1195 Evaluation rhoOilEff = 0.0;
1196 if ( std::abs(muGas.value() - muOil.value()) < cutOff ) {
1197 rhoOilEff = ( (1.0 - tlMixParamRho) * rhoOil) + (tlMixParamRho * rhoMixSolventGasOil);
1199 const Evaluation solventOilEffFraction = (muOilPow * (muOilEffPow - muSolventPow)) / (muOilEffPow * (muOilPow - muSolventPow));
1200 rhoOilEff = (rhoOil * solventOilEffFraction) + (rhoSolvent * (1.0 - solventOilEffFraction));
1203 Evaluation rhoSolventEff = 0.0;
1204 if ( std::abs( ( muSolventOilGasPow.value() - (muOilPow.value() * muGasPow.value()))) < cutOff ) {
1205 rhoSolventEff = ( (1.0 - tlMixParamRho) * rhoSolvent) + (tlMixParamRho * rhoMixSolventGasOil);
1207 const Evaluation sfraction_se = (muSolventOilGasPow - ( muOilPow * muGasPow * muSolventPow / muSolventEffPow) ) / ( muSolventOilGasPow - (muOilPow * muGasPow));
1208 rhoSolventEff = (rhoSolvent * sfraction_se) + (rhoGas * sgf * (1.0 - sfraction_se)) + (rhoOil * sof * (1.0 - sfraction_se));
1211 unsigned pvtRegionIdx = asImp_().pvtRegionIndex();
1213 const Evaluation bOilEff = rhoOilEff / (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs());
1214 const Evaluation bGasEff = rhoGasEff / (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) * fs.Rv());
1215 const Evaluation bSolventEff = rhoSolventEff / solventRefDensity();
1218 const auto& pmiscTable = SolventModule::pmisc(elemCtx, scvIdx, timeIdx);
1219 const Evaluation pmisc = pmiscTable.eval(po,
true);
1222 const Evaluation bo = fs.invB(oilPhaseIdx);
1223 const Evaluation bg = fs.invB(gasPhaseIdx);
1224 const Evaluation bs = solventInverseFormationVolumeFactor();
1227 fs.setInvB(oilPhaseIdx, pmisc * bOilEff + (1.0 - pmisc) * bo);
1228 fs.setInvB(gasPhaseIdx, pmisc * bGasEff + (1.0 - pmisc) * bg);
1229 solventInvFormationVolumeFactor_ = pmisc * bSolventEff + (1.0 - pmisc) * bs;
1232 fs.setDensity(oilPhaseIdx, fs.invB(oilPhaseIdx) * (FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) * fs.Rs()) );
1233 fs.setDensity(gasPhaseIdx, fs.invB(gasPhaseIdx) * (FluidSystem::referenceDensity(gasPhaseIdx, pvtRegionIdx) + FluidSystem::referenceDensity(oilPhaseIdx, pvtRegionIdx) * fs.Rv()) );
1234 solventDensity_ = solventInverseFormationVolumeFactor() * solventRefDensity();
1240 Evaluation& mobo = asImp_().mobility_[oilPhaseIdx];
1241 muOilEff = fs.invB(oilPhaseIdx) / (pmisc * bOilEff / muOilEff + (1.0 - pmisc) * bo / muOil);
1242 mobo *= muOil / muOilEff;
1244 Evaluation& mobg = asImp_().mobility_[gasPhaseIdx];
1245 muGasEff = fs.invB(gasPhaseIdx) / (pmisc * bGasEff / muGasEff + (1.0 - pmisc) * bg / muGas);
1246 mobg *= muGas / muGasEff;
1249 solventViscosity_ = solventInvFormationVolumeFactor_ / (pmisc * bSolventEff / muSolventEff + (1.0 - pmisc) * bs / muSolvent);
1253 Implementation& asImp_()
1254 {
return *
static_cast<Implementation*
>(
this); }
1256 Evaluation hydrocarbonSaturation_;
1257 Evaluation solventSaturation_;
1258 Evaluation solventDensity_;
1259 Evaluation solventViscosity_;
1260 Evaluation solventMobility_;
1261 Evaluation solventInvFormationVolumeFactor_;
1263 Scalar solventRefDensity_;
1266 template <
class TypeTag>
1269 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1270 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1271 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
1276 unsigned scvIdx OPM_UNUSED,
1277 unsigned timeIdx OPM_UNUSED)
1281 unsigned scvIdx OPM_UNUSED,
1282 unsigned timeIdx OPM_UNUSED)
1286 unsigned scvIdx OPM_UNUSED,
1287 unsigned timeIdx OPM_UNUSED)
1290 const Evaluation& solventSaturation()
const
1291 { OPM_THROW(std::runtime_error,
"solventSaturation() called but solvents are disabled"); }
1293 const Evaluation& solventDensity()
const
1294 { OPM_THROW(std::runtime_error,
"solventDensity() called but solvents are disabled"); }
1296 const Evaluation& solventViscosity()
const
1297 { OPM_THROW(std::runtime_error,
"solventViscosity() called but solvents are disabled"); }
1299 const Evaluation& solventMobility()
const
1300 { OPM_THROW(std::runtime_error,
"solventMobility() called but solvents are disabled"); }
1302 const Evaluation& solventInverseFormationVolumeFactor()
const
1303 { OPM_THROW(std::runtime_error,
"solventInverseFormationVolumeFactor() called but solvents are disabled"); }
1305 const Scalar& solventRefDensity()
const
1306 { OPM_THROW(std::runtime_error,
"solventRefDensity() called but solvents are disabled"); }
1316 template <
class TypeTag,
bool enableSolventV = GET_PROP_VALUE(TypeTag, EnableSolvent)>
1319 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) Implementation;
1321 typedef typename GET_PROP_TYPE(TypeTag, Scalar) Scalar;
1322 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1323 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1324 typedef typename GET_PROP_TYPE(TypeTag, IntensiveQuantities) IntensiveQuantities;
1325 typedef typename GET_PROP_TYPE(TypeTag, ExtensiveQuantities) ExtensiveQuantities;
1326 typedef typename GET_PROP_TYPE(TypeTag, FluidSystem) FluidSystem;
1327 typedef typename GET_PROP_TYPE(TypeTag, GridView) GridView;
1329 typedef Opm::MathToolbox<Evaluation> Toolbox;
1331 static constexpr
unsigned gasPhaseIdx = FluidSystem::gasPhaseIdx;
1332 static constexpr
int dimWorld = GridView::dimensionworld;
1334 typedef Dune::FieldVector<Scalar, dimWorld> DimVector;
1335 typedef Dune::FieldVector<Evaluation, dimWorld> DimEvalVector;
1342 template <
class Dummy =
bool>
1348 const auto& gradCalc = elemCtx.gradientCalculator();
1351 const auto& scvf = elemCtx.stencil(timeIdx).interiorFace(scvfIdx);
1352 const auto& faceNormal = scvf.normal();
1354 unsigned i = scvf.interiorIndex();
1355 unsigned j = scvf.exteriorIndex();
1358 DimEvalVector solventPGrad;
1360 gradCalc.calculateGradient(solventPGrad,
1364 Opm::Valgrind::CheckDefined(solventPGrad);
1370 const auto& gIn = elemCtx.problem().gravity(elemCtx, i, timeIdx);
1371 const auto& gEx = elemCtx.problem().gravity(elemCtx, j, timeIdx);
1373 const auto& intQuantsIn = elemCtx.intensiveQuantities(i, timeIdx);
1374 const auto& intQuantsEx = elemCtx.intensiveQuantities(j, timeIdx);
1376 const auto& posIn = elemCtx.pos(i, timeIdx);
1377 const auto& posEx = elemCtx.pos(j, timeIdx);
1378 const auto& posFace = scvf.integrationPos();
1381 DimVector distVecIn(posIn);
1382 DimVector distVecEx(posEx);
1383 DimVector distVecTotal(posEx);
1385 distVecIn -= posFace;
1386 distVecEx -= posFace;
1387 distVecTotal -= posIn;
1388 Scalar absDistTotalSquared = distVecTotal.two_norm2();
1391 auto rhoIn = intQuantsIn.solventDensity();
1392 auto pStatIn = - rhoIn*(gIn*distVecIn);
1396 Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1397 Scalar pStatEx = - rhoEx*(gEx*distVecEx);
1403 DimEvalVector f(distVecTotal);
1404 f *= (pStatEx - pStatIn)/absDistTotalSquared;
1407 for (
unsigned dimIdx = 0; dimIdx < dimWorld; ++dimIdx) {
1408 solventPGrad[dimIdx] += f[dimIdx];
1410 if (!Opm::isfinite(solventPGrad[dimIdx])) {
1411 OPM_THROW(Opm::NumericalProblem,
1412 "Non-finite potential gradient for solvent 'phase'");
1418 Evaluation solventPGradNormal = 0.0;
1419 for (
unsigned dimIdx = 0; dimIdx < faceNormal.size(); ++dimIdx)
1420 solventPGradNormal += solventPGrad[dimIdx]*faceNormal[dimIdx];
1422 if (solventPGradNormal > 0) {
1423 solventUpstreamDofIdx_ = j;
1424 solventDownstreamDofIdx_ = i;
1427 solventUpstreamDofIdx_ = i;
1428 solventDownstreamDofIdx_ = j;
1431 const auto& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1437 if (solventUpstreamDofIdx_ == i)
1438 solventVolumeFlux_ = solventPGradNormal*up.solventMobility();
1440 solventVolumeFlux_ = solventPGradNormal*Opm::scalarValue(up.solventMobility());
1447 template <
class Dummy =
bool>
1453 const ExtensiveQuantities& extQuants = asImp_();
1455 unsigned interiorDofIdx = extQuants.interiorIndex();
1456 unsigned exteriorDofIdx = extQuants.exteriorIndex();
1457 assert(interiorDofIdx != exteriorDofIdx);
1459 const auto& intQuantsIn = elemCtx.intensiveQuantities(interiorDofIdx, timeIdx);
1460 const auto& intQuantsEx = elemCtx.intensiveQuantities(exteriorDofIdx, timeIdx);
1462 unsigned I = elemCtx.globalSpaceIndex(interiorDofIdx, timeIdx);
1463 unsigned J = elemCtx.globalSpaceIndex(exteriorDofIdx, timeIdx);
1465 Scalar thpres = elemCtx.problem().thresholdPressure(I, J);
1466 Scalar trans = elemCtx.problem().transmissibility(elemCtx, interiorDofIdx, exteriorDofIdx);
1467 Scalar g = elemCtx.problem().gravity()[dimWorld - 1];
1469 Scalar zIn = elemCtx.problem().dofCenterDepth(elemCtx, interiorDofIdx, timeIdx);
1470 Scalar zEx = elemCtx.problem().dofCenterDepth(elemCtx, exteriorDofIdx, timeIdx);
1471 Scalar distZ = zIn - zEx;
1473 const Evaluation& rhoIn = intQuantsIn.solventDensity();
1474 Scalar rhoEx = Toolbox::value(intQuantsEx.solventDensity());
1475 const Evaluation& rhoAvg = rhoIn*0.5 + rhoEx*0.5;
1477 const Evaluation& pressureInterior = intQuantsIn.fluidState().pressure(gasPhaseIdx);
1478 Evaluation pressureExterior = Toolbox::value(intQuantsEx.fluidState().pressure(gasPhaseIdx));
1479 pressureExterior += distZ*g*rhoAvg;
1481 Evaluation pressureDiffSolvent = pressureExterior - pressureInterior;
1482 if (std::abs(Opm::scalarValue(pressureDiffSolvent)) > thpres) {
1483 if (pressureDiffSolvent < 0.0)
1484 pressureDiffSolvent += thpres;
1486 pressureDiffSolvent -= thpres;
1489 pressureDiffSolvent = 0.0;
1491 if (pressureDiffSolvent > 0.0) {
1492 solventUpstreamDofIdx_ = exteriorDofIdx;
1493 solventDownstreamDofIdx_ = interiorDofIdx;
1495 else if (pressureDiffSolvent < 0.0) {
1496 solventUpstreamDofIdx_ = interiorDofIdx;
1497 solventDownstreamDofIdx_ = exteriorDofIdx;
1503 solventUpstreamDofIdx_ = std::min(interiorDofIdx, exteriorDofIdx);
1504 solventDownstreamDofIdx_ = std::max(interiorDofIdx, exteriorDofIdx);
1505 solventVolumeFlux_ = 0.0;
1509 Scalar faceArea = elemCtx.stencil(timeIdx).interiorFace(scvfIdx).area();
1510 const IntensiveQuantities& up = elemCtx.intensiveQuantities(solventUpstreamDofIdx_, timeIdx);
1511 if (solventUpstreamDofIdx_ == interiorDofIdx)
1512 solventVolumeFlux_ =
1513 up.solventMobility()
1515 *pressureDiffSolvent;
1517 solventVolumeFlux_ =
1518 Opm::scalarValue(up.solventMobility())
1520 *pressureDiffSolvent;
1523 unsigned solventUpstreamIndex()
const
1524 {
return solventUpstreamDofIdx_; }
1526 unsigned solventDownstreamIndex()
const
1527 {
return solventDownstreamDofIdx_; }
1529 const Evaluation& solventVolumeFlux()
const
1530 {
return solventVolumeFlux_; }
1533 Implementation& asImp_()
1534 {
return *
static_cast<Implementation*
>(
this); }
1536 Evaluation solventVolumeFlux_;
1537 unsigned solventUpstreamDofIdx_;
1538 unsigned solventDownstreamDofIdx_;
1541 template <
class TypeTag>
1544 typedef typename GET_PROP_TYPE(TypeTag, ElementContext) ElementContext;
1545 typedef typename GET_PROP_TYPE(TypeTag, Evaluation) Evaluation;
1549 unsigned scvfIdx OPM_UNUSED,
1550 unsigned timeIdx OPM_UNUSED)
1554 unsigned scvfIdx OPM_UNUSED,
1555 unsigned timeIdx OPM_UNUSED)
1558 unsigned solventUpstreamIndex()
const
1559 { OPM_THROW(std::runtime_error,
"solventUpstreamIndex() called but solvents are disabled"); }
1561 unsigned solventDownstreamIndex()
const
1562 { OPM_THROW(std::runtime_error,
"solventDownstreamIndex() called but solvents are disabled"); }
1564 const Evaluation& solventVolumeFlux()
const
1565 { OPM_THROW(std::runtime_error,
"solventVolumeFlux() called but solvents are disabled"); }
static Scalar computeUpdateError(const PrimaryVariables &oldPv OPM_UNUSED, const EqVector &delta OPM_UNUSED)
Return how much a Newton-Raphson update is considered an error.
Definition: blackoilsolventmodules.hh:641
static void registerParameters()
Register all run-time parameters for the multi-phase VTK output module.
Definition: vtkblackoilsolventmodule.hh:97
void updateVolumeFluxTrans(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the volume flux of the polymer "phase" using the gas pressure potential diffe...
Definition: blackoilsolventmodules.hh:1449
Provides the volumetric quantities required for the equations needed by the solvents extension of the...
Definition: blackoilsolventmodules.hh:879
static void setSolventPvt(const SolventPvt &value)
Specify the solvent PVT of a all PVT regions.
Definition: blackoilsolventmodules.hh:478
#define GET_PROP_VALUE(TypeTag, PropTagName)
Access the value attribute of a property for a type tag.
Definition: propertysystem.hh:469
static void setMsfn(unsigned satRegionIdx, const TabulatedFunction &msfnKrsg, const TabulatedFunction &msfnKro)
Specify misicible relative permeability multipliers of a single region.
Definition: blackoilsolventmodules.hh:421
void solventPvtUpdate_(const ElementContext &elemCtx, unsigned scvIdx, unsigned timeIdx)
Update the intensive PVT properties needed to handle solvents from the primary variables.
Definition: blackoilsolventmodules.hh:1052
static void setSof2(unsigned satRegionIdx, const TabulatedFunction &sof2Krn)
Specify misicible hydrocabon relative permeability wrt water of a single region.
Definition: blackoilsolventmodules.hh:388
VTK output module for the black oil model's solvent related quantities.
Definition: vtkblackoilsolventmodule.hh:71
static void setMisc(unsigned miscRegionIdx, const TabulatedFunction &misc)
Misicibility function wrt solvent fraction of a single region.
Definition: blackoilsolventmodules.hh:399
static void setSorwmis(unsigned miscRegionIdx, const TabulatedFunction &sorwmis)
Misicibe residual oil saturation function wrt water saturation of a single region.
Definition: blackoilsolventmodules.hh:434
static void registerParameters()
Register all run-time parameters for the black-oil solvent module.
Definition: blackoilsolventmodules.hh:488
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: blackoilsolventmodules.hh:1317
static void assignPrimaryVars(PrimaryVariables &priVars, Scalar solventSaturation)
Assign the solvent specific primary variables to a PrimaryVariables object.
Definition: blackoilsolventmodules.hh:615
static void setSgcmis(unsigned miscRegionIdx, const TabulatedFunction &sgcwmis)
Misicibe critical gas saturation function wrt water saturation of a single region.
Definition: blackoilsolventmodules.hh:445
void setPhaseIndex(unsigned phaseIdx)
Set the index of the fluid phase for which the pressure should be returned.
Definition: quantitycallbacks.hh:108
This method contains all callback classes for quantities that are required by some extensive quantiti...
Contains the high level supplements required to extend the black oil model by solvents.
Definition: blackoilprimaryvariables.hh:47
Callback class for a phase pressure.
Definition: quantitycallbacks.hh:83
#define EWOMS_GET_PARAM(TypeTag, ParamType, ParamName)
Retrieve a runtime parameter.
Definition: parametersystem.hh:99
void updateVolumeFluxPerm(const ElementContext &elemCtx, unsigned scvfIdx, unsigned timeIdx)
Method which calculates the volume flux of the polymer "phase" using the pressure potential gradient ...
Definition: blackoilsolventmodules.hh:1344
static void setNumSatRegions(unsigned numRegions)
Specify the number of satuation regions.
Definition: blackoilsolventmodules.hh:364
static Scalar computeResidualError(const EqVector &resid)
Return how much a residual is considered an error.
Definition: blackoilsolventmodules.hh:653
static void setTlpmixpa(unsigned miscRegionIdx, const TabulatedFunction &tlPMixTable)
Todd-Longstaff mixing parameter multiplier wrt pressure of a single region.
Definition: blackoilsolventmodules.hh:469
static void setSsfn(unsigned satRegionIdx, const TabulatedFunction &ssfnKrg, const TabulatedFunction &ssfnKrs)
Specify the solvent saturation functions of a single region.
Definition: blackoilsolventmodules.hh:375
static void updatePrimaryVars(PrimaryVariables &newPv, const PrimaryVariables &oldPv, const EqVector &delta)
Do a Newton-Raphson update the primary variables of the solvents.
Definition: blackoilsolventmodules.hh:627
static void setTlmixpar(unsigned miscRegionIdx, const Scalar &tlMixParamViscosity, const Scalar &tlMixParamDensity)
Todd-Longstaff mixing parameters of a single region.
Definition: blackoilsolventmodules.hh:456
static void registerOutputModules(Model &model, Simulator &simulator)
Register all solvent specific VTK and ECL output modules.
Definition: blackoilsolventmodules.hh:500
void solventPreSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called before the saturation functions are doing their magic.
Definition: blackoilsolventmodules.hh:908
VTK output module for the black oil model's solvent related quantities.
static void setPmisc(unsigned miscRegionIdx, const TabulatedFunction &pmisc)
Misicibility function wrt pressure of a single region.
Definition: blackoilsolventmodules.hh:410
void solventPostSatFuncUpdate_(const ElementContext &elemCtx, unsigned dofIdx, unsigned timeIdx)
Called after the saturation functions have been doing their magic.
Definition: blackoilsolventmodules.hh:934