some updates for EclDefaultMaterial

it now uses shared pointers to the underlying twophase material law
parameter objects and it can be used with hysteresis.

opm/material/fluidmatrixinteractions/EclDefaultMaterial.hpp

@@ -290,27 +290,55 @@ public:
         typedef MathToolbox<Evaluation> Toolbox;
         typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
 
-        Scalar Swco = params.connateWaterSaturation();
+        Scalar Swco = params.Swl();
 
         Evaluation Sw = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(waterPhaseIdx));
-//        Evaluation So = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(oilPhaseIdx));
         Evaluation Sg = FsToolbox::template toLhs<Evaluation>(fluidState.saturation(gasPhaseIdx));
+
         Sw = Toolbox::min(1.0, Toolbox::max(Swco, Sw));
-        //So = Ewoms::Ad::min(1.0, Toolbox::max(0.0, So));
         Sg = Toolbox::min(1.0, Toolbox::max(0.0, Sg));
 
-        if (Toolbox::value(Sg) + Toolbox::value(Sw) - Swco < 1e-50)
-            return Toolbox::createConstant(1.0); // avoid division by zero
-        else {
-            const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sg + Sw);
-            const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), 1 - Sg - Sw + Swco);
+        Evaluation Sw_ow = Sg + Sw;
+        Evaluation So_go = 1 - Sw_ow + Swco;
+        const Evaluation& kro_ow = OilWaterMaterialLaw::twoPhaseSatKrn(params.oilWaterParams(), Sw_ow);
+        const Evaluation& kro_go = GasOilMaterialLaw::twoPhaseSatKrw(params.gasOilParams(), So_go);
 
+        Evaluation kro;
+        if (std::abs(Toolbox::value(Sg) + Toolbox::value(Sw) - Swco) < 1e-50)
+            kro = kro_ow; // avoid division by zero
+        else {
             const auto& weightOilWater = (Sw - Swco)/(Sg + Sw - Swco);
             const auto& weightGasOil = 1 - weightOilWater;
-
-            Evaluation kro = weightOilWater*kro_ow + weightGasOil*kro_go;
-            return Toolbox::min(1.0, Toolbox::max(0.0, kro));
+            kro = weightOilWater*kro_ow + weightGasOil*kro_go;
         }
+
+        return Toolbox::min(1.0, Toolbox::max(0.0, kro));
+    }
+
+    /*!
+     * \brief Update the hysteresis parameters after a time step.
+     *
+     * This assumes that the nested two-phase material laws are parameters for
+     * EclHysteresisLaw. If they are not, calling this methid will cause a compiler
+     * error. (But not calling it will still work.)
+     */
+    template <class FluidState>
+    static void updateHysteresis(Params &params, const FluidState &fluidState)
+    {
+        typedef MathToolbox<typename FluidState::Scalar> FsToolbox;
+
+        Scalar Swco = params.Swl();
+
+        Scalar Sw = FsToolbox::value(fluidState.saturation(waterPhaseIdx));
+        Scalar Sg = FsToolbox::value(fluidState.saturation(gasPhaseIdx));
+        Sw = std::min(1.0, std::max(Swco, Sw));
+        Sg = std::min(1.0, std::max(0.0, Sg));
+
+        Scalar Sw_ow = Sg + Sw;
+        Scalar So_go = 1 - Sw_ow + Swco;
+
+        params.oilWaterParams().update(/*pcSw=*/Sw, /*krwSw=*/1 - Sg, /*krnSw=*/Sw_ow);
+        params.gasOilParams().update(/*pcSw=*/1 - Sg, /*krwSw=*/So_go, /*krnSw=*/1 - Sg);
     }
 };
 } // namespace Opm

opm/material/fluidmatrixinteractions/EclDefaultMaterialParams.hpp

@@ -26,8 +26,8 @@
 #define OPM_ECL_DEFAULT_MATERIAL_PARAMS_HPP
 
 #include <type_traits>
-
 #include <cassert>
+#include <memory>
 
 namespace Opm {
 
@@ -72,24 +72,36 @@ public:
      * \brief The parameter object for the gas-oil twophase law.
      */
     const GasOilParams& gasOilParams() const
-    { assertFinalized_(); return gasOilParams_; }
+    { assertFinalized_(); return *gasOilParams_; }
+
+    /*!
+     * \brief The parameter object for the gas-oil twophase law.
+     */
+    GasOilParams& gasOilParams()
+    { assertFinalized_(); return *gasOilParams_; }
 
     /*!
      * \brief Set the parameter object for the gas-oil twophase law.
      */
-    void setGasOilParams(const GasOilParams& val)
+    void setGasOilParams(std::shared_ptr<GasOilParams> val)
     { gasOilParams_ = val; }
 
     /*!
      * \brief The parameter object for the oil-water twophase law.
      */
     const OilWaterParams& oilWaterParams() const
-    { assertFinalized_(); return oilWaterParams_; }
+    { assertFinalized_(); return *oilWaterParams_; }
+
+    /*!
+     * \brief The parameter object for the oil-water twophase law.
+     */
+    OilWaterParams& oilWaterParams()
+    { assertFinalized_(); return *oilWaterParams_; }
 
     /*!
      * \brief Set the parameter object for the oil-water twophase law.
      */
-    void setOilWaterParams(const OilWaterParams& val)
+    void setOilWaterParams(std::shared_ptr<OilWaterParams> val)
     { oilWaterParams_ = val; }
 
     /*!
@@ -103,14 +115,14 @@ public:
      * the rock's formation. For our application, this is basically a reduction of the
      * rock's porosity...
      */
-    void setConnateWaterSaturation(Scalar val)
-    { connateWaterSaturation_ = val; }
+    void setSwl(Scalar val)
+    { Swl_ = val; }
 
     /*!
      * \brief Return the saturation of "connate" water.
      */
-    Scalar connateWaterSaturation() const
-    { assertFinalized_(); return connateWaterSaturation_; }
+    Scalar Swl() const
+    { assertFinalized_(); return Swl_; }
 
 private:
 #ifndef NDEBUG
@@ -123,10 +135,10 @@ private:
     { }
 #endif
 
-    GasOilParams gasOilParams_;
-    OilWaterParams oilWaterParams_;
+    std::shared_ptr<GasOilParams> gasOilParams_;
+    std::shared_ptr<OilWaterParams> oilWaterParams_;
 
-    Scalar connateWaterSaturation_;
+    Scalar Swl_;
 };
 } // namespace Opm