Merge pull request #610 from jalvestad/restart_hyster

Add changes to correct problems/errors with Eclipse restarts from flow  with hysteresis and well crossflow

opm/output/eclipse/DoubHEAD.hpp

@@ -55,7 +55,8 @@ namespace Opm { namespace RestartIO {
         DoubHEAD& timeStamp(const TimeStamp& ts);
 
         DoubHEAD& drsdt(const Schedule&   sched,
-                        const std::size_t lookup_step);
+                        const std::size_t lookup_step,
+			const double      cnvT);
 
         const std::vector<double>& data() const
         {

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
@@ -115,6 +116,10 @@ namespace Opm { namespace RestartIO { namespace Helpers { namespace VectorItems
             BHPTarget      = 6, // Well's bottom hole pressure target
 
             DatumDepth     = 9, // Well's reference depth for BHP
+	    LiqRateTarget_2  = 33, //Well's liquid rate target/limit for a well on WCONINJH control or for a producer
+	    GasRateTarget_2  = 54, //Well's gas rate target/limit for a well on WCONINJH control or for producer
+	    BHPTarget_2	   = 55, //Well's bottom hole pressure target/limit
+
         };
     } // SWell
 

opm/parser/eclipse/EclipseState/Runspec.hpp

@@ -112,6 +112,49 @@ class WellSegmentDims {
     int nLatBranchMax;
 };
 
+class EclHysterConfig
+{
+public:
+    explicit EclHysterConfig(const Deck& deck);
+
+
+    /*!
+     * \brief Specify whether hysteresis is enabled or not.
+     */
+    //void setActive(bool yesno);
+
+    /*!
+     * \brief Returns whether hysteresis is enabled (active).
+     */
+    bool active() const;
+
+    /*!
+     * \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
+     */
+    int pcHysteresisModel() const;
+
+    /*!
+     * \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
+     */
+    int krHysteresisModel() const;
+
+private:
+    // enable hysteresis at all
+    bool activeHyst  { false };
+
+    // the capillary pressure and the relperm hysteresis models to be used
+    int pcHystMod { 0 };
+    int krHystMod { 0 };
+};
+
+
+
 class Runspec {
    public:
         explicit Runspec( const Deck& );
@@ -123,6 +166,7 @@ class Runspec {
         const Welldims& wellDimensions() const noexcept;
         const WellSegmentDims& wellSegmentDimensions() const noexcept;
         int eclPhaseMask( ) const noexcept;
+	const EclHysterConfig& hysterPar() const noexcept;
 
    private:
         Phases active_phases;
@@ -131,8 +175,10 @@ class Runspec {
         Welldims welldims;
         WellSegmentDims wsegdims;
         UDQParams udq_params;
+	EclHysterConfig 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.
@@ -443,11 +444,11 @@ namespace {
                 zero , zero , infty, infty, zero , dflt ,    //  12.. 17  ( 2)
                 infty, infty, infty, infty, infty, zero ,    //  18.. 23  ( 3)
                 one  , zero , zero , zero , zero , zero ,    //  24.. 29  ( 4)
-                zero , one  , zero , infty, zero , zero ,    //  30.. 35  ( 5)
+                zero , one  , zero , zero,  zero , zero ,    //  30.. 35  ( 5)
                 zero , zero , zero , zero , zero , zero ,    //  36.. 41  ( 6)
                 zero , zero , zero , zero , zero , zero ,    //  42.. 47  ( 7)
                 zero , zero , zero , zero , zero , zero ,    //  48.. 53  ( 8)
-                infty, zero , zero , zero , zero , zero ,    //  54.. 59  ( 9)
+                zero,  zero , zero , zero , zero , zero ,    //  54.. 59  ( 9)
                 zero , zero , zero , zero , zero , zero ,    //  60.. 65  (10)
                 zero , zero , zero , zero , zero , zero ,    //  66.. 71  (11)
                 zero , zero , zero , zero , zero , zero ,    //  72.. 77  (12)
@@ -510,11 +511,13 @@ namespace {
                 if ((pp.GasRate != 0.0) || (!predMode)) {
                     sWell[Ix::GasRateTarget] =
                         swprop(M::gas_surface_rate, pp.GasRate);
+		    sWell[Ix::GasRateTarget_2] = sWell[Ix::GasRateTarget];
                 }
 
                 if (pp.LiquidRate != 0.0 || (!predMode)) {
                     sWell[Ix::LiqRateTarget] =
                         swprop(M::liquid_surface_rate, pp.LiquidRate);
+		    sWell[Ix::LiqRateTarget_2] = sWell[Ix::LiqRateTarget];
                 }
                 else  {
                     sWell[Ix::LiqRateTarget] =
@@ -539,6 +542,7 @@ namespace {
 		sWell[Ix::BHPTarget] = pp.BHPLimit != 0.0
                     ? swprop(M::pressure, pp.BHPLimit)
                     : swprop(M::pressure, 1.0*::Opm::unit::atm);
+		sWell[Ix::BHPTarget_2] = sWell[Ix::BHPTarget];
             }
             else if (well.isInjector(sim_step)) {
                 const auto& ip = well.getInjectionProperties(sim_step);
@@ -552,9 +556,11 @@ namespace {
 		    }
 		    if (ip.injectorType == IT::WATER) {
 			sWell[Ix::WatRateTarget] = swprop(M::liquid_surface_rate, ip.surfaceInjectionRate);
-		    }
+			sWell[Ix::LiqRateTarget_2] = sWell[Ix::WatRateTarget];			    
+			}
 		    if (ip.injectorType == IT::GAS) {
 			sWell[Ix::GasRateTarget] = swprop(M::gas_surface_rate, ip.surfaceInjectionRate);
+			sWell[Ix::GasRateTarget_2] = sWell[Ix::GasRateTarget];
 		    }
                 }
 
@@ -569,6 +575,7 @@ namespace {
                 sWell[Ix::BHPTarget] = ip.hasInjectionControl(IP::BHP)
                     ? swprop(M::pressure, ip.BHPLimit)
                     : swprop(M::pressure, 1.0E05*::Opm::unit::psia);
+		sWell[Ix::BHPTarget_2] = sWell[Ix::BHPTarget];
             }
 
             sWell[Ix::DatumDepth] =

src/opm/output/eclipse/CreateDoubHead.cpp

@@ -78,11 +78,12 @@ createDoubHead(const EclipseState& es,
                const std::size_t   lookup_step,
                const double        simTime)
 {
+    const auto& usys = es.getDeckUnitSystem();
     const auto dh = DoubHEAD{}
-        .tuningParameters(sched.getTuning(), lookup_step,
-                          getTimeConv(es.getDeckUnitSystem()))
+        .tuningParameters(sched.getTuning(), lookup_step, 
+			  getTimeConv(usys))
         .timeStamp       (computeTimeStamp(sched, simTime))
-        .drsdt           (sched, lookup_step)
+        .drsdt           (sched, lookup_step, getTimeConv(usys))
         ;
 
     return dh.data();

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.active())
         ;
 
     return lh.data();

src/opm/output/eclipse/DoubHEAD.cpp

@@ -590,14 +590,15 @@ Opm::RestartIO::DoubHEAD::timeStamp(const TimeStamp& ts)
 
 Opm::RestartIO::DoubHEAD&
 Opm::RestartIO::DoubHEAD::drsdt(const Schedule&   sched,
-                                const std::size_t lookup_step)
+                                const std::size_t lookup_step,
+				const double      cnvT)
 {
     const auto& vappar =
         sched.getOilVaporizationProperties(lookup_step);
 
     this->data_[dRsdt] =
         (vappar.getType() == Opm::OilVaporizationEnum::DRDT)
-        ? vappar.getMaxDRSDT(0)
+        ? vappar.getMaxDRSDT(0)*cnvT
         : 1.0e+20;
 
     return *this;

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,87 @@ WellSegmentDims::WellSegmentDims(const Deck& deck) : WellSegmentDims()
     }
 }
 
+EclHysterConfig::EclHysterConfig(const Opm::Deck& deck) 
+    {
+
+        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")
+                activeHyst = true;
+        }
+
+        // check for the (deprecated) HYST keyword
+        if (deck.hasKeyword("HYST"))
+            activeHyst = true;
+
+        if (!activeHyst)
+	      return;
+
+        if (!deck.hasKeyword("EHYSTR"))
+            throw std::runtime_error("Enabling hysteresis via the HYST parameter for SATOPTS requires the "
+                                     "presence of the EHYSTR keyword");
+	    /*!
+	* \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
+	*/
+        const auto& ehystrKeyword = deck.getKeyword("EHYSTR");
+        if (deck.hasKeyword("NOHYKR"))
+            krHystMod = -1;
+        else {
+            krHystMod = ehystrKeyword.getRecord(0).getItem("relative_perm_hyst").get<int>(0);
+
+            if (krHystMod != 0 && krHystMod != 1)
+                throw std::runtime_error(
+                    "Only the Carlson relative permeability hysteresis 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.
+            /*!
+	* \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
+	*/
+        std::string whereFlag =
+            ehystrKeyword.getRecord(0).getItem("limiting_hyst_flag").getTrimmedString(0);
+        if (deck.hasKeyword("NOHYPC") || whereFlag == "KR")
+            pcHystMod = -1;
+        else {
+            // if capillary pressure hysteresis is enabled, Eclipse always uses the
+            // Killough model
+            pcHystMod = 0;
+
+            throw std::runtime_error("Capillary pressure hysteresis is not supported yet");
+        }
+    }
+
+
+bool EclHysterConfig::active() const 
+    { return activeHyst; }
+
+int EclHysterConfig::pcHysteresisModel() const
+    { return pcHystMod; }
+
+int EclHysterConfig::krHysteresisModel() const
+    { return krHystMod; }
+
 Runspec::Runspec( const Deck& deck ) :
     active_phases( Phases( deck.hasKeyword( "OIL" ),
                            deck.hasKeyword( "GAS" ),
@@ -123,7 +204,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 +230,11 @@ const WellSegmentDims& Runspec::wellSegmentDimensions() const noexcept
     return this->wsegdims;
 }
 
+const EclHysterConfig& 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,16 +29,18 @@ 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();
 
     BOOST_CHECK_EQUAL(v[  0], true);  //
     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
 }