From 1e64409f28f1ebd3138803d2a225685fdcc505de Mon Sep 17 00:00:00 2001 From: Tor Harald Sandve Date: Wed, 9 Dec 2015 11:30:16 +0100 Subject: [PATCH] Scaling of relative permeability endpoints by the miscibility function --- opm/autodiff/BlackoilSolventModel_impl.hpp | 42 ++++---- opm/autodiff/SolventPropsAdFromDeck.cpp | 109 +++++++++++++++++---- opm/autodiff/SolventPropsAdFromDeck.hpp | 16 +++ 3 files changed, 125 insertions(+), 42 deletions(-) diff --git a/opm/autodiff/BlackoilSolventModel_impl.hpp b/opm/autodiff/BlackoilSolventModel_impl.hpp index 8a212c0f3..479e7fe5d 100644 --- a/opm/autodiff/BlackoilSolventModel_impl.hpp +++ b/opm/autodiff/BlackoilSolventModel_impl.hpp @@ -533,52 +533,48 @@ namespace Opm { ? state.saturation[ pu.phase_pos[ Gas ] ] : zero); + Selector zero_selector(ss.value() + sg.value(), Selector::Zero); - ADB F_solvent = zero_selector.select(zero, ss / (ss + sg)); + const ADB F_solvent = zero_selector.select(zero, ss / (ss + sg)); if (is_miscible_ && canonicalPhaseIdx != Water ) { assert(active_[ Oil ]); assert(active_[ Gas ]); const ADB& so = state.saturation[ pu.phase_pos[ Oil ] ]; + const ADB& sw = (active_[ Water ] + ? state.saturation[ pu.phase_pos[ Water ] ] + : zero); const ADB misc = solvent_props_.miscibilityFunction(F_solvent, cells_); + const ADB sn = ss + so + sg; // adjust endpoints - //const int np = fluid_.numPhases(); - //V smin = V::Zero(np * nc); - //V smax = V::Constant(np*nc, 1.0); - //fluid_.getSaturationEndpoints(cells_, smin, smax); - //ADB sor = subset(smin, Span(nc, np, Oil)) * misc; //+ (ones - misc) * sorwmis; - //ADB sgc = subset(smin, Span(nc, np, Gas)) * misc; //+ (ones - misc) * sgcwmis; - ADB sor = V::Constant(nc, 0.0) * misc; - ADB sgc = V::Constant(nc, 0.0) * misc; + const V sgcr = fluid_.scaledCriticalGasSaturations(cells_); + const V sogcr = fluid_.scaledCriticalOilinGasSaturations(cells_); + const ADB sorwmis = solvent_props_.miscibleResidualOilSaturationFunction(sw, cells_); + const ADB sgcwmis = solvent_props_.miscibleCriticalGasSaturationFunction(sw, cells_); - //std::cout << sor.value().maxCoeff() << std::endl; - //std::cout << sgc.value().maxCoeff() << std::endl; + ADB sor = misc * sorwmis + (ones - misc) * sogcr; + ADB sgc = misc * sgcwmis + (ones - misc) * sgcr; - ADB sn_eff = sn - sor - sgc; - ADB ssg = ss + sg - sgc; + const ADB ssg = ss + sg - sgc; + const ADB sn_eff = sn - sor - sgc; - Selector negative_selector(ssg.value(), Selector::LessZero); - ssg = negative_selector.select(zero,ssg); + //std::cout << sn.value().minCoeff() << " " << sn.value().maxCoeff() < zeroSn_selector(sn_eff.value(), Selector::LessEqualZero); + Selector zeroSn_selector(sn_eff.value(), Selector::Zero); const ADB F_totalGas = zeroSn_selector.select( zero, ssg / sn_eff); kr_mod = (ones - misc) * kr_mod; if (canonicalPhaseIdx == Gas) { - const ADB mkrgt = F_totalGas * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_); + const ADB mkrgt = solvent_props_.miscibleSolventGasRelPermMultiplier(F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_); kr_mod += misc * mkrgt; } if (canonicalPhaseIdx == Oil) { - const ADB mkro = (ones - F_totalGas) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_); + const ADB mkro = solvent_props_.miscibleOilRelPermMultiplier(ones - F_totalGas, cells_) * solvent_props_.misicibleHydrocarbonWaterRelPerm(sn, cells_); kr_mod += misc * mkro; } diff --git a/opm/autodiff/SolventPropsAdFromDeck.cpp b/opm/autodiff/SolventPropsAdFromDeck.cpp index 4155c5e47..470bf3bc2 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.cpp +++ b/opm/autodiff/SolventPropsAdFromDeck.cpp @@ -141,6 +141,11 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, const std::vector& sn = sof2Table.getSoColumn(); const std::vector& krn = sof2Table.getKroColumn(); + for (size_t i = 0; i < sn.size(); ++i) { + std::cout << sn[i] << " " << krn[i] <(sn, krn); } @@ -148,6 +153,31 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, OPM_THROW(std::runtime_error, "SOF2 must be specified in MISCIBLE (SOLVENT) runs\n"); } + const TableContainer& miscTables = tables->getMiscTables(); + if (!miscTables.empty()) { + + int numRegions = miscTables.size(); + + if(numRegions > 1) { + OPM_THROW(std::runtime_error, "Only single table miscibility function supported for MISC"); + } + // resize the attributes of the object + misc_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::MiscTable& miscTable = miscTables.getTable(regionIdx); + + // Copy data + // solventFraction = Ss / (Ss + Sg); + const std::vector& solventFraction = miscTable.getSolventFractionColumn(); + const std::vector& misc = miscTable.getMiscibilityColumn(); + + misc_[regionIdx] = NonuniformTableLinear(solventFraction, misc); + + } + } else { + OPM_THROW(std::runtime_error, "MISC must be specified in MISCIBLE (SOLVENT) runs\n"); + } + // miscible relative permeability multipleiers const TableContainer& msfnTables = tables->getMsfnTables(); if (!msfnTables.empty()) { @@ -173,37 +203,49 @@ SolventPropsAdFromDeck::SolventPropsAdFromDeck(DeckConstPtr deck, mkro_[regionIdx] = NonuniformTableLinear(Ssg, kro); } - - } else { - OPM_THROW(std::runtime_error, "MSFN must be specified in MISCIBLE (SOLVENT) runs\n"); } - const TableContainer& miscTables = tables->getMiscTables(); - if (!miscTables.empty()) { + const TableContainer& sorwmisTables = tables->getSorwmisTables(); + if (!sorwmisTables.empty()) { - int numRegions = miscTables.size(); + int numRegions = sorwmisTables.size(); if(numRegions > 1) { - OPM_THROW(std::runtime_error, "Only single table miscibility function supported for MISC"); + OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SORWMIS"); } // resize the attributes of the object - misc_.resize(numRegions); + sorwmis_.resize(numRegions); for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { - const Opm::MiscTable& miscTable = miscTables.getTable(regionIdx); + const Opm::SorwmisTable& sorwmisTable = sorwmisTables.getTable(regionIdx); // Copy data - // solventFraction = Ss / (Ss + Sg); - const std::vector& solventFraction = miscTable.getSolventFractionColumn(); - const std::vector& misc = miscTable.getMiscibilityColumn(); - - misc_[regionIdx] = NonuniformTableLinear(solventFraction, misc); + const std::vector& sw = sorwmisTable.getWaterSaturationColumn(); + const std::vector& sorwmis = sorwmisTable.getMiscibleResidualOilColumn(); + sorwmis_[regionIdx] = NonuniformTableLinear(sw, sorwmis); } - - } else { - OPM_THROW(std::runtime_error, "MSFN must be specified in MISCIBLE (SOLVENT) runs\n"); } + const TableContainer& sgcwmisTables = tables->getSgcwmisTables(); + if (!sgcwmisTables.empty()) { + + int numRegions = sgcwmisTables.size(); + + if(numRegions > 1) { + OPM_THROW(std::runtime_error, "Only single table miscibility function supported for SGCWMIS"); + } + // resize the attributes of the object + sgcwmis_.resize(numRegions); + for (int regionIdx = 0; regionIdx < numRegions; ++regionIdx) { + const Opm::SgcwmisTable& sgcwmisTable = sgcwmisTables.getTable(regionIdx); + + // Copy data + const std::vector& sw = sgcwmisTable.getWaterSaturationColumn(); + const std::vector& sgcwmis = sgcwmisTable.getMiscibleResidualGasColumn(); + + sgcwmis_[regionIdx] = NonuniformTableLinear(sw, sgcwmis); + } + } } } @@ -266,21 +308,50 @@ ADB SolventPropsAdFromDeck::misicibleHydrocarbonWaterRelPerm(const ADB& Sn, ADB SolventPropsAdFromDeck::miscibleSolventGasRelPermMultiplier(const ADB& Ssg, const Cells& cells) const { - return SolventPropsAdFromDeck::makeAD(Ssg, cells, mkrsg_); + if (mkrsg_.size() > 0) { + return SolventPropsAdFromDeck::makeAD(Ssg, cells, mkrsg_); + } + // trivial function if not specified + return Ssg; } ADB SolventPropsAdFromDeck::miscibleOilRelPermMultiplier(const ADB& So, const Cells& cells) const { - return SolventPropsAdFromDeck::makeAD(So, cells, mkro_); + if (mkro_.size() > 0) { + return SolventPropsAdFromDeck::makeAD(So, cells, mkro_); + } + // trivial function if not specified + return So; } ADB SolventPropsAdFromDeck::miscibilityFunction(const ADB& solventFraction, const Cells& cells) const { + return SolventPropsAdFromDeck::makeAD(solventFraction, cells, misc_); } + +ADB SolventPropsAdFromDeck::miscibleCriticalGasSaturationFunction (const ADB& Sw, + const Cells& cells) const { + if (sgcwmis_.size()>0) { + return SolventPropsAdFromDeck::makeAD(Sw, cells, sgcwmis_); + } + // return zeros if not specified + return ADB::constant(V::Zero(Sw.size())); +} + + +ADB SolventPropsAdFromDeck::miscibleResidualOilSaturationFunction (const ADB& Sw, + const Cells& cells) const { + if (sorwmis_.size()>0) { + return SolventPropsAdFromDeck::makeAD(Sw, cells, sorwmis_); + } + // return zeros if not specified + return ADB::constant(V::Zero(Sw.size())); +} + ADB SolventPropsAdFromDeck::makeAD(const ADB& X_AD, const Cells& cells, std::vector> table) const { const int n = cells.size(); assert(Sn.value().size() == n); diff --git a/opm/autodiff/SolventPropsAdFromDeck.hpp b/opm/autodiff/SolventPropsAdFromDeck.hpp index ef409ffff..6fb349d9e 100644 --- a/opm/autodiff/SolventPropsAdFromDeck.hpp +++ b/opm/autodiff/SolventPropsAdFromDeck.hpp @@ -108,6 +108,20 @@ public: ADB miscibilityFunction (const ADB& solventFraction, const Cells& cells) const; + /// Miscible critical gas saturation function + /// \param[in] Sw Array of n water saturation values. + /// \param[in] cells Array of n cell indices to be associated with the saturation values. + /// \return Array of n miscible critical gas saturation values + ADB miscibleCriticalGasSaturationFunction (const ADB& Sw, + const Cells& cells) const; + + /// Miscible residual oil saturation function + /// \param[in] Sw Array of n water saturation values. + /// \param[in] cells Array of n cell indices to be associated with the saturation values. + /// \return Array of n miscible residual oil saturation values + ADB miscibleResidualOilSaturationFunction (const ADB& Sw, + const Cells& cells) const; + /// Solvent surface density /// \param[in] cells Array of n cell indices to be associated with the fraction values. /// \return Array of n solvent density values. @@ -129,6 +143,8 @@ private: std::vector > mkro_; std::vector > mkrsg_; std::vector > misc_; + std::vector > sorwmis_; + std::vector > sgcwmis_; }; } // namespace OPM