diff --git a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp index 8d5fa86453..4bb25dc5aa 100644 --- a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp +++ b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp @@ -707,12 +707,12 @@ std::vector RigFlowDiagSolverInterface::ca const Opm::ECLPhaseIndex queryPhaseIndex = queryPhaseArr[i]; const PvtCurve::Phase mapToPhase = mapToPhaseArr[i]; - std::vector graphArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getPvtCurve(rawCurveType, queryPhaseIndex, static_cast(activeCellIndex)); - for (Opm::FlowDiagnostics::Graph srcGraph : graphArr) + std::vector graphArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getPvtCurve(rawCurveType, queryPhaseIndex, static_cast(activeCellIndex)); + for (Opm::ECLPVT::PVTGraph srcGraph : graphArr) { - if (srcGraph.first.size() > 0) + if (srcGraph.press.size() > 0) { - retCurveArr.push_back({ mapToPhase, srcGraph.first, srcGraph.second }); + retCurveArr.push_back({ mapToPhase, srcGraph.press, srcGraph.value}); } } } diff --git a/ResInsightVersion.cmake b/ResInsightVersion.cmake index 8d5c10e930..cf0243eb50 100644 --- a/ResInsightVersion.cmake +++ b/ResInsightVersion.cmake @@ -14,7 +14,7 @@ set(ERT_GITHUB_SHA "2e36798b43daf18c112b91aa3febbf2fccd4a95f") set(OPM_FLOWDIAGNOSTICS_SHA "7e2be931d430796ed42efcfb5c6b67a8d5962f7f") # https://github.com/OPM/opm-flowdiagnostics-applications -set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "cf8a5d2381776a8887d61f825b40f03b1d1cc4b0") +set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "a773bcfc963705679ecc28f58048fc38936fbbc6") # https://github.com/OPM/opm-parser/blob/master/opm/parser/eclipse/Units/Units.hpp # This file was moved from opm-core to opm-parser october 2016 diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractPropCurves.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractPropCurves.cpp index 737819ee87..615aef4eb0 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractPropCurves.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractPropCurves.cpp @@ -22,6 +22,7 @@ #include #include +#include #include #include #include @@ -63,6 +64,34 @@ namespace { os.precision(oprec); } + template + void printGraph(OStream& os, + const std::string& name, + const std::vector& graphs) + { + const auto oprec = os.precision(16); + const auto oflags = os.setf(std::ios_base::scientific); + + auto k = 1; + for (const auto& graph : graphs) { + const auto& p = graph.press; + const auto& R = graph.mixRat; + const auto& f = graph.value; + + os << name << '{' << k << "} = [\n"; + + for (auto n = p.size(), i = 0*n; i < n; ++i) { + os << p[i] << ' ' << R[i] << ' ' << f[i] << '\n'; + } + + os << "];\n\n"; + k += 1; + } + + os.setf(oflags); + os.precision(oprec); + } + // ----------------------------------------------------------------- // Relative permeability @@ -274,6 +303,38 @@ namespace { printGraph(std::cout, "rsSat", graph); } + + // ----------------------------------------------------------------- + + std::unique_ptr + makeUnits(const std::string& unit, + const Opm::ECLInitFileData& init) + { + if ((unit == "si") || (unit == "SI") || (unit == "internal")) { + return {}; // No conversion needed. + } + + if ((unit == "metric") || (unit == "Metric") || (unit == "METRIC")) { + return Opm::ECLUnits::metricUnitConventions(); + } + + if ((unit == "field") || (unit == "Field") || (unit == "FIELD")) { + return Opm::ECLUnits::fieldUnitConventions(); + } + + if ((unit == "lab") || (unit == "Lab") || (unit == "LAB")) { + return Opm::ECLUnits::labUnitConventions(); + } + + if ((unit == "pvt-m") || (unit == "PVT-M") || (unit == "PVTM")) { + return Opm::ECLUnits::pvtmUnitConventions(); + } + + std::cerr << "Unit convention '" << unit << "' not recognized\n" + << "Using 'native' (input/serialised) conventions.\n"; + + return Opm::ECLUnits::serialisedUnitConventions(init); + } } // namespace Anonymous int main(int argc, char* argv[]) @@ -285,8 +346,16 @@ try { const auto rset = example::identifyResultSet(prm); const auto init = Opm::ECLInitFileData(rset.initFile()); const auto graph = Opm::ECLGraph::load(rset.gridFile(), init); - const auto sfunc = Opm::ECLSaturationFunc(graph, init, useEPS); - const auto pvtCC = Opm::ECLPVT::ECLPvtCurveCollection(graph, init); + + auto sfunc = Opm::ECLSaturationFunc(graph, init, useEPS); + auto pvtCC = Opm::ECLPVT::ECLPvtCurveCollection(graph, init); + + if (prm.has("unit")) { + auto units = makeUnits(prm.get("unit"), init); + + sfunc.setOutputUnits(units->clone()); + pvtCC.setOutputUnits(std::move(units)); + } // ----------------------------------------------------------------- // Relative permeability diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp index a8d3fcd306..bbe24a59c5 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.cpp @@ -29,6 +29,9 @@ #include #include +#include +#include +#include #include #include #include @@ -76,6 +79,21 @@ namespace { return tep; } + + bool validSaturation(const double s) + { + return (! (s < 0.0)) && (! (s > 1.0)); + } + + bool validSaturations(std::initializer_list sats) + { + return std::accumulate(std::begin(sats), + std::end (sats), true, + [](const bool result, const double s) -> bool + { + return result && validSaturation(s); + }); + } } // --------------------------------------------------------------------- @@ -85,10 +103,12 @@ namespace { class Opm::SatFunc::TwoPointScaling::Impl { public: - Impl(std::vector smin, - std::vector smax) - : smin_(std::move(smin)) - , smax_(std::move(smax)) + Impl(std::vector smin, + std::vector smax, + InvalidEndpointBehaviour handle_invalid) + : smin_ (std::move(smin)) + , smax_ (std::move(smax)) + , handle_invalid_(handle_invalid) { if (this->smin_.size() != this->smax_.size()) { throw std::invalid_argument { @@ -109,6 +129,11 @@ public: private: std::vector smin_; std::vector smax_; + + InvalidEndpointBehaviour handle_invalid_; + + void handleInvalidEndpoint(const SaturationAssoc& sp, + std::vector& effsat) const; }; std::vector @@ -127,6 +152,12 @@ Impl::eval(const TableEndPoints& tep, const auto sLO = this->smin_[cell]; const auto sHI = this->smax_[cell]; + if (! validSaturations({ sLO, sHI })) { + this->handleInvalidEndpoint(eval_pt, effsat); + + continue; + } + effsat.push_back(0.0); auto& s_eff = effsat.back(); @@ -192,6 +223,22 @@ Impl::reverse(const TableEndPoints& tep, return unscaledsat; } +void +Opm::SatFunc::TwoPointScaling::Impl:: +handleInvalidEndpoint(const SaturationAssoc& sp, + std::vector& effsat) const +{ + if (this->handle_invalid_ == InvalidEndpointBehaviour::UseUnscaled) { + // User requests that invalid scaling be treated as unscaled + // saturations. Pick that. + effsat.push_back(sp.sat); + return; + } + + // Nothing to do for IgnorePoint. In particular, we must not change the + // contents or the size of effsat. +} + // --------------------------------------------------------------------- // Class Opm::ThreePointScaling::Impl // --------------------------------------------------------------------- @@ -199,12 +246,14 @@ Impl::reverse(const TableEndPoints& tep, class Opm::SatFunc::ThreePointScaling::Impl { public: - Impl(std::vector smin , - std::vector sdisp, - std::vector smax ) - : smin_ (std::move(smin )) - , sdisp_(std::move(sdisp)) - , smax_ (std::move(smax )) + Impl(std::vector smin , + std::vector sdisp, + std::vector smax , + InvalidEndpointBehaviour handle_invalid) + : smin_ (std::move(smin )) + , sdisp_ (std::move(sdisp)) + , smax_ (std::move(smax )) + , handle_invalid_(handle_invalid) { if ((this->sdisp_.size() != this->smin_.size()) || (this->sdisp_.size() != this->smax_.size())) @@ -228,6 +277,11 @@ private: std::vector smin_; std::vector sdisp_; std::vector smax_; + + InvalidEndpointBehaviour handle_invalid_; + + void handleInvalidEndpoint(const SaturationAssoc& sp, + std::vector& effsat) const; }; std::vector @@ -241,13 +295,19 @@ Impl::eval(const TableEndPoints& tep, for (const auto& eval_pt : sp) { const auto cell = eval_pt.cell; - effsat.push_back(0.0); - auto& s_eff = effsat.back(); - const auto sLO = this->smin_ [cell]; const auto sR = this->sdisp_[cell]; const auto sHI = this->smax_ [cell]; + if (! validSaturations({ sLO, sR, sHI })) { + this->handleInvalidEndpoint(eval_pt, effsat); + + continue; + } + + effsat.push_back(0.0); + auto& s_eff = effsat.back(); + if (! (eval_pt.sat > sLO)) { // s <= sLO s_eff = tep.low; @@ -324,6 +384,22 @@ Impl::reverse(const TableEndPoints& tep, return unscaledsat; } +void +Opm::SatFunc::ThreePointScaling::Impl:: +handleInvalidEndpoint(const SaturationAssoc& sp, + std::vector& effsat) const +{ + if (this->handle_invalid_ == InvalidEndpointBehaviour::UseUnscaled) { + // User requests that invalid scaling be treated as unscaled + // saturations. Pick that. + effsat.push_back(sp.sat); + return; + } + + // Nothing to do for IgnorePoint. In particular, we must not change the + // contents or the size of effsat. +} + // --------------------------------------------------------------------- // EPS factory functions for two-point and three-point scaling options // --------------------------------------------------------------------- @@ -332,6 +408,7 @@ namespace Create { using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions; using RTEP = ::Opm::SatFunc::CreateEPS::RawTableEndPoints; using TEP = ::Opm::SatFunc::EPSEvalInterface::TableEndPoints; + using InvBeh = ::Opm::SatFunc::EPSEvalInterface::InvalidEndpointBehaviour; namespace TwoPoint { using EPS = ::Opm::SatFunc::TwoPointScaling; @@ -340,29 +417,35 @@ namespace Create { struct Kr { static EPSPtr G(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); static EPSPtr OG(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); static EPSPtr OW(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); static EPSPtr W(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); }; struct Pc { static EPSPtr GO(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); static EPSPtr OW(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); }; static EPSPtr @@ -381,19 +464,23 @@ namespace Create { struct Kr { static EPSPtr G(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); static EPSPtr OG(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); static EPSPtr OW(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); static EPSPtr W(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init); + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid); }; static EPSPtr @@ -410,7 +497,8 @@ namespace Create { // Implementation of Create::TwoPoint::scalingFunction() Create::TwoPoint::EPSPtr Create::TwoPoint::Kr::G(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto sgcr = G.rawLinearisedCellData(init, "SGCR"); auto sgu = G.rawLinearisedCellData(init, "SGU"); @@ -424,12 +512,17 @@ Create::TwoPoint::Kr::G(const ::Opm::ECLGraph& G, }; } - return EPSPtr { new EPS { std::move(sgcr), std::move(sgu) } }; + return EPSPtr { + new EPS { + std::move(sgcr), std::move(sgu), handle_invalid + } + }; } Create::TwoPoint::EPSPtr Create::TwoPoint::Kr::OG(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto sogcr = G.rawLinearisedCellData(init, "SOGCR"); @@ -475,12 +568,17 @@ Create::TwoPoint::Kr::OG(const ::Opm::ECLGraph& G, } } - return EPSPtr { new EPS { std::move(sogcr), std::move(smax) } }; + return EPSPtr { + new EPS { + std::move(sogcr), std::move(smax), handle_invalid + } + }; } Create::TwoPoint::EPSPtr Create::TwoPoint::Kr::OW(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto sowcr = G.rawLinearisedCellData(init, "SOWCR"); @@ -526,12 +624,17 @@ Create::TwoPoint::Kr::OW(const ::Opm::ECLGraph& G, } } - return EPSPtr { new EPS { std::move(sowcr), std::move(smax) } }; + return EPSPtr { + new EPS { + std::move(sowcr), std::move(smax), handle_invalid + } + }; } Create::TwoPoint::EPSPtr Create::TwoPoint::Kr::W(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto swcr = G.rawLinearisedCellData(init, "SWCR"); auto swu = G.rawLinearisedCellData(init, "SWU"); @@ -542,12 +645,17 @@ Create::TwoPoint::Kr::W(const ::Opm::ECLGraph& G, }; } - return EPSPtr { new EPS { std::move(swcr), std::move(swu) } }; + return EPSPtr { + new EPS { + std::move(swcr), std::move(swu), handle_invalid + } + }; } Create::TwoPoint::EPSPtr Create::TwoPoint::Pc::GO(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto sgl = G.rawLinearisedCellData(init, "SGL"); auto sgu = G.rawLinearisedCellData(init, "SGU"); @@ -561,12 +669,17 @@ Create::TwoPoint::Pc::GO(const ::Opm::ECLGraph& G, }; } - return EPSPtr { new EPS { std::move(sgl), std::move(sgu) } }; + return EPSPtr { + new EPS { + std::move(sgl), std::move(sgu), handle_invalid + } + }; } Create::TwoPoint::EPSPtr Create::TwoPoint::Pc::OW(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto swl = G.rawLinearisedCellData(init, "SWL"); auto swu = G.rawLinearisedCellData(init, "SWU"); @@ -580,7 +693,11 @@ Create::TwoPoint::Pc::OW(const ::Opm::ECLGraph& G, }; } - return EPSPtr { new EPS { std::move(swl), std::move(swu) } }; + return EPSPtr { + new EPS { + std::move(swl), std::move(swu), handle_invalid + } + }; } Create::TwoPoint::EPSPtr @@ -606,10 +723,10 @@ scalingFunction(const ::Opm::ECLGraph& G, } if (opt.thisPh == PhIdx::Aqua) { - return Create::TwoPoint::Kr::W(G, init); + return Create::TwoPoint::Kr::W(G, init, opt.handle_invalid); } - return Create::TwoPoint::Kr::OW(G, init); + return Create::TwoPoint::Kr::OW(G, init, opt.handle_invalid); } if (opt.subSys == SSys::OilGas) { @@ -621,10 +738,10 @@ scalingFunction(const ::Opm::ECLGraph& G, } if (opt.thisPh == PhIdx::Vapour) { - return Create::TwoPoint::Kr::G(G, init); + return Create::TwoPoint::Kr::G(G, init, opt.handle_invalid); } - return Create::TwoPoint::Kr::OG(G, init); + return Create::TwoPoint::Kr::OG(G, init, opt.handle_invalid); } } @@ -637,11 +754,11 @@ scalingFunction(const ::Opm::ECLGraph& G, } if (opt.thisPh == PhIdx::Vapour) { - return Create::TwoPoint::Pc::GO(G, init); + return Create::TwoPoint::Pc::GO(G, init, opt.handle_invalid); } if (opt.thisPh == PhIdx::Aqua) { - return Create::TwoPoint::Pc::OW(G, init); + return Create::TwoPoint::Pc::OW(G, init, opt.handle_invalid); } } @@ -726,7 +843,8 @@ Create::TwoPoint::unscaledEndPoints(const RTEP& ep, const EPSOpt& opt) Create::ThreePoint::EPSPtr Create::ThreePoint::Kr::G(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto sgcr = G.rawLinearisedCellData(init, "SGCR"); auto sgu = G.rawLinearisedCellData(init, "SGU"); @@ -777,13 +895,16 @@ Create::ThreePoint::Kr::G(const ::Opm::ECLGraph& G, } return EPSPtr { - new EPS { std::move(sgcr), std::move(sr), std::move(sgu) } + new EPS { + std::move(sgcr), std::move(sr), std::move(sgu), handle_invalid + } }; } Create::ThreePoint::EPSPtr Create::ThreePoint::Kr::OG(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto sogcr = G.rawLinearisedCellData(init, "SOGCR"); @@ -850,13 +971,16 @@ Create::ThreePoint::Kr::OG(const ::Opm::ECLGraph& G, } return EPSPtr { - new EPS { std::move(sogcr), std::move(sdisp), std::move(smax) } + new EPS { + std::move(sogcr), std::move(sdisp), std::move(smax), handle_invalid + } }; } Create::ThreePoint::EPSPtr Create::ThreePoint::Kr::OW(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto sowcr = G.rawLinearisedCellData(init, "SOWCR"); @@ -923,13 +1047,16 @@ Create::ThreePoint::Kr::OW(const ::Opm::ECLGraph& G, } return EPSPtr { - new EPS { std::move(sowcr), std::move(sdisp), std::move(smax) } + new EPS { + std::move(sowcr), std::move(sdisp), std::move(smax), handle_invalid + } }; } Create::ThreePoint::EPSPtr Create::ThreePoint::Kr::W(const ::Opm::ECLGraph& G, - const ::Opm::ECLInitFileData& init) + const ::Opm::ECLInitFileData& init, + const InvBeh handle_invalid) { auto swcr = G.rawLinearisedCellData(init, "SWCR"); auto swu = G.rawLinearisedCellData(init, "SWU"); @@ -979,7 +1106,9 @@ Create::ThreePoint::Kr::W(const ::Opm::ECLGraph& G, } return EPSPtr { - new EPS { std::move(swcr), std::move(sdisp), std::move(swu) } + new EPS { + std::move(swcr), std::move(sdisp), std::move(swu), handle_invalid + } }; } @@ -1008,10 +1137,10 @@ scalingFunction(const ::Opm::ECLGraph& G, } if (opt.thisPh == PhIdx::Aqua) { - return Create::ThreePoint::Kr::W(G, init); + return Create::ThreePoint::Kr::W(G, init, opt.handle_invalid); } - return Create::ThreePoint::Kr::OW(G, init); + return Create::ThreePoint::Kr::OW(G, init, opt.handle_invalid); } if (opt.subSys == SSys::OilGas) { @@ -1023,10 +1152,10 @@ scalingFunction(const ::Opm::ECLGraph& G, } if (opt.thisPh == PhIdx::Vapour) { - return Create::ThreePoint::Kr::G(G, init); + return Create::ThreePoint::Kr::G(G, init, opt.handle_invalid); } - return Create::ThreePoint::Kr::OG(G, init); + return Create::ThreePoint::Kr::OG(G, init, opt.handle_invalid); } // Invalid. @@ -1124,9 +1253,10 @@ Opm::SatFunc::EPSEvalInterface::~EPSEvalInterface() // Class Opm::SatFunc::TwoPointScaling Opm::SatFunc::TwoPointScaling:: -TwoPointScaling(std::vector smin, - std::vector smax) - : pImpl_(new Impl(std::move(smin), std::move(smax))) +TwoPointScaling(std::vector smin, + std::vector smax, + InvalidEndpointBehaviour handle_invalid) + : pImpl_(new Impl(std::move(smin), std::move(smax), handle_invalid)) {} Opm::SatFunc::TwoPointScaling::~TwoPointScaling() @@ -1182,10 +1312,14 @@ Opm::SatFunc::TwoPointScaling::clone() const // Class Opm::SatFunc::ThreePointScaling Opm::SatFunc::ThreePointScaling:: -ThreePointScaling(std::vector smin, - std::vector sdisp, - std::vector smax) - : pImpl_(new Impl(std::move(smin), std::move(sdisp), std::move(smax))) +ThreePointScaling(std::vector smin, + std::vector sdisp, + std::vector smax, + InvalidEndpointBehaviour handle_invalid) + : pImpl_(new Impl(std::move(smin) , + std::move(sdisp), + std::move(smax) , + handle_invalid)) {} Opm::SatFunc::ThreePointScaling::~ThreePointScaling() diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp index 26f8fb208f..acdb9e3cc4 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLEndPointScaling.hpp @@ -67,6 +67,19 @@ namespace Opm { namespace SatFunc { double sat; }; + /// Policy for how to handle an invalid end-point scaling (e.g., if + /// lower and/or upper scaled saturations have nonsensical values + /// like -1.0E+20). + enum class InvalidEndpointBehaviour { + /// Use the unscaled value for this point. + UseUnscaled, + + /// Ignore this scaling request (e.g., produce no graphical + /// output for a scaled saturation function when the scaled + /// end-points are meaningless). + IgnorePoint, + }; + /// Convenience type alias. using SaturationPoints = std::vector; @@ -128,8 +141,19 @@ namespace Opm { namespace SatFunc { /// \param[in] smin Left end points for a set of cells. /// /// \param[in] smax Right end points for a set of cells. - TwoPointScaling(std::vector smin, - std::vector smax); + /// + /// \param[in] handle_invalid How to treat scaling requests with + /// invalid scaled saturations. This can, for instance, happen + /// if the scaled saturations are present in the result set but + /// some (or all) cells have irreconcilable values (e.g., minimum + /// saturation greater than maximum saturation, smin < -1E+20, + /// smax < -1E+20). + /// + /// Default behaviour: Use unscaled saturation if this happens. + TwoPointScaling(std::vector smin, + std::vector smax, + const InvalidEndpointBehaviour handle_invalid + = InvalidEndpointBehaviour::UseUnscaled); /// Destructor. ~TwoPointScaling(); @@ -237,9 +261,20 @@ namespace Opm { namespace SatFunc { /// for a set of cells. /// /// \param[in] smax Right end points for a set of cells. - ThreePointScaling(std::vector smin, - std::vector sdisp, - std::vector smax); + /// + /// \param[in] handle_invalid How to treat scaling requests with + /// invalid scaled saturations. This can, for instance, happen + /// if the scaled saturations are present in the result set but + /// some (or all) cells have irreconcilable values (e.g., minimum + /// saturation greater than maximum saturation, smin < -1E+20, + /// smax < -1E+20). + /// + /// Default behaviour: Use unscaled saturation if this happens. + ThreePointScaling(std::vector smin, + std::vector sdisp, + std::vector smax, + const InvalidEndpointBehaviour handle_invalid + = InvalidEndpointBehaviour::UseUnscaled); /// Destructor. ~ThreePointScaling(); @@ -385,6 +420,13 @@ namespace Opm { namespace SatFunc { /// auto eps = CreateEPS::fromECLOutput(G, init, opt); /// \endcode ::Opm::ECLPhaseIndex thisPh; + + /// How to handle an invalid end-point scaling (e.g., if lower + /// and/or upper scaled saturations have nonsensical values like + /// -1.0E+20). + EPSEvalInterface::InvalidEndpointBehaviour handle_invalid { + EPSEvalInterface::InvalidEndpointBehaviour::UseUnscaled + }; }; /// Collection of raw saturation table end points. diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.cpp index 6eafb186d6..56bb6d075f 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.cpp @@ -25,6 +25,7 @@ #include #include +#include #include @@ -262,15 +263,25 @@ Opm::ECLPVT::PVDx::viscosity(const std::vector& p) const }); } -Opm::FlowDiagnostics::Graph +Opm::ECLPVT::PVTGraph Opm::ECLPVT::PVDx::getPvtCurve(const RawCurve curve) const { if (curve == RawCurve::SaturatedState) { // Not applicable to dry gas or dead oil. Return empty. - return FlowDiagnostics::Graph{}; + return {}; } - return extractRawPVTCurve(this->interp_, curve); + auto ret = PVTGraph{}; + + auto basic = extractRawPVTCurve(this->interp_, curve); + + // Dead oil/dry gas. Mixing ratio == 0. + ret.mixRat.assign(basic.first.size(), 0.0); + + ret.press = std::move(basic.first); + ret.value = std::move(basic.second); + + return ret; } // ===================================================================== diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.hpp index 5759e58b7c..ef56d044a2 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.hpp @@ -301,6 +301,22 @@ namespace Opm { namespace ECLPVT { SaturatedState, }; + /// Collection of Miscible PVT States and Function Values. + /// + /// Intended for graphical representation of miscible and immiscible PVT + /// functions. + struct PVTGraph { + /// Pressure values + std::vector press; + + /// Mixing ratio. Typically Rs or Rv. + std::vector mixRat; + + /// Function value. Typically Bg, Bo, mu_g or mu_o. Usually a copy + /// of the mixing ratio for the case of the saturated state curve. + std::vector value; + }; + template class DenseVector { public: @@ -485,8 +501,8 @@ namespace Opm { namespace ECLPVT { /// /// \param[in] curve PVT property curve descriptor /// - /// \return 2D graph for PVT property curve identified by - /// requests represented by \p func. + /// \return 2D graph for PVT property curve identified by requests + /// represented by \p func. Mixing ratio is a vector of zeros. /// /// Example: Retrieve formation volume factor curve. /// @@ -494,7 +510,7 @@ namespace Opm { namespace ECLPVT { /// const auto graph = /// pvdx.getPvtCurve(ECLPVT::RawCurve::FVF); /// \endcode - FlowDiagnostics::Graph getPvtCurve(const RawCurve curve) const; + PVTGraph getPvtCurve(const RawCurve curve) const; private: /// Extrapolation policy for property evaluator/interpolant. @@ -658,7 +674,8 @@ namespace Opm { namespace ECLPVT { /// \return Collection of 2D graphs for PVT property curve /// identified by requests represented by \p curve. One curve /// (vector element) for each tabulated value of the function's - /// primary lookup key. + /// primary lookup key. Primary key expanded and assigned to the + /// appropriate graph's mixing ratio vector (PVTO convention). /// /// Example: Retrieve formation volume factor curves. /// @@ -666,7 +683,7 @@ namespace Opm { namespace ECLPVT { /// const auto graph = /// pvtx.getPvtCurve(ECLPVT::RawCurve::FVF); /// \endcode - std::vector + std::vector getPvtCurve(const RawCurve curve) const { if (curve == RawCurve::SaturatedState) { @@ -828,35 +845,52 @@ namespace Opm { namespace ECLPVT { return result; } - std::vector + std::vector mainPvtCurve(const RawCurve curve) const { - auto ret = std::vector{}; + // Uses setup appropriate for PVTO. The PVTG case must be + // handled in the caller. + + auto ret = std::vector{}; ret.reserve(this->propInterp_.size()); + auto i = static_cast::size_type>(0); + for (const auto& interp : this->propInterp_) { - ret.push_back(extractRawPVTCurve(interp, curve)); + auto basic = extractRawPVTCurve(interp, curve); + + ret.emplace_back(); + auto& pvtgraph = ret.back(); + + pvtgraph.mixRat.assign(basic.first.size(), this->key_[i]); + + pvtgraph.press = std::move(basic.first); + pvtgraph.value = std::move(basic.second); + + i += 1; } return ret; } - std::vector saturatedStateCurve() const + std::vector saturatedStateCurve() const { - auto y = std::vector{}; - y.reserve(this->propInterp_.size()); + // Uses setup appropriate for PVTO. The PVTG case must be + // handled in the caller. + + auto curve = PVTGraph{}; + + curve.press.reserve(this->propInterp_.size()); + curve.mixRat = this->key_; + curve.value = this->key_; for (const auto& interp : this->propInterp_) { const auto& yi = interp.independentVariable(); - y.push_back(yi[0]); + curve.press.push_back(yi[0]); } - return { - FlowDiagnostics::Graph { - this->key_, std::move(y) - } - }; + return { curve }; } }; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.cpp index 3c8518cadb..4a8ac74c9b 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.cpp @@ -44,13 +44,8 @@ namespace { return pvtnum; } - std::vector emptyFDGraph() - { - return { Opm::FlowDiagnostics::Graph{} }; - } - template - std::vector + std::vector rawPvtCurve(const PVTInterp* pvt, const Opm::ECLPVT::RawCurve curve, const int regID) @@ -58,7 +53,7 @@ namespace { if (pvt == nullptr) { // Result set does not provide requisite tabulated properties. // Return empty. - return emptyFDGraph(); + return {}; } return pvt->getPvtCurve(curve, regID); @@ -144,123 +139,113 @@ namespace { return viscosity(pvt, regID, po, rs); } - std::vector - convertCurve(std::vector&& curves, - const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_x, - const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_y) + std::vector + convertCurve(std::vector&& curves, + const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_p, + const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_R, + const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_f) { for (auto& curve : curves) { - cvrt_x.appliedTo(curve.first); - cvrt_y.appliedTo(curve.second); + cvrt_p.appliedTo(curve.press); + cvrt_R.appliedTo(curve.mixRat); + cvrt_f.appliedTo(curve.value); } return curves; } - std::vector - convertFvfCurve(std::vector&& curve, - const Opm::ECLPhaseIndex phase, - const Opm::ECLUnits::UnitSystem& usysFrom, - const Opm::ECLUnits::UnitSystem& usysTo) + std::vector + convertFvfCurve(std::vector&& curve, + const Opm::ECLPhaseIndex phase, + const Opm::ECLUnits::UnitSystem& usysFrom, + const Opm::ECLUnits::UnitSystem& usysTo) { assert ((phase == Opm::ECLPhaseIndex::Liquid) || (phase == Opm::ECLPhaseIndex::Vapour)); + auto cvrt_p = Opm::ECLUnits::Convert::Pressure(); + cvrt_p.from(usysFrom).to(usysTo); + if (phase == Opm::ECLPhaseIndex::Liquid) { - // Oil FVF. First column is pressure, second column is Bo. - auto cvrt_x = Opm::ECLUnits::Convert::Pressure(); - cvrt_x.from(usysFrom).to(usysTo); + // Oil FVF. Mixing ratio is Rs, value is Bo. + auto cvrt_R = Opm::ECLUnits::Convert::DissolvedGasOilRatio(); + cvrt_R.from(usysFrom).to(usysTo); - auto cvrt_y = Opm::ECLUnits::Convert::OilFVF(); - cvrt_y.from(usysFrom).to(usysTo); + auto cvrt_f = Opm::ECLUnits::Convert::OilFVF(); + cvrt_f.from(usysFrom).to(usysTo); - return convertCurve(std::move(curve), cvrt_x, cvrt_y); + return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f); } - // Gas FVF. Need to distinguish miscible from immiscible cases. In - // the former, the first column is Rv (vapourised oil/gas ratio) and - // in the second case the first column is the gas pressure. - // - // Immiscible case identified by curve.size() <= 1. + // Gas FVF. Mixing ratio is Rv, value is Bg. + auto cvrt_R = Opm::ECLUnits::Convert::VaporisedOilGasRatio(); + cvrt_R.from(usysFrom).to(usysTo); - auto cvrt_y = Opm::ECLUnits::Convert::GasFVF(); - cvrt_y.from(usysFrom).to(usysTo); + auto cvrt_f = Opm::ECLUnits::Convert::GasFVF(); + cvrt_f.from(usysFrom).to(usysTo); - if (curve.size() <= 1) { - // Immiscible Gas FVF. First column is Pg. - auto cvrt_x = Opm::ECLUnits::Convert::Pressure(); - cvrt_x.from(usysFrom).to(usysTo); - - return convertCurve(std::move(curve), cvrt_x, cvrt_y); - } - - // Miscible Gas FVF. First column is Rv. - auto cvrt_x = Opm::ECLUnits::Convert::VaporisedOilGasRatio(); - cvrt_x.from(usysFrom).to(usysTo); - - return convertCurve(std::move(curve), cvrt_x, cvrt_y); + return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f); } - std::vector - convertViscosityCurve(std::vector&& curve, - const Opm::ECLPhaseIndex phase, - const Opm::ECLUnits::UnitSystem& usysFrom, - const Opm::ECLUnits::UnitSystem& usysTo) + std::vector + convertViscosityCurve(std::vector&& curve, + const Opm::ECLPhaseIndex phase, + const Opm::ECLUnits::UnitSystem& usysFrom, + const Opm::ECLUnits::UnitSystem& usysTo) { assert ((phase == Opm::ECLPhaseIndex::Liquid) || (phase == Opm::ECLPhaseIndex::Vapour)); - // This is the viscosity curve. Second column is always viscosity - // irrespective of phase or miscible/immiscible fluids. - auto cvrt_y = Opm::ECLUnits::Convert::Viscosity(); - cvrt_y.from(usysFrom).to(usysTo); + auto cvrt_p = Opm::ECLUnits::Convert::Pressure(); + cvrt_p.from(usysFrom).to(usysTo); - if ((phase == Opm::ECLPhaseIndex::Liquid) || (curve.size() <= 1)) { - // Graph is oil viscosity or immiscible gas viscosity. First - // column is pressure. - auto cvrt_x = Opm::ECLUnits::Convert::Pressure(); - cvrt_x.from(usysFrom).to(usysTo); + // This is the viscosity curve. Value is always viscosity + // irrespective of mixing ratios. + auto cvrt_f = Opm::ECLUnits::Convert::Viscosity(); + cvrt_f.from(usysFrom).to(usysTo); - return convertCurve(std::move(curve), cvrt_x, cvrt_y); - } - - // Miscible Gas viscosity. First column is Rv (vapourised oil/gas - // ratio). - auto cvrt_x = Opm::ECLUnits::Convert::VaporisedOilGasRatio(); - cvrt_x.from(usysFrom).to(usysTo); - - return convertCurve(std::move(curve), cvrt_x, cvrt_y); - } - - std::vector - convertSatStateCurve(std::vector&& curve, - const Opm::ECLPhaseIndex phase, - const Opm::ECLUnits::UnitSystem& usysFrom, - const Opm::ECLUnits::UnitSystem& usysTo) - { - assert ((phase == Opm::ECLPhaseIndex::Liquid) || - (phase == Opm::ECLPhaseIndex::Vapour)); - - // First column is pressure (Po or Pg). - auto cvrt_x = Opm::ECLUnits::Convert::Pressure(); - cvrt_x.from(usysFrom).to(usysTo); - - // Second column is Rs or Rv depending on 'phase'. if (phase == Opm::ECLPhaseIndex::Liquid) { - // Saturated state curve for miscible oil. Second column is Rs - // (dissolved gas/oil ratio). - auto cvrt_y = Opm::ECLUnits::Convert::DissolvedGasOilRatio(); - cvrt_y.from(usysFrom).to(usysTo); + // Oil viscosity. Mixing ratio is Rs. + auto cvrt_R = Opm::ECLUnits::Convert::DissolvedGasOilRatio(); + cvrt_R.from(usysFrom).to(usysTo); - return convertCurve(std::move(curve), cvrt_x, cvrt_y); + return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f); } - // Saturated state curve for miscible gas. Second column is Rv - // (vapourised oil/gas ratio). - auto cvrt_y = Opm::ECLUnits::Convert::VaporisedOilGasRatio(); - cvrt_y.from(usysFrom).to(usysTo); + // Gas viscosity. Mixing ratio is Rv. + auto cvrt_R = Opm::ECLUnits::Convert::VaporisedOilGasRatio(); + cvrt_R.from(usysFrom).to(usysTo); - return convertCurve(std::move(curve), cvrt_x, cvrt_y); + return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f); + } + + std::vector + convertSatStateCurve(std::vector&& curve, + const Opm::ECLPhaseIndex phase, + const Opm::ECLUnits::UnitSystem& usysFrom, + const Opm::ECLUnits::UnitSystem& usysTo) + { + assert ((phase == Opm::ECLPhaseIndex::Liquid) || + (phase == Opm::ECLPhaseIndex::Vapour)); + + auto cvrt_p = Opm::ECLUnits::Convert::Pressure(); + cvrt_p.from(usysFrom).to(usysTo); + + if (phase == Opm::ECLPhaseIndex::Liquid) { + // Saturated state curve for miscible oil. Mixing ratio (and + // function value) is Rs (dissolved gas/oil ratio). + auto cvrt_R = Opm::ECLUnits::Convert::DissolvedGasOilRatio(); + cvrt_R.from(usysFrom).to(usysTo); + + return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_R); + } + + // Saturated state curve for miscible gas. Second column (and + // function value) is Rv (vapourised oil/gas ratio). + auto cvrt_R = Opm::ECLUnits::Convert::VaporisedOilGasRatio(); + cvrt_R.from(usysFrom).to(usysTo); + + return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_R); } } @@ -281,7 +266,7 @@ setOutputUnits(std::unique_ptr usys) this->usys_output_ = std::move(usys); } -std::vector +std::vector Opm::ECLPVT::ECLPvtCurveCollection:: getPvtCurve(const RawCurve curve, const ECLPhaseIndex phase, @@ -290,7 +275,7 @@ getPvtCurve(const RawCurve curve, if (! this->isValidRequest(phase, activeCell)) { // Not a supported phase or cell index out of bounds. Not a valid // request so return empty. - return emptyFDGraph(); + return {}; } // PVTNUM is traditional one-based region identifier. Subtract one to @@ -442,11 +427,11 @@ isValidRequest(const ECLPhaseIndex phase, < this->pvtnum_.size(); } -std::vector +std::vector Opm::ECLPVT::ECLPvtCurveCollection:: -convertToOutputUnits(std::vector&& graph, - const RawCurve curve, - const ECLPhaseIndex phase) const +convertToOutputUnits(std::vector&& graph, + const RawCurve curve, + const ECLPhaseIndex phase) const { if (this->usys_output_ == nullptr) { // No defined system of units for outputs. Return unconverted (SI). diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.hpp index 857bf8c0bb..d95a3911cb 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.hpp @@ -103,7 +103,7 @@ namespace Opm { namespace ECLPVT { /// pvtCC.getPvtCurve(ECLPVT::RawCurve::Viscosity, /// ECLPhaseIndex::Vapour, 31415); /// \endcode - std::vector + std::vector getPvtCurve(const RawCurve curve, const ECLPhaseIndex phase, const int activeCell) const; @@ -256,10 +256,10 @@ namespace Opm { namespace ECLPVT { /// \param[in] phase Phase for which to compute the property curve. /// Must be \code ECLPhaseIndex::Vapour \endcode or \code /// ECLPhaseIndex::Liquid \endcode. - std::vector - convertToOutputUnits(std::vector&& graph, - const RawCurve curve, - const ECLPhaseIndex phase) const; + std::vector + convertToOutputUnits(std::vector&& graph, + const RawCurve curve, + const ECLPhaseIndex phase) const; }; }} // Opm::ECLPVT diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp index f65c822cbf..77d5db5942 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp @@ -98,7 +98,7 @@ public: viscosity(const std::vector& rv, const std::vector& pg) const = 0; - virtual std::vector + virtual std::vector getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0; virtual std::unique_ptr clone() const = 0; @@ -133,7 +133,7 @@ public: return this->interpolant_.viscosity(pg); } - virtual std::vector + virtual std::vector getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override { return { this->interpolant_.getPvtCurve(curve) }; @@ -190,10 +190,30 @@ public: return this->interp_.viscosity(key, x); } - virtual std::vector + virtual std::vector getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override { - return this->interp_.getPvtCurve(curve); + // Interpolant blindly assigns pressure values to mixing ratios and + // vice versa. Need to switch those. + + auto graphs = this->interp_.getPvtCurve(curve); + + if (curve == Opm::ECLPVT::RawCurve::SaturatedState) { + assert ((graphs.size() == 1) && "Logic Error"); + + graphs[0].press.swap(graphs[0].mixRat); + graphs[0].value = graphs[0].mixRat; // Just for consistency + } + else { + assert ((graphs.size() > 1) && "Logic Error"); + + // FVF or viscosity. Just swap Pg <-> Rv. + for (auto& graph : graphs) { + graph.press.swap(graph.mixRat); + } + } + + return graphs; } virtual std::unique_ptr clone() const override @@ -391,7 +411,7 @@ public: return this->rhoS_[region]; } - std::vector + std::vector getPvtCurve(const RegIdx region, const RawCurve curve) const; @@ -460,7 +480,7 @@ viscosity(const RegIdx region, return this->eval_[region]->viscosity(rv, pg); } -std::vector +std::vector Opm::ECLPVT::Gas::Impl:: getPvtCurve(const RegIdx region, const RawCurve curve) const @@ -543,7 +563,7 @@ double Opm::ECLPVT::Gas::surfaceMassDensity(const int region) const return this->pImpl_->surfaceMassDensity(region); } -std::vector +std::vector Opm::ECLPVT::Gas:: getPvtCurve(const RawCurve curve, const int region) const diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.hpp index e5215e8f77..f3d175eb25 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.hpp @@ -175,7 +175,7 @@ namespace Opm { namespace ECLPVT { /// const auto graph = /// pvtGas.getPvtCurve(ECLPVT::RawCurve::FVF, 0); /// \endcode - std::vector + std::vector getPvtCurve(const RawCurve curve, const int region) const; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp index fb32f147ec..5161784815 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp @@ -95,7 +95,7 @@ public: viscosity(const std::vector& rs, const std::vector& po) const = 0; - virtual std::vector + virtual std::vector getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0; virtual std::unique_ptr clone() const = 0; @@ -130,7 +130,7 @@ public: return this->interpolant_.viscosity(po); } - virtual std::vector + virtual std::vector getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override { return { this->interpolant_.getPvtCurve(curve) }; @@ -187,7 +187,7 @@ public: return this->interp_.viscosity(key, x); } - virtual std::vector + virtual std::vector getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override { return this->interp_.getPvtCurve(curve); @@ -406,7 +406,7 @@ public: return this->rhoS_[region]; } - std::vector + std::vector getPvtCurve(const RegIdx region, const RawCurve curve) const; @@ -475,7 +475,7 @@ viscosity(const RegIdx region, return this->eval_[region]->viscosity(rs, po); } -std::vector +std::vector Opm::ECLPVT::Oil::Impl:: getPvtCurve(const RegIdx region, const RawCurve curve) const @@ -558,7 +558,7 @@ double Opm::ECLPVT::Oil::surfaceMassDensity(const int region) const return this->pImpl_->surfaceMassDensity(region); } -std::vector +std::vector Opm::ECLPVT::Oil:: getPvtCurve(const RawCurve curve, const int region) const diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp index e58f01658b..78235de9b0 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp @@ -172,7 +172,7 @@ namespace Opm { namespace ECLPVT { /// const auto graph = /// pvtOil.getPvtCurve(ECLPVT::RawCurve::Viscosity, 3); /// \endcode - std::vector + std::vector getPvtCurve(const RawCurve curve, const int region) const; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp index 9f7c0a71df..3e486a3aa8 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.cpp @@ -25,6 +25,7 @@ #include #include +#include #include #include #include @@ -682,9 +683,12 @@ public: Impl(Impl&& rhs); Impl(const Impl& rhs); - void init(const ECLGraph& G, - const ECLInitFileData& init, - const bool useEPS); + void init(const ECLGraph& G, + const ECLInitFileData& init, + const bool useEPS, + const InvalidEPBehaviour handle_invalid); + + void setOutputUnits(std::unique_ptr usys); std::vector relperm(const ECLGraph& G, @@ -716,14 +720,16 @@ private: bool wat; }; - void define(const ECLGraph& G, - const ECLInitFileData& init, - const RawTEP& ep, - const bool use3PtScaling, - const ActPh& active) + void define(const ECLGraph& G, + const ECLInitFileData& init, + const RawTEP& ep, + const bool use3PtScaling, + const ActPh& active, + const InvalidEPBehaviour handle_invalid) { auto opt = Create::EPSOptions{}; - opt.use3PtScaling = use3PtScaling; + opt.use3PtScaling = use3PtScaling; + opt.handle_invalid = handle_invalid; if (active.oil) { this->create_oil_eps(G, init, ep, active, opt); @@ -1093,11 +1099,15 @@ private: ECLRegionMapping rmap_; - std::unique_ptr oil_; - std::unique_ptr gas_; - std::unique_ptr wat_; + std::unique_ptr usys_internal_{nullptr}; - std::unique_ptr eps_; + std::unique_ptr oil_{nullptr}; + std::unique_ptr gas_{nullptr}; + std::unique_ptr wat_{nullptr}; + + std::unique_ptr eps_{nullptr}; + + std::unique_ptr usys_output_{nullptr}; void initRelPermInterp(const EPSEvaluator::ActPh& active, const ECLInitFileData& init, @@ -1106,7 +1116,8 @@ private: void initEPS(const EPSEvaluator::ActPh& active, const bool use3PtScaling, const ECLGraph& G, - const ECLInitFileData& init); + const ECLInitFileData& init, + const InvalidEPBehaviour handle_invalid); std::vector kro(const ECLGraph& G, @@ -1224,25 +1235,56 @@ private: Opm::ECLSaturationFunc::Impl::Impl(const ECLGraph& G, const ECLInitFileData& init) - : satnum_(satnumVector(G, init)) - , rmap_ (satnum_) -{ -} + : satnum_ (satnumVector(G, init)) + , rmap_ (satnum_) + , usys_internal_(ECLUnits::internalUnitConventions()) +{} Opm::ECLSaturationFunc::Impl::Impl(Impl&& rhs) - : satnum_(std::move(rhs.satnum_)) - , rmap_ (std::move(rhs.rmap_)) - , oil_ (std::move(rhs.oil_ )) - , gas_ (std::move(rhs.gas_ )) - , wat_ (std::move(rhs.wat_ )) + : satnum_ (std::move(rhs.satnum_)) + , rmap_ (std::move(rhs.rmap_)) + , usys_internal_(std::move(rhs.usys_internal_)) + , oil_ (std::move(rhs.oil_)) + , gas_ (std::move(rhs.gas_)) + , wat_ (std::move(rhs.wat_)) + , eps_ (std::move(rhs.eps_)) + , usys_output_ (std::move(rhs.usys_output_)) {} +Opm::ECLSaturationFunc::Impl::Impl(const Impl& rhs) + : satnum_ (rhs.satnum_) + , rmap_ (rhs.rmap_) + , usys_internal_(rhs.usys_internal_->clone()) +{ + if (rhs.oil_) { + // Polymorphic object must use clone(). + this->oil_ = rhs.oil_->clone(); + } + + if (rhs.gas_) { + this->gas_.reset(new Gas::SatFunction(*rhs.gas_)); + } + + if (rhs.wat_) { + this->wat_.reset(new Water::SatFunction(*rhs.wat_)); + } + + if (rhs.eps_) { + this->eps_.reset(new EPSEvaluator(*rhs.eps_)); + } + + if (rhs.usys_output_) { + this->usys_output_ = rhs.usys_output_->clone(); + } +} + // --------------------------------------------------------------------- void -Opm::ECLSaturationFunc::Impl::init(const ECLGraph& G, - const ECLInitFileData& init, - const bool useEPS) +Opm::ECLSaturationFunc::Impl::init(const ECLGraph& G, + const ECLInitFileData& init, + const bool useEPS, + const InvalidEPBehaviour handle_invalid) { // Extract INTEHEAD from main grid const auto& ih = init.keywordData(INTEHEAD_KW); @@ -1263,7 +1305,7 @@ Opm::ECLSaturationFunc::Impl::init(const ECLGraph& G, lh[LOGIHEAD_ALT_ENDPOINT_SCALING_INDEX]); // Must be called *after* initRelPermInterp(). - this->initEPS(active, use3PtScaling, G, init); + this->initEPS(active, use3PtScaling, G, init, handle_invalid); } } } @@ -1322,36 +1364,25 @@ void Opm::ECLSaturationFunc:: Impl::initEPS(const EPSEvaluator::ActPh& active, const bool use3PtScaling, const ECLGraph& G, - const ECLInitFileData& init) + const ECLInitFileData& init, + const InvalidEPBehaviour handle_invalid) { const auto ep = this->extractRawTableEndPoints(active); this->eps_.reset(new EPSEvaluator()); - this->eps_->define(G, init, ep, use3PtScaling, active); + this->eps_->define(G, init, ep, use3PtScaling, active, handle_invalid); } // ##################################################################### -Opm::ECLSaturationFunc::Impl::Impl(const Impl& rhs) - : rmap_(rhs.rmap_) +void +Opm::ECLSaturationFunc::Impl:: +setOutputUnits(std::unique_ptr usys) { - if (rhs.oil_) { - // Polymorphic object must use clone(). - this->oil_ = rhs.oil_->clone(); - } - - if (rhs.gas_) { - this->gas_.reset(new Gas::SatFunction(*rhs.gas_)); - } - - if (rhs.wat_) { - this->wat_.reset(new Water::SatFunction(*rhs.wat_)); - } + this->usys_output_ = std::move(usys); } -// ##################################################################### - std::vector Opm::ECLSaturationFunc::Impl:: relperm(const ECLGraph& G, @@ -1717,7 +1748,13 @@ pcgoCurve(const ECLRegionMapping& rmap, // Region ID 'reg' is traditional, ECL-style one-based region ID // (SATNUM). Subtract one to create valid table index. - const auto pc = this->gas_->pcgo(regID - 1, sg_inp); + auto pc = this->gas_->pcgo(regID - 1, sg_inp); + + if (this->usys_output_ != nullptr) { + ::Opm::ECLUnits::Convert::Pressure() + .from(*this->usys_internal_) + .to (*this->usys_output_).appliedTo(pc); + } // FD::Graph == pair, vector> return FlowDiagnostics::Graph { @@ -1831,7 +1868,13 @@ pcowCurve(const ECLRegionMapping& rmap, // Region ID 'reg' is traditional, ECL-style one-based region ID // (SATNUM). Subtract one to create valid table index. - const auto pc = this->wat_->pcow(regID - 1, sw_inp); + auto pc = this->wat_->pcow(regID - 1, sw_inp); + + if (this->usys_output_ != nullptr) { + ::Opm::ECLUnits::Convert::Pressure() + .from(*this->usys_internal_) + .to (*this->usys_output_).appliedTo(pc); + } // FD::Graph == pair, vector> return FlowDiagnostics::Graph { @@ -1941,12 +1984,13 @@ extractRawTableEndPoints(const EPSEvaluator::ActPh& active) const // ===================================================================== Opm::ECLSaturationFunc:: -ECLSaturationFunc(const ECLGraph& G, - const ECLInitFileData& init, - const bool useEPS) +ECLSaturationFunc(const ECLGraph& G, + const ECLInitFileData& init, + const bool useEPS, + const InvalidEPBehaviour handle_invalid) : pImpl_(new Impl(G, init)) { - this->pImpl_->init(G, init, useEPS); + this->pImpl_->init(G, init, useEPS, handle_invalid); } Opm::ECLSaturationFunc::~ECLSaturationFunc() @@ -1976,6 +2020,13 @@ Opm::ECLSaturationFunc::operator=(const ECLSaturationFunc& rhs) return *this; } +void +Opm::ECLSaturationFunc:: +setOutputUnits(std::unique_ptr usys) +{ + this->pImpl_->setOutputUnits(std::move(usys)); +} + std::vector Opm::ECLSaturationFunc:: relperm(const ECLGraph& G, diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp index be31c7e2f4..ccf249d9f0 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLSaturationFunc.hpp @@ -21,7 +21,9 @@ #define OPM_ECLSATURATIONFUNC_HEADER_INCLUDED #include +#include #include +#include #include #include @@ -77,6 +79,9 @@ namespace Opm { ECLPhaseIndex thisPh; }; + using InvalidEPBehaviour = ::Opm::SatFunc:: + EPSEvalInterface::InvalidEndpointBehaviour; + /// Constructor /// /// \param[in] G Connected topology of current model's active cells. @@ -96,9 +101,18 @@ namespace Opm { /// /// Default value (\c true) means that effects of EPS are /// included if requisite data is present in the INIT result. - ECLSaturationFunc(const ECLGraph& G, - const ECLInitFileData& init, - const bool useEPS = true); + /// + /// \param[in] invalidIsUnscaled Whether or not treat invalid scaled + /// saturation end-points (e.g., SWL=-1.0E+20) as unscaled + /// saturations. True for "treat as unscaled", false for " + /// + /// Default value (\c true) means that invalid scalings are + /// treated as unscaled, false + ECLSaturationFunc(const ECLGraph& G, + const ECLInitFileData& init, + const bool useEPS = true, + const InvalidEPBehaviour handle_invalid + = InvalidEPBehaviour::UseUnscaled); /// Destructor. ~ECLSaturationFunc(); @@ -137,6 +151,21 @@ namespace Opm { /// \return \code *this \endcode. ECLSaturationFunc& operator=(const ECLSaturationFunc& rhs); + /// Define a collection of units of measure for output purposes. + /// + /// Capillary pressure curves produced by getSatFuncCurve() will be + /// reported in the pressure units of this system. If this function + /// is never called (or called with null pointer), then the output + /// units are implicitly set to the flow-diagnostics module's + /// internal units of measurement (meaning all properties and curves + /// will be reported in strict SI units). + /// + /// \param[in] usys Collection of units of measure for output + /// purposes. Typically the return value from one of the \code + /// *UnitConvention() \endcode functions of the \code ECLUnits + /// \endcode namespace. + void setOutputUnits(std::unique_ptr usys); + /// Compute relative permeability values in all active cells for a /// single phase. ///