This adds bubble/dew point pressure output to legacy flow

opm/autodiff/BlackoilModelBase.hpp

@@ -105,6 +105,9 @@ namespace Opm {
 
             std::vector<double> soMax; // Maximum oil saturation
 
+            std::vector<double> Pb; // Bubble point pressure
+            std::vector<double> Pd; // Dew point pressure
+
             //Hysteresis parameters
             std::vector<double> krnswdc_ow;
             std::vector<double> krnswdc_go;

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -464,6 +464,8 @@ typedef Eigen::Array<double,
         , rsSat(ADB::null())
         , rvSat(ADB::null())
         , soMax()
+        , Pb()
+        , Pd()
         , krnswdc_ow()
         , krnswdc_go()
         , pcswmdc_ow()
@@ -670,6 +672,13 @@ typedef Eigen::Array<double,
                     sd_.soMax = fluid_.satOilMax();
                     fluid_.getGasOilHystParams(sd_.krnswdc_go, sd_.pcswmdc_go, cells_);
                     fluid_.getOilWaterHystParams(sd_.krnswdc_ow, sd_.pcswmdc_ow, cells_);
+
+                    sd_.Pb = fluid_.bubblePointPressure(cells_,
+                            state.temperature.value(),
+                            state.rs.value());
+                    sd_.Pd = fluid_.dewPointPressure(cells_,
+                            state.temperature.value(),
+                            state.rv.value());
                 }
             }
             else {

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -942,16 +942,66 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         }
     }
 
+    // Get max oil saturation vector
     const std::vector<double>& BlackoilPropsAdFromDeck::satOilMax() const
     {
         return satOilMax_;
     }
 
+    // Set max oil saturation vector
     void BlackoilPropsAdFromDeck::setSatOilMax(const std::vector<double>& max_sat) {
         assert(satOilMax_.size() == max_sat.size());
         satOilMax_ = max_sat;
     }
 
+    /// Bubble point pressures
+    std::vector<double> BlackoilPropsAdFromDeck::bubblePointPressure(const Cells& cells,
+            const V& T,
+            const V& rs) const
+    {
+        if (!phase_usage_.phase_used[Gas]) {
+            OPM_THROW(std::runtime_error, "Cannot call bubblePointPressure(): gas phase not active.");
+        }
+        const int n = cells.size();
+        std::vector<double> Pb(n, 0.0);
+
+        for (int i = 0; i < n; ++i) {
+            unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
+            try {
+                Pb[i] = FluidSystem::oilPvt().saturationPressure(pvtRegionIdx, T[i], rs[i]);
+            }
+            catch (const NumericalProblem&) {
+                // Ignore
+            }
+        }
+
+        return Pb;
+    }
+
+    /// Dew point pressures
+    std::vector<double> BlackoilPropsAdFromDeck::dewPointPressure(const Cells& cells,
+            const V& T,
+            const V& rv) const
+    {
+        if (!phase_usage_.phase_used[Gas]) {
+            OPM_THROW(std::runtime_error, "Cannot call dewPointPressure(): gas phase not active.");
+        }
+        const int n = cells.size();
+        std::vector<double> Pd(n, 0.0);
+
+        for (int i = 0; i < n; ++i) {
+            unsigned pvtRegionIdx = cellPvtRegionIdx_[cells[i]];
+            try {
+                Pd[i] = FluidSystem::gasPvt().saturationPressure(pvtRegionIdx, T[i], rv[i]);
+            }
+            catch (const NumericalProblem&) {
+                // Ignore
+            }
+        }
+
+        return Pd;
+    }
+
     /// Set capillary pressure scaling according to pressure diff. and initial water saturation.
     /// \param[in]  saturation Array of n*numPhases saturation values.
     /// \param[in]  pc         Array of n*numPhases capillary pressure values.

opm/autodiff/BlackoilPropsAdFromDeck.hpp

@@ -387,17 +387,30 @@ namespace Opm
                                    std::vector<double>& krnswdc,
                                    const std::vector<int>& cells) const;
 
-        /// Update for max oil saturation.                  
+        /// Update for max oil saturation.
         /// \param[in] saturation Saturations for all phases
         void updateSatOilMax(const std::vector<double>& saturation);
 
-        /// Returns the max oil saturation
+        /// Returns the max oil saturation vector
         const std::vector<double>& satOilMax() const;
 
-        /// Set max oil saturation (for restarting)
+        /// Force set max oil saturation (used for restarting)
-        /// \param[in] max_sat Max oil saturations. Note that this is *only* oil saturations
+        /// \param[in] max_sat Max oil saturations.
+        /// Note that this is a vector of *only* oil saturations (no other phases)
+        /// @see The similar function updateSatOilMax(const std::vector<double>& saturation)
+        /// @see satOilMax()
         void setSatOilMax(const std::vector<double>& max_sat);
 
+        /// Returns the bubble point pressures
+        std::vector<double> bubblePointPressure(const Cells& cells,
+                const V& T,
+                const V& rs) const;
+
+        /// Returns the dew point pressures
+        std::vector<double> dewPointPressure(const Cells& cells,
+                const V& T,
+                const V& rv) const;
+
         /// Set capillary pressure scaling according to pressure diff. and initial water saturation.
         /// \param[in]  saturation Array of n*numPhases saturation values.
         /// \param[in]  pc         Array of n*numPhases capillary pressure values.

opm/autodiff/Compat.cpp

@@ -175,27 +175,27 @@ void solutionToSim( const data::Solution& sol,
         state.getCellData("SSOL") = sol.data("SSOL");
     }
 
-    if ( sol.has("SOMAX")) {
+    if ( sol.has("SOMAX" ) ) {
         state.registerCellData("SOMAX", 1);
         state.getCellData("SOMAX") = sol.data("SOMAX");
     }
 
-    if ( sol.has("PCSWM_OW")) {
+    if ( sol.has("PCSWM_OW" ) ) {
         state.registerCellData("PCSWMDC_OW", 1);
         state.getCellData("PCSWMDC_OW") = sol.data("PCSWM_OW");
     }
 
-    if ( sol.has("KRNSW_OW")) {
+    if ( sol.has("KRNSW_OW" ) ) {
         state.registerCellData("KRNSWMDC_OW", 1);
         state.getCellData("KRNSWMDC_OW") = sol.data("KRNSW_OW");
     }
 
-    if ( sol.has("PCSWM_GO")) {
+    if ( sol.has("PCSWM_GO" ) ) {
         state.registerCellData("PCSWMDC_GO", 1);
         state.getCellData("PCSWMDC_GO") = sol.data("PCSWM_GO");
     }
 
-    if ( sol.has("KRNSW_GO")) {
+    if ( sol.has("KRNSW_GO" ) ) {
         state.registerCellData("KRNSWMDC_GO", 1);
         state.getCellData("KRNSWMDC_GO") = sol.data("KRNSW_GO");
     }

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -546,10 +546,12 @@ namespace Opm
             addToSimData( simData, "RVSAT", sd.rvSat );
 
             addToSimData( simData, "SOMAX", sd.soMax );
-            addToSimData( simData, "PCSWMDC_OW", sd.pcswmdc_ow);
+            addToSimData( simData, "PBUB", sd.Pb );
-            addToSimData( simData, "KRNSWMDC_OW", sd.krnswdc_ow);
+            addToSimData( simData, "PDEW", sd.Pd );
-            addToSimData( simData, "PCSWMDC_GO", sd.pcswmdc_go);
+            addToSimData( simData, "PCSWMDC_OW", sd.pcswmdc_ow );
-            addToSimData( simData, "KRNSWMDC_GO", sd.krnswdc_go);
+            addToSimData( simData, "KRNSWMDC_OW", sd.krnswdc_ow );
+            addToSimData( simData, "PCSWMDC_GO", sd.pcswmdc_go );
+            addToSimData( simData, "KRNSWMDC_GO", sd.krnswdc_go );
 
             return simData;
         }
@@ -772,7 +774,7 @@ namespace Opm
                                   data::TargetType::RESTART_AUXILIARY);
                 }
                 else if (log) {
-                    Opm::OpmLog::warning("Output of bubble point pressure requested but not available in this simulator. Ignoring.");
+                    Opm::OpmLog::warning("Bubble point pressure unavailable", "Output of bubble point pressure requested but not available in this simulator. Ignoring.");
                 }
 
                 if (sd.hasCellData("PDEW")) {
@@ -782,7 +784,7 @@ namespace Opm
                                   data::TargetType::RESTART_AUXILIARY);
                 }
                 else if (log) {
-                    Opm::OpmLog::warning("Output of dew point pressure requested but not available in this simulator. Ignoring.");
+                    Opm::OpmLog::warning("Dew point pressure unavailable", "Output of dew point pressure requested but not available in this simulator. Ignoring.");
                 }
 
             }

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp

@@ -84,6 +84,8 @@ namespace Opm {
                 , rsSat(ADB::null())
                 , rvSat(ADB::null())
                 , soMax() // FIXME: Not handled properly
+                , Pb() //FIXME: Not handled properly
+                , Pd() //FIXME: Not handled properly
                 , krnswdc_ow() // FIXME: Not handled properly
                 , krnswdc_go() // FIXME: Not handled properly
                 , pcswmdc_ow() // FIXME: Not handled properly
@@ -99,6 +101,8 @@ namespace Opm {
             ADB rsSat;
             ADB rvSat;
             std::vector<double> soMax;
+            std::vector<double> Pb;
+            std::vector<double> Pd;
             std::vector<double> krnswdc_ow;
             std::vector<double> krnswdc_go;
             std::vector<double> pcswmdc_ow;