From 2f23dbc47af540fb3588bd7ec1ddea4a19af9a0f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Fri, 19 Dec 2025 12:36:05 +0100 Subject: [PATCH] Guard Against Missing Phases When Extracting Curves This commit adds certain checks for active phases before attempting to construct curve interpolants and retrieving curve points from those interpolants. Notably, we make sure that we don't return water-related curves (e.g, krw, krow, pcow) in a gas/oil system, nor gas-related curves in an oil/water system. Finally, add initial support for gas/water systems. The Pcgw capillary pressure curve is, however, not fully finalised. --- examples/extractPropCurves.cpp | 42 ++++++- opm/utility/ECLSaturationFunc.cpp | 193 ++++++++++++++++++++++++------ opm/utility/ECLSaturationFunc.hpp | 3 + 3 files changed, 201 insertions(+), 37 deletions(-) diff --git a/examples/extractPropCurves.cpp b/examples/extractPropCurves.cpp index ebe944b..57463cc 100644 --- a/examples/extractPropCurves.cpp +++ b/examples/extractPropCurves.cpp @@ -244,6 +244,32 @@ namespace { printGraph(std::cout, "Sg", curveName(curveSet, "Pcgo"), graph); } + void pcgw(const Opm::ECLSaturationFunc& sfunc, + const Opm::ECLInitFileData& init, + const std::string& gridID, + const int activeCell, + const Opm::ECLSaturationFunc::SatFuncScaling& scaling, + const CurveSet curveSet) + { + using RC = Opm::ECLSaturationFunc::RawCurve; + + auto func = std::vector{}; + func.reserve(1); + + // Request pcgw (gas/water capillary pressure (Pg-Pw) in G/W system) + func.push_back(RC{ + RC::Function::CapPress, + RC::SubSystem::GasWater, + Opm::ECLPhaseIndex::Vapour, + curveSet, + }); + + const auto graph = sfunc + .getSatFuncCurve(func, init, gridID, activeCell, scaling); + + printGraph(std::cout, "Sg", curveName(curveSet, "Pcgw"), graph); + } + void pcow(const Opm::ECLSaturationFunc& sfunc, const Opm::ECLInitFileData& init, const std::string& gridID, @@ -456,7 +482,8 @@ try { const auto gridID = prm.getDefault("gridName", std::string{}); const auto cellID = getActiveCell(init, rset, prm); - const auto pvtNum = init.keywordData("PVTNUM", gridID)[cellID]; + const auto pvtNum = init.haveKeywordData("PVTNUM") + ? init .keywordData("PVTNUM", gridID)[cellID] : 1; const auto haveHyst = init.haveHysteresis(); auto sfunc = Opm::ECLSaturationFunc(init); @@ -543,6 +570,19 @@ try { } } + if (prm.getDefault("pcgw", false)) { + pcgw(sfunc, init, gridID, cellID, scaling, CurveSet::Drainage); + } + + if (prm.getDefault("ipcgw", false)) { + if (haveHyst) { + pcgw(sfunc, init, gridID, cellID, scaling, CurveSet::Imbibition); + } + else { + std::cerr << "IPcgw not available in this result set\n"; + } + } + if (prm.getDefault("pcow", false)) { pcow(sfunc, init, gridID, cellID, scaling, CurveSet::Drainage); } diff --git a/opm/utility/ECLSaturationFunc.cpp b/opm/utility/ECLSaturationFunc.cpp index 26eb370..ccd96f2 100644 --- a/opm/utility/ECLSaturationFunc.cpp +++ b/opm/utility/ECLSaturationFunc.cpp @@ -152,6 +152,21 @@ namespace { return (scaling.enable & T::Vertical) != 0; } + bool hasTabulatedGasSaturationFunctions(const std::vector& tabdims) + { + return tabdims[TABDIMS_NSSGFN_ITEM] * tabdims[TABDIMS_NTSGFN_ITEM] > 0; + } + + bool hasTabulatedOilSaturationFunctions(const std::vector& tabdims) + { + return tabdims[TABDIMS_NSSOFN_ITEM] * tabdims[TABDIMS_NTSOFN_ITEM] > 0; + } + + bool hasTabulatedWaterSaturationFunctions(const std::vector& tabdims) + { + return tabdims[TABDIMS_NSSWFN_ITEM] * tabdims[TABDIMS_NTSWFN_ITEM] > 0; + } + int saturationRegion(const ::Opm::ECLInitFileData& init, const ::Opm::ECLSaturationFunc::RawCurve& curve, const std::string& gridID, @@ -234,8 +249,8 @@ namespace { } std::vector - pcgo(const std::size_t regID, - const std::vector& sg) const + pcg(const std::size_t regID, + const std::vector& sg) const { const auto t = this->table(regID); const auto c = this->pccol(); @@ -1353,7 +1368,7 @@ class Opm::ECLSaturationFunc::Impl vertFuncVal(G, init, ep, opt, [&host] (const int regID, const double sat) -> double { - return host.gas_->pcgo(regID, { sat })[0]; + return host.gas_->pcg(regID, { sat })[0]; }); eps.vertscaling = Create::Vertical:: @@ -1508,6 +1523,24 @@ class Opm::ECLSaturationFunc::Impl const ECLRestartData& rstrt, const SatFuncScaling& scaling = SatFuncScaling{}) const; + FlowDiagnostics::Graph + gasOilCurve(const RawCurve& func, + const SatFunc::CreateEPS::CurveSet curveSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const; + + FlowDiagnostics::Graph + gasWaterCurve(const RawCurve& func, + const SatFunc::CreateEPS::CurveSet curveSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const; + FlowDiagnostics::Graph krgCurve(const SatFunc::CreateEPS::CurveSet curveSet, const ECLInitFileData& init, @@ -1524,6 +1557,22 @@ class Opm::ECLSaturationFunc::Impl const int satnum, const SatFuncScaling& scaling) const; + FlowDiagnostics::Graph + pcgwCurve(const SatFunc::CreateEPS::CurveSet curveSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const; + + FlowDiagnostics::Graph + pcgCurve(const SatFunc::CreateEPS::CurveSet curveSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const; + std::vector krw(const ECLGraph& G, const ECLRestartData& rstrt, @@ -1718,31 +1767,31 @@ Impl::initRelPermInterp(const EPSEvaluator::ActPh& active, const auto& tabdims = init.keywordData("TABDIMS"); const auto& tab = init.keywordData("TAB"); - if (active.gas) { - this->gas_.reset(new Gas::SatFunction(tabdims, tab, usys)); + if (active.gas && hasTabulatedGasSaturationFunctions(tabdims)) { + this->gas_ = std::make_unique(tabdims, tab, usys); this->setNumTables(this->gas_->numTables()); } - if (active.wat) { - this->wat_.reset(new Water::SatFunction(tabdims, tab, usys)); + if (active.wat && hasTabulatedWaterSaturationFunctions(tabdims)) { + this->wat_ = std::make_unique(tabdims, tab, usys); this->setNumTables(this->wat_->numTables()); } - if (active.oil) { + if (active.oil && hasTabulatedOilSaturationFunctions(tabdims)) { if (! isThreePh) { using KrModel = Oil::TwoPhase; - if (active.gas) { + if (this->gas_ != nullptr) { const auto subsys = KrModel::SubSys::OilGas; - this->oil_.reset(new KrModel(subsys, tabdims, tab)); + this->oil_ = std::make_unique(subsys, tabdims, tab); } - else if (active.wat) { + else if (this->wat_ != nullptr) { const auto subsys = KrModel::SubSys::OilWater; - this->oil_.reset(new KrModel(subsys, tabdims, tab)); + this->oil_ = std::make_unique(subsys, tabdims, tab); } else { throw std::invalid_argument { @@ -1751,13 +1800,15 @@ Impl::initRelPermInterp(const EPSEvaluator::ActPh& active, } } - if (isThreePh) { + if (isThreePh && (this->wat_ != nullptr)) { using KrModel = Oil::ECLStdThreePhase; - this->oil_.reset(new KrModel(tabdims, tab, this->wat_->swco())); + this->oil_ = std::make_unique(tabdims, tab, this->wat_->swco()); } - this->setNumTables(this->oil_->numTables()); + if (this->oil_ != nullptr) { + this->setNumTables(this->oil_->numTables()); + } } } @@ -1853,11 +1904,8 @@ getSatFuncCurve(const std::vector& func, switch (fi.thisPh) { case ECLPhaseIndex::Aqua: { - if (! (this->wat_ && (fi.subsys == RawCurve::SubSystem::OilWater))) - { - // Water is not an active phase, or we're being asked to - // extract a saturation function curve for water in the O/G - // subsystem. No such curve. + if (!this->wat_) { + // Water is not an active phase. No such curve. graph.emplace_back(); continue; } @@ -1878,6 +1926,15 @@ getSatFuncCurve(const std::vector& func, continue; } + if (((fi.subsys == RawCurve::SubSystem::OilGas) && !this->gas_) || + ((fi.subsys == RawCurve::SubSystem::OilWater) && !this->wat_)) + { + // Caller requests relative permeability for oil in an inactive + // two-phase system. No such curve. + graph.emplace_back(); + continue; + } + const auto kro = this->kroCurve (fi.subsys, crvSet, init, gridID, activeCell, satnum, scaling); @@ -1893,20 +1950,21 @@ getSatFuncCurve(const std::vector& func, case ECLPhaseIndex::Vapour: { - if (! (this->gas_ && (fi.subsys == RawCurve::SubSystem::OilGas))) - { - // Gas is not an active phase, or we're being asked to - // extract a saturation function curve for gas in the O/W - // subsystem. No such curve. + if (!this->gas_) { graph.emplace_back(); continue; } - auto crv = (fi.curve == RawCurve::Function::RelPerm) - ? this->krgCurve (crvSet, init, gridID, activeCell, satnum, scaling) - : this->pcgoCurve(crvSet, init, gridID, activeCell, satnum, scaling); + if (fi.subsys == RawCurve::SubSystem::OilGas) { + graph.push_back(this->gasOilCurve(fi, crvSet, init, gridID, + activeCell, satnum, scaling)); + } - graph.push_back(std::move(crv)); + if (fi.subsys == RawCurve::SubSystem::GasWater) { + graph.push_back(this->gasWaterCurve(fi, crvSet, init, gridID, + activeCell, satnum, scaling)); + + } } break; @@ -2127,6 +2185,36 @@ krg(const ECLGraph& G, return kr; } +Opm::FlowDiagnostics::Graph +Opm::ECLSaturationFunc::Impl:: +gasOilCurve(const RawCurve& func, + const SatFunc::CreateEPS::CurveSet crvSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const +{ + return (func.curve == RawCurve::Function::RelPerm) + ? this->krgCurve(crvSet, init, gridID, activeCell, satnum, scaling) + : this->pcgoCurve(crvSet, init, gridID, activeCell, satnum, scaling); +} + +Opm::FlowDiagnostics::Graph +Opm::ECLSaturationFunc::Impl:: +gasWaterCurve(const RawCurve& func, + const SatFunc::CreateEPS::CurveSet crvSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const +{ + return (func.curve == RawCurve::Function::RelPerm) + ? this->krgCurve(crvSet, init, gridID, activeCell, satnum, scaling) + : this->pcgwCurve(crvSet, init, gridID, activeCell, satnum, scaling); +} + Opm::FlowDiagnostics::Graph Opm::ECLSaturationFunc::Impl:: krgCurve(const SatFunc::CreateEPS::CurveSet curveSet, @@ -2231,10 +2319,43 @@ pcgoCurve(const SatFunc::CreateEPS::CurveSet curveSet, const int satnum, const SatFuncScaling& scaling) const { - if (! this->gas_) { + if ((this->gas_ == nullptr) || + ((this->oil_ == nullptr) && (this->wat_ != nullptr))) + { + // No Pcgo curve if gas not active or oil not active but water is. return FlowDiagnostics::Graph{}; } + return this->pcgCurve(curveSet, init, gridID, activeCell, satnum, scaling); +} + +Opm::FlowDiagnostics::Graph +Opm::ECLSaturationFunc::Impl:: +pcgwCurve(const SatFunc::CreateEPS::CurveSet curveSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const +{ + if (!this->gas_ || !this->wat_ || this->oil_) { + // G/W system not applicable if oil active or if one or + // both of the gas and water phases are *not* active. + return FlowDiagnostics::Graph{}; + } + + return this->pcgCurve(curveSet, init, gridID, activeCell, satnum, scaling); +} + +Opm::FlowDiagnostics::Graph +Opm::ECLSaturationFunc::Impl:: +pcgCurve(const SatFunc::CreateEPS::CurveSet curveSet, + const ECLInitFileData& init, + const std::string& gridID, + const std::size_t activeCell, + const int satnum, + const SatFuncScaling& scaling) const +{ auto opt = SatFunc::CreateEPS::EPSOptions{}; opt.use3PtScaling = this->use3PtScaling_; @@ -2284,14 +2405,14 @@ pcgoCurve(const SatFunc::CreateEPS::CurveSet curveSet, const auto nPoints = abscissas.size(); // Subtract one from SATNUM to create valid table index. - auto pc = this->gas_->pcgo(satnum - 1, sg_lookup); + auto pc = this->gas_->pcg(satnum - 1, sg_lookup); if (enableVerticalSFScaling(scaling)) { const auto fvals = SatFunc::CreateEPS::Vertical:: unscaledFunctionValues(init, gridID, satnum, *this->tep_, opt, [this](const int satreg, const double s) { - return this->gas_->pcgo(satreg, std::vector{s}).front(); + return this->gas_->pcg(satreg, std::vector{s}).front(); }); auto eps = SatFunc::CreateEPS::Vertical:: @@ -2467,7 +2588,7 @@ pcowCurve(const SatFunc::CreateEPS::CurveSet curveSet, const int satnum, const SatFuncScaling& scaling) const { - if (! this->wat_) { + if (!this->wat_ || !this->oil_) { return FlowDiagnostics::Graph{}; } @@ -2654,7 +2775,7 @@ extractRawTableEndPoints(const EPSEvaluator::ActPh& active) this->tep_ = std::make_unique(); - if (active.oil) { + if (active.oil && (this->oil_ != nullptr)) { this->tep_->conn.oil = this->oil_->soco(); this->tep_->crit.oil_in_gas = this->oil_->sogcr(); this->tep_->crit.oil_in_water = this->oil_->sowcr(); @@ -2667,7 +2788,7 @@ extractRawTableEndPoints(const EPSEvaluator::ActPh& active) this->tep_->smax.oil = zero; } - if (active.gas) { + if (active.gas && (this->gas_ != nullptr)) { this->tep_->conn.gas = this->gas_->sgco(); this->tep_->crit.gas = this->gas_->sgcr(); this->tep_->smax.gas = this->gas_->sgmax(); @@ -2678,7 +2799,7 @@ extractRawTableEndPoints(const EPSEvaluator::ActPh& active) this->tep_->smax.gas = zero; } - if (active.wat) { + if (active.wat && (this->wat_ != nullptr)) { this->tep_->conn.water = this->wat_->swco(); this->tep_->crit.water = this->wat_->swcr(); this->tep_->smax.water = this->wat_->swmax(); diff --git a/opm/utility/ECLSaturationFunc.hpp b/opm/utility/ECLSaturationFunc.hpp index c30b8a5..599bdf9 100644 --- a/opm/utility/ECLSaturationFunc.hpp +++ b/opm/utility/ECLSaturationFunc.hpp @@ -94,6 +94,9 @@ namespace Opm { /// Oil-Water subsystem OilWater, + + /// Gas-Water sybsystem + GasWater, }; /// Which process does this request reference.