Fixed hysteresis input/output in flow_legacy

opm/autodiff/BlackoilModelBase.hpp

@@ -105,6 +105,12 @@ namespace Opm {
 
             std::vector<double> soMax; // Maximum oil saturation
 
+            //Hysteresis parameters
+            std::vector<double> krnswdc_ow;
+            std::vector<double> krnswdc_go;
+            std::vector<double> pcswmdc_ow;
+            std::vector<double> pcswmdc_go;
+
             std::array<V, fipValues> fip;
         };
 

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -464,6 +464,10 @@ typedef Eigen::Array<double,
         , rsSat(ADB::null())
         , rvSat(ADB::null())
         , soMax()
+        , krnswdc_ow()
+        , krnswdc_go()
+        , pcswmdc_ow()
+        , pcswmdc_go()
         , fip()
     {
     }
@@ -664,6 +668,8 @@ typedef Eigen::Array<double,
                         state.rv = sd_.rvSat;
                     }
                     sd_.soMax = fluid_.satOilMax();
+                    fluid_.getGasOilHystParams(sd_.krnswdc_go, sd_.pcswmdc_go, cells_);
+                    fluid_.getOilWaterHystParams(sd_.krnswdc_ow, sd_.pcswmdc_ow, cells_);
                 }
             }
             else {

opm/autodiff/BlackoilModelEbos.hpp

@@ -1290,15 +1290,17 @@ namespace Opm {
 
                 somax[cellIdx] = ebosSimulator().model().maxOilSaturation(cellIdx);
 
-                const auto matLawManager = ebosSimulator().problem().materialLawManager();
+                const auto& matLawManager = ebosSimulator().problem().materialLawManager();
-                matLawManager->oilWaterHysteresisParams(
+                if (matLawManager->enableHysteresis()) {
-                        pcSwMdc_ow[cellIdx],
+                    matLawManager->oilWaterHysteresisParams(
-                        krnSwMdc_ow[cellIdx],
+                            pcSwMdc_ow[cellIdx],
-                        cellIdx);
+                            krnSwMdc_ow[cellIdx],
-                matLawManager->gasOilHysteresisParams(
+                            cellIdx);
-                        pcSwMdc_go[cellIdx],
+                    matLawManager->gasOilHysteresisParams(
-                        krnSwMdc_go[cellIdx],
+                            pcSwMdc_go[cellIdx],
-                        cellIdx);
+                            krnSwMdc_go[cellIdx],
+                            cellIdx);
+                }
 
                 if (aqua_active) {
                     saturation[ satIdx + aqua_pos ] = fs.saturation(FluidSystem::waterPhaseIdx).value();

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -876,6 +876,58 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         satprops_->updateSatHyst(n, cells.data(), saturation.data());
     }
 
+    /// Set gas-oil hysteresis parameters
+    /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+    /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+    void BlackoilPropsAdFromDeck::setGasOilHystParams(const std::vector<double>& pcswmdc,
+                             const std::vector<double>& krnswdc,
+                             const std::vector<int>& cells)
+    {
+        const int n = cells.size();
+        assert(pcswmdc.size() == n);
+        assert(krnswdc.size() == n);
+        satprops_->setGasOilHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
+    }
+
+    /// Get gas-oil hysteresis parameters
+    /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+    /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+    void BlackoilPropsAdFromDeck::getGasOilHystParams(std::vector<double>& pcswmdc,
+                             std::vector<double>& krnswdc,
+                             const std::vector<int>& cells) const
+    {
+        const int n = cells.size();
+        pcswmdc.resize(n);
+        krnswdc.resize(n);
+        satprops_->getGasOilHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
+    }
+
+    /// Set oil-water hysteresis parameters
+    /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+    /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+    void BlackoilPropsAdFromDeck::setOilWaterHystParams(const std::vector<double>& pcswmdc,
+                               const std::vector<double>& krnswdc,
+                               const std::vector<int>& cells)
+    {
+        const int n = cells.size();
+        assert(pcswmdc.size() == n);
+        assert(krnswdc.size() == n);
+        satprops_->setOilWaterHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
+    }
+
+    /// Get oil-water hysteresis parameters
+    /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+    /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+    void BlackoilPropsAdFromDeck::getOilWaterHystParams(std::vector<double>& pcswmdc,
+                               std::vector<double>& krnswdc,
+                               const std::vector<int>& cells) const
+    {
+        const int n = cells.size();
+        pcswmdc.resize(n);
+        krnswdc.resize(n);
+        satprops_->getOilWaterHystParams(n, cells.data(), pcswmdc.data(), krnswdc.data());
+    }
+
     /// Update for max oil saturation.
     void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
     {
@@ -895,6 +947,11 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         return satOilMax_;
     }
 
+    void BlackoilPropsAdFromDeck::setSatOilMax(const std::vector<double>& max_sat) {
+        assert(satOilMax_.size() == max_sat.size());
+        satOilMax_ = max_sat;
+    }
+
     /// 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

@@ -359,12 +359,45 @@ namespace Opm
         void updateSatHyst(const std::vector<double>& saturation,
                            const std::vector<int>& cells);
 
+        /// Set gas-oil hysteresis parameters
+        /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+        /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+        void setGasOilHystParams(const std::vector<double>& pcswmdc,
+                                 const std::vector<double>& krnswdc,
+                                 const std::vector<int>& cells);
+
+        /// Get gas-oil hysteresis parameters
+        /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+        /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+        void getGasOilHystParams(std::vector<double>& pcswmdc,
+                                 std::vector<double>& krnswdc,
+                                 const std::vector<int>& cells) const;
+
+        /// Set oil-water hysteresis parameters
+        /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+        /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+        void setOilWaterHystParams(const std::vector<double>& pcswmdc,
+                                   const std::vector<double>& krnswdc,
+                                   const std::vector<int>& cells);
+
+        /// Set oil-water hysteresis parameters
+        /// \param[in]  pcswmdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::pcSwMdc(...))
+        /// \param[in]  krnswdc  Vector of hysteresis parameters (@see EclHysteresisTwoPhaseLawParams::krnSwMdc(...))
+        void getOilWaterHystParams(std::vector<double>& pcswmdc,
+                                   std::vector<double>& krnswdc,
+                                   const std::vector<int>& cells) const;
+
         /// Update for max oil saturation.                  
+        /// \param[in] saturation Saturations for all phases
         void updateSatOilMax(const std::vector<double>& saturation);
 
         /// Returns the max oil saturation
         const std::vector<double>& satOilMax() const;
 
+        /// Set max oil saturation (for restarting)
+        /// \param[in] max_sat Max oil saturations. Note that this is *only* oil saturations
+        void setSatOilMax(const std::vector<double>& max_sat);
+
         /// 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/SimulatorBase.hpp

@@ -186,6 +186,7 @@ namespace Opm
                                    const WellState& well_state,
                                    DynamicListEconLimited& list_econ_limited) const;
 
+        void initHysteresisParams(ReservoirState& state);
 
         // Data.
         typedef RateConverter::

opm/autodiff/SimulatorBase_impl.hpp

@@ -94,6 +94,7 @@ namespace Opm
             // This is a restart, populate WellState and ReservoirState state objects from restart file
             output_writer_.initFromRestartFile(props_.phaseUsage(), grid_, state, prev_well_state, extra);
             initHydroCarbonState(state, props_.phaseUsage(), Opm::UgGridHelpers::numCells(grid_), has_disgas_, has_vapoil_);
+            initHysteresisParams(state);
         }
 
         // Create timers and file for writing timing info.
@@ -844,4 +845,25 @@ namespace Opm
                                                           well_state, list_econ_limited);
     }
 
+    template <class Implementation>
+    void
+    SimulatorBase<Implementation>::
+    initHysteresisParams(ReservoirState& state)
+    {
+        typedef std::vector<double> VectorType;
+
+        const VectorType& somax = state.getCellData( "SOMAX" );
+
+        VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
+        VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
+
+        VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
+        VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
+
+        props_.setSatOilMax(somax);
+        props_.setOilWaterHystParams(pcSwMdc_ow, krnSwMdc_ow, allcells_);
+        props_.setGasOilHystParams(pcSwMdc_go, krnSwMdc_go, allcells_);
+    }
+
+
 } // namespace Opm

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -797,22 +797,25 @@ protected:
             ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx);
         }
 
-        VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
+        if (ebosSimulator_.problem().materialLawManager()->enableHysteresis()) {
-        VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
-
-        VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
-        VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
-
-        for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
             auto matLawManager = ebosSimulator_.problem().materialLawManager();
-            matLawManager->setOilWaterHysteresisParams(
+
-                    pcSwMdc_ow[cellIdx],
+            VectorType& pcSwMdc_ow = state.getCellData( "PCSWMDC_OW" );
-                    krnSwMdc_ow[cellIdx],
+            VectorType& krnSwMdc_ow = state.getCellData( "KRNSWMDC_OW" );
-                    cellIdx);
+
-            matLawManager->setGasOilHysteresisParams(
+            VectorType& pcSwMdc_go = state.getCellData( "PCSWMDC_GO" );
-                    pcSwMdc_go[cellIdx],
+            VectorType& krnSwMdc_go = state.getCellData( "KRNSWMDC_GO" );
-                    krnSwMdc_go[cellIdx],
+
-                    cellIdx);
+            for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
+                matLawManager->setOilWaterHysteresisParams(
+                        pcSwMdc_ow[cellIdx],
+                        krnSwMdc_ow[cellIdx],
+                        cellIdx);
+                matLawManager->setGasOilHysteresisParams(
+                        pcSwMdc_go[cellIdx],
+                        krnSwMdc_go[cellIdx],
+                        cellIdx);
+            }
         }
     }
 

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -556,6 +556,10 @@ namespace Opm
             addToSimData( simData, "RVSAT", sd.rvSat );
 
             addToSimData( simData, "SOMAX", sd.soMax );
+            addToSimData( simData, "PCSWMDC_OW", sd.pcswmdc_ow);
+            addToSimData( simData, "KRNSWMDC_OW", sd.krnswdc_ow);
+            addToSimData( simData, "PCSWMDC_GO", sd.pcswmdc_go);
+            addToSimData( simData, "KRNSWMDC_GO", sd.krnswdc_go);
 
             return simData;
         }

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp

@@ -82,7 +82,11 @@ namespace Opm {
                 : rq(num_phases)
                 , rsSat(ADB::null())
                 , rvSat(ADB::null())
-                , soMax()
+                , soMax() // FIXME: Not handled properly
+                , krnswdc_ow() // FIXME: Not handled properly
+                , krnswdc_go() // FIXME: Not handled properly
+                , pcswmdc_ow() // FIXME: Not handled properly
+                , pcswmdc_go() // FIXME: Not handled properly
                 , fip()
             {
             }
@@ -94,6 +98,10 @@ namespace Opm {
             ADB rsSat;
             ADB rvSat;
             std::vector<double> soMax;
+            std::vector<double> krnswdc_ow;
+            std::vector<double> krnswdc_go;
+            std::vector<double> pcswmdc_ow;
+            std::vector<double> pcswmdc_go;
             std::array<V, fipValues> fip;
         };