Added changes to extend the eclipse restart file export to handle hysteresis option in eclipse

opm/output/eclipse/LogiHEAD.hpp

@@ -38,7 +38,9 @@ namespace Opm { namespace RestartIO {
 
         LogiHEAD& variousParam(const bool e300_radial,
                                const bool e100_radial,
-                               const int  nswlmx);
+                               const int  nswlmx,
+			       const bool enableHyster
+			      );
 
         const std::vector<bool>& data() const
         {

opm/output/eclipse/VectorItems/well.hpp

@@ -41,6 +41,7 @@ namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems
             VFPTab = 11, // ID (one-based) of well's current VFP table
 
             item18 = 17, // Unknown
+	    XFlow  = 22,
             item25 = 24, // Unknown
             item32 = 31, // Unknown
             item48 = 47, // Unknown

opm/parser/eclipse/EclipseState/Runspec.hpp

@@ -112,6 +112,69 @@ private:
     int nLatBranchMax;
 };
 
+class EclHysteresisConfig
+{
+public:
+      EclHysteresisConfig();
+      explicit EclHysteresisConfig(const Deck& deck);
+
+
+    /*!
+     * \brief Specify whether hysteresis is enabled or not.
+     */
+    void setEnableHysteresis(bool yesno);
+
+    /*!
+     * \brief Returns whether hysteresis is enabled.
+     */
+    const bool enableHysteresis() const;
+
+    /*!
+     * \brief Set the type of the hysteresis model which is used for capillary pressure.
+     *
+     * -1: capillary pressure hysteresis is disabled
+     * 0: use the Killough model for capillary pressure hysteresis
+     */
+    void setPcHysteresisModel(int value);
+
+    /*!
+     * \brief Return the type of the hysteresis model which is used for capillary pressure.
+     *
+     * -1: capillary pressure hysteresis is disabled
+     * 0: use the Killough model for capillary pressure hysteresis
+     */
+    const int pcHysteresisModel() const;
+
+    /*!
+     * \brief Set the type of the hysteresis model which is used for relative permeability.
+     *
+     * -1: relperm hysteresis is disabled
+     * 0: use the Carlson model for relative permeability hysteresis of the non-wetting
+     *    phase and the drainage curve for the relperm of the wetting phase
+     * 1: use the Carlson model for relative permeability hysteresis of the non-wetting
+     *    phase and the imbibition curve for the relperm of the wetting phase
+     */
+    void setKrHysteresisModel(int value);
+
+    /*!
+     * \brief Return the type of the hysteresis model which is used for relative permeability.
+     *
+     * -1: relperm hysteresis is disabled
+     * 0: use the Carlson model for relative permeability hysteresis
+     */
+    const int krHysteresisModel() const;
+
+private:
+    // enable hysteresis at all
+    bool enableHysteresis_;
+
+    // the capillary pressure and the relperm hysteresis models to be used
+    int pcHysteresisModel_;
+    int krHysteresisModel_;
+};
+
+
+
 class Runspec {
    public:
         explicit Runspec( const Deck& );
@@ -123,6 +186,7 @@ class Runspec {
         const Welldims& wellDimensions() const noexcept;
         const WellSegmentDims& wellSegmentDimensions() const noexcept;
         int eclPhaseMask( ) const noexcept;
+	const EclHysteresisConfig& hysterPar() const noexcept;
 
    private:
         Phases active_phases;
@@ -131,8 +195,10 @@ class Runspec {
         Welldims welldims;
         WellSegmentDims wsegdims;
         UDQParams udq_params;
+	EclHysteresisConfig hystpar;
 };
 
+
 }
 
 #endif // OPM_RUNSPEC_HPP

src/opm/output/eclipse/AggregateWellData.cpp

@@ -340,6 +340,7 @@ namespace {
             iWell[Ix::WType]  = wellType  (well, sim_step);
             iWell[Ix::WCtrl]  = ctrlMode  (well, sim_step);
             iWell[Ix::VFPTab] = wellVFPTab(well, sim_step);
+	    iWell[Ix::XFlow]  = well.getAllowCrossFlow() ? 1 : 0;
 
             // The following items aren't fully characterised yet, but
             // needed for restart of M2.  Will need further refinement.

src/opm/output/eclipse/CreateLogiHead.cpp

@@ -36,9 +36,10 @@ createLogiHead(const EclipseState& es)
 {
     const auto& rspec = es.runspec();
     const auto& wsd   = rspec.wellSegmentDimensions();
+    const auto& hystPar = rspec.hysterPar();
 
     const auto lh = LogiHEAD{}
-        .variousParam(false, false, wsd.maxSegmentedWells())
+        .variousParam(false, false, wsd.maxSegmentedWells(), hystPar.enableHysteresis())
         ;
 
     return lh.data();

src/opm/output/eclipse/LogiHEAD.cpp

@@ -10,7 +10,7 @@ lh_002	=	2	,		//	FALSE
 lh_003	=	3	,		//	FALSE		Flag set to FALSE for a non-radial model, TRUE for a radial model (ECLIPSE 300 and other simulators)
 lh_004	=	4	,		//	FALSE		Flag set to FALSE for a non-radial model, TRUE for a radial model (ECLIPSE 100)
 lh_005	=	5	,		//	FALSE
-lh_006	=	6	,		//	FALSE
+lh_006	=	6	,		//	FALSE 		Flag set to FALSE when no hysteresis, TRUE when hysteresis option is used
 lh_007	=	7	,		//	FALSE
 lh_008	=	8	,		//	FALSE
 lh_009	=	9	,		//	FALSE
@@ -142,12 +142,13 @@ Opm::RestartIO::LogiHEAD::LogiHEAD()
 
 Opm::RestartIO::LogiHEAD&
 Opm::RestartIO::LogiHEAD::
-variousParam(const bool e300_radial, const bool e100_radial, const int nswlmx)
+variousParam(const bool e300_radial, const bool e100_radial, const int nswlmx, const bool enableHyster)
 {
     this -> data_[lh_000] = true;
     this -> data_[lh_001] = true;
     this -> data_[lh_003] = e300_radial;
     this -> data_[lh_004] = e100_radial;
+    this -> data_[lh_006] = enableHyster;
     //this -> data_[lh_016] = true;
     //this -> data_[lh_018] = true;
     //this -> data_[lh_031] = true;

src/opm/parser/eclipse/EclipseState/Runspec.cpp

@@ -111,6 +111,90 @@ WellSegmentDims::WellSegmentDims(const Deck& deck) : WellSegmentDims()
     }
 }
 
+EclHysteresisConfig::EclHysteresisConfig() :
+    enableHysteresis_( false ),
+    pcHysteresisModel_( 0 ),
+    krHysteresisModel_( 0 )
+    {}
+
+
+EclHysteresisConfig::EclHysteresisConfig(const Opm::Deck& deck) : EclHysteresisConfig()
+    {
+        enableHysteresis_ = false;
+
+        if (!deck.hasKeyword("SATOPTS"))
+            return;
+
+        const auto& satoptsItem = deck.getKeyword("SATOPTS").getRecord(0).getItem(0);
+        for (unsigned i = 0; i < satoptsItem.size(); ++i) {
+            std::string satoptsValue = satoptsItem.get< std::string >(0);
+            std::transform(satoptsValue.begin(),
+                           satoptsValue.end(),
+                           satoptsValue.begin(),
+                           ::toupper);
+
+            if (satoptsValue == "HYSTER")
+                enableHysteresis_ = true;
+        }
+
+        // check for the (deprecated) HYST keyword
+        if (deck.hasKeyword("HYST"))
+            enableHysteresis_ = true;
+
+        if (!enableHysteresis_)
+            return;
+
+        if (!deck.hasKeyword("EHYSTR"))
+            throw std::runtime_error("Enabling hysteresis via the HYST parameter for SATOPTS requires the "
+                                     "presence of the EHYSTR keyword");
+
+        const auto& ehystrKeyword = deck.getKeyword("EHYSTR");
+        if (deck.hasKeyword("NOHYKR"))
+            krHysteresisModel_ = -1;
+        else {
+            krHysteresisModel_ = ehystrKeyword.getRecord(0).getItem("relative_perm_hyst").get<int>(0);
+
+            if (krHysteresisModel_ != 0 && krHysteresisModel_ != 1)
+                throw std::runtime_error(
+                    "Only the Carlson relative permeability hystersis models (indicated by '0' or "
+                    "'1' for the second item of the 'EHYSTR' keyword) are supported");
+        }
+
+        // this is slightly screwed: it is possible to specify contradicting hysteresis
+        // models with HYPC/NOHYPC and the fifth item of EHYSTR. Let's ignore that for
+        // now.
+        std::string whereFlag =
+            ehystrKeyword.getRecord(0).getItem("limiting_hyst_flag").getTrimmedString(0);
+        if (deck.hasKeyword("NOHYPC") || whereFlag == "KR")
+            pcHysteresisModel_ = -1;
+        else {
+            // if capillary pressure hysteresis is enabled, Eclipse always uses the
+            // Killough model
+            pcHysteresisModel_ = 0;
+
+            throw std::runtime_error("Capillary pressure hysteresis is not supported yet");
+        }
+    }
+    
+void EclHysteresisConfig::setEnableHysteresis(bool yesno)
+    { enableHysteresis_ = yesno; }
+
+const bool EclHysteresisConfig::enableHysteresis() const 
+    { return enableHysteresis_; }
+    
+void EclHysteresisConfig::setPcHysteresisModel(int value)
+    { pcHysteresisModel_ = value; }
+    
+const int EclHysteresisConfig::pcHysteresisModel() const
+    { return pcHysteresisModel_; }
+    
+void EclHysteresisConfig::setKrHysteresisModel(int value)
+    { krHysteresisModel_ = value; }
+    
+const int EclHysteresisConfig::krHysteresisModel() const
+    { return krHysteresisModel_; }
+
+    
 Runspec::Runspec( const Deck& deck ) :
     active_phases( Phases( deck.hasKeyword( "OIL" ),
                            deck.hasKeyword( "GAS" ),
@@ -123,7 +207,8 @@ Runspec::Runspec( const Deck& deck ) :
     endscale( deck ),
     welldims( deck ),
     wsegdims( deck ),
-    udq_params( deck )
+    udq_params( deck ),
+    hystpar( deck )
 {}
 
 const Phases& Runspec::phases() const noexcept {
@@ -148,6 +233,11 @@ const WellSegmentDims& Runspec::wellSegmentDimensions() const noexcept
     return this->wsegdims;
 }
 
+const EclHysteresisConfig& Runspec::hysterPar() const noexcept
+{
+    return this->hystpar;
+}
+
 /*
   Returns an integer in the range 0...7 which can be used to indicate
   available phases in Eclipse restart and init files.

tests/test_LogiHEAD.cpp

@@ -29,9 +29,10 @@ BOOST_AUTO_TEST_CASE(Radial_Settings_and_Init)
 {
     const auto e300_radial = false;
     const auto e100_radial = true;
+    const auto enableHyster = true;
 
     const auto lh = Opm::RestartIO::LogiHEAD{}
-        .variousParam(e300_radial, e100_radial, 4);
+        .variousParam(e300_radial, e100_radial, 4, enableHyster);
 
     const auto& v = lh.data();
 
@@ -39,6 +40,7 @@ BOOST_AUTO_TEST_CASE(Radial_Settings_and_Init)
     BOOST_CHECK_EQUAL(v[  1], true);  //
     BOOST_CHECK_EQUAL(v[  3], false); // E300 Radial
     BOOST_CHECK_EQUAL(v[  4], true);  // E100 Radial
+    BOOST_CHECK_EQUAL(v[  6], true);  // enableHyster
     BOOST_CHECK_EQUAL(v[ 75], true);  // MS Well Simulation Case
 }