From 5efe0f705e83fdbb904957a6ed5c79fd0de72fd1 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Jacob=20St=C3=B8ren?= Date: Wed, 22 Nov 2017 15:28:38 +0100 Subject: [PATCH] #2174 Upgrade opm-flowdiagnostics-applications to 75b333335f6cd055d3130d460c6d87444fb7aed4 to get more of the PVT/RelPerm functionality. --- ResInsightVersion.cmake | 2 +- .../examples/exampleSetup.hpp | 2 +- .../examples/extractPropCurves.cpp | 116 ++++- .../opm/utility/ECLCaseUtilities.cpp | 19 + .../opm/utility/ECLEndPointScaling.cpp | 6 + .../opm/utility/ECLPvtCommon.cpp | 35 +- .../opm/utility/ECLPvtCommon.hpp | 91 ++++ .../opm/utility/ECLPvtCurveCollection.cpp | 72 ++- .../opm/utility/ECLPvtCurveCollection.hpp | 43 +- .../opm/utility/ECLPvtGas.cpp | 137 ++++- .../opm/utility/ECLPvtGas.hpp | 19 +- .../opm/utility/ECLPvtOil.cpp | 137 ++++- .../opm/utility/ECLPvtOil.hpp | 19 +- .../opm/utility/ECLSaturationFunc.cpp | 485 +++++++++++++----- .../opm/utility/ECLSaturationFunc.hpp | 3 + 15 files changed, 950 insertions(+), 236 deletions(-) diff --git a/ResInsightVersion.cmake b/ResInsightVersion.cmake index d626d63d25..730726692e 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 "e6ae223fbeb0c9bd83c63a5fb6f57e4a3ad2ec18") +set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "75b333335f6cd055d3130d460c6d87444fb7aed4") # 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/exampleSetup.hpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp index 07641c775b..55d3240cbb 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/exampleSetup.hpp @@ -247,7 +247,7 @@ namespace example { } { - auto wsol = Opm::ECLWellSolution{}; + auto wsol = Opm::ECLWellSolution{-1.0, false}; well_fluxes = wsol.solution(*restart, graph.activeGrids()); } 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 a64a5792f8..737819ee87 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 @@ -37,24 +37,28 @@ namespace { template - void printGraph(OStream& os, - const std::string& name, - const Opm::FlowDiagnostics::Graph& graph) + void printGraph(OStream& os, + const std::string& name, + const std::vector& graphs) { - const auto& x = graph.first; - const auto& y = graph.second; - const auto oprec = os.precision(16); const auto oflags = os.setf(std::ios_base::scientific); - os << name << " = [\n"; + auto k = 1; + for (const auto& graph : graphs) { + const auto& x = graph.first; + const auto& y = graph.second; - for (auto n = x.size(), i = 0*n; i < n; ++i) { - os << x[i] << ' ' << y[i] << '\n'; + os << name << '{' << k << "} = [\n"; + + for (auto n = x.size(), i = 0*n; i < n; ++i) { + os << x[i] << ' ' << y[i] << '\n'; + } + + os << "];\n\n"; + k += 1; } - os << "];\n\n"; - os.setf(oflags); os.precision(oprec); } @@ -81,7 +85,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, useEPS); - printGraph(std::cout, "krg", graph[0]); + printGraph(std::cout, "krg", graph); } void krog(const Opm::ECLSaturationFunc& sfunc, @@ -103,7 +107,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, useEPS); - printGraph(std::cout, "krog", graph[0]); + printGraph(std::cout, "krog", graph); } void krow(const Opm::ECLSaturationFunc& sfunc, @@ -125,7 +129,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, useEPS); - printGraph(std::cout, "krow", graph[0]); + printGraph(std::cout, "krow", graph); } void krw(const Opm::ECLSaturationFunc& sfunc, @@ -147,7 +151,54 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, useEPS); - printGraph(std::cout, "krw", graph[0]); + printGraph(std::cout, "krw", graph); + } + + // ----------------------------------------------------------------- + // Capillary pressure + + void pcgo(const Opm::ECLSaturationFunc& sfunc, + const int activeCell, + const bool useEPS) + { + using RC = Opm::ECLSaturationFunc::RawCurve; + + auto func = std::vector{}; + func.reserve(1); + + // Request pcgo (gas/oil capillary pressure (Pg-Po) in G/O system) + func.push_back(RC{ + RC::Function::CapPress, + RC::SubSystem::OilGas, + Opm::ECLPhaseIndex::Vapour + }); + + const auto graph = + sfunc.getSatFuncCurve(func, activeCell, useEPS); + + printGraph(std::cout, "pcgo", graph); + } + + void pcow(const Opm::ECLSaturationFunc& sfunc, + const int activeCell, + const bool useEPS) + { + using RC = Opm::ECLSaturationFunc::RawCurve; + + auto func = std::vector{}; + func.reserve(1); + + // Request pcow (oil/water capillary pressure (Po-Pw) in O/W system) + func.push_back(RC{ + RC::Function::CapPress, + RC::SubSystem::OilWater, + Opm::ECLPhaseIndex::Aqua + }); + + const auto graph = + sfunc.getSatFuncCurve(func, activeCell, useEPS); + + printGraph(std::cout, "pcow", graph); } // ----------------------------------------------------------------- @@ -196,6 +247,33 @@ namespace { printGraph(std::cout, "mu_o", graph); } + + // ----------------------------------------------------------------- + // Saturated states (RvSat(Pg) and RsSat(Po)) + + void rvSat(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, + const int activeCell) + { + using RC = Opm::ECLPVT::RawCurve; + using PI = Opm::ECLPhaseIndex; + + const auto graph = pvtCurves + .getPvtCurve(RC::SaturatedState, PI::Vapour, activeCell); + + printGraph(std::cout, "rvSat", graph); + } + + void rsSat(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, + const int activeCell) + { + using RC = Opm::ECLPVT::RawCurve; + using PI = Opm::ECLPhaseIndex; + + const auto graph = pvtCurves + .getPvtCurve(RC::SaturatedState, PI::Liquid, activeCell); + + printGraph(std::cout, "rsSat", graph); + } } // namespace Anonymous int main(int argc, char* argv[]) @@ -218,6 +296,11 @@ try { if (prm.getDefault("krow", false)) { krow(sfunc, cellID, useEPS); } if (prm.getDefault("krw" , false)) { krw (sfunc, cellID, useEPS); } + // ----------------------------------------------------------------- + // Capillary pressure + if (prm.getDefault("pcgo", false)) { pcgo(sfunc, cellID, useEPS); } + if (prm.getDefault("pcow", false)) { pcow(sfunc, cellID, useEPS); } + // ----------------------------------------------------------------- // PVT Curves @@ -225,6 +308,9 @@ try { if (prm.getDefault("mu_g", false)) { mu_g(pvtCC, cellID); } if (prm.getDefault("Bo" , false)) { Bo (pvtCC, cellID); } if (prm.getDefault("mu_o", false)) { mu_o(pvtCC, cellID); } + + if (prm.getDefault("rvSat", false)) { rvSat(pvtCC, cellID); } + if (prm.getDefault("rsSat", false)) { rsSat(pvtCC, cellID); } } catch (const std::exception& e) { std::cerr << "Caught Exception: " << e.what() << '\n'; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLCaseUtilities.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLCaseUtilities.cpp index 1c6f00f3dd..8700f01176 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLCaseUtilities.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLCaseUtilities.cpp @@ -1,3 +1,22 @@ +/* + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + #include #include 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 cea4bacffa..a8d3fcd306 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 @@ -989,7 +989,10 @@ scalingFunction(const ::Opm::ECLGraph& G, const ::Opm::ECLInitFileData& init, const ::Opm::SatFunc::CreateEPS::EPSOptions& opt) { +#if !defined(NDEBUG) using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; +#endif // !defined(NDEBUG) + using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; using PhIdx = ::Opm::ECLPhaseIndex; @@ -1035,7 +1038,10 @@ scalingFunction(const ::Opm::ECLGraph& G, std::vector Create::ThreePoint::unscaledEndPoints(const RTEP& ep, const EPSOpt& opt) { +#if !defined(NDEBUG) using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; +#endif // !defined(NDEBUG) + using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; using PhIdx = ::Opm::ECLPhaseIndex; 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 5df4b6c941..7c7083d5c6 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 @@ -265,40 +265,7 @@ Opm::ECLPVT::PVDx::viscosity(const std::vector& p) const Opm::FlowDiagnostics::Graph Opm::ECLPVT::PVDx::getPvtCurve(const RawCurve curve) const { - assert ((curve == RawCurve::FVF) || - (curve == RawCurve::Viscosity)); - - const auto colID = (curve == RawCurve::FVF) - ? std::size_t{0} : std::size_t{1}; - - auto x = this->interp_.independentVariable(); - auto y = this->interp_.resultVariable(colID); - - assert ((x.size() == y.size()) && "Setup Error"); - - // Post-process ordinates according to which curve is requested. - if (curve == RawCurve::FVF) { - // y == 1/B. Convert to proper FVF. - for (auto& yi : y) { - yi = 1.0 / yi; - } - } - else { - // y == 1/(B*mu). Extract viscosity term through the usual - // conversion formula: - // - // (1 / B) / (1 / (B*mu)). - const auto b = this->interp_.resultVariable(0); // 1/B - - assert ((b.size() == y.size()) && "Setup Error"); - - for (auto n = y.size(), i = 0*n; i < n; ++i) { - y[i] = b[i] / y[i]; - } - } - - // Graph == pair, vector> - return FlowDiagnostics::Graph { std::move(x), std::move(y) }; + return extractRawPVTCurve(this->interp_, curve); } // ===================================================================== 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 17cd9637f4..c054754335 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 @@ -28,6 +28,8 @@ #include #include +#include +#include #include #include #include @@ -294,6 +296,9 @@ namespace Opm { namespace ECLPVT { /// Viscosity Viscosity, + + /// Enveloping curve for saturated state (wet gas or live oil) + SaturatedState, }; template @@ -390,6 +395,47 @@ namespace Opm { namespace ECLPVT { return u -= v; } + template + FlowDiagnostics::Graph + extractRawPVTCurve(const Interp1D::PiecewisePolynomial::Linear& interpolant, + const RawCurve curve) + { + assert ((curve == RawCurve::FVF) || + (curve == RawCurve::Viscosity)); + + const auto colID = (curve == RawCurve::FVF) + ? std::size_t{0} : std::size_t{1}; + + auto x = interpolant.independentVariable(); + auto y = interpolant.resultVariable(colID); + + assert ((x.size() == y.size()) && "Setup Error"); + + // Post-process ordinates according to which curve is requested. + if (curve == RawCurve::FVF) { + // y == 1/B. Convert to proper FVF. + for (auto& yi : y) { + yi = 1.0 / yi; + } + } + else { + // y == 1/(B*mu). Extract viscosity term through the usual + // conversion formula: + // + // (1 / B) / (1 / (B*mu)). + const auto b = interpolant.resultVariable(0); // 1/B + + assert ((b.size() == y.size()) && "Setup Error"); + + for (auto n = y.size(), i = 0*n; i < n; ++i) { + y[i] = b[i] / y[i]; + } + } + + // Graph == pair, vector> + return FlowDiagnostics::Graph { std::move(x), std::move(y) }; + } + /// Evaluate pressure-dependent properties (formation volume factor, /// viscosity &c) for dead oil (PVDO) or dry gas (PVDG) from tabulated /// functions as represented in an ECL result set (ECLInitData). @@ -602,6 +648,51 @@ namespace Opm { namespace ECLPVT { }); } + /// Retrieve 2D graph representation PVT property function. + /// + /// \param[in] curve PVT property curve descriptor + /// + /// \return Collection of 2D graphs for PVT property curve + /// identified by requests represented by \p func. + /// + /// \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. + /// + /// Example: Retrieve formation volume factor curves. + /// + /// \code + /// const auto graph = + /// pvtx.getPvtCurve(ECLPVT::RawCurve::FVF); + /// \endcode + std::vector + getPvtCurve(const RawCurve curve) const + { + auto ret = std::vector{}; + ret.reserve(this->propInterp_.size()); + + for (const auto& interp : this->propInterp_) { + ret.push_back(extractRawPVTCurve(interp, curve)); + } + + return ret; + } + + std::vector getSaturatedPoints() const + { + auto y = std::vector{}; + y.reserve(this->propInterp_.size()); + + for (const auto& interp : this->propInterp_) { + const auto& yi = interp.independentVariable(); + + y.push_back(yi[0]); + } + + return y; + } + private: using InnerEvalPoint = typename std::decay< decltype(std::declval().classifyPoint(0.0)) 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 38d7aa93e4..7bc5a391b5 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 @@ -1,9 +1,31 @@ +/* + Copyright 2017 Statoil ASA. + + This file is part of the Open Porous Media Project (OPM). + + OPM is free software: you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation, either version 3 of the License, or + (at your option) any later version. + + OPM is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with OPM. If not, see . +*/ + #include #include +#include #include +#include + namespace { std::vector pvtnumVector(const ::Opm::ECLGraph& G, @@ -19,6 +41,35 @@ namespace { return pvtnum; } + + std::unique_ptr + createUnitSystem(const ::Opm::ECLInitFileData& init) + { + const auto& ih = init.keywordData(INTEHEAD_KW); + + return ::Opm::ECLUnits::createUnitSystem(ih[ INTEHEAD_UNIT_INDEX ]); + } + + std::vector emptyFDGraph() + { + return { Opm::FlowDiagnostics::Graph{} }; + } + + template + std::vector + rawPvtCurve(PvtPtr& pvt, // ref to shared_ptr<> + const Opm::ECLPVT::RawCurve curve, + const int regID, + const Opm::ECLUnits::UnitSystem& usys) + { + if (pvt != nullptr) { + return pvt->getPvtCurve(curve, regID, usys); + } + + // Result set does not provide requisite tabulated properties. + // Return empty. + return emptyFDGraph(); + } } Opm::ECLPVT::ECLPvtCurveCollection:: @@ -27,9 +78,10 @@ ECLPvtCurveCollection(const ECLGraph& G, : pvtnum_(pvtnumVector(G, init)) , gas_ (CreateGasPVTInterpolant::fromECLOutput(init)) // u_p<> -> s_p<> , oil_ (CreateOilPVTInterpolant::fromECLOutput(init)) // u_p<> -> s_p<> + , usys_ (createUnitSystem(init)) // u_p<> -> s_p<> {} -Opm::FlowDiagnostics::Graph +std::vector Opm::ECLPVT::ECLPvtCurveCollection:: getPvtCurve(const RawCurve curve, const ECLPhaseIndex phase, @@ -38,14 +90,14 @@ getPvtCurve(const RawCurve curve, if (phase == ECLPhaseIndex::Aqua) { // Not supported at this time. // Return empty. - return FlowDiagnostics::Graph{}; + return emptyFDGraph(); } if (static_castpvtnum_.size())>(activeCell) >= this->pvtnum_.size()) { // Active cell index out of bounds. Return empty. - return FlowDiagnostics::Graph{}; + return emptyFDGraph(); } // PVTNUM is traditional one-based region identifier. Subtract one to @@ -53,16 +105,10 @@ getPvtCurve(const RawCurve curve, const auto regID = this->pvtnum_[activeCell] - 1; if (phase == ECLPhaseIndex::Liquid) { - // return this->oil_->getPvtCurve(curve, regID); - return FlowDiagnostics::Graph{}; + // Caller requests oil properties. + return rawPvtCurve(this->oil_, curve, regID, *this->usys_); } - if (this->gas_) { - return this->gas_->getPvtCurve(curve, regID); - } - else { - // Result set does not provide tabulated gas properties. Return - // empty. - return FlowDiagnostics::Graph{}; - } + // Call requests gas properties. + return rawPvtCurve(this->gas_, curve, regID, *this->usys_); } 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 dfd761afcd..a3a7fd1e5e 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 @@ -25,6 +25,7 @@ #include #include #include +#include #include #include @@ -45,10 +46,47 @@ namespace Opm { namespace ECLPVT { class ECLPvtCurveCollection { public: + /// Constructor + /// + /// \param[in] G Connected topology of current model's active cells. + /// Needed to linearise region mapping (e.g., SATNUM) that is + /// distributed on local grids to all of the model's active cells + /// (\code member function G.rawLinearisedCellData() \endcode). + /// + /// \param[in] init Container of tabulated PVT functions for all PVT + /// regions in the model \p G. ECLPvtCurveCollection(const ECLGraph& G, const ECLInitFileData& init); - FlowDiagnostics::Graph + /// Retrieve 2D graph representation of Phase PVT property function + /// in a specific active cell. + /// + /// \param[in] curve PVT property curve descriptor + /// + /// \param[in] phase Phase for which to compute extract graph + /// representation of PVT property function. + /// + /// \param[in] activeCell Index of particular active cell in model.. + /// + /// \return Collection of 2D graphs for PVT property curve + /// identified by requests represented by \p curve, \p phase and + /// \p activeCell. One curve (vector element) for each tabulated + /// node of the primary look-up key. Single curve (i.e., a + /// single element vector) in the case of dry gas (no vaporised + /// oil) or dead oil (no dissolved gas). + /// + /// No curves for water or dead oil with constant compressibility + /// (i.e., keyword 'PVCDO' in the input deck). + /// + /// Example: Retrieve collection of gas viscosity curves pertaining + /// to model's active cell 31415. + /// + /// \code + /// const auto curves = + /// pvtCC.getPvtCurve(ECLPVT::RawCurve::Viscosity, + /// ECLPhaseIndex::Vapour, 31415); + /// \endcode + std::vector getPvtCurve(const RawCurve curve, const ECLPhaseIndex phase, const int activeCell) const; @@ -62,6 +100,9 @@ namespace Opm { namespace ECLPVT { /// Oil PVT property evaluator. std::shared_ptr oil_; + + /// Unit handling (SI -> result-set convention) + std::shared_ptr usys_; }; }} // 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 05ba66abc5..06be02035f 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 @@ -25,10 +25,13 @@ #include #include +#include + #include #include #include #include +#include #include #include #include @@ -95,8 +98,9 @@ public: viscosity(const std::vector& rv, const std::vector& pg) const = 0; - virtual Opm::FlowDiagnostics::Graph - getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0; + virtual std::vector + getPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const = 0; virtual std::unique_ptr clone() const = 0; }; @@ -130,10 +134,31 @@ public: return this->interpolant_.viscosity(pg); } - virtual Opm::FlowDiagnostics::Graph - getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override + virtual std::vector + getPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const override { - return this->interpolant_.getPvtCurve(curve); + if (curve == Opm::ECLPVT::RawCurve::SaturatedState) { + // Not applicable to dry gas. Return empty. + return { Opm::FlowDiagnostics::Graph{} }; + } + + auto pvtcurve = this->interpolant_.getPvtCurve(curve); + + const auto x_unit = usys.pressure(); + const auto y_unit = (curve == ::Opm::ECLPVT::RawCurve::FVF) + ? (usys.reservoirVolume() / usys.surfaceVolumeGas()) + : usys.viscosity(); + + auto& x = pvtcurve.first; + auto& y = pvtcurve.second; + + for (auto n = x.size(), i = 0*n; i < n; ++i) { + x[i] = ::Opm::unit::convert::to(x[i], x_unit); + y[i] = ::Opm::unit::convert::to(y[i], y_unit); + } + + return { std::move(pvtcurve) }; } virtual std::unique_ptr clone() const override @@ -158,7 +183,8 @@ public: WetGas(std::vector key, std::vector propInterp) - : interp_(std::move(key), std::move(propInterp)) + : key_ (key) // Copy + , interp_(std::move(key), std::move(propInterp)) {} virtual std::vector @@ -187,13 +213,15 @@ public: return this->interp_.viscosity(key, x); } - virtual Opm::FlowDiagnostics::Graph - getPvtCurve(const Opm::ECLPVT::RawCurve /* curve */) const override + virtual std::vector + getPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const override { - throw std::runtime_error { - "Property Evaluator for Wet Gas Does not " - "Support Retrieving Raw Curves (Blame BSKA)" - }; + if (curve == ::Opm::ECLPVT::RawCurve::SaturatedState) { + return this->saturatedStateCurve(usys); + } + + return this->mainPvtCurve(curve, usys); } virtual std::unique_ptr clone() const override @@ -204,9 +232,68 @@ public: private: using TableInterpolant = ::Opm::ECLPVT::PVTx; - TableInterpolant interp_; + std::vector key_; + TableInterpolant interp_; + + std::vector + mainPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const; + + std::vector + saturatedStateCurve(const Opm::ECLUnits::UnitSystem& usys) const; }; +std::vector +WetGas::mainPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const +{ + auto curves = this->interp_.getPvtCurve(curve); + + const auto x_unit = usys.vaporisedOilGasRat(); + const auto y_unit = (curve == ::Opm::ECLPVT::RawCurve::FVF) + ? (usys.reservoirVolume() / usys.surfaceVolumeGas()) + : usys.viscosity(); + + for (auto& crv : curves) { + auto& x = crv.first; + auto& y = crv.second; + + for (auto n = x.size(), i = 0*n; i < n; ++i) { + x[i] = ::Opm::unit::convert::to(x[i], x_unit); + y[i] = ::Opm::unit::convert::to(y[i], y_unit); + } + } + + return curves; +} + +std::vector +WetGas::saturatedStateCurve(const Opm::ECLUnits::UnitSystem& usys) const +{ + const auto key_unit = usys.pressure(); + const auto rv_unit = usys.vaporisedOilGasRat(); + + auto rv = this->interp_.getSaturatedPoints(); + auto p = this->key_; + + if (rv.size() != p.size()) { + throw std::invalid_argument { + "Inconsistent table sizes of saturated gas function (RvSat(p))" + }; + } + + for (auto n = rv.size(), i = 0*n; i < n; ++i) { + p [i] = ::Opm::unit::convert::to(p [i], key_unit); + rv[i] = ::Opm::unit::convert::to(rv[i], rv_unit); + } + + auto graph = Opm::FlowDiagnostics::Graph { + std::move(p), std::move(rv) + }; + + return { std::move(graph) }; +} + // ##################################################################### namespace { @@ -391,9 +478,10 @@ public: return this->rhoS_[region]; } - FlowDiagnostics::Graph - getPvtCurve(const RegIdx region, - const RawCurve curve) const; + std::vector + getPvtCurve(const RegIdx region, + const RawCurve curve, + const ECLUnits::UnitSystem& usys) const; private: std::vector eval_; @@ -460,14 +548,15 @@ viscosity(const RegIdx region, return this->eval_[region]->viscosity(rv, pg); } -Opm::FlowDiagnostics::Graph +std::vector Opm::ECLPVT::Gas::Impl:: -getPvtCurve(const RegIdx region, - const RawCurve curve) const +getPvtCurve(const RegIdx region, + const RawCurve curve, + const ECLUnits::UnitSystem& usys) const { this->validateRegIdx(region); - return this->eval_[region]->getPvtCurve(curve); + return this->eval_[region]->getPvtCurve(curve, usys); } void @@ -543,11 +632,13 @@ double Opm::ECLPVT::Gas::surfaceMassDensity(const int region) const return this->pImpl_->surfaceMassDensity(region); } -Opm::FlowDiagnostics::Graph +std::vector Opm::ECLPVT::Gas:: -getPvtCurve(const RawCurve curve, const int region) const +getPvtCurve(const RawCurve curve, + const int region, + const ECLUnits::UnitSystem& usys) const { - return this->pImpl_->getPvtCurve(region, curve); + return this->pImpl_->getPvtCurve(region, curve, usys); } // ===================================================================== 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 df7f4dec39..b14a26838e 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 @@ -22,6 +22,7 @@ #include #include +#include #include #include @@ -161,8 +162,16 @@ namespace Opm { namespace ECLPVT { /// \param[in] region Region ID. Non-negative integer typically /// derived from the PVTNUM mapping vector. /// - /// \return 2D graph for PVT property curve identified by - /// requests represented by \p func and \p region. + /// \param[in] usys Unit system. Collection of units to which the + /// raw, tabulated PVT curve data will be converted. Usually + /// created by function \code ECLUnits::createUnitSystem() + /// \endcode or a similar facility. + /// + /// \return Collection of 2D graphs for PVT property curve + /// identified by requests represented by \p func and \p region. + /// One curve (vector element) for each pressure node. Single + /// curve (i.e., a single vector element) in the case of dry gas + /// (no vaporised oil). /// /// Example: Retrieve gas formation volume factor curve in PVT /// region 0 (zero based, i.e., cells for which PVTNUM==1). @@ -171,8 +180,10 @@ namespace Opm { namespace ECLPVT { /// const auto graph = /// pvtGas.getPvtCurve(ECLPVT::RawCurve::FVF, 0); /// \endcode - FlowDiagnostics::Graph - getPvtCurve(const RawCurve curve, const int region) const; + std::vector + getPvtCurve(const RawCurve curve, + const int region, + const ECLUnits::UnitSystem& usys) const; private: /// Implementation class. 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 bdeb488048..3924baad95 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 @@ -24,10 +24,13 @@ #include #include +#include + #include #include #include #include +#include #include #include #include @@ -92,8 +95,9 @@ public: viscosity(const std::vector& rs, const std::vector& po) const = 0; - virtual Opm::FlowDiagnostics::Graph - getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0; + virtual std::vector + getPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const = 0; virtual std::unique_ptr clone() const = 0; }; @@ -127,10 +131,31 @@ public: return this->interpolant_.viscosity(po); } - virtual Opm::FlowDiagnostics::Graph - getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override + virtual std::vector + getPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const override { - return this->interpolant_.getPvtCurve(curve); + if (curve == Opm::ECLPVT::RawCurve::SaturatedState) { + // Not applicable to dead oil. Return empty. + return { Opm::FlowDiagnostics::Graph{} }; + } + + auto pvtcurve = this->interpolant_.getPvtCurve(curve); + + const auto x_unit = usys.pressure(); + const auto y_unit = (curve == ::Opm::ECLPVT::RawCurve::FVF) + ? (usys.reservoirVolume() / usys.surfaceVolumeLiquid()) + : usys.viscosity(); + + auto& x = pvtcurve.first; + auto& y = pvtcurve.second; + + for (auto n = x.size(), i = 0*n; i < n; ++i) { + x[i] = ::Opm::unit::convert::to(x[i], x_unit); + y[i] = ::Opm::unit::convert::to(y[i], y_unit); + } + + return { std::move(pvtcurve) }; } virtual std::unique_ptr clone() const override @@ -155,7 +180,8 @@ public: LiveOil(std::vector key, std::vector propInterp) - : interp_(std::move(key), std::move(propInterp)) + : key_ (key) // Copy + , interp_(std::move(key), std::move(propInterp)) {} virtual std::vector @@ -184,13 +210,15 @@ public: return this->interp_.viscosity(key, x); } - virtual Opm::FlowDiagnostics::Graph - getPvtCurve(const Opm::ECLPVT::RawCurve /* curve */) const override + virtual std::vector + getPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const override { - throw std::runtime_error { - "Property Evaluator for Live Oil Does not " - "Support Retrieving Raw Curves (Blame BSKA)" - }; + if (curve == ::Opm::ECLPVT::RawCurve::SaturatedState) { + return this->saturatedStateCurve(usys); + } + + return this->mainPvtCurve(curve, usys); } virtual std::unique_ptr clone() const override @@ -201,9 +229,68 @@ public: private: using TableInterpolant = ::Opm::ECLPVT::PVTx; - TableInterpolant interp_; + std::vector key_; + TableInterpolant interp_; + + std::vector + mainPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const; + + std::vector + saturatedStateCurve(const Opm::ECLUnits::UnitSystem& usys) const; }; +std::vector +LiveOil::mainPvtCurve(const Opm::ECLPVT::RawCurve curve, + const Opm::ECLUnits::UnitSystem& usys) const +{ + auto curves = this->interp_.getPvtCurve(curve); + + const auto x_unit = usys.pressure(); + const auto y_unit = (curve == ::Opm::ECLPVT::RawCurve::FVF) + ? (usys.reservoirVolume() / usys.surfaceVolumeLiquid()) + : usys.viscosity(); + + for (auto& crv : curves) { + auto& x = crv.first; + auto& y = crv.second; + + for (auto n = x.size(), i = 0*n; i < n; ++i) { + x[i] = ::Opm::unit::convert::to(x[i], x_unit); + y[i] = ::Opm::unit::convert::to(y[i], y_unit); + } + } + + return curves; +} + +std::vector +LiveOil::saturatedStateCurve(const Opm::ECLUnits::UnitSystem& usys) const +{ + const auto press_unit = usys.pressure(); + const auto rs_unit = usys.dissolvedGasOilRat(); + + auto rs = this->key_; + auto p = this->interp_.getSaturatedPoints(); + + if (rs.size() != p.size()) { + throw std::invalid_argument { + "Inconsistent table sizes of saturated gas function (RsSat(p))" + }; + } + + for (auto n = rs.size(), i = 0*n; i < n; ++i) { + p [i] = ::Opm::unit::convert::to(p [i], press_unit); + rs[i] = ::Opm::unit::convert::to(rs[i], rs_unit); + } + + auto graph = Opm::FlowDiagnostics::Graph { + std::move(p), std::move(rs) + }; + + return { std::move(graph) }; +} + // ##################################################################### namespace { @@ -388,9 +475,10 @@ public: return this->rhoS_[region]; } - FlowDiagnostics::Graph - getPvtCurve(const RegIdx region, - const RawCurve curve) const; + std::vector + getPvtCurve(const RegIdx region, + const RawCurve curve, + const ECLUnits::UnitSystem& usys) const; private: std::vector eval_; @@ -457,14 +545,15 @@ viscosity(const RegIdx region, return this->eval_[region]->viscosity(rs, po); } -Opm::FlowDiagnostics::Graph +std::vector Opm::ECLPVT::Oil::Impl:: -getPvtCurve(const RegIdx region, - const RawCurve curve) const +getPvtCurve(const RegIdx region, + const RawCurve curve, + const ECLUnits::UnitSystem& usys) const { this->validateRegIdx(region); - return this->eval_[region]->getPvtCurve(curve); + return this->eval_[region]->getPvtCurve(curve, usys); } void @@ -540,11 +629,13 @@ double Opm::ECLPVT::Oil::surfaceMassDensity(const int region) const return this->pImpl_->surfaceMassDensity(region); } -Opm::FlowDiagnostics::Graph +std::vector Opm::ECLPVT::Oil:: -getPvtCurve(const RawCurve curve, const int region) const +getPvtCurve(const RawCurve curve, + const int region, + const ECLUnits::UnitSystem& usys) const { - return this->pImpl_->getPvtCurve(region, curve); + return this->pImpl_->getPvtCurve(region, curve, usys); } // ===================================================================== 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 ac2a91b4bf..5e0614f6d2 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 @@ -21,6 +21,7 @@ #define OPM_ECLPVTOIL_HEADER_INCLUDED #include +#include #include #include @@ -158,8 +159,16 @@ namespace Opm { namespace ECLPVT { /// \param[in] region Region ID. Non-negative integer typically /// derived from the PVTNUM mapping vector. /// - /// \return 2D graph for PVT property curve identified by - /// requests represented by \p func and \p region. + /// \param[in] usys Unit system. Collection of units to which the + /// raw, tabulated PVT curve data will be converted. Usually + /// created by function \code ECLUnits::createUnitSystem() + /// \endcode or a similar facility. + /// + /// \return Collection of 2D graphs for PVT property curve + /// identified by requests represented by \p func and \p region. + /// One curve (vector element) for each dissolved gas/oil ratio + /// node. Single curve (i.e., a single vector element) in the + /// case of dead oil (no dissolved gas). /// /// Example: Retrieve oil viscosity curve in PVT region 3 (zero /// based, i.e., those cells for which PVTNUM==4). @@ -168,8 +177,10 @@ namespace Opm { namespace ECLPVT { /// const auto graph = /// pvtOil.getPvtCurve(ECLPVT::RawCurve::Viscosity, 3); /// \endcode - FlowDiagnostics::Graph - getPvtCurve(const RawCurve curve, const int region) const; + std::vector + getPvtCurve(const RawCurve curve, + const int region, + const ECLUnits::UnitSystem& usys) const; private: /// Implementation class. 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 cec6fbc4c0..9f7c0a71df 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 @@ -94,11 +94,34 @@ namespace { return satnum; } + + Opm::FlowDiagnostics::Graph + transformOilCurve(const Opm::FlowDiagnostics::Graph& curve, + const double So_offset) + { + auto Sx = std::vector{}; Sx.reserve(curve.first.size()); + { + const auto& So = curve.first; + + std::transform(So.rbegin(), So.rend(), std::back_inserter(Sx), + [So_offset](const double so) + { + return So_offset - so; + }); + } + + auto y = std::vector{ + curve.second.rbegin(), + curve.second.rend() + }; + + return { std::move(Sx), std::move(y) }; + } } // Anonymous // ===================================================================== -namespace Relperm { +namespace { namespace Gas { namespace Details { Opm::ECLPropTableRawData @@ -109,12 +132,12 @@ namespace Relperm { unitConverter(const int usys); } - class KrFunction + class SatFunction { public: - KrFunction(const std::vector& tabdims, - const std::vector& tab, - const int usys) + SatFunction(const std::vector& tabdims, + const std::vector& tab, + const int usys) : func_(Details::tableData(tabdims, tab), Details::unitConverter(usys)) {} @@ -144,6 +167,16 @@ namespace Relperm { return this->func_.interpolate(t, c, sg); } + std::vector + pcgo(const std::size_t regID, + const std::vector& sg) const + { + const auto t = this->table(regID); + const auto c = this->pccol(); + + return this->func_.interpolate(t, c, sg); + } + const std::vector& saturationPoints(const std::size_t regID) const { @@ -163,6 +196,11 @@ namespace Relperm { { return ::Opm::SatFuncInterpolant::ResultColumn{0}; } + + Opm::SatFuncInterpolant::ResultColumn pccol() const + { + return ::Opm::SatFuncInterpolant::ResultColumn{1}; + } }; } // namespace Gas @@ -425,12 +463,12 @@ namespace Relperm { unitConverter(const int usys); } // namespace Details - class KrFunction + class SatFunction { public: - KrFunction(const std::vector& tabdims, - const std::vector& tab, - const int usys) + SatFunction(const std::vector& tabdims, + const std::vector& tab, + const int usys) : func_(Details::tableData(tabdims, tab), Details::unitConverter(usys)) {} @@ -460,6 +498,16 @@ namespace Relperm { return this->func_.interpolate(t, c, sw); } + std::vector + pcow(const std::size_t regID, + const std::vector& sw) const + { + const auto t = this->table(regID); + const auto c = this->pccol(); + + return this->func_.interpolate(t, c, sw); + } + const std::vector& saturationPoints(const std::size_t regID) const { @@ -479,13 +527,18 @@ namespace Relperm { { return ::Opm::SatFuncInterpolant::ResultColumn{0}; } + + Opm::SatFuncInterpolant::ResultColumn pccol() const + { + return ::Opm::SatFuncInterpolant::ResultColumn{1}; + } }; } // namespace Water -} // namespace Relperm +} // Anonymous Opm::ECLPropTableRawData -Relperm::Gas::Details::tableData(const std::vector& tabdims, - const std::vector& tab) +Gas::Details::tableData(const std::vector& tabdims, + const std::vector& tab) { auto t = Opm::ECLPropTableRawData{}; @@ -505,7 +558,7 @@ Relperm::Gas::Details::tableData(const std::vector& tabdims, } Opm::SatFuncInterpolant::ConvertUnits -Relperm::Gas::Details::unitConverter(const int usys) +Gas::Details::unitConverter(const int usys) { using CU = Opm::SatFuncInterpolant::ConvertUnits; using Cvrt = CU::Converter; @@ -530,9 +583,9 @@ Relperm::Gas::Details::unitConverter(const int usys) // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Opm::ECLPropTableRawData -Relperm::Oil::Details::tableData(const std::vector& tabdims, - const bool isTwoP, - const std::vector& tab) +Oil::Details::tableData(const std::vector& tabdims, + const bool isTwoP, + const std::vector& tab) { auto t = Opm::ECLPropTableRawData{}; @@ -556,7 +609,7 @@ Relperm::Oil::Details::tableData(const std::vector& tabdims, } Opm::SatFuncInterpolant::ConvertUnits -Relperm::Oil::Details::unitConverter(const bool isTwoP) +Oil::Details::unitConverter(const bool isTwoP) { using CU = Opm::SatFuncInterpolant::ConvertUnits; using Cvrt = CU::Converter; @@ -575,8 +628,8 @@ Relperm::Oil::Details::unitConverter(const bool isTwoP) // - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Opm::ECLPropTableRawData -Relperm::Water::Details::tableData(const std::vector& tabdims, - const std::vector& tab) +Water::Details::tableData(const std::vector& tabdims, + const std::vector& tab) { auto t = Opm::ECLPropTableRawData{}; @@ -596,7 +649,7 @@ Relperm::Water::Details::tableData(const std::vector& tabdims, } Opm::SatFuncInterpolant::ConvertUnits -Relperm::Water::Details::unitConverter(const int usys) +Water::Details::unitConverter(const int usys) { using CU = Opm::SatFuncInterpolant::ConvertUnits; using Cvrt = CU::Converter; @@ -667,12 +720,10 @@ private: const ECLInitFileData& init, const RawTEP& ep, const bool use3PtScaling, - const FuncCat curve, const ActPh& active) { auto opt = Create::EPSOptions{}; opt.use3PtScaling = use3PtScaling; - opt.curve = curve; if (active.oil) { this->create_oil_eps(G, init, ep, active, opt); @@ -687,56 +738,81 @@ private: } } - void scaleOG(const ECLRegionMapping& rmap, - std::vector& so) const + void scaleKrOG(const ECLRegionMapping& rmap, + std::vector& so) const { this->scale(this->oil_in_og_, rmap, so); } - void scaleOW(const ECLRegionMapping& rmap, - std::vector& so) const + void scaleKrOW(const ECLRegionMapping& rmap, + std::vector& so) const { this->scale(this->oil_in_ow_, rmap, so); } - void scaleGas(const ECLRegionMapping& rmap, - std::vector& sg) const + void scaleKrGas(const ECLRegionMapping& rmap, + std::vector& sg) const { - this->scale(this->gas_, rmap, sg); + this->scale(this->gas_.kr, rmap, sg); } - void scaleWat(const ECLRegionMapping& rmap, - std::vector& sw) const + void scaleKrWat(const ECLRegionMapping& rmap, + std::vector& sw) const { - this->scale(this->wat_, rmap, sw); + this->scale(this->wat_.kr, rmap, sw); } - void reverseScaleOG(const ECLRegionMapping& rmap, - std::vector& so) const + void scalePcGO(const ECLRegionMapping& rmap, + std::vector& sg) const + { + this->scale(this->gas_.pc, rmap, sg); + } + + void scalePcOW(const ECLRegionMapping& rmap, + std::vector& sw) const + { + this->scale(this->wat_.pc, rmap, sw); + } + + void reverseScaleKrOG(const ECLRegionMapping& rmap, + std::vector& so) const { this->reverseScale(this->oil_in_og_, rmap, so); } - void reverseScaleOW(const ECLRegionMapping& rmap, - std::vector& so) const + void reverseScaleKrOW(const ECLRegionMapping& rmap, + std::vector& so) const { this->reverseScale(this->oil_in_ow_, rmap, so); } - void reverseScaleGas(const ECLRegionMapping& rmap, - std::vector& sg) const + void reverseScaleKrGas(const ECLRegionMapping& rmap, + std::vector& sg) const { - this->reverseScale(this->gas_, rmap, sg); + this->reverseScale(this->gas_.kr, rmap, sg); } - void reverseScaleWat(const ECLRegionMapping& rmap, - std::vector& sw) const + void reverseScaleKrWat(const ECLRegionMapping& rmap, + std::vector& sw) const { - this->reverseScale(this->wat_, rmap, sw); + this->reverseScale(this->wat_.kr, rmap, sw); + } + + void reverseScalePcGO(const ECLRegionMapping& rmap, + std::vector& sg) const + { + this->reverseScale(this->gas_.pc, rmap, sg); + } + + void reverseScalePcOW(const ECLRegionMapping& rmap, + std::vector& sw) const + { + this->reverseScale(this->wat_.pc, rmap, sw); } private: using Create = ::Opm::SatFunc::CreateEPS; + using FCat = ::Opm::SatFunc::CreateEPS::FunctionCategory; using SSys = ::Opm::SatFunc::CreateEPS::SubSystem; using PhIdx = ::Opm::ECLPhaseIndex; @@ -791,10 +867,15 @@ private: EndPtsPtr tep; }; - EPS oil_in_og_; - EPS oil_in_ow_; - EPS gas_; - EPS wat_; + struct FullEPS { + EPS kr; + EPS pc; + }; + + EPS oil_in_og_; + EPS oil_in_ow_; + FullEPS gas_; + FullEPS wat_; void scale(const EPS& eps, const ECLRegionMapping& rmap, @@ -903,6 +984,7 @@ private: Create::EPSOptions& opt) { opt.thisPh = PhIdx::Liquid; + opt.curve = FCat::Relperm; // Oil EPS only applies to Kr if (active.gas) { opt.subSys = SSys::OilGas; @@ -931,10 +1013,33 @@ private: opt.thisPh = PhIdx::Vapour; opt.subSys = SSys::OilGas; - this->gas_.scaling = - Create::fromECLOutput(G, init, opt); + // Scaling for Gas relative permeabilty + { + opt.curve = FCat::Relperm; - this->gas_.tep = this->endPoints(ep, opt); + this->gas_.kr.scaling = + Create::fromECLOutput(G, init, opt); + + this->gas_.kr.tep = this->endPoints(ep, opt); + } + + // Scaling for Gas/Oil capillary pressure + { + // Save setting for Alternative/3-pt scaling. + // Capillary pressure uses two-point scaling exclusively. + const auto use3Pt = opt.use3PtScaling; + opt.use3PtScaling = false; + + opt.curve = FCat::CapPress; + + this->gas_.pc.scaling = + Create::fromECLOutput(G, init, opt); + + this->gas_.pc.tep = this->endPoints(ep, opt); + + // Restore original setting for Alternative scaling option. + opt.use3PtScaling = use3Pt; + } } void create_wat_eps(const ECLGraph& G, @@ -945,10 +1050,34 @@ private: opt.thisPh = PhIdx::Aqua; opt.subSys = SSys::OilWater; - this->wat_.scaling = - Create::fromECLOutput(G, init, opt); + // Scaling for Water relative permeabilty + { + opt.curve = FCat::Relperm; + + this->wat_.kr.scaling = + Create::fromECLOutput(G, init, opt); + + this->wat_.kr.tep = this->endPoints(ep, opt); + } + + // Scaling for Oil/Water capillary pressure + { + // Save setting for Alternative/3-pt scaling. + // Capillary pressure uses two-point scaling exclusively. + const auto use3Pt = opt.use3PtScaling; + opt.use3PtScaling = false; + + opt.curve = FCat::CapPress; + + this->wat_.pc.scaling = + Create::fromECLOutput(G, init, opt); + + this->wat_.pc.tep = this->endPoints(ep, opt); + + // Restore original setting for Alternative scaling option. + opt.use3PtScaling = use3Pt; + } - this->wat_.tep = this->endPoints(ep, opt); } EndPtsPtr @@ -964,9 +1093,9 @@ private: ECLRegionMapping rmap_; - std::unique_ptr oil_; - std::unique_ptr gas_; - std::unique_ptr wat_; + std::unique_ptr oil_; + std::unique_ptr gas_; + std::unique_ptr wat_; std::unique_ptr eps_; @@ -1002,6 +1131,12 @@ private: const std::vector& sg, const bool useEPS) const; + FlowDiagnostics::Graph + pcgoCurve(const ECLRegionMapping& rmap, + const std::size_t regID, + const std::vector& sg, + const bool useEPS) const; + std::vector krw(const ECLGraph& G, const ECLRestartData& rstrt, @@ -1013,18 +1148,32 @@ private: const std::vector& sw, const bool useEPS) const; - void scaleGasSat(const ECLRegionMapping& rmap, - const bool useEPS, - std::vector& sg) const; + FlowDiagnostics::Graph + pcowCurve(const ECLRegionMapping& rmap, + const std::size_t regID, + const std::vector& sw, + const bool useEPS) const; - void scaleOilSat(const ECLRegionMapping& rmap, - const bool useEPS, - std::vector& so_g, - std::vector& so_w) const; - - void scaleWaterSat(const ECLRegionMapping& rmap, + void scaleKrGasSat(const ECLRegionMapping& rmap, const bool useEPS, - std::vector& sw) const; + std::vector& sg) const; + + void scalePcGasSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& sg) const; + + void scaleKrOilSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& so_g, + std::vector& so_w) const; + + void scaleKrWaterSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& sw) const; + + void scalePcWaterSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& sw) const; void uniqueReverseScaleSat(std::vector& s) const; @@ -1133,16 +1282,16 @@ Impl::initRelPermInterp(const EPSEvaluator::ActPh& active, const auto& tab = init.keywordData("TAB"); if (active.gas) { - this->gas_.reset(new Relperm::Gas::KrFunction(tabdims, tab, usys)); + this->gas_.reset(new Gas::SatFunction(tabdims, tab, usys)); } if (active.wat) { - this->wat_.reset(new Relperm::Water::KrFunction(tabdims, tab, usys)); + this->wat_.reset(new Water::SatFunction(tabdims, tab, usys)); } if (active.oil) { if (! isThreePh) { - using KrModel = Relperm::Oil::TwoPhase; + using KrModel = Oil::TwoPhase; if (active.gas) { const auto subsys = KrModel::SubSys::OilGas; @@ -1162,7 +1311,7 @@ Impl::initRelPermInterp(const EPSEvaluator::ActPh& active, } if (isThreePh) { - using KrModel = Relperm::Oil::ECLStdThreePhase; + using KrModel = Oil::ECLStdThreePhase; this->oil_.reset(new KrModel(tabdims, tab, this->wat_->swco())); } @@ -1177,12 +1326,9 @@ Impl::initEPS(const EPSEvaluator::ActPh& active, { const auto ep = this->extractRawTableEndPoints(active); - const auto curve = ::Opm::SatFunc::CreateEPS:: - FunctionCategory::Relperm; - this->eps_.reset(new EPSEvaluator()); - this->eps_->define(G, init, ep, use3PtScaling, curve, active); + this->eps_->define(G, init, ep, use3PtScaling, active); } // ##################################################################### @@ -1196,11 +1342,11 @@ Opm::ECLSaturationFunc::Impl::Impl(const Impl& rhs) } if (rhs.gas_) { - this->gas_.reset(new Relperm::Gas::KrFunction(*rhs.gas_)); + this->gas_.reset(new Gas::SatFunction(*rhs.gas_)); } if (rhs.wat_) { - this->wat_.reset(new Relperm::Water::KrFunction(*rhs.wat_)); + this->wat_.reset(new Water::SatFunction(*rhs.wat_)); } } @@ -1253,8 +1399,8 @@ getSatFuncCurve(const std::vector& func, if (! (this->wat_ && (fi.subsys == RawCurve::SubSystem::OilWater))) { // Water is not an active phase, or we're being asked to - // extract the relative permeability curve for water in the - // O/G subsystem. No such curve. + // extract a saturation function curve for water in the O/G + // subsystem. No such curve. graph.emplace_back(); continue; } @@ -1272,10 +1418,13 @@ getSatFuncCurve(const std::vector& func, .reset(new ECLRegionMapping(this->satnum_, subset)); } - auto krw = this->krwCurve(*wat_assoc.second, regID, - wat_assoc.first, useEPS); + auto crv = (fi.curve == RawCurve::Function::RelPerm) + ? this->krwCurve (*wat_assoc.second, regID, + wat_assoc.first, useEPS) + : this->pcowCurve(*wat_assoc.second, regID, + wat_assoc.first, useEPS); - graph.push_back(std::move(krw)); + graph.push_back(std::move(crv)); } break; @@ -1300,10 +1449,15 @@ getSatFuncCurve(const std::vector& func, .reset(new ECLRegionMapping(this->satnum_, subset)); } - auto kro = this->kroCurve(*oil_assoc.second, regID, - fi.subsys, oil_assoc.first, useEPS); + const auto kro = + this->kroCurve(*oil_assoc.second, regID, + fi.subsys, oil_assoc.first, useEPS); - graph.push_back(std::move(kro)); + const auto So_off = (fi.subsys == RawCurve::SubSystem::OilGas) + ? oil_assoc.first.back() // Sg = Max{So} - So in G/O system + : 1.0; // Sw = 1.0 - So in O/W system + + graph.push_back(transformOilCurve(kro, So_off)); } break; @@ -1312,8 +1466,8 @@ getSatFuncCurve(const std::vector& func, if (! (this->gas_ && (fi.subsys == RawCurve::SubSystem::OilGas))) { // Gas is not an active phase, or we're being asked to - // extract the relative permeability curve for gas in the - // O/W subsystem. No such curve. + // extract a saturation function curve for gas in the O/W + // subsystem. No such curve. graph.emplace_back(); continue; } @@ -1331,10 +1485,13 @@ getSatFuncCurve(const std::vector& func, .reset(new ECLRegionMapping(this->satnum_, subset)); } - auto krg = this->krgCurve(*gas_assoc.second, regID, - gas_assoc.first, useEPS); + auto crv = (fi.curve == RawCurve::Function::RelPerm) + ? this->krgCurve (*gas_assoc.second, regID, + gas_assoc.first, useEPS) + : this->pcgoCurve(*gas_assoc.second, regID, + gas_assoc.first, useEPS); - graph.push_back(std::move(krg)); + graph.push_back(std::move(crv)); } break; @@ -1364,7 +1521,7 @@ kro(const ECLGraph& G, auto so_g = oil_saturation(sg, sw, G, rstrt); auto so_w = so_g; - this->scaleOilSat(this->rmap_, useEPS, so_g, so_w); + this->scaleKrOilSat(this->rmap_, useEPS, so_g, so_w); // Allocate result. Member function scatterRegionResult() depends on // having an allocated result vector into which to write the values from @@ -1377,20 +1534,20 @@ kro(const ECLGraph& G, (const int reg, const ECLRegionMapping& rmap) { - const auto So_g = Relperm::Oil::KrFunction::SOil { + const auto So_g = Oil::KrFunction::SOil { this->gatherRegionSubset(reg, rmap, so_g) }; - const auto So_w = Relperm::Oil::KrFunction::SOil { + const auto So_w = Oil::KrFunction::SOil { this->gatherRegionSubset(reg, rmap, so_w) }; - const auto Sg = Relperm::Oil::KrFunction::SGas { + const auto Sg = Oil::KrFunction::SGas { // Empty in case of Oil/Water system this->gatherRegionSubset(reg, rmap, sg) }; - const auto Sw = Relperm::Oil::KrFunction::SWat { + const auto Sw = Oil::KrFunction::SWat { // Empty in case of Oil/Gas system this->gatherRegionSubset(reg, rmap, sw) }; @@ -1421,8 +1578,8 @@ kroCurve(const ECLRegionMapping& rmap, auto so_g = so; auto so_w = so_g; if (useEPS && this->eps_) { - this->eps_->reverseScaleOG(rmap, so_g); - this->eps_->reverseScaleOW(rmap, so_w); + this->eps_->reverseScaleKrOG(rmap, so_g); + this->eps_->reverseScaleKrOW(rmap, so_w); this->uniqueReverseScaleSat(so_g); this->uniqueReverseScaleSat(so_w); @@ -1438,9 +1595,9 @@ kroCurve(const ECLRegionMapping& rmap, so_g.resize(so.size(), 1.0); so_w.resize(so.size(), 1.0); - this->scaleOilSat(rmap, useEPS, so_g, so_w); + this->scaleKrOilSat(rmap, useEPS, so_g, so_w); - const auto so_inp = Relperm::Oil::KrFunction::SOil { + const auto so_inp = Oil::KrFunction::SOil { (subsys == RawCurve::SubSystem::OilGas) ? so_g : so_w }; @@ -1470,7 +1627,7 @@ krg(const ECLGraph& G, auto sg = G.rawLinearisedCellData(rstrt, "SGAS"); if (useEPS && this->eps_) { - this->eps_->scaleGas(this->rmap_, sg); + this->eps_->scaleKrGas(this->rmap_, sg); } // Allocate result. Member function scatterRegionResult() depends on @@ -1510,7 +1667,7 @@ krgCurve(const ECLRegionMapping& rmap, auto sg_inp = sg; if (useEPS && this->eps_) { - this->eps_->reverseScaleGas(rmap, sg_inp); + this->eps_->reverseScaleKrGas(rmap, sg_inp); this->uniqueReverseScaleSat(sg_inp); } @@ -1520,7 +1677,7 @@ krgCurve(const ECLRegionMapping& rmap, // Re-expand input arrays to match size requirements of EPS evaluator. sg_inp.resize(sg.size(), 1.0); - this->scaleGasSat(rmap, useEPS, sg_inp); + this->scaleKrGasSat(rmap, useEPS, sg_inp); // Region ID 'reg' is traditional, ECL-style one-based region ID // (SATNUM). Subtract one to create valid table index. @@ -1533,6 +1690,42 @@ krgCurve(const ECLRegionMapping& rmap, }; } +Opm::FlowDiagnostics::Graph +Opm::ECLSaturationFunc::Impl:: +pcgoCurve(const ECLRegionMapping& rmap, + const std::size_t regID, + const std::vector& sg, + const bool useEPS) const +{ + if (! this->gas_) { + return FlowDiagnostics::Graph{}; + } + + auto sg_inp = sg; + if (useEPS && this->eps_) { + this->eps_->reverseScalePcGO(rmap, sg_inp); + this->uniqueReverseScaleSat(sg_inp); + } + + auto abscissas = sg_inp; + const auto nPoints = abscissas.size(); + + // Re-expand input arrays to match size requirements of EPS evaluator. + sg_inp.resize(sg.size(), 1.0); + + this->scalePcGasSat(rmap, useEPS, sg_inp); + + // 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); + + // FD::Graph == pair, vector> + return FlowDiagnostics::Graph { + std::move(abscissas), + { std::begin(pc), std::begin(pc) + nPoints } + }; +} + std::vector Opm::ECLSaturationFunc::Impl:: krw(const ECLGraph& G, @@ -1548,7 +1741,7 @@ krw(const ECLGraph& G, auto sw = G.rawLinearisedCellData(rstrt, "SWAT"); if (useEPS && this->eps_) { - this->eps_->scaleWat(this->rmap_, sw); + this->eps_->scaleKrWat(this->rmap_, sw); } // Allocate result. Member function scatterRegionResult() depends on @@ -1588,7 +1781,7 @@ krwCurve(const ECLRegionMapping& rmap, auto sw_inp = sw; if (useEPS && this->eps_) { - this->eps_->reverseScaleWat(rmap, sw_inp); + this->eps_->reverseScaleKrWat(rmap, sw_inp); this->uniqueReverseScaleSat(sw_inp); } @@ -1598,7 +1791,7 @@ krwCurve(const ECLRegionMapping& rmap, // Re-expand input arrays to match size requirements of EPS evaluator. sw_inp.resize(sw.size(), 1.0); - this->scaleWaterSat(rmap, useEPS, sw_inp); + this->scaleKrWaterSat(rmap, useEPS, sw_inp); // Region ID 'reg' is traditional, ECL-style one-based region ID // (SATNUM). Subtract one to create valid table index. @@ -1611,41 +1804,99 @@ krwCurve(const ECLRegionMapping& rmap, }; } +Opm::FlowDiagnostics::Graph +Opm::ECLSaturationFunc::Impl:: +pcowCurve(const ECLRegionMapping& rmap, + const std::size_t regID, + const std::vector& sw, + const bool useEPS) const +{ + if (! this->wat_) { + return FlowDiagnostics::Graph{}; + } + + auto sw_inp = sw; + if (useEPS && this->eps_) { + this->eps_->reverseScalePcOW(rmap, sw_inp); + this->uniqueReverseScaleSat(sw_inp); + } + + auto abscissas = sw_inp; + const auto nPoints = abscissas.size(); + + // Re-expand input arrays to match size requirements of EPS evaluator. + sw_inp.resize(sw.size(), 1.0); + + this->scalePcWaterSat(rmap, useEPS, sw_inp); + + // 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); + + // FD::Graph == pair, vector> + return FlowDiagnostics::Graph { + std::move(abscissas), + { std::begin(pc), std::begin(pc) + nPoints } + }; +} + void Opm::ECLSaturationFunc::Impl:: -scaleGasSat(const ECLRegionMapping& rmap, - const bool useEPS, - std::vector& sg) const +scaleKrGasSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& sg) const { if (useEPS && this->eps_) { - this->eps_->scaleGas(rmap, sg); + this->eps_->scaleKrGas(rmap, sg); } } void Opm::ECLSaturationFunc::Impl:: -scaleOilSat(const ECLRegionMapping& rmap, - const bool useEPS, - std::vector& so_g, - std::vector& so_w) const +scalePcGasSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& sg) const +{ + if (useEPS && this->eps_) { + this->eps_->scalePcGO(rmap, sg); + } +} + +void +Opm::ECLSaturationFunc::Impl:: +scaleKrOilSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& so_g, + std::vector& so_w) const { if (useEPS && this->eps_) { // Independent scaling of So in O/G and O/W sub-systems of an O/G/W // run. Performs duplicate work in a two-phase case. Need to take // action if this becomes a bottleneck. - this->eps_->scaleOG(rmap, so_g); - this->eps_->scaleOW(rmap, so_w); + this->eps_->scaleKrOG(rmap, so_g); + this->eps_->scaleKrOW(rmap, so_w); } } void Opm::ECLSaturationFunc::Impl:: -scaleWaterSat(const ECLRegionMapping& rmap, - const bool useEPS, - std::vector& sg) const +scaleKrWaterSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& sg) const { if (useEPS && this->eps_) { - this->eps_->scaleWat(rmap, sg); + this->eps_->scaleKrWat(rmap, sg); + } +} + +void +Opm::ECLSaturationFunc::Impl:: +scalePcWaterSat(const ECLRegionMapping& rmap, + const bool useEPS, + std::vector& sg) const +{ + if (useEPS && this->eps_) { + this->eps_->scalePcOW(rmap, sg); } } 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 72e265152b..be31c7e2f4 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 @@ -52,6 +52,9 @@ namespace Opm { enum class Function { /// Relative permeability functions RelPerm, + + /// Capillary pressure + CapPress, }; /// Which one-dimensional sub-system does this request reference.