#2174 Upgrade opm-flowdiagnostics-applications to 75b333335f6cd055d3130d460c6d87444fb7aed4

to get more of the PVT/RelPerm functionality.

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

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());
             }
 

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/examples/extractPropCurves.cpp

@@ -37,24 +37,28 @@
 
 namespace {
     template <class OStream>
-    void printGraph(OStream&                           os,
+    void printGraph(OStream&                                        os,
-                    const std::string&                 name,
+                    const std::string&                              name,
-                    const Opm::FlowDiagnostics::Graph& graph)
+                    const std::vector<Opm::FlowDiagnostics::Graph>& 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 << name << '{' << k << "} = [\n";
-            os << x[i] << ' ' << y[i] << '\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<RC>{};
+        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<RC>{};
+        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';

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 <http://www.gnu.org/licenses/>.
+*/
+
 #include <opm/utility/ECLCaseUtilities.hpp>
 
 #include <exception>

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::TEP>
 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;
 

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.cpp

@@ -265,40 +265,7 @@ Opm::ECLPVT::PVDx::viscosity(const std::vector<double>& p) const
 Opm::FlowDiagnostics::Graph
 Opm::ECLPVT::PVDx::getPvtCurve(const RawCurve curve) const
 {
-    assert ((curve == RawCurve::FVF) ||
+    return extractRawPVTCurve(this->interp_, curve);
-            (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<double>, vector<double>>
-    return FlowDiagnostics::Graph { std::move(x), std::move(y) };
 }
 
 // =====================================================================

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCommon.hpp

@@ -28,6 +28,8 @@
 #include <opm/utility/ECLUnitHandling.hpp>
 
 #include <algorithm>
+#include <cassert>
+#include <cstddef>
 #include <array>
 #include <functional>
 #include <initializer_list>
@@ -294,6 +296,9 @@ namespace Opm { namespace ECLPVT {
 
         /// Viscosity
         Viscosity,
+
+        /// Enveloping curve for saturated state (wet gas or live oil)
+        SaturatedState,
     };
 
     template <std::size_t N>
@@ -390,6 +395,47 @@ namespace Opm { namespace ECLPVT {
         return u -= v;
     }
 
+    template <class Extrapolation, bool IsAscendingRange>
+    FlowDiagnostics::Graph
+    extractRawPVTCurve(const Interp1D::PiecewisePolynomial::Linear<Extrapolation, IsAscendingRange>& 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<double>, vector<double>>
+        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<FlowDiagnostics::Graph>
+        getPvtCurve(const RawCurve curve) const
+        {
+            auto ret = std::vector<FlowDiagnostics::Graph>{};
+            ret.reserve(this->propInterp_.size());
+
+            for (const auto& interp : this->propInterp_) {
+                ret.push_back(extractRawPVTCurve(interp, curve));
+            }
+
+            return ret;
+        }
+
+        std::vector<double> getSaturatedPoints() const
+        {
+            auto y = std::vector<double>{};
+            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<SubtableInterpolant>().classifyPoint(0.0))

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 <http://www.gnu.org/licenses/>.
+*/
+
 #include <opm/utility/ECLPvtCurveCollection.hpp>
 
 #include <opm/utility/ECLResultData.hpp>
 
+#include <initializer_list>
 #include <vector>
 
+#include <ert/ecl/ecl_kw_magic.h>
+
 namespace {
     std::vector<int>
     pvtnumVector(const ::Opm::ECLGraph&        G,
@@ -19,6 +41,35 @@ namespace {
 
         return pvtnum;
     }
+
+    std::unique_ptr<const Opm::ECLUnits::UnitSystem>
+    createUnitSystem(const ::Opm::ECLInitFileData& init)
+    {
+        const auto& ih = init.keywordData<int>(INTEHEAD_KW);
+
+        return ::Opm::ECLUnits::createUnitSystem(ih[ INTEHEAD_UNIT_INDEX ]);
+    }
+
+    std::vector<Opm::FlowDiagnostics::Graph> emptyFDGraph()
+    {
+        return { Opm::FlowDiagnostics::Graph{} };
+    }
+
+    template <class PvtPtr>
+    std::vector<Opm::FlowDiagnostics::Graph>
+    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::FlowDiagnostics::Graph>
 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_cast<decltype(this->pvtnum_.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);
+        // Caller requests oil properties.
-        return FlowDiagnostics::Graph{};
+        return rawPvtCurve(this->oil_, curve, regID, *this->usys_);
     }
 
-    if (this->gas_) {
+    // Call requests gas properties.
-        return this->gas_->getPvtCurve(curve, regID);
+    return rawPvtCurve(this->gas_, curve, regID, *this->usys_);
-    }
-    else {
-        // Result set does not provide tabulated gas properties.  Return
-        // empty.
-        return FlowDiagnostics::Graph{};
-    }
 }

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtCurveCollection.hpp

@@ -25,6 +25,7 @@
 #include <opm/utility/ECLPvtCommon.hpp>
 #include <opm/utility/ECLPvtGas.hpp>
 #include <opm/utility/ECLPvtOil.hpp>
+#include <opm/utility/ECLUnitHandling.hpp>
 
 #include <memory>
 #include <vector>
@@ -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<FlowDiagnostics::Graph>
         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> oil_;
+
+        /// Unit handling (SI -> result-set convention)
+        std::shared_ptr<const ECLUnits::UnitSystem> usys_;
     };
 
 }} // Opm::ECLPVT

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.cpp

@@ -25,10 +25,13 @@
 #include <opm/utility/ECLResultData.hpp>
 #include <opm/utility/ECLUnitHandling.hpp>
 
+#include <opm/parser/eclipse/Units/Units.hpp>
+
 #include <cassert>
 #include <cmath>
 #include <exception>
 #include <functional>
+#include <initializer_list>
 #include <memory>
 #include <stdexcept>
 #include <string>
@@ -95,8 +98,9 @@ public:
     viscosity(const std::vector<double>& rv,
               const std::vector<double>& pg) const = 0;
 
-    virtual Opm::FlowDiagnostics::Graph
+    virtual std::vector<Opm::FlowDiagnostics::Graph>
-    getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0;
+    getPvtCurve(const Opm::ECLPVT::RawCurve      curve,
+                const Opm::ECLUnits::UnitSystem& usys) const = 0;
 
     virtual std::unique_ptr<PVxGBase> clone() const = 0;
 };
@@ -130,10 +134,31 @@ public:
         return this->interpolant_.viscosity(pg);
     }
 
-    virtual Opm::FlowDiagnostics::Graph
+    virtual std::vector<Opm::FlowDiagnostics::Graph>
-    getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override
+    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<PVxGBase> clone() const override
@@ -158,7 +183,8 @@ public:
 
     WetGas(std::vector<double>              key,
            std::vector<SubtableInterpolant> propInterp)
-        : interp_(std::move(key), std::move(propInterp))
+        : key_   (key)          // Copy
+        , interp_(std::move(key), std::move(propInterp))
     {}
 
     virtual std::vector<double>
@@ -187,13 +213,15 @@ public:
         return this->interp_.viscosity(key, x);
     }
 
-    virtual Opm::FlowDiagnostics::Graph
+    virtual std::vector<Opm::FlowDiagnostics::Graph>
-    getPvtCurve(const Opm::ECLPVT::RawCurve /* curve */) const override
+    getPvtCurve(const Opm::ECLPVT::RawCurve      curve,
+                const Opm::ECLUnits::UnitSystem& usys) const override
     {
-        throw std::runtime_error {
+        if (curve == ::Opm::ECLPVT::RawCurve::SaturatedState) {
-            "Property Evaluator for Wet Gas Does not "
+            return this->saturatedStateCurve(usys);
-            "Support Retrieving Raw Curves (Blame BSKA)"
+        }
-        };
+
+        return this->mainPvtCurve(curve, usys);
     }
 
     virtual std::unique_ptr<PVxGBase> clone() const override
@@ -204,9 +232,68 @@ public:
 private:
     using TableInterpolant = ::Opm::ECLPVT::PVTx<SubtableInterpolant>;
 
-    TableInterpolant interp_;
+    std::vector<double> key_;
+    TableInterpolant    interp_;
+
+    std::vector<Opm::FlowDiagnostics::Graph>
+    mainPvtCurve(const Opm::ECLPVT::RawCurve      curve,
+                 const Opm::ECLUnits::UnitSystem& usys) const;
+
+    std::vector<Opm::FlowDiagnostics::Graph>
+    saturatedStateCurve(const Opm::ECLUnits::UnitSystem& usys) const;
 };
 
+std::vector<Opm::FlowDiagnostics::Graph>
+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<Opm::FlowDiagnostics::Graph>
+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
+    std::vector<FlowDiagnostics::Graph>
-    getPvtCurve(const RegIdx   region,
+    getPvtCurve(const RegIdx                region,
-                const RawCurve curve) const;
+                const RawCurve              curve,
+                const ECLUnits::UnitSystem& usys) const;
 
 private:
     std::vector<EvalPtr> eval_;
@@ -460,14 +548,15 @@ viscosity(const RegIdx               region,
     return this->eval_[region]->viscosity(rv, pg);
 }
 
-Opm::FlowDiagnostics::Graph
+std::vector<Opm::FlowDiagnostics::Graph>
 Opm::ECLPVT::Gas::Impl::
-getPvtCurve(const RegIdx   region,
+getPvtCurve(const RegIdx                region,
-            const RawCurve curve) const
+            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::FlowDiagnostics::Graph>
 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);
 }
 
 // =====================================================================

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtGas.hpp

@@ -22,6 +22,7 @@
 
 #include <opm/flowdiagnostics/DerivedQuantities.hpp>
 #include <opm/utility/ECLPvtCommon.hpp>
+#include <opm/utility/ECLUnitHandling.hpp>
 
 #include <memory>
 #include <vector>
@@ -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
+        /// \param[in] usys Unit system.  Collection of units to which the
-        ///    requests represented by \p func and \p region.
+        ///    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
+        std::vector<FlowDiagnostics::Graph>
-        getPvtCurve(const RawCurve curve, const int region) const;
+        getPvtCurve(const RawCurve              curve,
+                    const int                   region,
+                    const ECLUnits::UnitSystem& usys) const;
 
     private:
         /// Implementation class.

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.cpp

@@ -24,10 +24,13 @@
 #include <opm/utility/ECLResultData.hpp>
 #include <opm/utility/ECLUnitHandling.hpp>
 
+#include <opm/parser/eclipse/Units/Units.hpp>
+
 #include <cassert>
 #include <cmath>
 #include <exception>
 #include <functional>
+#include <initializer_list>
 #include <memory>
 #include <stdexcept>
 #include <string>
@@ -92,8 +95,9 @@ public:
     viscosity(const std::vector<double>& rs,
               const std::vector<double>& po) const = 0;
 
-    virtual Opm::FlowDiagnostics::Graph
+    virtual std::vector<Opm::FlowDiagnostics::Graph>
-    getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0;
+    getPvtCurve(const Opm::ECLPVT::RawCurve      curve,
+                const Opm::ECLUnits::UnitSystem& usys) const = 0;
 
     virtual std::unique_ptr<PVxOBase> clone() const = 0;
 };
@@ -127,10 +131,31 @@ public:
         return this->interpolant_.viscosity(po);
     }
 
-    virtual Opm::FlowDiagnostics::Graph
+    virtual std::vector<Opm::FlowDiagnostics::Graph>
-    getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override
+    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<PVxOBase> clone() const override
@@ -155,7 +180,8 @@ public:
 
     LiveOil(std::vector<double>              key,
             std::vector<SubtableInterpolant> propInterp)
-        : interp_(std::move(key), std::move(propInterp))
+        : key_   (key)          // Copy
+        , interp_(std::move(key), std::move(propInterp))
     {}
 
     virtual std::vector<double>
@@ -184,13 +210,15 @@ public:
         return this->interp_.viscosity(key, x);
     }
 
-    virtual Opm::FlowDiagnostics::Graph
+    virtual std::vector<Opm::FlowDiagnostics::Graph>
-    getPvtCurve(const Opm::ECLPVT::RawCurve /* curve */) const override
+    getPvtCurve(const Opm::ECLPVT::RawCurve      curve,
+                const Opm::ECLUnits::UnitSystem& usys) const override
     {
-        throw std::runtime_error {
+        if (curve == ::Opm::ECLPVT::RawCurve::SaturatedState) {
-            "Property Evaluator for Live Oil Does not "
+            return this->saturatedStateCurve(usys);
-            "Support Retrieving Raw Curves (Blame BSKA)"
+        }
-        };
+
+        return this->mainPvtCurve(curve, usys);
     }
 
     virtual std::unique_ptr<PVxOBase> clone() const override
@@ -201,9 +229,68 @@ public:
 private:
     using TableInterpolant = ::Opm::ECLPVT::PVTx<SubtableInterpolant>;
 
-    TableInterpolant interp_;
+    std::vector<double> key_;
+    TableInterpolant    interp_;
+
+    std::vector<Opm::FlowDiagnostics::Graph>
+    mainPvtCurve(const Opm::ECLPVT::RawCurve      curve,
+                 const Opm::ECLUnits::UnitSystem& usys) const;
+
+    std::vector<Opm::FlowDiagnostics::Graph>
+    saturatedStateCurve(const Opm::ECLUnits::UnitSystem& usys) const;
 };
 
+std::vector<Opm::FlowDiagnostics::Graph>
+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<Opm::FlowDiagnostics::Graph>
+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
+    std::vector<FlowDiagnostics::Graph>
-    getPvtCurve(const RegIdx   region,
+    getPvtCurve(const RegIdx                region,
-                const RawCurve curve) const;
+                const RawCurve              curve,
+                const ECLUnits::UnitSystem& usys) const;
 
 private:
     std::vector<EvalPtr> eval_;
@@ -457,14 +545,15 @@ viscosity(const RegIdx               region,
     return this->eval_[region]->viscosity(rs, po);
 }
 
-Opm::FlowDiagnostics::Graph
+std::vector<Opm::FlowDiagnostics::Graph>
 Opm::ECLPVT::Oil::Impl::
-getPvtCurve(const RegIdx   region,
+getPvtCurve(const RegIdx                region,
-            const RawCurve curve) const
+            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::FlowDiagnostics::Graph>
 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);
 }
 
 // =====================================================================

ThirdParty/custom-opm-flowdiag-app/opm-flowdiagnostics-applications/opm/utility/ECLPvtOil.hpp

@@ -21,6 +21,7 @@
 #define OPM_ECLPVTOIL_HEADER_INCLUDED
 
 #include <opm/utility/ECLPvtCommon.hpp>
+#include <opm/utility/ECLUnitHandling.hpp>
 
 #include <memory>
 #include <vector>
@@ -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
+        /// \param[in] usys Unit system.  Collection of units to which the
-        ///    requests represented by \p func and \p region.
+        ///    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
+        std::vector<FlowDiagnostics::Graph>
-        getPvtCurve(const RawCurve curve, const int region) const;
+        getPvtCurve(const RawCurve              curve,
+                    const int                   region,
+                    const ECLUnits::UnitSystem& usys) const;
 
     private:
         /// Implementation class.

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<double>{};  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<double>{
+            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<int>&    tabdims,
+            SatFunction(const std::vector<int>&    tabdims,
-                       const std::vector<double>& tab,
+                        const std::vector<double>& tab,
-                       const int                  usys)
+                        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<double>
+            pcgo(const std::size_t          regID,
+                 const std::vector<double>& sg) const
+            {
+                const auto t = this->table(regID);
+                const auto c = this->pccol();
+
+                return this->func_.interpolate(t, c, sg);
+            }
+
             const std::vector<double>&
             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<int>&    tabdims,
+            SatFunction(const std::vector<int>&    tabdims,
-                       const std::vector<double>& tab,
+                        const std::vector<double>& tab,
-                       const int                  usys)
+                        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<double>
+            pcow(const std::size_t          regID,
+                 const std::vector<double>& sw) const
+            {
+                const auto t = this->table(regID);
+                const auto c = this->pccol();
+
+                return this->func_.interpolate(t, c, sw);
+            }
+
             const std::vector<double>&
             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<int>&    tabdims,
+Gas::Details::tableData(const std::vector<int>&    tabdims,
-                                 const std::vector<double>& tab)
+                        const std::vector<double>& tab)
 {
     auto t = Opm::ECLPropTableRawData{};
 
@@ -505,7 +558,7 @@ Relperm::Gas::Details::tableData(const std::vector<int>&    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<int>&    tabdims,
+Oil::Details::tableData(const std::vector<int>&    tabdims,
-                                 const bool                 isTwoP,
+                        const bool                 isTwoP,
-                                 const std::vector<double>& tab)
+                        const std::vector<double>& tab)
 {
     auto t = Opm::ECLPropTableRawData{};
 
@@ -556,7 +609,7 @@ Relperm::Oil::Details::tableData(const std::vector<int>&    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<int>&    tabdims,
+Water::Details::tableData(const std::vector<int>&    tabdims,
-                                   const std::vector<double>& tab)
+                          const std::vector<double>& tab)
 {
     auto t = Opm::ECLPropTableRawData{};
 
@@ -596,7 +649,7 @@ Relperm::Water::Details::tableData(const std::vector<int>&    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,
+        void scaleKrOG(const ECLRegionMapping& rmap,
-                     std::vector<double>&    so) const
+                       std::vector<double>&    so) const
         {
             this->scale(this->oil_in_og_, rmap, so);
         }
 
-        void scaleOW(const ECLRegionMapping& rmap,
+        void scaleKrOW(const ECLRegionMapping& rmap,
-                     std::vector<double>&    so) const
+                       std::vector<double>&    so) const
         {
             this->scale(this->oil_in_ow_, rmap, so);
         }
 
-        void scaleGas(const ECLRegionMapping& rmap,
+        void scaleKrGas(const ECLRegionMapping& rmap,
-                      std::vector<double>&    sg) const
+                        std::vector<double>&    sg) const
         {
-            this->scale(this->gas_, rmap, sg);
+            this->scale(this->gas_.kr, rmap, sg);
         }
 
-        void scaleWat(const ECLRegionMapping& rmap,
+        void scaleKrWat(const ECLRegionMapping& rmap,
-                      std::vector<double>&    sw) const
+                        std::vector<double>&    sw) const
         {
-            this->scale(this->wat_, rmap, sw);
+            this->scale(this->wat_.kr, rmap, sw);
         }
 
-        void reverseScaleOG(const ECLRegionMapping& rmap,
+        void scalePcGO(const ECLRegionMapping& rmap,
-                            std::vector<double>&    so) const
+                       std::vector<double>&    sg) const
+        {
+            this->scale(this->gas_.pc, rmap, sg);
+        }
+
+        void scalePcOW(const ECLRegionMapping& rmap,
+                       std::vector<double>&    sw) const
+        {
+            this->scale(this->wat_.pc, rmap, sw);
+        }
+
+        void reverseScaleKrOG(const ECLRegionMapping& rmap,
+                              std::vector<double>&    so) const
         {
             this->reverseScale(this->oil_in_og_, rmap, so);
         }
 
-        void reverseScaleOW(const ECLRegionMapping& rmap,
+        void reverseScaleKrOW(const ECLRegionMapping& rmap,
-                            std::vector<double>&    so) const
+                              std::vector<double>&    so) const
         {
             this->reverseScale(this->oil_in_ow_, rmap, so);
         }
 
-        void reverseScaleGas(const ECLRegionMapping& rmap,
+        void reverseScaleKrGas(const ECLRegionMapping& rmap,
-                             std::vector<double>&    sg) const
+                               std::vector<double>&    sg) const
         {
-            this->reverseScale(this->gas_, rmap, sg);
+            this->reverseScale(this->gas_.kr, rmap, sg);
         }
 
-        void reverseScaleWat(const ECLRegionMapping& rmap,
+        void reverseScaleKrWat(const ECLRegionMapping& rmap,
-                             std::vector<double>&    sw) const
+                               std::vector<double>&    sw) const
         {
-            this->reverseScale(this->wat_, rmap, sw);
+            this->reverseScale(this->wat_.kr, rmap, sw);
+        }
+
+        void reverseScalePcGO(const ECLRegionMapping& rmap,
+                              std::vector<double>&    sg) const
+        {
+            this->reverseScale(this->gas_.pc, rmap, sg);
+        }
+
+        void reverseScalePcOW(const ECLRegionMapping& rmap,
+                              std::vector<double>&    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_;
+        struct FullEPS {
-        EPS oil_in_ow_;
+            EPS kr;
-        EPS gas_;
+            EPS pc;
-        EPS wat_;
+        };
+
+        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 =
+            // Scaling for Gas relative permeabilty
-                Create::fromECLOutput(G, init, opt);
+            {
+                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 =
+            // Scaling for Water relative permeabilty
-                Create::fromECLOutput(G, init, opt);
+            {
+                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<Relperm::Oil::KrFunction>   oil_;
+    std::unique_ptr<Oil::KrFunction>    oil_;
-    std::unique_ptr<Relperm::Gas::KrFunction>   gas_;
+    std::unique_ptr<Gas::SatFunction>   gas_;
-    std::unique_ptr<Relperm::Water::KrFunction> wat_;
+    std::unique_ptr<Water::SatFunction> wat_;
 
     std::unique_ptr<EPSEvaluator> eps_;
 
@@ -1002,6 +1131,12 @@ private:
              const std::vector<double>& sg,
              const bool                 useEPS) const;
 
+    FlowDiagnostics::Graph
+    pcgoCurve(const ECLRegionMapping&    rmap,
+              const std::size_t          regID,
+              const std::vector<double>& sg,
+              const bool                 useEPS) const;
+
     std::vector<double>
     krw(const ECLGraph&       G,
         const ECLRestartData& rstrt,
@@ -1013,18 +1148,32 @@ private:
              const std::vector<double>& sw,
              const bool                 useEPS) const;
 
-    void scaleGasSat(const ECLRegionMapping& rmap,
+    FlowDiagnostics::Graph
-                     const bool              useEPS,
+    pcowCurve(const ECLRegionMapping&    rmap,
-                     std::vector<double>&    sg) const;
+              const std::size_t          regID,
+              const std::vector<double>& sw,
+              const bool                 useEPS) const;
 
-    void scaleOilSat(const ECLRegionMapping& rmap,
+    void scaleKrGasSat(const ECLRegionMapping& rmap,
-                     const bool              useEPS,
-                     std::vector<double>&    so_g,
-                     std::vector<double>&    so_w) const;
-
-    void scaleWaterSat(const ECLRegionMapping& rmap,
                        const bool              useEPS,
-                       std::vector<double>&    sw) const;
+                       std::vector<double>&    sg) const;
+
+    void scalePcGasSat(const ECLRegionMapping& rmap,
+                       const bool              useEPS,
+                       std::vector<double>&    sg) const;
+
+    void scaleKrOilSat(const ECLRegionMapping& rmap,
+                       const bool              useEPS,
+                       std::vector<double>&    so_g,
+                       std::vector<double>&    so_w) const;
+
+    void scaleKrWaterSat(const ECLRegionMapping& rmap,
+                         const bool              useEPS,
+                         std::vector<double>&    sw) const;
+
+    void scalePcWaterSat(const ECLRegionMapping& rmap,
+                         const bool              useEPS,
+                         std::vector<double>&    sw) const;
 
     void uniqueReverseScaleSat(std::vector<double>& s) const;
 
@@ -1133,16 +1282,16 @@ Impl::initRelPermInterp(const EPSEvaluator::ActPh& active,
     const auto& tab     = init.keywordData<double>("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<RawCurve>& 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
+                // extract a saturation function curve for water in the O/G
-                // O/G subsystem.  No such curve.
+                // subsystem.  No such curve.
                 graph.emplace_back();
                 continue;
             }
@@ -1272,10 +1418,13 @@ getSatFuncCurve(const std::vector<RawCurve>& func,
                     .reset(new ECLRegionMapping(this->satnum_, subset));
             }
 
-            auto krw = this->krwCurve(*wat_assoc.second, regID,
+            auto crv = (fi.curve == RawCurve::Function::RelPerm)
-                                      wat_assoc.first, useEPS);
+                ? 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<RawCurve>& func,
                     .reset(new ECLRegionMapping(this->satnum_, subset));
             }
 
-            auto kro = this->kroCurve(*oil_assoc.second, regID,
+            const auto kro =
-                                      fi.subsys, oil_assoc.first, useEPS);
+                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<RawCurve>& 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
+                // extract a saturation function curve for gas in the O/W
-                // O/W subsystem.  No such curve.
+                // subsystem.  No such curve.
                 graph.emplace_back();
                 continue;
             }
@@ -1331,10 +1485,13 @@ getSatFuncCurve(const std::vector<RawCurve>& func,
                     .reset(new ECLRegionMapping(this->satnum_, subset));
             }
 
-            auto krg = this->krgCurve(*gas_assoc.second, regID,
+            auto crv = (fi.curve == RawCurve::Function::RelPerm)
-                                      gas_assoc.first, useEPS);
+                ? 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_->reverseScaleKrOG(rmap, so_g);
-        this->eps_->reverseScaleOW(rmap, so_w);
+        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<double>(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<double>& 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<double>, vector<double>>
+    return FlowDiagnostics::Graph {
+        std::move(abscissas),
+        { std::begin(pc), std::begin(pc) + nPoints }
+    };
+}
+
 std::vector<double>
 Opm::ECLSaturationFunc::Impl::
 krw(const ECLGraph&       G,
@@ -1548,7 +1741,7 @@ krw(const ECLGraph&       G,
     auto sw = G.rawLinearisedCellData<double>(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<double>& 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<double>, vector<double>>
+    return FlowDiagnostics::Graph {
+        std::move(abscissas),
+        { std::begin(pc), std::begin(pc) + nPoints }
+    };
+}
+
 void
 Opm::ECLSaturationFunc::Impl::
-scaleGasSat(const ECLRegionMapping& rmap,
+scaleKrGasSat(const ECLRegionMapping& rmap,
-            const bool              useEPS,
+              const bool              useEPS,
-            std::vector<double>&    sg) const
+              std::vector<double>&    sg) const
 {
     if (useEPS && this->eps_) {
-        this->eps_->scaleGas(rmap, sg);
+        this->eps_->scaleKrGas(rmap, sg);
     }
 }
 
 void
 Opm::ECLSaturationFunc::Impl::
-scaleOilSat(const ECLRegionMapping& rmap,
+scalePcGasSat(const ECLRegionMapping& rmap,
-            const bool              useEPS,
+              const bool              useEPS,
-            std::vector<double>&    so_g,
+              std::vector<double>&    sg) const
-            std::vector<double>&    so_w) const
+{
+    if (useEPS && this->eps_) {
+        this->eps_->scalePcGO(rmap, sg);
+    }
+}
+
+void
+Opm::ECLSaturationFunc::Impl::
+scaleKrOilSat(const ECLRegionMapping& rmap,
+              const bool              useEPS,
+              std::vector<double>&    so_g,
+              std::vector<double>&    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_->scaleKrOG(rmap, so_g);
-        this->eps_->scaleOW(rmap, so_w);
+        this->eps_->scaleKrOW(rmap, so_w);
     }
 }
 
 void
 Opm::ECLSaturationFunc::Impl::
-scaleWaterSat(const ECLRegionMapping& rmap,
+scaleKrWaterSat(const ECLRegionMapping& rmap,
-              const bool              useEPS,
+                const bool              useEPS,
-              std::vector<double>&    sg) const
+                std::vector<double>&    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<double>&    sg) const
+{
+    if (useEPS && this->eps_) {
+        this->eps_->scalePcOW(rmap, sg);
     }
 }
 

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.