Added calculation of PVT related dynamic properties and display of these in the PVT FVF and Viscosity plots as a point marker. Started using flow diagnostics library's output unit conversion for PVT curves to align output units with the case units. #1990

ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.cpp

@@ -121,7 +121,7 @@ void RigFlowDiagTimeStepResult::addResult(const RigFlowDiagResultAddress& resAdd
 class RigOpmFlowDiagStaticData : public cvf::Object
 {
 public:
-    RigOpmFlowDiagStaticData(const std::string& grid, const std::string& init)
+    RigOpmFlowDiagStaticData(const std::string& grid, const std::string& init, RiaEclipseUnitTools::UnitSystem caseUnitSystem)
     {
         Opm::ECLInitFileData initData(init);
 
@@ -136,9 +136,15 @@ public:
         {
             m_eclPvtCurveCollection.reset(new Opm::ECLPVT::ECLPvtCurveCollection(*m_eclGraph, initData));
 
-            // Tried to set output units, but currently causes crashes when getting PVT curve data
+            // Try and set output unit system to the same system as the eclipse case system
-            //m_eclPvtCurveCollection->setOutputUnits(Opm::ECLUnits::metricUnitConventions());
+            std::unique_ptr<const Opm::ECLUnits::UnitSystem> eclUnitSystem;
-            //m_eclPvtCurveCollection->setOutputUnits(Opm::ECLUnits::fieldUnitConventions());
+            if      (caseUnitSystem == RiaEclipseUnitTools::UNITS_METRIC) eclUnitSystem = Opm::ECLUnits::metricUnitConventions();
+            else if (caseUnitSystem == RiaEclipseUnitTools::UNITS_FIELD)  eclUnitSystem = Opm::ECLUnits::fieldUnitConventions();
+            else if (caseUnitSystem == RiaEclipseUnitTools::UNITS_LAB)    eclUnitSystem = Opm::ECLUnits::labUnitConventions();
+            if (eclUnitSystem)
+            {
+                m_eclPvtCurveCollection->setOutputUnits(eclUnitSystem->clone());
+            }
         }
         catch (...)
         {
@@ -146,6 +152,7 @@ public:
         }
     }
 
+public:
     std::unique_ptr<Opm::ECLGraph>                  m_eclGraph;
     std::vector<double>                             m_poreVolume;
     std::unique_ptr<Opm::FlowDiagnostics::Toolbox>  m_fldToolbox;
@@ -565,7 +572,10 @@ bool RigFlowDiagSolverInterface::ensureStaticDataObjectInstanceCreated()
         std::string initFileName = getInitFileName();
         if (initFileName.empty()) return false;
 
-        m_opmFlowDiagStaticData = new RigOpmFlowDiagStaticData(gridFileName.toStdString(), initFileName);
+        const RigEclipseCaseData* eclipseCaseData = m_eclipseCase->eclipseCaseData();
+        RiaEclipseUnitTools::UnitSystem caseUnitSystem = eclipseCaseData ? eclipseCaseData->unitsType() : RiaEclipseUnitTools::UNITS_UNKNOWN;
+
+        m_opmFlowDiagStaticData = new RigOpmFlowDiagStaticData(gridFileName.toStdString(), initFileName, caseUnitSystem);
     }
 
     return m_opmFlowDiagStaticData.notNull() ? true : false;
@@ -615,7 +625,7 @@ RigFlowDiagSolverInterface::FlowCharacteristicsResultFrame RigFlowDiagSolverInte
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-std::vector<RigFlowDiagSolverInterface::RelPermCurve> RigFlowDiagSolverInterface::calculateRelPermCurvesForActiveCell(size_t activeCellIndex)
+std::vector<RigFlowDiagSolverInterface::RelPermCurve> RigFlowDiagSolverInterface::calculateRelPermCurves(size_t activeCellIndex)
 {
     std::vector<RelPermCurve> retCurveArr;
 
@@ -670,7 +680,7 @@ std::vector<RigFlowDiagSolverInterface::RelPermCurve> RigFlowDiagSolverInterface
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-std::vector<RigFlowDiagSolverInterface::PvtCurve> RigFlowDiagSolverInterface::calculatePvtCurvesForActiveCell(PvtCurveType pvtCurveType, size_t activeCellIndex)
+std::vector<RigFlowDiagSolverInterface::PvtCurve> RigFlowDiagSolverInterface::calculatePvtCurves(PvtCurveType pvtCurveType, size_t activeCellIndex)
 {
     std::vector<PvtCurve> retCurveArr;
 
@@ -710,6 +720,82 @@ std::vector<RigFlowDiagSolverInterface::PvtCurve> RigFlowDiagSolverInterface::ca
     return retCurveArr;
 }
 
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+bool RigFlowDiagSolverInterface::calculatePvtDynamicPropertiesFvf(size_t activeCellIndex, double pressure, double rs, double rv, double* bo, double* bg)
+{
+    if (bo) *bo = HUGE_VAL;
+    if (bg) *bg = HUGE_VAL;
+
+    if (!ensureStaticDataObjectInstanceCreated())
+    {
+        return false;
+    }
+
+    // Bo
+    {
+        std::vector<double> phasePress = { pressure };
+        std::vector<double> mixRatio = { rs };
+        std::vector<double> valArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getDynamicPropertyNative(Opm::ECLPVT::RawCurve::FVF, Opm::ECLPhaseIndex::Liquid, static_cast<int>(activeCellIndex), phasePress, mixRatio);
+        if (valArr.size() > 0)
+        {
+            *bo = valArr[0];
+        }
+    }
+
+    // Bg
+    {
+        std::vector<double> phasePress = { pressure };
+        std::vector<double> mixRatio = { rv };
+        std::vector<double> valArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getDynamicPropertyNative(Opm::ECLPVT::RawCurve::FVF, Opm::ECLPhaseIndex::Vapour, static_cast<int>(activeCellIndex), phasePress, mixRatio);
+        if (valArr.size() > 0)
+        {
+            *bg = valArr[0];
+        }
+    }
+
+    return true;
+}
+
+//--------------------------------------------------------------------------------------------------
+/// 
+//--------------------------------------------------------------------------------------------------
+bool RigFlowDiagSolverInterface::calculatePvtDynamicPropertiesViscosity(size_t activeCellIndex, double pressure, double rs, double rv, double* mu_o, double* mu_g)
+{
+    if (mu_o) *mu_o = HUGE_VAL;
+    if (mu_g) *mu_g = HUGE_VAL;
+
+    if (!ensureStaticDataObjectInstanceCreated())
+    {
+        return false;
+    }
+
+    // mu_o
+    {
+        std::vector<double> phasePress = { pressure };
+        std::vector<double> mixRatio = { rs };
+        std::vector<double> valArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getDynamicPropertyNative(Opm::ECLPVT::RawCurve::Viscosity, Opm::ECLPhaseIndex::Liquid, static_cast<int>(activeCellIndex), phasePress, mixRatio);
+        if (valArr.size() > 0)
+        {
+            *mu_o = valArr[0];
+        }
+    }
+
+    // mu_o
+    {
+        std::vector<double> phasePress = { pressure };
+        std::vector<double> mixRatio = { rv };
+        std::vector<double> valArr = m_opmFlowDiagStaticData->m_eclPvtCurveCollection->getDynamicPropertyNative(Opm::ECLPVT::RawCurve::Viscosity, Opm::ECLPhaseIndex::Vapour, static_cast<int>(activeCellIndex), phasePress, mixRatio);
+        if (valArr.size() > 0)
+        {
+            *mu_g = valArr[0];
+        }
+    }
+
+    return true;
+}
+
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------

ApplicationCode/ReservoirDataModel/RigFlowDiagSolverInterface.h

@@ -116,8 +116,10 @@ public:
                                                                 const std::vector<size_t>& selected_cell_indices,
                                                                 double max_pv_fraction);
 
-    std::vector<RelPermCurve>      calculateRelPermCurvesForActiveCell(size_t activeCellIndex);
+    std::vector<RelPermCurve>      calculateRelPermCurves(size_t activeCellIndex);
-    std::vector<PvtCurve>          calculatePvtCurvesForActiveCell(PvtCurveType pvtCurveType, size_t activeCellIndex);
+    std::vector<PvtCurve>          calculatePvtCurves(PvtCurveType pvtCurveType, size_t activeCellIndex);
+    bool                           calculatePvtDynamicPropertiesFvf(size_t activeCellIndex, double pressure, double rs, double rv, double* bo, double* bg);
+    bool                           calculatePvtDynamicPropertiesViscosity(size_t activeCellIndex, double pressure, double rs, double rv, double* mu_o, double* mu_g);
 
 private:
     std::string getInitFileName() const;

ApplicationCode/UserInterface/RiuPvtPlotPanel.cpp

@@ -132,12 +132,14 @@ void RiuPvtPlotPanel::setPlotDefaults(QwtPlot* plot)
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-void RiuPvtPlotPanel::setPlotData(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& fvfCurveArr, const std::vector<RigFlowDiagSolverInterface::PvtCurve>& viscosityCurveArr, double pressure)
+void RiuPvtPlotPanel::setPlotData(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& fvfCurveArr, const std::vector<RigFlowDiagSolverInterface::PvtCurve>& viscosityCurveArr, FvfDynProps fvfDynProps, ViscosityDynProps viscosityDynProps, double pressure)
 {
     //cvf::Trace::show("RiuPvtPlotPanel::setPlotData()");
 
     m_allFvfCurvesArr = fvfCurveArr;
     m_allViscosityCurvesArr = viscosityCurveArr;
+    m_fvfDynProps = fvfDynProps;
+    m_viscosityDynProps = viscosityDynProps;
     m_pressure = pressure;
 
     plotUiSelectedCurves();
@@ -157,6 +159,8 @@ void RiuPvtPlotPanel::clearPlot()
 
     m_allFvfCurvesArr.clear();
     m_allViscosityCurvesArr.clear();
+    m_fvfDynProps = FvfDynProps();
+    m_viscosityDynProps = ViscosityDynProps();
     m_pressure = HUGE_VAL;
 
     plotUiSelectedCurves();
@@ -199,19 +203,31 @@ void RiuPvtPlotPanel::plotUiSelectedCurves()
     }
 
     QString phaseString = "";
-    if      (phaseToPlot == RigFlowDiagSolverInterface::PvtCurve::GAS)   phaseString = "Gas ";
+    double fvfPointMarkerYValue = HUGE_VAL;
-    else if (phaseToPlot == RigFlowDiagSolverInterface::PvtCurve::OIL)   phaseString = "Oil ";
+    double viscosityPointMarkerYValue = HUGE_VAL;
+    if (phaseToPlot == RigFlowDiagSolverInterface::PvtCurve::GAS)
+    {
+        phaseString = "Gas ";
+        fvfPointMarkerYValue = m_fvfDynProps.bg;
+        viscosityPointMarkerYValue = m_viscosityDynProps.mu_g;
+    }
+    else if (phaseToPlot == RigFlowDiagSolverInterface::PvtCurve::OIL)
+    {
+        phaseString = "Oil ";
+        fvfPointMarkerYValue = m_fvfDynProps.bo;
+        viscosityPointMarkerYValue = m_viscosityDynProps.mu_o;
+    }
 
     {
         const QString plotTitle = phaseString + "Formation Volume Factor";
         const QString yAxisTitle = phaseString + "Formation Volume Factor";
-        plotCurvesInQwt(selectedFvfCurves, m_pressure, plotTitle, yAxisTitle, m_fvfPlot, &m_fvfPlotMarkers);
+        plotCurvesInQwt(selectedFvfCurves, m_pressure, fvfPointMarkerYValue, plotTitle, yAxisTitle, m_fvfPlot, &m_fvfPlotMarkers);
     }
     
     {
         const QString plotTitle = phaseString + "Viscosity";
         const QString yAxisTitle = phaseString + "Viscosity";
-        plotCurvesInQwt(selectedViscosityCurves, m_pressure, plotTitle, yAxisTitle, m_viscosityPlot, &m_viscosityPlotMarkers);
+        plotCurvesInQwt(selectedViscosityCurves, m_pressure, viscosityPointMarkerYValue, plotTitle, yAxisTitle, m_viscosityPlot, &m_viscosityPlotMarkers);
     }
 
     
@@ -220,7 +236,7 @@ void RiuPvtPlotPanel::plotUiSelectedCurves()
 //--------------------------------------------------------------------------------------------------
 /// 
 //--------------------------------------------------------------------------------------------------
-void RiuPvtPlotPanel::plotCurvesInQwt(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& curveArr, double pressure, QString plotTitle, QString yAxisTitle, QwtPlot* plot, std::vector<QwtPlotMarker*>* myPlotMarkers)
+void RiuPvtPlotPanel::plotCurvesInQwt(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& curveArr, double pressure, double pointMarkerYValue, QString plotTitle, QString yAxisTitle, QwtPlot* plot, std::vector<QwtPlotMarker*>* myPlotMarkers)
 {
     plot->detachItems(QwtPlotItem::Rtti_PlotCurve);
 
@@ -268,7 +284,7 @@ void RiuPvtPlotPanel::plotCurvesInQwt(const std::vector<RigFlowDiagSolverInterfa
         QwtPlotMarker* lineMarker = new QwtPlotMarker;
         lineMarker->setXValue(pressure);
         lineMarker->setLineStyle(QwtPlotMarker::VLine);
-        lineMarker->setLinePen(QPen(Qt::black, 1, Qt::DashLine));
+        lineMarker->setLinePen(QPen(QColor(128,0,255), 1, Qt::DashLine));
         lineMarker->setLabel(QString("PRESSURE"));
         lineMarker->setLabelAlignment(Qt::AlignTop | Qt::AlignRight);
         lineMarker->setLabelOrientation(Qt::Vertical);
@@ -276,6 +292,20 @@ void RiuPvtPlotPanel::plotCurvesInQwt(const std::vector<RigFlowDiagSolverInterfa
         myPlotMarkers->push_back(lineMarker);
     }
 
+    if (pressure != HUGE_VAL && pointMarkerYValue != HUGE_VAL)
+    {
+        QwtPlotMarker* pointMarker = new QwtPlotMarker;
+        pointMarker->setValue(pressure, pointMarkerYValue);
+
+        QwtSymbol* symbol = new QwtSymbol(QwtSymbol::Ellipse);
+        symbol->setSize(13, 13);
+        symbol->setPen(QPen(QColor(128, 128, 255), 2));
+        symbol->setBrush(Qt::NoBrush);
+        pointMarker->setSymbol(symbol);
+        pointMarker->attach(plot);
+        myPlotMarkers->push_back(pointMarker);
+    }
+
     plot->setTitle(plotTitle);
 
     plot->setAxisTitle(QwtPlot::xBottom, "Pressure");

ApplicationCode/UserInterface/RiuPvtPlotPanel.h

@@ -22,6 +22,7 @@
 
 #include <QWidget>
 
+#include <cmath>
 #include <memory>
 
 class RiuPvtPlotUpdater;
@@ -40,18 +41,31 @@ class RiuPvtPlotPanel : public QWidget
 {
     Q_OBJECT
 
+public:
+    struct FvfDynProps
+    {
+        double bo = HUGE_VAL;
+        double bg = HUGE_VAL;
+    };
+
+    struct ViscosityDynProps
+    {
+        double mu_o = HUGE_VAL;
+        double mu_g = HUGE_VAL;
+    };
+
 public:
     RiuPvtPlotPanel(QDockWidget* parent);
     virtual ~RiuPvtPlotPanel();
 
-    void                setPlotData(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& fvfCurveArr, const std::vector<RigFlowDiagSolverInterface::PvtCurve>& viscosityCurveArr, double pressure);
+    void                setPlotData(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& fvfCurveArr, const std::vector<RigFlowDiagSolverInterface::PvtCurve>& viscosityCurveArr, FvfDynProps fvfDynProps, ViscosityDynProps viscosityDynProps, double pressure);
     void                clearPlot();
     RiuPvtPlotUpdater*  plotUpdater();
 
 private:
     void            plotUiSelectedCurves();
     static void     setPlotDefaults(QwtPlot* plot);
-    static void     plotCurvesInQwt(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& curveArr, double pressure, QString plotTitle, QString yAxisTitle, QwtPlot* plot, std::vector<QwtPlotMarker*>* myPlotMarkers);
+    static void     plotCurvesInQwt(const std::vector<RigFlowDiagSolverInterface::PvtCurve>& curveArr, double pressure, double pointMarkerYValue, QString plotTitle, QString yAxisTitle, QwtPlot* plot, std::vector<QwtPlotMarker*>* myPlotMarkers);
 
 private slots:
     void            slotPhaseComboCurrentIndexChanged(int);
@@ -59,6 +73,8 @@ private slots:
 private:
     std::vector<RigFlowDiagSolverInterface::PvtCurve>   m_allFvfCurvesArr;
     std::vector<RigFlowDiagSolverInterface::PvtCurve>   m_allViscosityCurvesArr;
+    FvfDynProps                                         m_fvfDynProps;
+    ViscosityDynProps                                   m_viscosityDynProps;
     double                                              m_pressure;
 
     QComboBox*                                          m_phaseComboBox;

ApplicationCode/UserInterface/RiuPvtPlotUpdater.cpp

@@ -134,8 +134,8 @@ bool RiuPvtPlotUpdater::queryDataAndUpdatePlot(const RimEclipseView& eclipseView
         {
             //cvf::Trace::show("Update PVT plot for active cell index: %d", static_cast<int>(activeCellIndex));
 
-            std::vector<RigFlowDiagSolverInterface::PvtCurve> fvfCurveArr = eclipseResultCase->flowDiagSolverInterface()->calculatePvtCurvesForActiveCell(RigFlowDiagSolverInterface::PVT_CT_FVF, activeCellIndex);
+            std::vector<RigFlowDiagSolverInterface::PvtCurve> fvfCurveArr = eclipseResultCase->flowDiagSolverInterface()->calculatePvtCurves(RigFlowDiagSolverInterface::PVT_CT_FVF, activeCellIndex);
-            std::vector<RigFlowDiagSolverInterface::PvtCurve> viscosityCurveArr = eclipseResultCase->flowDiagSolverInterface()->calculatePvtCurvesForActiveCell(RigFlowDiagSolverInterface::PVT_CT_VISCOSITY, activeCellIndex);
+            std::vector<RigFlowDiagSolverInterface::PvtCurve> viscosityCurveArr = eclipseResultCase->flowDiagSolverInterface()->calculatePvtCurves(RigFlowDiagSolverInterface::PVT_CT_VISCOSITY, activeCellIndex);
 
             const size_t timeStepIndex = static_cast<size_t>(eclipseView.currentTimeStep());
 
@@ -153,7 +153,13 @@ bool RiuPvtPlotUpdater::queryDataAndUpdatePlot(const RimEclipseView& eclipseView
             const double cellPressure = pressureAccessor.notNull() ? pressureAccessor->cellScalar(gridLocalCellIndex) : HUGE_VAL;
             //cvf::Trace::show("cellRS = %f  cellRV = %f  cellPressure = %f", cellRS, cellRV, cellPressure);
 
-            plotPanel->setPlotData(fvfCurveArr, viscosityCurveArr, cellPressure);
+            RiuPvtPlotPanel::FvfDynProps fvfDynProps;
+            eclipseResultCase->flowDiagSolverInterface()->calculatePvtDynamicPropertiesFvf(activeCellIndex, cellPressure, cellRS, cellRV, &fvfDynProps.bo, &fvfDynProps.bg);
+
+            RiuPvtPlotPanel::ViscosityDynProps viscosityDynProps;
+            eclipseResultCase->flowDiagSolverInterface()->calculatePvtDynamicPropertiesViscosity(activeCellIndex, cellPressure, cellRS, cellRV, &viscosityDynProps.mu_o, &viscosityDynProps.mu_g);
+
+            plotPanel->setPlotData(fvfCurveArr, viscosityCurveArr, fvfDynProps, viscosityDynProps, cellPressure);
 
             return true;
         }

ApplicationCode/UserInterface/RiuRelativePermeabilityPlotUpdater.cpp

@@ -134,7 +134,7 @@ bool RiuRelativePermeabilityPlotUpdater::queryDataAndUpdatePlot(const RimEclipse
         {
             //cvf::Trace::show("Updating RelPerm plot for active cell index: %d", static_cast<int>(activeCellIndex));
 
-            std::vector<RigFlowDiagSolverInterface::RelPermCurve> relPermCurveArr = eclipseResultCase->flowDiagSolverInterface()->calculateRelPermCurvesForActiveCell(activeCellIndex);
+            std::vector<RigFlowDiagSolverInterface::RelPermCurve> relPermCurveArr = eclipseResultCase->flowDiagSolverInterface()->calculateRelPermCurves(activeCellIndex);
 
             // Make sure we load the results that we'll query below
             RigCaseCellResultsData* cellResultsData = eclipseCaseData->results(RiaDefines::MATRIX_MODEL);