From dc8931f78290e8c6463c77642c4d4bea9947e10c Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Tue, 23 Oct 2018 15:58:30 +0200 Subject: [PATCH 1/3] Update Flow Diagnostics Application Library Fixes incorrect horizontal and vertical end-point scaling of model's saturation functions. API Change: No longer supports user-selected behaviour for treating scaled end-points with a sentinel value (-1.0E+20). That option was introduced due to incomplete understanding of the semantics of the sentinel value. Now that we understand the meaning (use actual, unscaled end-point value from input table), we no longer need the option. Update the calling code in RigFlowDiagSolverInterface.cpp accordingly. --- .../RigFlowDiagSolverInterface.cpp | 1 - .../custom-opm-flowdiag-app/CMakeLists.txt | 1 - .../examples/extractPropCurves.cpp | 67 +--- .../opm/utility/ECLEndPointScaling.cpp | 223 ++++++------ .../opm/utility/ECLEndPointScaling.hpp | 20 -- .../opm/utility/ECLPvtGas.cpp | 2 + .../opm/utility/ECLPvtOil.cpp | 2 + .../opm/utility/ECLResultData.cpp | 24 +- .../opm/utility/ECLSaturationFunc.cpp | 3 +- .../opm/utility/ECLSaturationFunc.hpp | 9 +- .../opm/parser/eclipse/Units/Units.hpp | 326 ------------------ 11 files changed, 135 insertions(+), 543 deletions(-) delete mode 100644 ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp diff --git a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp index c943688b2e..50d3d9bad8 100644 --- a/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp +++ b/ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp @@ -682,7 +682,6 @@ std::vector RigFlowDiagSolverInterface if (!useEps) { scaling.enable = static_cast(0); } - scaling.invalid = Opm::SatFunc::EPSEvalInterface::InvalidEndpointBehaviour::IgnorePoint; std::vector graphArr = m_opmFlowDiagStaticData->m_eclSaturationFunc->getSatFuncCurve(satFuncRequests, static_cast(activeCellIndex), scaling); for (size_t i = 0; i < graphArr.size(); i++) { diff --git a/ThirdParty/custom-opm-flowdiag-app/CMakeLists.txt b/ThirdParty/custom-opm-flowdiag-app/CMakeLists.txt index c37bd33e02..d54ae136a9 100644 --- a/ThirdParty/custom-opm-flowdiag-app/CMakeLists.txt +++ b/ThirdParty/custom-opm-flowdiag-app/CMakeLists.txt @@ -5,7 +5,6 @@ project (custom-opm-flowdiag-app) include_directories( ../custom-opm-flowdiagnostics/opm-flowdiagnostics opm-flowdiagnostics-applications - opmCore ) include (opm-flowdiagnostics-applications/CMakeLists_files.cmake) 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 cbfbfed881..acf60e4dd4 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 @@ -54,13 +54,13 @@ namespace { const auto& x = graph.first; const auto& y = graph.second; - os << name << '{' << k << "} = [\n"; + os << name << '{' << k << "} = extendTab([\n"; for (auto n = x.size(), i = 0*n; i < n; ++i) { os << x[i] << ' ' << y[i] << '\n'; } - os << "];\n\n"; + os << "]);\n\n"; k += 1; } @@ -118,7 +118,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, scaling); - printGraph(std::cout, "krg", graph); + printGraph(std::cout, "crv.krg", graph); } void krog(const Opm::ECLSaturationFunc& sfunc, @@ -140,7 +140,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, scaling); - printGraph(std::cout, "krog", graph); + printGraph(std::cout, "crv.krog", graph); } void krow(const Opm::ECLSaturationFunc& sfunc, @@ -162,7 +162,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, scaling); - printGraph(std::cout, "krow", graph); + printGraph(std::cout, "crv.krow", graph); } void krw(const Opm::ECLSaturationFunc& sfunc, @@ -184,7 +184,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, scaling); - printGraph(std::cout, "krw", graph); + printGraph(std::cout, "crv.krw", graph); } // ----------------------------------------------------------------- @@ -209,7 +209,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, scaling); - printGraph(std::cout, "pcgo", graph); + printGraph(std::cout, "crv.pcgo", graph); } void pcow(const Opm::ECLSaturationFunc& sfunc, @@ -231,7 +231,7 @@ namespace { const auto graph = sfunc.getSatFuncCurve(func, activeCell, scaling); - printGraph(std::cout, "pcow", graph); + printGraph(std::cout, "crv.pcow", graph); } // ----------------------------------------------------------------- @@ -245,7 +245,7 @@ namespace { const auto graph = pvtCurves .getPvtCurve(RC::FVF, Opm::ECLPhaseIndex::Vapour, activeCell); - printGraph(std::cout, "Bg", graph); + printGraph(std::cout, "crv.Bg", graph); } void mu_g(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, @@ -256,7 +256,7 @@ namespace { const auto graph = pvtCurves .getPvtCurve(RC::Viscosity, Opm::ECLPhaseIndex::Vapour, activeCell); - printGraph(std::cout, "mu_g", graph); + printGraph(std::cout, "crv.mu_g", graph); } void Bo(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, @@ -267,7 +267,7 @@ namespace { const auto graph = pvtCurves .getPvtCurve(RC::FVF, Opm::ECLPhaseIndex::Liquid, activeCell); - printGraph(std::cout, "Bo", graph); + printGraph(std::cout, "crv.Bo", graph); } void mu_o(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, @@ -278,7 +278,7 @@ namespace { const auto graph = pvtCurves .getPvtCurve(RC::Viscosity, Opm::ECLPhaseIndex::Liquid, activeCell); - printGraph(std::cout, "mu_o", graph); + printGraph(std::cout, "crv.mu_o", graph); } // ----------------------------------------------------------------- @@ -293,7 +293,7 @@ namespace { const auto graph = pvtCurves .getPvtCurve(RC::SaturatedState, PI::Vapour, activeCell); - printGraph(std::cout, "rvSat", graph); + printGraph(std::cout, "crv.rvSat", graph); } void rsSat(const Opm::ECLPVT::ECLPvtCurveCollection& pvtCurves, @@ -305,7 +305,7 @@ namespace { const auto graph = pvtCurves .getPvtCurve(RC::SaturatedState, PI::Liquid, activeCell); - printGraph(std::cout, "rsSat", graph); + printGraph(std::cout, "crv.rsSat", graph); } // ----------------------------------------------------------------- @@ -340,37 +340,6 @@ namespace { return Opm::ECLUnits::serialisedUnitConventions(init); } - Opm::SatFunc::EPSEvalInterface::InvalidEndpointBehaviour - handleInvalid(const std::string& behaviour) - { - using IEB = Opm::SatFunc:: - EPSEvalInterface::InvalidEndpointBehaviour; - - if ((behaviour == "ignore") || - (behaviour == "ignore-point") || - (behaviour == "ignore_point") || - (behaviour == "ignorepoint")) - { - return IEB::IgnorePoint; - } - - return IEB::UseUnscaled; - } - - auto handleInvalid(const Opm::ParameterGroup& prm) - -> decltype(handleInvalid("ignore")) - { - for (const auto* param : { "handle_invalid" , - "hInv", "handleInv" }) - { - if (prm.has(param)) { - return handleInvalid(prm.get(param)); - } - } - - return handleInvalid("ignore"); - } - int getActiveCell(const Opm::ECLGraph& G, const Opm::ParameterGroup& prm) { @@ -424,8 +393,7 @@ namespace { prm.getDefault("useEPS", std::string{"off"}); auto scaling = Opm::ECLSaturationFunc::SatFuncScaling{}; - scaling.enable = static_cast(0); - scaling.invalid = handleInvalid(prm); + scaling.enable = static_cast(0); if (std::regex_search(useEPS, horiz)) { scaling.enable |= T::Horizontal; @@ -460,6 +428,8 @@ try { pvtCC.setOutputUnits(std::move(units)); } + std::cout << "function crv = pcurves()\n"; + // ----------------------------------------------------------------- // Relative permeability @@ -470,7 +440,8 @@ try { // ----------------------------------------------------------------- // Capillary pressure - if (prm.getDefault("pcgo", false)) { pcgo(sfunc, cellID, scaling); } + if (prm.getDefault("pcog", false) || // Alias pcog -> pcgo + prm.getDefault("pcgo", false)) { pcgo(sfunc, cellID, scaling); } if (prm.getDefault("pcow", false)) { pcow(sfunc, cellID, scaling); } // ----------------------------------------------------------------- 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 47b5024db6..6b16d2c88b 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 @@ -1,4 +1,5 @@ /* + Copyright 2018 Equinor ASA. Copyright 2017 Statoil ASA. This file is part of the Open Porous Media Project (OPM). @@ -34,8 +35,8 @@ #include #include #include -#include #include +#include #include #include #include @@ -217,21 +218,6 @@ namespace { return (std::abs(s) < 1.0e+20) ? s : dflt; } - 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); - }); - } - bool haveScaledRelPermAtCritSat(const ::Opm::ECLGraph& G, const ::Opm::ECLInitFileData& init, @@ -263,9 +249,7 @@ class Opm::SatFunc::TwoPointScaling::Impl public: Impl(std::vector smin, std::vector smax) - : smin_ (std::move(smin)) - , smax_ (std::move(smax)) - , handle_invalid_(InvalidEndpointBehaviour::UseUnscaled) + : smin_(std::move(smin)), smax_(std::move(smax)) { if (this->smin_.size() != this->smax_.size()) { throw std::invalid_argument { @@ -287,11 +271,6 @@ private: std::vector smin_; std::vector smax_; - InvalidEndpointBehaviour handle_invalid_; - - void handleInvalidEndpoint(const SaturationAssoc& sp, - std::vector& effsat) const; - double sMin(const std::vector::size_type cell, const TableEndPoints& tep) const { @@ -327,12 +306,6 @@ Impl::eval(const TableEndPoints& tep, const auto sLO = this->sMin(cell, tep); const auto sHI = this->sMax(cell, tep); - if (! validSaturations({ sLO, sHI })) { - this->handleInvalidEndpoint(eval_pt, effsat); - - continue; - } - effsat.push_back(0.0); auto& s_eff = effsat.back(); @@ -372,12 +345,6 @@ Impl::reverse(const TableEndPoints& tep, const auto sLO = this->sMin(cell, tep); const auto sHI = this->sMax(cell, tep); - if (! validSaturations({ sLO, sHI })) { - this->handleInvalidEndpoint(eval_pt, unscaledsat); - - continue; - } - unscaledsat.push_back(0.0); auto& s_unsc = unscaledsat.back(); @@ -404,26 +371,6 @@ 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; - } - - if (this->handle_invalid_ == InvalidEndpointBehaviour::IgnorePoint) { - // User requests that invalid scaling be ignored. Signal invalid - // scaled saturation to caller as NaN. - effsat.push_back(std::nan("")); - return; - } -} - // --------------------------------------------------------------------- // Class Opm::SatFunc::PureVerticalScaling::Impl // --------------------------------------------------------------------- @@ -433,6 +380,9 @@ class Opm::SatFunc::PureVerticalScaling::Impl public: explicit Impl(std::vector fmax) : fmax_(std::move(fmax)) + , inf_ (std::numeric_limits::has_infinity + ? std::numeric_limits::infinity() + : std::numeric_limits::max()) {} std::vector @@ -442,6 +392,16 @@ public: private: std::vector fmax_; + + double inf_; + + std::vector + zeroMaxVal(const SaturationPoints& sp) const; + + std::vector + nonZeroMaxVal(const SaturationPoints& sp, + const double maxVal, + std::vector val) const; }; std::vector @@ -452,10 +412,40 @@ vertScale(const FunctionValues& f, { assert (sp.size() == val.size() && "Internal Error in Vertical Scaling"); - auto ret = std::move(val); - const auto maxVal = f.max.val; + if (! (std::abs(maxVal) > 0.0)) { + return this->zeroMaxVal(sp); + } + + return this->nonZeroMaxVal(sp, maxVal, std::move(val)); +} + +std::vector +Opm::SatFunc::PureVerticalScaling::Impl:: +zeroMaxVal(const SaturationPoints& sp) const +{ + auto ret = std::vector(sp.size(), 0.0); + + for (auto n = sp.size(), i = 0*n; i < n; ++i) { + const auto fmax = this->fmax_[ sp[i].cell ]; + + ret[i] = (std::abs(fmax) > 0.0) + ? (std::signbit(fmax) ? -this->inf_ : this->inf_) + : 0.0; + } + + return ret; +} + +std::vector +Opm::SatFunc::PureVerticalScaling::Impl:: +nonZeroMaxVal(const SaturationPoints& sp, + const double maxVal, + std::vector val) const +{ + auto ret = std::move(val); + for (auto n = sp.size(), i = 0*n; i < n; ++i) { ret[i] *= this->fmax_[ sp[i].cell ] / maxVal; } @@ -513,24 +503,29 @@ vertScale(const FunctionValues& f, const auto c = sp[i].cell; const auto s = sp[i].sat; - const auto sr = this->sdisp_[c]; + const auto sr = std::min(this->sdisp_[c], this->smax_[c]); const auto fr = this->fdisp_[c]; const auto sm = this->smax_ [c]; const auto fm = this->fmax_ [c]; - if (! (s > sr)) { - // s <= sr: Pure vertical scaling in left interval. + if (s < sr) { + // s < sr: Pure vertical scaling in left interval. y *= fr / fdisp; } else if (sepfv) { - // s > sr; normal case: Kr(Smax) > Kr(Sr) + // s \in [sr, sm), sm > sr; normal case: Kr(Smax) > Kr(Sr). + // + // Linear function between (sr,fr) and (sm,fm) in terms of + // function value 'y'. const auto t = (y - fdisp) / (fmax - fdisp); y = fr + t*(fm - fr); } else if (s < sm) { - // s > sr; special case: Kr(Smax) == Kr(Sr) in table. Use - // linear function between (sr,fr) and (sm,fm). + // s \in [sr, sm), sm > sr; special case: Kr(Smax) == Kr(Sr). + // + // Use linear function between (sr,fr) and (sm,fm) in terms of + // saturation value 's'. const auto t = (s - sr) / (sm - sr); y = fr + t*(fm - fr); @@ -554,10 +549,7 @@ public: Impl(std::vector smin , std::vector sdisp, std::vector smax ) - : smin_ (std::move(smin )) - , sdisp_ (std::move(sdisp)) - , smax_ (std::move(smax )) - , handle_invalid_(InvalidEndpointBehaviour::UseUnscaled) + : smin_(std::move(smin)), sdisp_(std::move(sdisp)), smax_(std::move(smax)) { if ((this->sdisp_.size() != this->smin_.size()) || (this->sdisp_.size() != this->smax_.size())) @@ -582,11 +574,6 @@ private: std::vector sdisp_; std::vector smax_; - InvalidEndpointBehaviour handle_invalid_; - - void handleInvalidEndpoint(const SaturationAssoc& sp, - std::vector& effsat) const; - double sMin(const std::vector::size_type cell, const TableEndPoints& tep) const { @@ -630,12 +617,6 @@ Impl::eval(const TableEndPoints& tep, const auto sR = this->sDisp(cell, tep); const auto sHI = this->sMax (cell, tep); - if (! validSaturations({ sLO, sR, sHI })) { - this->handleInvalidEndpoint(eval_pt, effsat); - - continue; - } - effsat.push_back(0.0); auto& s_eff = effsat.back(); @@ -643,22 +624,26 @@ Impl::eval(const TableEndPoints& tep, // s <= sLO s_eff = tep.low; } - else if (! (eval_pt.sat < sHI)) { - // s >= sHI - s_eff = tep.high; - } - else if (eval_pt.sat < sR) { - // s \in (sLO, sR) + else if (eval_pt.sat < std::min(sR, sHI)) { + // s in scaled interval [sLO, sR) + // Map to tabulated saturation in [tep.low, tep.disp) const auto t = (eval_pt.sat - sLO) / (sR - sLO); s_eff = tep.low + t*(tep.disp - tep.low); } - else { - // s \in (sR, sHI) + else if (eval_pt.sat < sHI) { + // s in scaled interval [sR, sHI) + // Map to tabulated saturation in [tep.disp, tep.high) + assert (sHI > sR); + const auto t = (eval_pt.sat - sR) / (sHI - sR); s_eff = tep.disp + t*(tep.high - tep.disp); } + else { + // s >= sHI + s_eff = tep.high; + } } return effsat; @@ -679,12 +664,6 @@ Impl::reverse(const TableEndPoints& tep, const auto sR = this->sDisp(cell, tep); const auto sHI = this->sMax (cell, tep); - if (! validSaturations({ sLO, sR, sHI })) { - this->handleInvalidEndpoint(eval_pt, unscaledsat); - - continue; - } - unscaledsat.push_back(0.0); auto& s_unsc = unscaledsat.back(); @@ -693,54 +672,45 @@ Impl::reverse(const TableEndPoints& tep, // Map to Minimum Input Saturation in cell (sLO). s_unsc = sLO; } - else if (! (eval_pt.sat < tep.high)) { - // s >= maximum tabulated saturation. - // Map to Maximum Input Saturation in cell (sHI). - s_unsc = sHI; - } else if (eval_pt.sat < tep.disp) { - // s in tabulated interval (tep.low, tep.disp) + // s in tabulated interval [tep.low, tep.disp) // Map to Input Saturation in (sLO, sR) const auto t = (eval_pt.sat - tep.low) / (tep.disp - tep.low); - s_unsc = sLO + t*(sR - sLO); + s_unsc = std::min(sLO + t*(sR - sLO), sHI); } - else { - // s in tabulated interval (tep.disp, tep.high) - // Map to Input Saturation in (sR, sHI) + else if (eval_pt.sat < tep.high) { + // s in tabulated interval [tep.disp, tep.high) + // Map to Input Saturation in [sR, sHI) + assert (tep.high > tep.disp); + const auto t = (eval_pt.sat - tep.disp) / (tep.high - tep.disp); - s_unsc = sR + t*(sHI - sR); + s_unsc = std::min(sR + t*std::max(sHI - sR, 0.0), sHI); + } + else { + // s >= maximum tabulated saturation. + // + // Map to Maximum Input Saturation in cell (sHI) if maximum + // tabulated saturation is strictly greater than critical + // displacing saturation--otherwise map to critical displacing + // saturation. + // + // Needed to handle cases in which \code tep.disp==tep.high + // \endcode but scaled versions of these might differ, i.e. when + // sR < sHI, but the corresponding saturation points in the + // underlying input table coincide. + s_unsc = (tep.high > tep.disp) ? sHI : sR; } } 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; - } - - if (this->handle_invalid_ == InvalidEndpointBehaviour::IgnorePoint) { - // User requests that invalid scaling be ignored. Signal invalid - // scaled saturation to caller as NaN. - effsat.push_back(std::nan("")); - return; - } -} - // --------------------------------------------------------------------- // EPS factory functions for two-point and three-point scaling options // --------------------------------------------------------------------- @@ -749,7 +719,6 @@ 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; @@ -1742,7 +1711,7 @@ KrGO::twoPointMethod(const ::Opm::ECLGraph& G, auto t = std::vector(sogcr.size(), 0.0); { const auto& sgco = tep.conn.gas; - const auto& swco = tep.conn.gas; + const auto& swco = tep.conn.water; const auto& sgcr = tep.crit.gas; for (auto n = sgcr.size(), i = 0*n; i < n; ++i) { @@ -1790,7 +1759,7 @@ KrOW::twoPointMethod(const ::Opm::ECLGraph& G, auto t = std::vector(sowcr.size(), 0.0); { const auto& sgco = tep.conn.gas; - const auto& swco = tep.conn.gas; + const auto& swco = tep.conn.water; const auto& swcr = tep.crit.water; for (auto n = swcr.size(), i = 0*n; i < n; ++i) { 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 384ce66c97..4fbc84aba2 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,19 +67,6 @@ 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; @@ -574,13 +561,6 @@ 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/ECLPvtGas.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp index fef7dd5140..c712476ee3 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 @@ -90,6 +90,8 @@ namespace { class PVxGBase { public: + virtual ~PVxGBase() {} + virtual std::vector formationVolumeFactor(const std::vector& rv, const std::vector& pg) const = 0; 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 362ae167f8..7ace98f6c7 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 @@ -129,6 +129,8 @@ namespace { class PVxOBase { public: + virtual ~PVxOBase() {} + virtual std::vector formationVolumeFactor(const std::vector& rs, const std::vector& po) const = 0; diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp index d2e977cd5f..d0effa9510 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLResultData.cpp @@ -473,13 +473,13 @@ namespace { explicit InitFileSections(const ecl_file_type* init); struct Section { - Section(const ecl_file_view_type* blk) + Section(ecl_file_view_type* blk) : block (blk) , first_kw(Details::firstBlockKeyword(block)) {} - const ecl_file_view_type* block; - std::string first_kw; + ecl_file_view_type* block; + std::string first_kw; }; const std::vector
& sections() const @@ -502,7 +502,7 @@ namespace { } private: - const ecl_file_view_type* init_; + ecl_file_view_type* init_; std::vector
sect_; }; @@ -796,7 +796,7 @@ ECLImpl::InitFileSections::InitFileSections(const ecl_file_type* init) ? sectID - 1 : 0*sectID; // Main section 'sectID': [ start_kw, LGRSGONE ] - const auto* sect = + auto* sect = ecl_file_view_add_blockview2(this->init_, start_kw, end_kw, start_kw_occurrence); @@ -808,12 +808,12 @@ ECLImpl::InitFileSections::InitFileSections(const ecl_file_type* init) const auto occurrence = 0; // Main grid sub-section of 'sectID': [ start_kw, LGR ] - const auto* main_grid_sect = + auto* main_grid_sect = ecl_file_view_add_blockview2(sect, firstkw.c_str(), LGR_KW, occurrence); // LGR sub-section of 'sectID': [ LGR, LGRSGONE ] - const auto* lgr_sect = + auto* lgr_sect = ecl_file_view_add_blockview2(sect, LGR_KW, end_kw, occurrence); @@ -965,7 +965,7 @@ private: /// Saved original active view from host (prior to restricting view /// to single grid ID). - const ecl_file_view_type* save_; + ecl_file_view_type* save_; }; /// Casename prefix. Mostly to implement copy ctor. @@ -986,7 +986,7 @@ private: std::unique_ptr gridIDCache_; /// Current active result-set view. - mutable const ecl_file_view_type* activeBlock_{ nullptr }; + mutable ecl_file_view_type* activeBlock_{ nullptr }; /// Support for passing \code *this \endcode to ERT functions that /// require an \c ecl_file_type, particularly the function that selects @@ -1276,7 +1276,7 @@ private: /// Sections of the INIT result set. ECLImpl::InitFileSections sections_; - mutable const ecl_file_view_type* activeBlock_{ nullptr }; + mutable ecl_file_view_type* activeBlock_{ nullptr }; /// Negative look-up cache for haveKeywordData() queries. mutable MissingKW missing_kw_; @@ -1415,6 +1415,10 @@ lookup(const std::string& vector, const std::string& gridName) const } } + // Status of 'vector' unknown for 'gridName'. Actually look for the + // vector in gridName's sections (main grid if gridName.empty(), + // otherwise named local grid). + if (gridName.empty()) { return this->lookupMainGrid(key); } 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 d7c8c3a72b..0465752cc2 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 @@ -758,8 +758,7 @@ private: const ActPh& active) { auto opt = Create::EPSOptions{}; - opt.use3PtScaling = use3PtScaling; - opt.handle_invalid = SatFuncScaling::IEB::UseUnscaled; + opt.use3PtScaling = use3PtScaling; if (active.oil) { this->create_oil_eps(host, G, init, ep, active, opt); 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 ff9926d0aa..7b94f868c6 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 @@ -85,19 +85,12 @@ namespace Opm { Vertical = 1 << 1u, }; - using IEB = SatFunc::EPSEvalInterface::InvalidEndpointBehaviour; - SatFuncScaling() - : enable (Type::Horizontal | Type::Vertical) - , invalid(IEB::UseUnscaled) + : enable(Type::Horizontal | Type::Vertical) {} // Default: Use both Horizontal and Vertical if specified. unsigned char enable; - - // Default: Use unscaled values in case of invalid scaled - // saturations occurring in the result set. - IEB invalid; }; /// Constructor diff --git a/ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp b/ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp deleted file mode 100644 index cc1dde7e66..0000000000 --- a/ThirdParty/custom-opm-flowdiag-app/opmCore/opm/parser/eclipse/Units/Units.hpp +++ /dev/null @@ -1,326 +0,0 @@ -//=========================================================================== -// -// File: Units.hpp -// -// Created: Thu Jul 2 09:19:08 2009 -// -// Author(s): Halvor M Nilsen -// -// $Date$ -// -// $Revision$ -// -//=========================================================================== - -/* -Copyright 2009, 2010, 2011, 2012 SINTEF ICT, Applied Mathematics. -Copyright 2009, 2010, 2011, 2012 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 . -*/ - -#ifndef OPM_UNITS_HEADER -#define OPM_UNITS_HEADER - -/** -The unit sets emplyed in ECLIPSE, in particular the FIELD units, -are quite inconsistent. Ideally one should choose units for a set -of base quantities like Mass,Time and Length and then derive the -units for e.g. pressure and flowrate in a consisten manner. However -that is not the case; for instance in the metric system we have: - -[Length] = meters -[time] = days -[mass] = kg - -This should give: - -[Pressure] = [mass] / ([length] * [time]^2) = kg / (m * days * days) - -Instead pressure is given in Bars. When it comes to FIELD units the -number of such examples is long. -*/ -namespace Opm { -namespace prefix - /// Conversion prefix for units. -{ -constexpr const double micro = 1.0e-6; /**< Unit prefix [\f$\mu\f$] */ -constexpr const double milli = 1.0e-3; /**< Unit prefix [m] */ -constexpr const double centi = 1.0e-2; /**< Non-standard unit prefix [c] */ -constexpr const double deci = 1.0e-1; /**< Non-standard unit prefix [d] */ -constexpr const double kilo = 1.0e3; /**< Unit prefix [k] */ -constexpr const double mega = 1.0e6; /**< Unit prefix [M] */ -constexpr const double giga = 1.0e9; /**< Unit prefix [G] */ -} // namespace prefix - -namespace unit - /// Definition of various units. - /// All the units are defined in terms of international standard - /// units (SI). Example of use: We define a variable \c k which - /// gives a permeability. We want to set \c k to \f$1\,mD\f$. - /// \code - /// using namespace Opm::unit - /// double k = 0.001*darcy; - /// \endcode - /// We can also use one of the prefixes defined in Opm::prefix - /// \code - /// using namespace Opm::unit - /// using namespace Opm::prefix - /// double k = 1.0*milli*darcy; - /// \endcode -{ -///\name Common powers -/// @{ -constexpr double square(double v) { return v * v; } -constexpr double cubic (double v) { return v * v * v; } -/// @} - -// -------------------------------------------------------------- -// Basic (fundamental) units and conversions -// -------------------------------------------------------------- - -/// \name Length -/// @{ -constexpr const double meter = 1; -constexpr const double inch = 2.54 * prefix::centi*meter; -constexpr const double feet = 12 * inch; -/// @} - -/// \name Time -/// @{ -constexpr const double second = 1; -constexpr const double minute = 60 * second; -constexpr const double hour = 60 * minute; -constexpr const double day = 24 * hour; -constexpr const double year = 365 * day; -/// @} - -/// \name Volume -/// @{ -constexpr const double gallon = 231 * cubic(inch); -constexpr const double stb = 42 * gallon; -constexpr const double liter = 1 * cubic(prefix::deci*meter); -/// @} - -/// \name Mass -/// @{ -constexpr const double kilogram = 1; -constexpr const double gram = 1.0e-3 * kilogram; -// http://en.wikipedia.org/wiki/Pound_(mass)#Avoirdupois_pound -constexpr const double pound = 0.45359237 * kilogram; -/// @} - -// -------------------------------------------------------------- -// Standardised constants -// -------------------------------------------------------------- - -/// \name Standardised constant -/// @{ -constexpr const double gravity = 9.80665 * meter/square(second); -/// @} - -// -------------------------------------------------------------- -// Derived units and conversions -// -------------------------------------------------------------- - -/// \name Force -/// @{ -constexpr const double Newton = kilogram*meter / square(second); // == 1 -constexpr const double dyne = 1e-5*Newton; -constexpr const double lbf = pound * gravity; // Pound-force - /// @} - - /// \name Pressure - /// @{ -constexpr const double Pascal = Newton / square(meter); // == 1 -constexpr const double barsa = 100000 * Pascal; -constexpr const double atm = 101325 * Pascal; -constexpr const double psia = lbf / square(inch); -/// @} - -/// \name Temperature. This one is more complicated -/// because the unit systems used by Eclipse (i.e. degrees -/// Celsius and degrees Fahrenheit require to add or -/// subtract an offset for the conversion between from/to -/// Kelvin -/// @{ -constexpr const double degCelsius = 1.0; // scaling factor °C -> K -constexpr const double degCelsiusOffset = 273.15; // offset for the °C -> K conversion - -constexpr const double degFahrenheit = 5.0/9; // scaling factor °F -> K -constexpr const double degFahrenheitOffset = 255.37; // offset for the °C -> K conversion - /// @} - - /// \name Viscosity - /// @{ -constexpr const double Pas = Pascal * second; // == 1 -constexpr const double Poise = prefix::deci*Pas; -/// @} - -namespace perm_details { -constexpr const double p_grad = atm / (prefix::centi*meter); -constexpr const double area = square(prefix::centi*meter); -constexpr const double flux = cubic (prefix::centi*meter) / second; -constexpr const double velocity = flux / area; -constexpr const double visc = prefix::centi*Poise; -constexpr const double darcy = (velocity * visc) / p_grad; -// == 1e-7 [m^2] / 101325 -// == 9.869232667160130e-13 [m^2] -} -/// \name Permeability -/// @{ -/// -/// A porous medium with a permeability of 1 darcy permits a flow (flux) -/// of \f$1\,\mathit{cm}^3/s\f$ of a fluid with viscosity -/// \f$1\,\mathit{cP}\f$ (\f$1\,mPa\cdot s\f$) under a pressure gradient -/// of \f$1\,\mathit{atm}/\mathit{cm}\f$ acting across an area of -/// \f$1\,\mathit{cm}^2\f$. -/// -constexpr const double darcy = perm_details::darcy; -/// @} - -/** -* Unit conversion routines. -*/ -namespace convert { -/** -* Convert from external units of measurements to equivalent -* internal units of measurements. Note: The internal units of -* measurements are *ALWAYS*, and exclusively, SI. -* -* Example: Convert a double @c kx, containing a permeability value -* in units of milli-darcy (mD) to the equivalent value in SI units -* (i.e., \f$m^2\f$). -* \code -* using namespace Opm::unit; -* using namespace Opm::prefix; -* convert::from(kx, milli*darcy); -* \endcode -* -* @param[in] q Physical quantity. -* @param[in] unit Physical unit of measurement. -* @return Value of @c q in equivalent SI units of measurements. -*/ -constexpr double from(const double q, const double unit) -{ - return q * unit; -} - -/** -* Convert from internal units of measurements to equivalent -* external units of measurements. Note: The internal units of -* measurements are *ALWAYS*, and exclusively, SI. -* -* Example: Convert a std::vector p, containing -* pressure values in the SI unit Pascal (i.e., unit::Pascal) to the -* equivalent values in Psi (unit::psia). -* \code -* using namespace Opm::unit; -* std::transform(p.begin(), p.end(), p.begin(), -* boost::bind(convert::to, _1, psia)); -* \endcode -* -* @param[in] q Physical quantity, measured in SI units. -* @param[in] unit Physical unit of measurement. -* @return Value of @c q in unit unit. -*/ -constexpr double to(const double q, const double unit) -{ - return q / unit; -} -} // namespace convert -} - -namespace Metric { -using namespace prefix; -using namespace unit; -constexpr const double Pressure = barsa; -constexpr const double Temperature = degCelsius; -constexpr const double TemperatureOffset = degCelsiusOffset; -constexpr const double AbsoluteTemperature = degCelsius; // actually [K], but the these two are identical -constexpr const double Length = meter; -constexpr const double Time = day; -constexpr const double Mass = kilogram; -constexpr const double Permeability = milli*darcy; -constexpr const double Transmissibility = centi*Poise*cubic(meter)/(day*barsa); -constexpr const double LiquidSurfaceVolume = cubic(meter); -constexpr const double GasSurfaceVolume = cubic(meter); -constexpr const double ReservoirVolume = cubic(meter); -constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume; -constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume; -constexpr const double Density = kilogram/cubic(meter); -constexpr const double PolymerDensity = kilogram/cubic(meter); -constexpr const double Salinity = kilogram/cubic(meter); -constexpr const double Viscosity = centi*Poise; -constexpr const double Timestep = day; -constexpr const double SurfaceTension = dyne/(centi*meter); -} - - -namespace Field { -using namespace prefix; -using namespace unit; -constexpr const double Pressure = psia; -constexpr const double Temperature = degFahrenheit; -constexpr const double TemperatureOffset = degFahrenheitOffset; -constexpr const double AbsoluteTemperature = degFahrenheit; // actually [°R], but the these two are identical -constexpr const double Length = feet; -constexpr const double Time = day; -constexpr const double Mass = pound; -constexpr const double Permeability = milli*darcy; -constexpr const double Transmissibility = centi*Poise*stb/(day*psia); -constexpr const double LiquidSurfaceVolume = stb; -constexpr const double GasSurfaceVolume = 1000*cubic(feet); -constexpr const double ReservoirVolume = stb; -constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume; -constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume; -constexpr const double Density = pound/cubic(feet); -constexpr const double PolymerDensity = pound/stb; -constexpr const double Salinity = pound/stb; -constexpr const double Viscosity = centi*Poise; -constexpr const double Timestep = day; -constexpr const double SurfaceTension = dyne/(centi*meter); -} - - -namespace Lab { -using namespace prefix; -using namespace unit; -constexpr const double Pressure = atm; -constexpr const double Temperature = degCelsius; -constexpr const double TemperatureOffset = degCelsiusOffset; -constexpr const double AbsoluteTemperature = degCelsius; // actually [K], but the these two are identical -constexpr const double Length = centi*meter; -constexpr const double Time = hour; -constexpr const double Mass = gram; -constexpr const double Permeability = milli*darcy; -constexpr const double Transmissibility = centi*Poise*cubic(centi*meter)/(hour*atm); -constexpr const double LiquidSurfaceVolume = cubic(centi*meter); -constexpr const double GasSurfaceVolume = cubic(centi*meter); -constexpr const double ReservoirVolume = cubic(centi*meter); -constexpr const double GasDissolutionFactor = GasSurfaceVolume/LiquidSurfaceVolume; -constexpr const double OilDissolutionFactor = LiquidSurfaceVolume/GasSurfaceVolume; -constexpr const double Density = gram/cubic(centi*meter); -constexpr const double PolymerDensity = gram/cubic(centi*meter); -constexpr const double Salinity = gram/cubic(centi*meter); -constexpr const double Viscosity = centi*Poise; -constexpr const double Timestep = hour; -constexpr const double SurfaceTension = dyne/(centi*meter); -} - -} - -#endif // OPM_UNITS_HEADER From d6862265b7a4e4df2f24183d964dca01eae983f2 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 24 Oct 2018 13:03:59 +0200 Subject: [PATCH 2/3] Remark That 3pt EPS Changes Function Shape in [Sr, Smax] The three-point vertical scaling option will usually change the shape of the relative permeability function in the saturation interval between the scaled critical displacing saturation and the maximum phase saturation. In particular, we usually replace a general, non-linear function with a linear approximation when applying three-point vertical scaling. Add a comment to that effect. --- .../opm/utility/ECLEndPointScaling.cpp | 8 ++++++-- 1 file changed, 6 insertions(+), 2 deletions(-) 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 6b16d2c88b..91522e8c59 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 @@ -516,7 +516,9 @@ vertScale(const FunctionValues& f, // s \in [sr, sm), sm > sr; normal case: Kr(Smax) > Kr(Sr). // // Linear function between (sr,fr) and (sm,fm) in terms of - // function value 'y'. + // function value 'y'. This usually alters the shape of the + // relative permeability function in this interval (e.g., + // roughly quadratic to linear). const auto t = (y - fdisp) / (fmax - fdisp); y = fr + t*(fm - fr); @@ -525,7 +527,9 @@ vertScale(const FunctionValues& f, // s \in [sr, sm), sm > sr; special case: Kr(Smax) == Kr(Sr). // // Use linear function between (sr,fr) and (sm,fm) in terms of - // saturation value 's'. + // saturation value 's'. This usually alters the shape of the + // relative permeability function in this interval (e.g., + // roughly quadratic to linear). const auto t = (s - sr) / (sm - sr); y = fr + t*(fm - fr); From 405f6426f3b3c890705121dbb78edf4d48f4716b Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Wed, 24 Oct 2018 13:07:20 +0200 Subject: [PATCH 3/3] Bring EPS Unit Test Up To Date WRT Changing Behaviour A previous commit changed the reverse mapping behaviour in the case of coincident tabulated scaled displacing saturation and maximum phase saturation in the two-phase systems. Bring the corresponding unit test in line with the changed behaviour (now maps high saturation values to scaled critical displacing saturation Sr rather than Smax in this situation). Pointy Hat: [at]bska Thanks to: [at]atgeirr --- .../tests/test_eclendpointscaling.cpp | 15 +++++++++++---- 1 file changed, 11 insertions(+), 4 deletions(-) diff --git a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp index ef4a92c7a9..3a65935ec4 100644 --- a/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp +++ b/ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/tests/test_eclendpointscaling.cpp @@ -972,6 +972,13 @@ BOOST_AUTO_TEST_CASE (KrOW_3pt) } // Reverse: Lookup/Table so -> Scaled SO. + // + // Note that "tep.disp" is equal to "tep.high" in the input table. In + // this particular case, the reverse scaling will map high saturation + // values to 'Sr' rather than the maximum mobile So (= 0.85) in order to + // enable distinguishing scaled critical saturation from scaled maximum + // saturation for purpose of converting to "regular" phase saturations + // (i.e., SW or SG). { const auto SO = std::vector { 0.20000000000000001, // so = 0 @@ -993,10 +1000,10 @@ BOOST_AUTO_TEST_CASE (KrOW_3pt) 0.6583231470163351, // so = 0.75 0.68888209801089018, // so = 0.80000000000000004 0.71944104900544503, // so = 0.84999999999999998 - 0.84999999999999998, // so = 0.90000000000000002 - 0.84999999999999998, // so = 0.94999999999999996 - 0.84999999999999998, // so = 0.99990000000000001 - 0.84999999999999998, // so = 1 + 0.75, // so = 0.90000000000000002 + 0.75, // so = 0.94999999999999996 + 0.75, // so = 0.99990000000000001 + 0.75, // so = 1 }; const auto so = interp.saturationPoints(Interp::InTable{0});