Upgraded opm-flowdiagnostics-applications to a773bcfc963705679ecc28f58048fc38936fbbc6

ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp

@@ -707,12 +707,12 @@ std::vector<RigFlowDiagSolverInterface::PvtCurve> RigFlowDiagSolverInterface::ca
         const Opm::ECLPhaseIndex queryPhaseIndex = queryPhaseArr[i];
         const PvtCurve::Phase mapToPhase = mapToPhaseArr[i];
 
-        std::vector<Opm::FlowDiagnostics::Graph> graphArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getPvtCurve(rawCurveType, queryPhaseIndex, static_cast<int>(activeCellIndex));
+        std::vector<Opm::ECLPVT::PVTGraph> graphArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getPvtCurve(rawCurveType, queryPhaseIndex, static_cast<int>(activeCellIndex));
-        for (Opm::FlowDiagnostics::Graph srcGraph : graphArr)
+        for (Opm::ECLPVT::PVTGraph srcGraph : graphArr)
         {
-            if (srcGraph.first.size() > 0)
+            if (srcGraph.press.size() > 0)
             {
-                retCurveArr.push_back({ mapToPhase, srcGraph.first, srcGraph.second });
+                retCurveArr.push_back({ mapToPhase, srcGraph.press, srcGraph.value});
             }
         }
     }

ResInsightVersion.cmake

@@ -14,7 +14,7 @@ set(ERT_GITHUB_SHA "2e36798b43daf18c112b91aa3febbf2fccd4a95f")
 set(OPM_FLOWDIAGNOSTICS_SHA "7e2be931d430796ed42efcfb5c6b67a8d5962f7f")
 
 # https://github.com/OPM/opm-flowdiagnostics-applications
-set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "cf8a5d2381776a8887d61f825b40f03b1d1cc4b0")
+set(OPM_FLOWDIAGNOSTICS_APPLICATIONS_SHA "a773bcfc963705679ecc28f58048fc38936fbbc6")
 
 # https://github.com/OPM/opm-parser/blob/master/opm/parser/eclipse/Units/Units.hpp
 # This file was moved from opm-core to opm-parser october 2016

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

@@ -22,6 +22,7 @@
 
 #include <opm/utility/ECLCaseUtilities.hpp>
 #include <opm/utility/ECLPhaseIndex.hpp>
+#include <opm/utility/ECLPropertyUnitConversion.hpp>
 #include <opm/utility/ECLPvtCommon.hpp>
 #include <opm/utility/ECLPvtCurveCollection.hpp>
 #include <opm/utility/ECLResultData.hpp>
@@ -63,6 +64,34 @@ namespace {
         os.precision(oprec);
     }
 
+    template <class OStream>
+    void printGraph(OStream&                                  os,
+                    const std::string&                        name,
+                    const std::vector<Opm::ECLPVT::PVTGraph>& graphs)
+    {
+        const auto oprec  = os.precision(16);
+        const auto oflags = os.setf(std::ios_base::scientific);
+
+        auto k = 1;
+        for (const auto& graph : graphs) {
+            const auto& p = graph.press;
+            const auto& R = graph.mixRat;
+            const auto& f = graph.value;
+
+            os << name << '{' << k << "} = [\n";
+
+            for (auto n = p.size(), i = 0*n; i < n; ++i) {
+                os << p[i] << ' ' << R[i] << ' ' << f[i] << '\n';
+            }
+
+            os << "];\n\n";
+            k += 1;
+        }
+
+        os.setf(oflags);
+        os.precision(oprec);
+    }
+
     // -----------------------------------------------------------------
     // Relative permeability
 
@@ -274,6 +303,38 @@ namespace {
 
         printGraph(std::cout, "rsSat", graph);
     }
+
+    // -----------------------------------------------------------------
+
+    std::unique_ptr<const Opm::ECLUnits::UnitSystem>
+    makeUnits(const std::string&          unit,
+              const Opm::ECLInitFileData& init)
+    {
+        if ((unit == "si") || (unit == "SI") || (unit == "internal")) {
+            return {};          // No conversion needed.
+        }
+
+        if ((unit == "metric") || (unit == "Metric") || (unit == "METRIC")) {
+            return Opm::ECLUnits::metricUnitConventions();
+        }
+
+        if ((unit == "field") || (unit == "Field") || (unit == "FIELD")) {
+            return Opm::ECLUnits::fieldUnitConventions();
+        }
+
+        if ((unit == "lab") || (unit == "Lab") || (unit == "LAB")) {
+            return Opm::ECLUnits::labUnitConventions();
+        }
+
+        if ((unit == "pvt-m") || (unit == "PVT-M") || (unit == "PVTM")) {
+            return Opm::ECLUnits::pvtmUnitConventions();
+        }
+
+        std::cerr << "Unit convention '" << unit << "' not recognized\n"
+                  << "Using 'native' (input/serialised) conventions.\n";
+
+        return Opm::ECLUnits::serialisedUnitConventions(init);
+    }
 } // namespace Anonymous
 
 int main(int argc, char* argv[])
@@ -285,8 +346,16 @@ try {
     const auto rset  = example::identifyResultSet(prm);
     const auto init  = Opm::ECLInitFileData(rset.initFile());
     const auto graph = Opm::ECLGraph::load(rset.gridFile(), init);
-    const auto sfunc = Opm::ECLSaturationFunc(graph, init, useEPS);
+
-    const auto pvtCC = Opm::ECLPVT::ECLPvtCurveCollection(graph, init);
+    auto sfunc = Opm::ECLSaturationFunc(graph, init, useEPS);
+    auto pvtCC = Opm::ECLPVT::ECLPvtCurveCollection(graph, init);
+
+    if (prm.has("unit")) {
+        auto units = makeUnits(prm.get<std::string>("unit"), init);
+        
+        sfunc.setOutputUnits(units->clone());
+        pvtCC.setOutputUnits(std::move(units));
+    }
 
     // -----------------------------------------------------------------
     // Relative permeability

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

@@ -29,6 +29,9 @@
 
 #include <cassert>
 #include <exception>
+#include <initializer_list>
+#include <iterator>
+#include <numeric>
 #include <stdexcept>
 #include <utility>
 #include <vector>
@@ -76,6 +79,21 @@ namespace {
 
         return tep;
     }
+
+    bool validSaturation(const double s)
+    {
+        return (! (s < 0.0)) && (! (s > 1.0));
+    }
+
+    bool validSaturations(std::initializer_list<double> sats)
+    {
+        return std::accumulate(std::begin(sats),
+                               std::end  (sats), true,
+            [](const bool result, const double s) -> bool
+        {
+            return result && validSaturation(s);
+        });
+    }
 }
 
 // ---------------------------------------------------------------------
@@ -85,10 +103,12 @@ namespace {
 class Opm::SatFunc::TwoPointScaling::Impl
 {
 public:
-    Impl(std::vector<double> smin,
+    Impl(std::vector<double>      smin,
-         std::vector<double> smax)
+         std::vector<double>      smax,
-        : smin_(std::move(smin))
+         InvalidEndpointBehaviour handle_invalid)
-        , smax_(std::move(smax))
+        : smin_          (std::move(smin))
+        , smax_          (std::move(smax))
+        , handle_invalid_(handle_invalid)
     {
         if (this->smin_.size() != this->smax_.size()) {
             throw std::invalid_argument {
@@ -109,6 +129,11 @@ public:
 private:
     std::vector<double> smin_;
     std::vector<double> smax_;
+
+    InvalidEndpointBehaviour handle_invalid_;
+
+    void handleInvalidEndpoint(const SaturationAssoc& sp,
+                               std::vector<double>&   effsat) const;
 };
 
 std::vector<double>
@@ -127,6 +152,12 @@ Impl::eval(const TableEndPoints&   tep,
         const auto sLO = this->smin_[cell];
         const auto sHI = this->smax_[cell];
 
+        if (! validSaturations({ sLO, sHI })) {
+            this->handleInvalidEndpoint(eval_pt, effsat);
+
+            continue;
+        }
+
         effsat.push_back(0.0);
         auto& s_eff = effsat.back();
 
@@ -192,6 +223,22 @@ Impl::reverse(const TableEndPoints&   tep,
     return unscaledsat;
 }
 
+void
+Opm::SatFunc::TwoPointScaling::Impl::
+handleInvalidEndpoint(const SaturationAssoc& sp,
+                      std::vector<double>&   effsat) const
+{
+    if (this->handle_invalid_ == InvalidEndpointBehaviour::UseUnscaled) {
+        // User requests that invalid scaling be treated as unscaled
+        // saturations.  Pick that.
+        effsat.push_back(sp.sat);
+        return;
+    }
+
+    // Nothing to do for IgnorePoint.  In particular, we must not change the
+    // contents or the size of effsat.
+}
+
 // ---------------------------------------------------------------------
 // Class Opm::ThreePointScaling::Impl
 // ---------------------------------------------------------------------
@@ -199,12 +246,14 @@ Impl::reverse(const TableEndPoints&   tep,
 class Opm::SatFunc::ThreePointScaling::Impl
 {
 public:
-    Impl(std::vector<double> smin ,
+    Impl(std::vector<double>      smin ,
-         std::vector<double> sdisp,
+         std::vector<double>      sdisp,
-         std::vector<double> smax )
+         std::vector<double>      smax ,
-        : smin_ (std::move(smin ))
+         InvalidEndpointBehaviour handle_invalid)
-        , sdisp_(std::move(sdisp))
+        : smin_          (std::move(smin ))
-        , smax_ (std::move(smax ))
+        , sdisp_         (std::move(sdisp))
+        , smax_          (std::move(smax ))
+        , handle_invalid_(handle_invalid)
     {
         if ((this->sdisp_.size() != this->smin_.size()) ||
             (this->sdisp_.size() != this->smax_.size()))
@@ -228,6 +277,11 @@ private:
     std::vector<double> smin_;
     std::vector<double> sdisp_;
     std::vector<double> smax_;
+
+    InvalidEndpointBehaviour handle_invalid_;
+
+    void handleInvalidEndpoint(const SaturationAssoc& sp,
+                               std::vector<double>&   effsat) const;
 };
 
 std::vector<double>
@@ -241,13 +295,19 @@ Impl::eval(const TableEndPoints&   tep,
     for (const auto& eval_pt : sp) {
         const auto cell = eval_pt.cell;
 
-        effsat.push_back(0.0);
-        auto& s_eff = effsat.back();
-
         const auto sLO = this->smin_ [cell];
         const auto sR  = this->sdisp_[cell];
         const auto sHI = this->smax_ [cell];
 
+        if (! validSaturations({ sLO, sR, sHI })) {
+            this->handleInvalidEndpoint(eval_pt, effsat);
+
+            continue;
+        }
+
+        effsat.push_back(0.0);
+        auto& s_eff = effsat.back();
+
         if (! (eval_pt.sat > sLO)) {
             // s <= sLO
             s_eff = tep.low;
@@ -324,6 +384,22 @@ Impl::reverse(const TableEndPoints&   tep,
     return unscaledsat;
 }
 
+void
+Opm::SatFunc::ThreePointScaling::Impl::
+handleInvalidEndpoint(const SaturationAssoc& sp,
+                      std::vector<double>&   effsat) const
+{
+    if (this->handle_invalid_ == InvalidEndpointBehaviour::UseUnscaled) {
+        // User requests that invalid scaling be treated as unscaled
+        // saturations.  Pick that.
+        effsat.push_back(sp.sat);
+        return;
+    }
+
+    // Nothing to do for IgnorePoint.  In particular, we must not change the
+    // contents or the size of effsat.
+}
+
 // ---------------------------------------------------------------------
 // EPS factory functions for two-point and three-point scaling options
 // ---------------------------------------------------------------------
@@ -332,6 +408,7 @@ namespace Create {
     using EPSOpt = ::Opm::SatFunc::CreateEPS::EPSOptions;
     using RTEP   = ::Opm::SatFunc::CreateEPS::RawTableEndPoints;
     using TEP    = ::Opm::SatFunc::EPSEvalInterface::TableEndPoints;
+    using InvBeh = ::Opm::SatFunc::EPSEvalInterface::InvalidEndpointBehaviour;
 
     namespace TwoPoint {
         using EPS    = ::Opm::SatFunc::TwoPointScaling;
@@ -340,29 +417,35 @@ namespace Create {
         struct Kr {
             static EPSPtr
             G(const ::Opm::ECLGraph&        G,
-              const ::Opm::ECLInitFileData& init);
+              const ::Opm::ECLInitFileData& init,
+              const InvBeh                  handle_invalid);
 
             static EPSPtr
             OG(const ::Opm::ECLGraph&        G,
-               const ::Opm::ECLInitFileData& init);
+               const ::Opm::ECLInitFileData& init,
+               const InvBeh                  handle_invalid);
 
             static EPSPtr
             OW(const ::Opm::ECLGraph&        G,
-               const ::Opm::ECLInitFileData& init);
+               const ::Opm::ECLInitFileData& init,
+               const InvBeh                  handle_invalid);
 
             static EPSPtr
             W(const ::Opm::ECLGraph&        G,
-              const ::Opm::ECLInitFileData& init);
+              const ::Opm::ECLInitFileData& init,
+              const InvBeh                  handle_invalid);
         };
 
         struct Pc {
             static EPSPtr
             GO(const ::Opm::ECLGraph&        G,
-               const ::Opm::ECLInitFileData& init);
+               const ::Opm::ECLInitFileData& init,
+               const InvBeh                  handle_invalid);
 
             static EPSPtr
             OW(const ::Opm::ECLGraph&        G,
-               const ::Opm::ECLInitFileData& init);
+               const ::Opm::ECLInitFileData& init,
+               const InvBeh                  handle_invalid);
         };
 
         static EPSPtr
@@ -381,19 +464,23 @@ namespace Create {
         struct Kr {
             static EPSPtr
             G(const ::Opm::ECLGraph&        G,
-              const ::Opm::ECLInitFileData& init);
+              const ::Opm::ECLInitFileData& init,
+              const InvBeh                  handle_invalid);
 
             static EPSPtr
             OG(const ::Opm::ECLGraph&        G,
-               const ::Opm::ECLInitFileData& init);
+               const ::Opm::ECLInitFileData& init,
+               const InvBeh                  handle_invalid);
 
             static EPSPtr
             OW(const ::Opm::ECLGraph&        G,
-               const ::Opm::ECLInitFileData& init);
+               const ::Opm::ECLInitFileData& init,
+               const InvBeh                  handle_invalid);
 
             static EPSPtr
             W(const ::Opm::ECLGraph&        G,
-              const ::Opm::ECLInitFileData& init);
+              const ::Opm::ECLInitFileData& init,
+              const InvBeh                  handle_invalid);
         };
 
         static EPSPtr
@@ -410,7 +497,8 @@ namespace Create {
 // Implementation of Create::TwoPoint::scalingFunction()
 Create::TwoPoint::EPSPtr
 Create::TwoPoint::Kr::G(const ::Opm::ECLGraph&        G,
-                        const ::Opm::ECLInitFileData& init)
+                        const ::Opm::ECLInitFileData& init,
+                        const InvBeh                  handle_invalid)
 {
     auto sgcr = G.rawLinearisedCellData<double>(init, "SGCR");
     auto sgu  = G.rawLinearisedCellData<double>(init, "SGU");
@@ -424,12 +512,17 @@ Create::TwoPoint::Kr::G(const ::Opm::ECLGraph&        G,
         };
     }
 
-    return EPSPtr { new EPS { std::move(sgcr), std::move(sgu) } };
+    return EPSPtr {
+        new EPS {
+            std::move(sgcr), std::move(sgu), handle_invalid
+        }
+    };
 }
 
 Create::TwoPoint::EPSPtr
 Create::TwoPoint::Kr::OG(const ::Opm::ECLGraph&        G,
-                         const ::Opm::ECLInitFileData& init)
+                         const ::Opm::ECLInitFileData& init,
+                         const InvBeh                  handle_invalid)
 {
     auto sogcr = G.rawLinearisedCellData<double>(init, "SOGCR");
 
@@ -475,12 +568,17 @@ Create::TwoPoint::Kr::OG(const ::Opm::ECLGraph&        G,
         }
     }
 
-    return EPSPtr { new EPS { std::move(sogcr), std::move(smax) } };
+    return EPSPtr {
+        new EPS {
+            std::move(sogcr), std::move(smax), handle_invalid
+        }
+    };
 }
 
 Create::TwoPoint::EPSPtr
 Create::TwoPoint::Kr::OW(const ::Opm::ECLGraph&        G,
-                         const ::Opm::ECLInitFileData& init)
+                         const ::Opm::ECLInitFileData& init,
+                         const InvBeh                  handle_invalid)
 {
     auto sowcr = G.rawLinearisedCellData<double>(init, "SOWCR");
 
@@ -526,12 +624,17 @@ Create::TwoPoint::Kr::OW(const ::Opm::ECLGraph&        G,
         }
     }
 
-    return EPSPtr { new EPS { std::move(sowcr), std::move(smax) } };
+    return EPSPtr {
+        new EPS {
+            std::move(sowcr), std::move(smax), handle_invalid
+        }
+    };
 }
 
 Create::TwoPoint::EPSPtr
 Create::TwoPoint::Kr::W(const ::Opm::ECLGraph&        G,
-                        const ::Opm::ECLInitFileData& init)
+                        const ::Opm::ECLInitFileData& init,
+                        const InvBeh                  handle_invalid)
 {
     auto swcr = G.rawLinearisedCellData<double>(init, "SWCR");
     auto swu  = G.rawLinearisedCellData<double>(init, "SWU");
@@ -542,12 +645,17 @@ Create::TwoPoint::Kr::W(const ::Opm::ECLGraph&        G,
         };
     }
 
-    return EPSPtr { new EPS { std::move(swcr), std::move(swu) } };
+    return EPSPtr {
+        new EPS {
+            std::move(swcr), std::move(swu), handle_invalid
+        }
+    };
 }
 
 Create::TwoPoint::EPSPtr
 Create::TwoPoint::Pc::GO(const ::Opm::ECLGraph&        G,
-                         const ::Opm::ECLInitFileData& init)
+                         const ::Opm::ECLInitFileData& init,
+                         const InvBeh                  handle_invalid)
 {
     auto sgl = G.rawLinearisedCellData<double>(init, "SGL");
     auto sgu = G.rawLinearisedCellData<double>(init, "SGU");
@@ -561,12 +669,17 @@ Create::TwoPoint::Pc::GO(const ::Opm::ECLGraph&        G,
         };
     }
 
-    return EPSPtr { new EPS { std::move(sgl), std::move(sgu) } };
+    return EPSPtr {
+        new EPS {
+            std::move(sgl), std::move(sgu), handle_invalid
+        }
+    };
 }
 
 Create::TwoPoint::EPSPtr
 Create::TwoPoint::Pc::OW(const ::Opm::ECLGraph&        G,
-                         const ::Opm::ECLInitFileData& init)
+                         const ::Opm::ECLInitFileData& init,
+                         const InvBeh                  handle_invalid)
 {
     auto swl = G.rawLinearisedCellData<double>(init, "SWL");
     auto swu = G.rawLinearisedCellData<double>(init, "SWU");
@@ -580,7 +693,11 @@ Create::TwoPoint::Pc::OW(const ::Opm::ECLGraph&        G,
         };
     }
 
-    return EPSPtr { new EPS { std::move(swl), std::move(swu) } };
+    return EPSPtr {
+        new EPS {
+            std::move(swl), std::move(swu), handle_invalid
+        }
+    };
 }
 
 Create::TwoPoint::EPSPtr
@@ -606,10 +723,10 @@ scalingFunction(const ::Opm::ECLGraph&                       G,
             }
 
             if (opt.thisPh == PhIdx::Aqua) {
-                return Create::TwoPoint::Kr::W(G, init);
+                return Create::TwoPoint::Kr::W(G, init, opt.handle_invalid);
             }
 
-            return Create::TwoPoint::Kr::OW(G, init);
+            return Create::TwoPoint::Kr::OW(G, init, opt.handle_invalid);
         }
 
         if (opt.subSys == SSys::OilGas) {
@@ -621,10 +738,10 @@ scalingFunction(const ::Opm::ECLGraph&                       G,
             }
 
             if (opt.thisPh == PhIdx::Vapour) {
-                return Create::TwoPoint::Kr::G(G, init);
+                return Create::TwoPoint::Kr::G(G, init, opt.handle_invalid);
             }
 
-            return Create::TwoPoint::Kr::OG(G, init);
+            return Create::TwoPoint::Kr::OG(G, init, opt.handle_invalid);
         }
     }
 
@@ -637,11 +754,11 @@ scalingFunction(const ::Opm::ECLGraph&                       G,
         }
 
         if (opt.thisPh == PhIdx::Vapour) {
-            return Create::TwoPoint::Pc::GO(G, init);
+            return Create::TwoPoint::Pc::GO(G, init, opt.handle_invalid);
         }
 
         if (opt.thisPh == PhIdx::Aqua) {
-            return Create::TwoPoint::Pc::OW(G, init);
+            return Create::TwoPoint::Pc::OW(G, init, opt.handle_invalid);
         }
     }
 
@@ -726,7 +843,8 @@ Create::TwoPoint::unscaledEndPoints(const RTEP& ep, const EPSOpt& opt)
 
 Create::ThreePoint::EPSPtr
 Create::ThreePoint::Kr::G(const ::Opm::ECLGraph&        G,
-                          const ::Opm::ECLInitFileData& init)
+                          const ::Opm::ECLInitFileData& init,
+                          const InvBeh                  handle_invalid)
 {
     auto sgcr = G.rawLinearisedCellData<double>(init, "SGCR");
     auto sgu  = G.rawLinearisedCellData<double>(init, "SGU");
@@ -777,13 +895,16 @@ Create::ThreePoint::Kr::G(const ::Opm::ECLGraph&        G,
     }
 
     return EPSPtr {
-        new EPS { std::move(sgcr), std::move(sr), std::move(sgu) }
+        new EPS {
+            std::move(sgcr), std::move(sr), std::move(sgu), handle_invalid
+        }
     };
 }
 
 Create::ThreePoint::EPSPtr
 Create::ThreePoint::Kr::OG(const ::Opm::ECLGraph&        G,
-                           const ::Opm::ECLInitFileData& init)
+                           const ::Opm::ECLInitFileData& init,
+                           const InvBeh                  handle_invalid)
 {
     auto sogcr = G.rawLinearisedCellData<double>(init, "SOGCR");
 
@@ -850,13 +971,16 @@ Create::ThreePoint::Kr::OG(const ::Opm::ECLGraph&        G,
     }
 
     return EPSPtr {
-        new EPS { std::move(sogcr), std::move(sdisp), std::move(smax) }
+        new EPS {
+            std::move(sogcr), std::move(sdisp), std::move(smax), handle_invalid
+        }
     };
 }
 
 Create::ThreePoint::EPSPtr
 Create::ThreePoint::Kr::OW(const ::Opm::ECLGraph&        G,
-                           const ::Opm::ECLInitFileData& init)
+                           const ::Opm::ECLInitFileData& init,
+                           const InvBeh                  handle_invalid)
 {
     auto sowcr = G.rawLinearisedCellData<double>(init, "SOWCR");
 
@@ -923,13 +1047,16 @@ Create::ThreePoint::Kr::OW(const ::Opm::ECLGraph&        G,
     }
 
     return EPSPtr {
-        new EPS { std::move(sowcr), std::move(sdisp), std::move(smax) }
+        new EPS {
+            std::move(sowcr), std::move(sdisp), std::move(smax), handle_invalid
+        }
     };
 }
 
 Create::ThreePoint::EPSPtr
 Create::ThreePoint::Kr::W(const ::Opm::ECLGraph&        G,
-                          const ::Opm::ECLInitFileData& init)
+                          const ::Opm::ECLInitFileData& init,
+                          const InvBeh                  handle_invalid)
 {
     auto swcr = G.rawLinearisedCellData<double>(init, "SWCR");
     auto swu  = G.rawLinearisedCellData<double>(init, "SWU");
@@ -979,7 +1106,9 @@ Create::ThreePoint::Kr::W(const ::Opm::ECLGraph&        G,
     }
 
     return EPSPtr {
-        new EPS { std::move(swcr), std::move(sdisp), std::move(swu) }
+        new EPS {
+            std::move(swcr), std::move(sdisp), std::move(swu), handle_invalid
+        }
     };
 }
 
@@ -1008,10 +1137,10 @@ scalingFunction(const ::Opm::ECLGraph&                       G,
         }
 
         if (opt.thisPh == PhIdx::Aqua) {
-            return Create::ThreePoint::Kr::W(G, init);
+            return Create::ThreePoint::Kr::W(G, init, opt.handle_invalid);
         }
 
-        return Create::ThreePoint::Kr::OW(G, init);
+        return Create::ThreePoint::Kr::OW(G, init, opt.handle_invalid);
     }
 
     if (opt.subSys == SSys::OilGas) {
@@ -1023,10 +1152,10 @@ scalingFunction(const ::Opm::ECLGraph&                       G,
         }
 
         if (opt.thisPh == PhIdx::Vapour) {
-            return Create::ThreePoint::Kr::G(G, init);
+            return Create::ThreePoint::Kr::G(G, init, opt.handle_invalid);
         }
 
-        return Create::ThreePoint::Kr::OG(G, init);
+        return Create::ThreePoint::Kr::OG(G, init, opt.handle_invalid);
     }
 
     // Invalid.
@@ -1124,9 +1253,10 @@ Opm::SatFunc::EPSEvalInterface::~EPSEvalInterface()
 
 // Class Opm::SatFunc::TwoPointScaling
 Opm::SatFunc::TwoPointScaling::
-TwoPointScaling(std::vector<double> smin,
+TwoPointScaling(std::vector<double>      smin,
-                std::vector<double> smax)
+                std::vector<double>      smax,
-    : pImpl_(new Impl(std::move(smin), std::move(smax)))
+                InvalidEndpointBehaviour handle_invalid)
+    : pImpl_(new Impl(std::move(smin), std::move(smax), handle_invalid))
 {}
 
 Opm::SatFunc::TwoPointScaling::~TwoPointScaling()
@@ -1182,10 +1312,14 @@ Opm::SatFunc::TwoPointScaling::clone() const
 
 // Class Opm::SatFunc::ThreePointScaling
 Opm::SatFunc::ThreePointScaling::
-ThreePointScaling(std::vector<double> smin,
+ThreePointScaling(std::vector<double>      smin,
-                  std::vector<double> sdisp,
+                  std::vector<double>      sdisp,
-                  std::vector<double> smax)
+                  std::vector<double>      smax,
-    : pImpl_(new Impl(std::move(smin), std::move(sdisp), std::move(smax)))
+                  InvalidEndpointBehaviour handle_invalid)
+    : pImpl_(new Impl(std::move(smin) ,
+                      std::move(sdisp),
+                      std::move(smax) ,
+                      handle_invalid))
 {}
 
 Opm::SatFunc::ThreePointScaling::~ThreePointScaling()

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

@@ -67,6 +67,19 @@ namespace Opm { namespace SatFunc {
             double sat;
         };
 
+        /// Policy for how to handle an invalid end-point scaling (e.g., if
+        /// lower and/or upper scaled saturations have nonsensical values
+        /// like -1.0E+20).
+        enum class InvalidEndpointBehaviour {
+            /// Use the unscaled value for this point.
+            UseUnscaled,
+
+            /// Ignore this scaling request (e.g., produce no graphical
+            /// output for a scaled saturation function when the scaled
+            /// end-points are meaningless).
+            IgnorePoint,
+        };
+
         /// Convenience type alias.
         using SaturationPoints = std::vector<SaturationAssoc>;
 
@@ -128,8 +141,19 @@ namespace Opm { namespace SatFunc {
         /// \param[in] smin Left end points for a set of cells.
         ///
         /// \param[in] smax Right end points for a set of cells.
-        TwoPointScaling(std::vector<double> smin,
+        ///
-                        std::vector<double> smax);
+        /// \param[in] handle_invalid How to treat scaling requests with
+        ///    invalid scaled saturations.  This can, for instance, happen
+        ///    if the scaled saturations are present in the result set but
+        ///    some (or all) cells have irreconcilable values (e.g., minimum
+        ///    saturation greater than maximum saturation, smin < -1E+20,
+        ///    smax < -1E+20).
+        ///
+        ///    Default behaviour: Use unscaled saturation if this happens.
+        TwoPointScaling(std::vector<double>            smin,
+                        std::vector<double>            smax,
+                        const InvalidEndpointBehaviour handle_invalid
+                        = InvalidEndpointBehaviour::UseUnscaled);
 
         /// Destructor.
         ~TwoPointScaling();
@@ -237,9 +261,20 @@ namespace Opm { namespace SatFunc {
         ///    for a set of cells.
         ///
         /// \param[in] smax Right end points for a set of cells.
-        ThreePointScaling(std::vector<double> smin,
+        ///
-                          std::vector<double> sdisp,
+        /// \param[in] handle_invalid How to treat scaling requests with
-                          std::vector<double> smax);
+        ///    invalid scaled saturations.  This can, for instance, happen
+        ///    if the scaled saturations are present in the result set but
+        ///    some (or all) cells have irreconcilable values (e.g., minimum
+        ///    saturation greater than maximum saturation, smin < -1E+20,
+        ///    smax < -1E+20).
+        ///
+        ///    Default behaviour: Use unscaled saturation if this happens.
+        ThreePointScaling(std::vector<double>            smin,
+                          std::vector<double>            sdisp,
+                          std::vector<double>            smax,
+                          const InvalidEndpointBehaviour handle_invalid
+                          = InvalidEndpointBehaviour::UseUnscaled);
 
         /// Destructor.
         ~ThreePointScaling();
@@ -385,6 +420,13 @@ namespace Opm { namespace SatFunc {
             ///   auto eps = CreateEPS::fromECLOutput(G, init, opt);
             /// \endcode
             ::Opm::ECLPhaseIndex thisPh;
+
+            /// How to handle an invalid end-point scaling (e.g., if lower
+            /// and/or upper scaled saturations have nonsensical values like
+            /// -1.0E+20).
+            EPSEvalInterface::InvalidEndpointBehaviour handle_invalid {
+                EPSEvalInterface::InvalidEndpointBehaviour::UseUnscaled
+            };
         };
 
         /// Collection of raw saturation table end points.

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

@@ -25,6 +25,7 @@
 #include <opm/parser/eclipse/Units/Units.hpp>
 
 #include <functional>
+#include <utility>
 
 #include <ert/ecl/ecl_kw_magic.h>
 
@@ -262,15 +263,25 @@ Opm::ECLPVT::PVDx::viscosity(const std::vector<double>& p) const
     });
 }
 
-Opm::FlowDiagnostics::Graph
+Opm::ECLPVT::PVTGraph
 Opm::ECLPVT::PVDx::getPvtCurve(const RawCurve curve) const
 {
     if (curve == RawCurve::SaturatedState) {
         // Not applicable to dry gas or dead oil.  Return empty.
-        return FlowDiagnostics::Graph{};
+        return {};
     }
 
-    return extractRawPVTCurve(this->interp_, curve);
+    auto ret = PVTGraph{};
+
+    auto basic = extractRawPVTCurve(this->interp_, curve);
+
+    // Dead oil/dry gas.  Mixing ratio == 0.
+    ret.mixRat.assign(basic.first.size(), 0.0);
+
+    ret.press = std::move(basic.first);
+    ret.value = std::move(basic.second);
+
+    return ret;
 }
 
 // =====================================================================

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

@@ -301,6 +301,22 @@ namespace Opm { namespace ECLPVT {
         SaturatedState,
     };
 
+    /// Collection of Miscible PVT States and Function Values.
+    ///
+    /// Intended for graphical representation of miscible and immiscible PVT
+    /// functions.
+    struct PVTGraph {
+        /// Pressure values
+        std::vector<double> press;
+
+        /// Mixing ratio.  Typically Rs or Rv.
+        std::vector<double> mixRat;
+
+        /// Function value.  Typically Bg, Bo, mu_g or mu_o.  Usually a copy
+        /// of the mixing ratio for the case of the saturated state curve.
+        std::vector<double> value;
+    };
+
     template <std::size_t N>
     class DenseVector {
     public:
@@ -485,8 +501,8 @@ namespace Opm { namespace ECLPVT {
         ///
         /// \param[in] curve PVT property curve descriptor
         ///
-        /// \return 2D graph for PVT property curve identified by
+        /// \return 2D graph for PVT property curve identified by requests
-        ///    requests represented by \p func.
+        ///    represented by \p func.  Mixing ratio is a vector of zeros.
         ///
         /// Example: Retrieve formation volume factor curve.
         ///
@@ -494,7 +510,7 @@ namespace Opm { namespace ECLPVT {
         ///       const auto graph =
         ///           pvdx.getPvtCurve(ECLPVT::RawCurve::FVF);
         ///    \endcode
-        FlowDiagnostics::Graph getPvtCurve(const RawCurve curve) const;
+        PVTGraph getPvtCurve(const RawCurve curve) const;
 
     private:
         /// Extrapolation policy for property evaluator/interpolant.
@@ -658,7 +674,8 @@ namespace Opm { namespace ECLPVT {
         /// \return Collection of 2D graphs for PVT property curve
         ///    identified by requests represented by \p curve.  One curve
         ///    (vector element) for each tabulated value of the function's
-        ///    primary lookup key.
+        ///    primary lookup key.  Primary key expanded and assigned to the
+        ///    appropriate graph's mixing ratio vector (PVTO convention).
         ///
         /// Example: Retrieve formation volume factor curves.
         ///
@@ -666,7 +683,7 @@ namespace Opm { namespace ECLPVT {
         ///       const auto graph =
         ///           pvtx.getPvtCurve(ECLPVT::RawCurve::FVF);
         ///    \endcode
-        std::vector<FlowDiagnostics::Graph>
+        std::vector<PVTGraph>
         getPvtCurve(const RawCurve curve) const
         {
             if (curve == RawCurve::SaturatedState) {
@@ -828,35 +845,52 @@ namespace Opm { namespace ECLPVT {
             return result;
         }
 
-        std::vector<FlowDiagnostics::Graph>
+        std::vector<PVTGraph>
         mainPvtCurve(const RawCurve curve) const
         {
-            auto ret = std::vector<FlowDiagnostics::Graph>{};
+            // Uses setup appropriate for PVTO.  The PVTG case must be
+            // handled in the caller.
+
+            auto ret = std::vector<PVTGraph>{};
             ret.reserve(this->propInterp_.size());
 
+            auto i = static_cast<std::vector<double>::size_type>(0);
+
             for (const auto& interp : this->propInterp_) {
-                ret.push_back(extractRawPVTCurve(interp, curve));
+                auto basic = extractRawPVTCurve(interp, curve);
+
+                ret.emplace_back();
+                auto& pvtgraph = ret.back();
+
+                pvtgraph.mixRat.assign(basic.first.size(), this->key_[i]);
+
+                pvtgraph.press = std::move(basic.first);
+                pvtgraph.value = std::move(basic.second);
+
+                i += 1;
             }
 
             return ret;
         }
 
-        std::vector<FlowDiagnostics::Graph> saturatedStateCurve() const
+        std::vector<PVTGraph> saturatedStateCurve() const
         {
-            auto y = std::vector<double>{};
+            // Uses setup appropriate for PVTO.  The PVTG case must be
-            y.reserve(this->propInterp_.size());
+            // handled in the caller.
+
+            auto curve = PVTGraph{};
+
+            curve.press.reserve(this->propInterp_.size());
+            curve.mixRat = this->key_;
+            curve.value  = this->key_;
 
             for (const auto& interp : this->propInterp_) {
                 const auto& yi = interp.independentVariable();
 
-                y.push_back(yi[0]);
+                curve.press.push_back(yi[0]);
             }
 
-            return {
+            return { curve };
-                FlowDiagnostics::Graph {
-                    this->key_, std::move(y)
-                }
-            };
         }
     };
 

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

@@ -44,13 +44,8 @@ namespace {
         return pvtnum;
     }
 
-    std::vector<Opm::FlowDiagnostics::Graph> emptyFDGraph()
-    {
-        return { Opm::FlowDiagnostics::Graph{} };
-    }
-
     template <class PVTInterp>
-    std::vector<Opm::FlowDiagnostics::Graph>
+    std::vector<Opm::ECLPVT::PVTGraph>
     rawPvtCurve(const PVTInterp*            pvt,
                 const Opm::ECLPVT::RawCurve curve,
                 const int                   regID)
@@ -58,7 +53,7 @@ namespace {
         if (pvt == nullptr) {
             // Result set does not provide requisite tabulated properties.
             // Return empty.
-            return emptyFDGraph();
+            return {};
         }
 
         return pvt->getPvtCurve(curve, regID);
@@ -144,123 +139,113 @@ namespace {
         return viscosity(pvt, regID, po, rs);
     }
 
-    std::vector<Opm::FlowDiagnostics::Graph>
+    std::vector<Opm::ECLPVT::PVTGraph>
-    convertCurve(std::vector<Opm::FlowDiagnostics::Graph>&&      curves,
+    convertCurve(std::vector<Opm::ECLPVT::PVTGraph>&&            curves,
-                 const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_x,
+                 const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_p,
-                 const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_y)
+                 const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_R,
+                 const Opm::ECLUnits::Convert::PhysicalQuantity& cvrt_f)
     {
         for (auto& curve : curves) {
-            cvrt_x.appliedTo(curve.first);
+            cvrt_p.appliedTo(curve.press);
-            cvrt_y.appliedTo(curve.second);
+            cvrt_R.appliedTo(curve.mixRat);
+            cvrt_f.appliedTo(curve.value);
         }
 
         return curves;
     }
 
-    std::vector<Opm::FlowDiagnostics::Graph>
+    std::vector<Opm::ECLPVT::PVTGraph>
-    convertFvfCurve(std::vector<Opm::FlowDiagnostics::Graph>&& curve,
+    convertFvfCurve(std::vector<Opm::ECLPVT::PVTGraph>&& curve,
-                    const Opm::ECLPhaseIndex                   phase,
+                    const Opm::ECLPhaseIndex             phase,
-                    const Opm::ECLUnits::UnitSystem&           usysFrom,
+                    const Opm::ECLUnits::UnitSystem&     usysFrom,
-                    const Opm::ECLUnits::UnitSystem&           usysTo)
+                    const Opm::ECLUnits::UnitSystem&     usysTo)
     {
         assert ((phase == Opm::ECLPhaseIndex::Liquid) ||
                 (phase == Opm::ECLPhaseIndex::Vapour));
 
+        auto cvrt_p = Opm::ECLUnits::Convert::Pressure();
+        cvrt_p.from(usysFrom).to(usysTo);
+
         if (phase == Opm::ECLPhaseIndex::Liquid) {
-             // Oil FVF.  First column is pressure, second column is Bo.
+             // Oil FVF.  Mixing ratio is Rs, value is Bo.
-            auto cvrt_x = Opm::ECLUnits::Convert::Pressure();
+            auto cvrt_R = Opm::ECLUnits::Convert::DissolvedGasOilRatio();
-            cvrt_x.from(usysFrom).to(usysTo);
+            cvrt_R.from(usysFrom).to(usysTo);
 
-            auto cvrt_y = Opm::ECLUnits::Convert::OilFVF();
+            auto cvrt_f = Opm::ECLUnits::Convert::OilFVF();
-            cvrt_y.from(usysFrom).to(usysTo);
+            cvrt_f.from(usysFrom).to(usysTo);
 
-            return convertCurve(std::move(curve), cvrt_x, cvrt_y);
+            return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f);
         }
 
-        // Gas FVF.  Need to distinguish miscible from immiscible cases.  In
+        // Gas FVF.  Mixing ratio is Rv, value is Bg.
-        // the former, the first column is Rv (vapourised oil/gas ratio) and
+        auto cvrt_R = Opm::ECLUnits::Convert::VaporisedOilGasRatio();
-        // in the second case the first column is the gas pressure.
+        cvrt_R.from(usysFrom).to(usysTo);
-        //
-        // Immiscible case identified by curve.size() <= 1.
 
-        auto cvrt_y = Opm::ECLUnits::Convert::GasFVF();
+        auto cvrt_f = Opm::ECLUnits::Convert::GasFVF();
-        cvrt_y.from(usysFrom).to(usysTo);
+        cvrt_f.from(usysFrom).to(usysTo);
 
-        if (curve.size() <= 1) {
+        return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f);
-            // Immiscible Gas FVF.  First column is Pg.
-            auto cvrt_x = Opm::ECLUnits::Convert::Pressure();
-            cvrt_x.from(usysFrom).to(usysTo);
-
-            return convertCurve(std::move(curve), cvrt_x, cvrt_y);
-        }
-
-        // Miscible Gas FVF.  First column is Rv.
-        auto cvrt_x = Opm::ECLUnits::Convert::VaporisedOilGasRatio();
-        cvrt_x.from(usysFrom).to(usysTo);
-
-        return convertCurve(std::move(curve), cvrt_x, cvrt_y);
     }
 
-    std::vector<Opm::FlowDiagnostics::Graph>
+    std::vector<Opm::ECLPVT::PVTGraph>
-    convertViscosityCurve(std::vector<Opm::FlowDiagnostics::Graph>&& curve,
+    convertViscosityCurve(std::vector<Opm::ECLPVT::PVTGraph>&& curve,
-                          const Opm::ECLPhaseIndex                   phase,
+                          const Opm::ECLPhaseIndex             phase,
-                          const Opm::ECLUnits::UnitSystem&           usysFrom,
+                          const Opm::ECLUnits::UnitSystem&     usysFrom,
-                          const Opm::ECLUnits::UnitSystem&           usysTo)
+                          const Opm::ECLUnits::UnitSystem&     usysTo)
     {
         assert ((phase == Opm::ECLPhaseIndex::Liquid) ||
                 (phase == Opm::ECLPhaseIndex::Vapour));
 
-        // This is the viscosity curve.  Second column is always viscosity
+        auto cvrt_p = Opm::ECLUnits::Convert::Pressure();
-        // irrespective of phase or miscible/immiscible fluids.
+        cvrt_p.from(usysFrom).to(usysTo);
-        auto cvrt_y = Opm::ECLUnits::Convert::Viscosity();
-        cvrt_y.from(usysFrom).to(usysTo);
 
-        if ((phase == Opm::ECLPhaseIndex::Liquid) || (curve.size() <= 1)) {
+        // This is the viscosity curve.  Value is always viscosity
-             // Graph is oil viscosity or immiscible gas viscosity.  First
+        // irrespective of mixing ratios.
-             // column is pressure.
+        auto cvrt_f = Opm::ECLUnits::Convert::Viscosity();
-            auto cvrt_x = Opm::ECLUnits::Convert::Pressure();
+        cvrt_f.from(usysFrom).to(usysTo);
-            cvrt_x.from(usysFrom).to(usysTo);
 
-            return convertCurve(std::move(curve), cvrt_x, cvrt_y);
-        }
-
-        // Miscible Gas viscosity.  First column is Rv (vapourised oil/gas
-        // ratio).
-        auto cvrt_x = Opm::ECLUnits::Convert::VaporisedOilGasRatio();
-        cvrt_x.from(usysFrom).to(usysTo);
-
-        return convertCurve(std::move(curve), cvrt_x, cvrt_y);
-    }
-
-    std::vector<Opm::FlowDiagnostics::Graph>
-    convertSatStateCurve(std::vector<Opm::FlowDiagnostics::Graph>&& curve,
-                         const Opm::ECLPhaseIndex                   phase,
-                         const Opm::ECLUnits::UnitSystem&           usysFrom,
-                         const Opm::ECLUnits::UnitSystem&           usysTo)
-    {
-        assert ((phase == Opm::ECLPhaseIndex::Liquid) ||
-                (phase == Opm::ECLPhaseIndex::Vapour));
-
-        // First column is pressure (Po or Pg).
-        auto cvrt_x = Opm::ECLUnits::Convert::Pressure();
-        cvrt_x.from(usysFrom).to(usysTo);
-
-        // Second column is Rs or Rv depending on 'phase'.
         if (phase == Opm::ECLPhaseIndex::Liquid) {
-            // Saturated state curve for miscible oil.  Second column is Rs
+             // Oil viscosity.  Mixing ratio is Rs.
-            // (dissolved gas/oil ratio).
+            auto cvrt_R = Opm::ECLUnits::Convert::DissolvedGasOilRatio();
-            auto cvrt_y = Opm::ECLUnits::Convert::DissolvedGasOilRatio();
+            cvrt_R.from(usysFrom).to(usysTo);
-            cvrt_y.from(usysFrom).to(usysTo);
 
-            return convertCurve(std::move(curve), cvrt_x, cvrt_y);
+            return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f);
         }
 
-        // Saturated state curve for miscible gas.  Second column is Rv
+        // Gas viscosity.  Mixing ratio is Rv.
-        // (vapourised oil/gas ratio).
+        auto cvrt_R = Opm::ECLUnits::Convert::VaporisedOilGasRatio();
-        auto cvrt_y = Opm::ECLUnits::Convert::VaporisedOilGasRatio();
+        cvrt_R.from(usysFrom).to(usysTo);
-        cvrt_y.from(usysFrom).to(usysTo);
 
-        return convertCurve(std::move(curve), cvrt_x, cvrt_y);
+        return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_f);
+    }
+
+    std::vector<Opm::ECLPVT::PVTGraph>
+    convertSatStateCurve(std::vector<Opm::ECLPVT::PVTGraph>&& curve,
+                         const Opm::ECLPhaseIndex             phase,
+                         const Opm::ECLUnits::UnitSystem&     usysFrom,
+                         const Opm::ECLUnits::UnitSystem&     usysTo)
+    {
+        assert ((phase == Opm::ECLPhaseIndex::Liquid) ||
+                (phase == Opm::ECLPhaseIndex::Vapour));
+
+        auto cvrt_p = Opm::ECLUnits::Convert::Pressure();
+        cvrt_p.from(usysFrom).to(usysTo);
+
+        if (phase == Opm::ECLPhaseIndex::Liquid) {
+            // Saturated state curve for miscible oil.  Mixing ratio (and
+            // function value) is Rs (dissolved gas/oil ratio).
+            auto cvrt_R = Opm::ECLUnits::Convert::DissolvedGasOilRatio();
+            cvrt_R.from(usysFrom).to(usysTo);
+
+            return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_R);
+        }
+
+        // Saturated state curve for miscible gas.  Second column (and
+        // function value) is Rv (vapourised oil/gas ratio).
+        auto cvrt_R = Opm::ECLUnits::Convert::VaporisedOilGasRatio();
+        cvrt_R.from(usysFrom).to(usysTo);
+
+        return convertCurve(std::move(curve), cvrt_p, cvrt_R, cvrt_R);
     }
 }
 
@@ -281,7 +266,7 @@ setOutputUnits(std::unique_ptr<const ECLUnits::UnitSystem> usys)
     this->usys_output_ = std::move(usys);
 }
 
-std::vector<Opm::FlowDiagnostics::Graph>
+std::vector<Opm::ECLPVT::PVTGraph>
 Opm::ECLPVT::ECLPvtCurveCollection::
 getPvtCurve(const RawCurve      curve,
             const ECLPhaseIndex phase,
@@ -290,7 +275,7 @@ getPvtCurve(const RawCurve      curve,
     if (! this->isValidRequest(phase, activeCell)) {
         // Not a supported phase or cell index out of bounds.  Not a valid
         // request so return empty.
-        return emptyFDGraph();
+        return {};
     }
 
     // PVTNUM is traditional one-based region identifier.  Subtract one to
@@ -442,11 +427,11 @@ isValidRequest(const ECLPhaseIndex phase,
         <  this->pvtnum_.size();
 }
 
-std::vector<Opm::FlowDiagnostics::Graph>
+std::vector<Opm::ECLPVT::PVTGraph>
 Opm::ECLPVT::ECLPvtCurveCollection::
-convertToOutputUnits(std::vector<FlowDiagnostics::Graph>&& graph,
+convertToOutputUnits(std::vector<PVTGraph>&& graph,
-                     const RawCurve                        curve,
+                     const RawCurve          curve,
-                     const ECLPhaseIndex                   phase) const
+                     const ECLPhaseIndex     phase) const
 {
     if (this->usys_output_ == nullptr) {
         // No defined system of units for outputs.  Return unconverted (SI).

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

@@ -103,7 +103,7 @@ namespace Opm { namespace ECLPVT {
         ///           pvtCC.getPvtCurve(ECLPVT::RawCurve::Viscosity,
         ///                             ECLPhaseIndex::Vapour, 31415);
         ///    \endcode
-        std::vector<FlowDiagnostics::Graph>
+        std::vector<PVTGraph>
         getPvtCurve(const RawCurve      curve,
                     const ECLPhaseIndex phase,
                     const int           activeCell) const;
@@ -256,10 +256,10 @@ namespace Opm { namespace ECLPVT {
         /// \param[in] phase Phase for which to compute the property curve.
         ///    Must be \code ECLPhaseIndex::Vapour \endcode or \code
         ///    ECLPhaseIndex::Liquid \endcode.
-        std::vector<FlowDiagnostics::Graph>
+        std::vector<PVTGraph>
-        convertToOutputUnits(std::vector<FlowDiagnostics::Graph>&& graph,
+        convertToOutputUnits(std::vector<PVTGraph>&& graph,
-                             const RawCurve                        curve,
+                             const RawCurve          curve,
-                             const ECLPhaseIndex                   phase) const;
+                             const ECLPhaseIndex     phase) const;
     };
 
 }} // Opm::ECLPVT

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

@@ -98,7 +98,7 @@ public:
     viscosity(const std::vector<double>& rv,
               const std::vector<double>& pg) const = 0;
 
-    virtual std::vector<Opm::FlowDiagnostics::Graph>
+    virtual std::vector<Opm::ECLPVT::PVTGraph>
     getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0;
 
     virtual std::unique_ptr<PVxGBase> clone() const = 0;
@@ -133,7 +133,7 @@ public:
         return this->interpolant_.viscosity(pg);
     }
 
-    virtual std::vector<Opm::FlowDiagnostics::Graph>
+    virtual std::vector<Opm::ECLPVT::PVTGraph>
     getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override
     {
         return { this->interpolant_.getPvtCurve(curve) };
@@ -190,10 +190,30 @@ public:
         return this->interp_.viscosity(key, x);
     }
 
-    virtual std::vector<Opm::FlowDiagnostics::Graph>
+    virtual std::vector<Opm::ECLPVT::PVTGraph>
     getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override
     {
-        return this->interp_.getPvtCurve(curve);
+        // Interpolant blindly assigns pressure values to mixing ratios and
+        // vice versa.  Need to switch those.
+
+        auto graphs = this->interp_.getPvtCurve(curve);
+
+        if (curve == Opm::ECLPVT::RawCurve::SaturatedState) {
+            assert ((graphs.size() == 1) && "Logic Error");
+
+            graphs[0].press.swap(graphs[0].mixRat);
+            graphs[0].value = graphs[0].mixRat; // Just for consistency
+        }
+        else {
+            assert ((graphs.size() > 1) && "Logic Error");
+
+            // FVF or viscosity.  Just swap Pg <-> Rv.
+            for (auto& graph : graphs) {
+                graph.press.swap(graph.mixRat);
+            }
+        }
+
+        return graphs;
     }
 
     virtual std::unique_ptr<PVxGBase> clone() const override
@@ -391,7 +411,7 @@ public:
         return this->rhoS_[region];
     }
 
-    std::vector<FlowDiagnostics::Graph>
+    std::vector<PVTGraph>
     getPvtCurve(const RegIdx   region,
                 const RawCurve curve) const;
 
@@ -460,7 +480,7 @@ viscosity(const RegIdx               region,
     return this->eval_[region]->viscosity(rv, pg);
 }
 
-std::vector<Opm::FlowDiagnostics::Graph>
+std::vector<Opm::ECLPVT::PVTGraph>
 Opm::ECLPVT::Gas::Impl::
 getPvtCurve(const RegIdx   region,
             const RawCurve curve) const
@@ -543,7 +563,7 @@ double Opm::ECLPVT::Gas::surfaceMassDensity(const int region) const
     return this->pImpl_->surfaceMassDensity(region);
 }
 
-std::vector<Opm::FlowDiagnostics::Graph>
+std::vector<Opm::ECLPVT::PVTGraph>
 Opm::ECLPVT::Gas::
 getPvtCurve(const RawCurve curve,
             const int      region) const

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

@@ -175,7 +175,7 @@ namespace Opm { namespace ECLPVT {
         ///       const auto graph =
         ///           pvtGas.getPvtCurve(ECLPVT::RawCurve::FVF, 0);
         ///    \endcode
-        std::vector<FlowDiagnostics::Graph>
+        std::vector<PVTGraph>
         getPvtCurve(const RawCurve curve,
                     const int      region) const;
 

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

@@ -95,7 +95,7 @@ public:
     viscosity(const std::vector<double>& rs,
               const std::vector<double>& po) const = 0;
 
-    virtual std::vector<Opm::FlowDiagnostics::Graph>
+    virtual std::vector<Opm::ECLPVT::PVTGraph>
     getPvtCurve(const Opm::ECLPVT::RawCurve curve) const = 0;
 
     virtual std::unique_ptr<PVxOBase> clone() const = 0;
@@ -130,7 +130,7 @@ public:
         return this->interpolant_.viscosity(po);
     }
 
-    virtual std::vector<Opm::FlowDiagnostics::Graph>
+    virtual std::vector<Opm::ECLPVT::PVTGraph>
     getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override
     {
         return { this->interpolant_.getPvtCurve(curve) };
@@ -187,7 +187,7 @@ public:
         return this->interp_.viscosity(key, x);
     }
 
-    virtual std::vector<Opm::FlowDiagnostics::Graph>
+    virtual std::vector<Opm::ECLPVT::PVTGraph>
     getPvtCurve(const Opm::ECLPVT::RawCurve curve) const override
     {
         return this->interp_.getPvtCurve(curve);
@@ -406,7 +406,7 @@ public:
         return this->rhoS_[region];
     }
 
-    std::vector<FlowDiagnostics::Graph>
+    std::vector<PVTGraph>
     getPvtCurve(const RegIdx   region,
                 const RawCurve curve) const;
 
@@ -475,7 +475,7 @@ viscosity(const RegIdx               region,
     return this->eval_[region]->viscosity(rs, po);
 }
 
-std::vector<Opm::FlowDiagnostics::Graph>
+std::vector<Opm::ECLPVT::PVTGraph>
 Opm::ECLPVT::Oil::Impl::
 getPvtCurve(const RegIdx   region,
             const RawCurve curve) const
@@ -558,7 +558,7 @@ double Opm::ECLPVT::Oil::surfaceMassDensity(const int region) const
     return this->pImpl_->surfaceMassDensity(region);
 }
 
-std::vector<Opm::FlowDiagnostics::Graph>
+std::vector<Opm::ECLPVT::PVTGraph>
 Opm::ECLPVT::Oil::
 getPvtCurve(const RawCurve curve,
             const int      region) const

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

@@ -172,7 +172,7 @@ namespace Opm { namespace ECLPVT {
         ///       const auto graph =
         ///           pvtOil.getPvtCurve(ECLPVT::RawCurve::Viscosity, 3);
         ///    \endcode
-        std::vector<FlowDiagnostics::Graph>
+        std::vector<PVTGraph>
         getPvtCurve(const RawCurve curve,
                     const int      region) const;
 

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

@@ -25,6 +25,7 @@
 
 #include <opm/utility/ECLEndPointScaling.hpp>
 #include <opm/utility/ECLGraph.hpp>
+#include <opm/utility/ECLPropertyUnitConversion.hpp>
 #include <opm/utility/ECLPropTable.hpp>
 #include <opm/utility/ECLRegionMapping.hpp>
 #include <opm/utility/ECLResultData.hpp>
@@ -682,9 +683,12 @@ public:
     Impl(Impl&& rhs);
     Impl(const Impl& rhs);
 
-    void init(const ECLGraph&        G,
+    void init(const ECLGraph&          G,
-              const ECLInitFileData& init,
+              const ECLInitFileData&   init,
-              const bool             useEPS);
+              const bool               useEPS,
+              const InvalidEPBehaviour handle_invalid);
+
+    void setOutputUnits(std::unique_ptr<const ECLUnits::UnitSystem> usys);
 
     std::vector<double>
     relperm(const ECLGraph&             G,
@@ -716,14 +720,16 @@ private:
             bool wat;
         };
 
-        void define(const ECLGraph&        G,
+        void define(const ECLGraph&          G,
-                    const ECLInitFileData& init,
+                    const ECLInitFileData&   init,
-                    const RawTEP&          ep,
+                    const RawTEP&            ep,
-                    const bool             use3PtScaling,
+                    const bool               use3PtScaling,
-                    const ActPh&           active)
+                    const ActPh&             active,
+                    const InvalidEPBehaviour handle_invalid)
         {
             auto opt = Create::EPSOptions{};
-            opt.use3PtScaling = use3PtScaling;
+            opt.use3PtScaling  = use3PtScaling;
+            opt.handle_invalid = handle_invalid;
 
             if (active.oil) {
                 this->create_oil_eps(G, init, ep, active, opt);
@@ -1093,11 +1099,15 @@ private:
 
     ECLRegionMapping rmap_;
 
-    std::unique_ptr<Oil::KrFunction>    oil_;
+    std::unique_ptr<const ECLUnits::UnitSystem> usys_internal_{nullptr};
-    std::unique_ptr<Gas::SatFunction>   gas_;
-    std::unique_ptr<Water::SatFunction> wat_;
 
-    std::unique_ptr<EPSEvaluator> eps_;
+    std::unique_ptr<Oil::KrFunction>    oil_{nullptr};
+    std::unique_ptr<Gas::SatFunction>   gas_{nullptr};
+    std::unique_ptr<Water::SatFunction> wat_{nullptr};
+
+    std::unique_ptr<EPSEvaluator> eps_{nullptr};
+
+    std::unique_ptr<const ECLUnits::UnitSystem> usys_output_{nullptr};
 
     void initRelPermInterp(const EPSEvaluator::ActPh& active,
                            const ECLInitFileData&     init,
@@ -1106,7 +1116,8 @@ private:
     void initEPS(const EPSEvaluator::ActPh& active,
                  const bool                 use3PtScaling,
                  const ECLGraph&            G,
-                 const ECLInitFileData&     init);
+                 const ECLInitFileData&     init,
+                 const InvalidEPBehaviour   handle_invalid);
 
     std::vector<double>
     kro(const ECLGraph&       G,
@@ -1224,25 +1235,56 @@ private:
 
 Opm::ECLSaturationFunc::Impl::Impl(const ECLGraph&        G,
                                    const ECLInitFileData& init)
-    : satnum_(satnumVector(G, init))
+    : satnum_       (satnumVector(G, init))
-    , rmap_  (satnum_)
+    , rmap_         (satnum_)
-{
+    , usys_internal_(ECLUnits::internalUnitConventions())
-}
+{}
 
 Opm::ECLSaturationFunc::Impl::Impl(Impl&& rhs)
-    : satnum_(std::move(rhs.satnum_))
+    : satnum_       (std::move(rhs.satnum_))
-    , rmap_  (std::move(rhs.rmap_))
+    , rmap_         (std::move(rhs.rmap_))
-    , oil_   (std::move(rhs.oil_ ))
+    , usys_internal_(std::move(rhs.usys_internal_))
-    , gas_   (std::move(rhs.gas_ ))
+    , oil_          (std::move(rhs.oil_))
-    , wat_   (std::move(rhs.wat_ ))
+    , gas_          (std::move(rhs.gas_))
+    , wat_          (std::move(rhs.wat_))
+    , eps_          (std::move(rhs.eps_))
+    , usys_output_  (std::move(rhs.usys_output_))
 {}
 
+Opm::ECLSaturationFunc::Impl::Impl(const Impl& rhs)
+    : satnum_       (rhs.satnum_)
+    , rmap_         (rhs.rmap_)
+    , usys_internal_(rhs.usys_internal_->clone())
+{
+    if (rhs.oil_) {
+        // Polymorphic object must use clone().
+        this->oil_ = rhs.oil_->clone();
+    }
+
+    if (rhs.gas_) {
+        this->gas_.reset(new Gas::SatFunction(*rhs.gas_));
+    }
+
+    if (rhs.wat_) {
+        this->wat_.reset(new Water::SatFunction(*rhs.wat_));
+    }
+
+    if (rhs.eps_) {
+        this->eps_.reset(new EPSEvaluator(*rhs.eps_));
+    }
+
+    if (rhs.usys_output_) {
+        this->usys_output_ = rhs.usys_output_->clone();
+    }
+}
+
 // ---------------------------------------------------------------------
 
 void
-Opm::ECLSaturationFunc::Impl::init(const ECLGraph&        G,
+Opm::ECLSaturationFunc::Impl::init(const ECLGraph&          G,
-                                   const ECLInitFileData& init,
+                                   const ECLInitFileData&   init,
-                                   const bool             useEPS)
+                                   const bool               useEPS,
+                                   const InvalidEPBehaviour handle_invalid)
 {
     // Extract INTEHEAD from main grid
     const auto& ih   = init.keywordData<int>(INTEHEAD_KW);
@@ -1263,7 +1305,7 @@ Opm::ECLSaturationFunc::Impl::init(const ECLGraph&        G,
                 lh[LOGIHEAD_ALT_ENDPOINT_SCALING_INDEX]);
 
             // Must be called *after* initRelPermInterp().
-            this->initEPS(active, use3PtScaling, G, init);
+            this->initEPS(active, use3PtScaling, G, init, handle_invalid);
         }
     }
 }
@@ -1322,36 +1364,25 @@ void Opm::ECLSaturationFunc::
 Impl::initEPS(const EPSEvaluator::ActPh& active,
               const bool                 use3PtScaling,
               const ECLGraph&            G,
-              const ECLInitFileData&     init)
+              const ECLInitFileData&     init,
+              const InvalidEPBehaviour   handle_invalid)
 {
     const auto ep = this->extractRawTableEndPoints(active);
 
     this->eps_.reset(new EPSEvaluator());
 
-    this->eps_->define(G, init, ep, use3PtScaling, active);
+    this->eps_->define(G, init, ep, use3PtScaling, active, handle_invalid);
 }
 
 // #####################################################################
 
-Opm::ECLSaturationFunc::Impl::Impl(const Impl& rhs)
+void
-    : rmap_(rhs.rmap_)
+Opm::ECLSaturationFunc::Impl::
+setOutputUnits(std::unique_ptr<const ECLUnits::UnitSystem> usys)
 {
-    if (rhs.oil_) {
+    this->usys_output_ = std::move(usys);
-        // Polymorphic object must use clone().
-        this->oil_ = rhs.oil_->clone();
-    }
-
-    if (rhs.gas_) {
-        this->gas_.reset(new Gas::SatFunction(*rhs.gas_));
-    }
-
-    if (rhs.wat_) {
-        this->wat_.reset(new Water::SatFunction(*rhs.wat_));
-    }
 }
 
-// #####################################################################
-
 std::vector<double>
 Opm::ECLSaturationFunc::Impl::
 relperm(const ECLGraph&       G,
@@ -1717,7 +1748,13 @@ pcgoCurve(const ECLRegionMapping&    rmap,
 
     // Region ID 'reg' is traditional, ECL-style one-based region ID
     // (SATNUM).  Subtract one to create valid table index.
-    const auto pc = this->gas_->pcgo(regID - 1, sg_inp);
+    auto pc = this->gas_->pcgo(regID - 1, sg_inp);
+
+    if (this->usys_output_ != nullptr) {
+        ::Opm::ECLUnits::Convert::Pressure()
+            .from(*this->usys_internal_)
+            .to  (*this->usys_output_).appliedTo(pc);
+    }
 
     // FD::Graph == pair<vector<double>, vector<double>>
     return FlowDiagnostics::Graph {
@@ -1831,7 +1868,13 @@ pcowCurve(const ECLRegionMapping&    rmap,
 
     // Region ID 'reg' is traditional, ECL-style one-based region ID
     // (SATNUM).  Subtract one to create valid table index.
-    const auto pc = this->wat_->pcow(regID - 1, sw_inp);
+    auto pc = this->wat_->pcow(regID - 1, sw_inp);
+
+    if (this->usys_output_ != nullptr) {
+        ::Opm::ECLUnits::Convert::Pressure()
+            .from(*this->usys_internal_)
+            .to  (*this->usys_output_).appliedTo(pc);
+    }
 
     // FD::Graph == pair<vector<double>, vector<double>>
     return FlowDiagnostics::Graph {
@@ -1941,12 +1984,13 @@ extractRawTableEndPoints(const EPSEvaluator::ActPh& active) const
 // =====================================================================
 
 Opm::ECLSaturationFunc::
-ECLSaturationFunc(const ECLGraph&        G,
+ECLSaturationFunc(const ECLGraph&          G,
-                  const ECLInitFileData& init,
+                  const ECLInitFileData&   init,
-                  const bool             useEPS)
+                  const bool               useEPS,
+                  const InvalidEPBehaviour handle_invalid)
     : pImpl_(new Impl(G, init))
 {
-    this->pImpl_->init(G, init, useEPS);
+    this->pImpl_->init(G, init, useEPS, handle_invalid);
 }
 
 Opm::ECLSaturationFunc::~ECLSaturationFunc()
@@ -1976,6 +2020,13 @@ Opm::ECLSaturationFunc::operator=(const ECLSaturationFunc& rhs)
     return *this;
 }
 
+void
+Opm::ECLSaturationFunc::
+setOutputUnits(std::unique_ptr<const ECLUnits::UnitSystem> usys)
+{
+    this->pImpl_->setOutputUnits(std::move(usys));
+}
+
 std::vector<double>
 Opm::ECLSaturationFunc::
 relperm(const ECLGraph&       G,

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

@@ -21,7 +21,9 @@
 #define OPM_ECLSATURATIONFUNC_HEADER_INCLUDED
 
 #include <opm/flowdiagnostics/DerivedQuantities.hpp>
+#include <opm/utility/ECLEndPointScaling.hpp>
 #include <opm/utility/ECLPhaseIndex.hpp>
+#include <opm/utility/ECLUnitHandling.hpp>
 
 #include <memory>
 #include <vector>
@@ -77,6 +79,9 @@ namespace Opm {
             ECLPhaseIndex thisPh;
         };
 
+        using InvalidEPBehaviour = ::Opm::SatFunc::
+            EPSEvalInterface::InvalidEndpointBehaviour;
+
         /// Constructor
         ///
         /// \param[in] G Connected topology of current model's active cells.
@@ -96,9 +101,18 @@ namespace Opm {
         ///
         ///    Default value (\c true) means that effects of EPS are
         ///    included if requisite data is present in the INIT result.
-        ECLSaturationFunc(const ECLGraph&        G,
+        ///
-                          const ECLInitFileData& init,
+        /// \param[in] invalidIsUnscaled Whether or not treat invalid scaled
-                          const bool             useEPS = true);
+        ///    saturation end-points (e.g., SWL=-1.0E+20) as unscaled
+        ///    saturations.  True for "treat as unscaled", false for "
+        ///
+        ///    Default value (\c true) means that invalid scalings are
+        ///    treated as unscaled, false
+        ECLSaturationFunc(const ECLGraph&          G,
+                          const ECLInitFileData&   init,
+                          const bool               useEPS = true,
+                          const InvalidEPBehaviour handle_invalid
+                          = InvalidEPBehaviour::UseUnscaled);
 
         /// Destructor.
         ~ECLSaturationFunc();
@@ -137,6 +151,21 @@ namespace Opm {
         /// \return \code *this \endcode.
         ECLSaturationFunc& operator=(const ECLSaturationFunc& rhs);
 
+        /// Define a collection of units of measure for output purposes.
+        ///
+        /// Capillary pressure curves produced by getSatFuncCurve() will be
+        /// reported in the pressure units of this system.  If this function
+        /// is never called (or called with null pointer), then the output
+        /// units are implicitly set to the flow-diagnostics module's
+        /// internal units of measurement (meaning all properties and curves
+        /// will be reported in strict SI units).
+        ///
+        /// \param[in] usys Collection of units of measure for output
+        ///    purposes.  Typically the return value from one of the \code
+        ///    *UnitConvention() \endcode functions of the \code ECLUnits
+        ///    \endcode namespace.
+        void setOutputUnits(std::unique_ptr<const ECLUnits::UnitSystem> usys);
+
         /// Compute relative permeability values in all active cells for a
         /// single phase.
         ///