SWATINIT: Initialisation and capillary pressure scaling.

opm/core/props/BlackoilPropertiesBasic.cpp

@@ -239,5 +239,16 @@ namespace Opm
     }
 
 
+    /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+    /// \param[in]     cell   Cell index. 
+    /// \param[in]     pcow   P_oil - P_water.
+    /// \param[in/out] swat   Water saturation. / Possibly modified Water saturation.      
+    void BlackoilPropertiesBasic::swatInitScaling(const int cell, 
+                                                  const double pcow, 
+                                                  double & swat)
+    {
+        OPM_THROW(std::runtime_error, "BlackoilPropertiesBasic::swatInitScaling()  --  not implemented.");
+    }
+
 } // namespace Opm
 

opm/core/props/BlackoilPropertiesBasic.hpp

@@ -175,6 +175,16 @@ namespace Opm
                               double* smin,
                               double* smax) const;
 
+
+        /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+        /// \param[in]     cell   Cell index.
+        /// \param[in]     pcow   P_oil - P_water.
+        /// \param[in/out] swat   Water saturation. / Possibly modified Water saturation.      
+        virtual  void swatInitScaling(const int cell,
+                                      const double pcow, 
+                                      double & swat);
+
+
     private:
         RockBasic rock_;
 	PvtPropertiesBasic pvt_;

opm/core/props/BlackoilPropertiesFromDeck.cpp

@@ -305,5 +305,16 @@ namespace Opm
     }
  
  
+    /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+    /// \param[in]     cell   Cell index. 
+    /// \param[in]     pcow   P_oil - P_water.
+    /// \param[in/out] swat   Water saturation. / Possibly modified Water saturation.
+    void BlackoilPropertiesFromDeck::swatInitScaling(const int cell,
+                                                     const double pcow, 
+                                                     double & swat)
+    {
+        satprops_->swatInitScaling(cell, pcow, swat);
+    }
+
 } // namespace Opm
 

opm/core/props/BlackoilPropertiesFromDeck.hpp

@@ -216,6 +216,15 @@ namespace Opm
                               double* smin,
                               double* smax) const;
                       
+                     
+        /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+        /// \param[in]     cell   Cell index. 
+        /// \param[in]     pcow   P_oil - P_water.
+        /// \param[in/out] swat   Water saturation. / Possibly modified Water saturation.
+        virtual  void swatInitScaling(const int cell,
+                                      const double pcow, 
+                                      double & swat);
+
     private:
         int getTableIndex_(const int* pvtTableIdx, int cellIdx) const
         {

opm/core/props/BlackoilPropertiesInterface.hpp

@@ -160,6 +160,16 @@ namespace Opm
                               const int* cells,
                               double* smin,
                               double* smax) const = 0;
+
+
+        /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+        /// \param[in]     cell  Cell index.
+        /// \param[in]     pcow  P_oil - P_water.
+        /// \param[in/out] swat  Water saturation. / Possibly modified Water saturation.      
+        virtual  void swatInitScaling(const int cell, 
+                                      const double pcow, 
+                                      double & swat) = 0;
+
     };
 
 

opm/core/props/satfunc/SaturationPropsFromDeck.hpp

@@ -132,6 +132,14 @@ namespace Opm
                            const int* cells,
                            const double* s);
 
+        /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+        /// \param[in]     cell  Cell index. 
+        /// \param[in]     pcow  P_oil - P_water.
+        /// \param[in/out] swat  Water saturation. / Possibly modified Water saturation.        
+        void swatInitScaling(const int cell, 
+                             const double pcow, 
+                             double & swat);
+
     private:
         PhaseUsage phase_usage_;
         std::vector<SatFuncSet> satfuncset_;

opm/core/props/satfunc/SaturationPropsFromDeck_impl.hpp

@@ -412,6 +412,39 @@ namespace Opm
         } 
     }
 
+
+    /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+    /// \param[in]     cell  Cell index. 
+    /// \param[in]     pcow  P_oil - P_water.
+    /// \param[in/out] swat  Water saturation. / Possibly modified Water saturation.    
+    template <class SatFuncSet>
+    void SaturationPropsFromDeck<SatFuncSet>::swatInitScaling(const int cell, 
+                                                              const double pcow, 
+                                                              double & swat)
+    {
+        if (phase_usage_.phase_used[Aqua]) {
+            // TODO: Mixed wettability systems - see ecl kw OPTIONS switch 74
+            if (swat <= eps_transf_[cell].wat.smin) {
+                swat = eps_transf_[cell].wat.smin;
+            } else if (pcow < 1.0e-8) {
+                swat = eps_transf_[cell].wat.smax;
+            } else {
+                const int wpos = phase_usage_.phase_pos[BlackoilPhases::Aqua];
+                const int np = phase_usage_.num_phases;
+                double s[np];
+                s[wpos] = swat; 
+                double pc[np];
+                funcForCell(cell).evalPc(s, pc, &(eps_transf_[cell]));          
+                if (pc[wpos] > 1.0e-8) {
+                    eps_transf_[cell].wat.pcFactor *= pcow/pc[wpos];
+                }
+            }
+        } else {
+            OPM_THROW(std::runtime_error, "swatInitScaling: no water phase! ");
+        }
+    }
+
+
     // Map the cell number to the correct function set.
     template <class SatFuncSet>
     const typename SaturationPropsFromDeck<SatFuncSet>::Funcs&

opm/core/props/satfunc/SaturationPropsInterface.hpp

@@ -81,6 +81,13 @@ namespace Opm
                                    const int* cells,
                                    const double* s) = 0;
  
+        /// Update capillary pressure scaling according to pressure diff. and initial water saturation.
+        /// \param[in]     cell  Cell index.
+        /// \param[in]     pcow  P_oil - P_water.
+        /// \param[in/out] swat  Water saturation. / Possibly modified Water saturation.      
+        virtual void swatInitScaling(const int cell, 
+                                     const double pcow, 
+                                     double & swat) = 0;
 
     };
 

opm/core/simulator/initStateEquil.hpp

@@ -150,8 +150,9 @@ namespace Opm
         std::vector< std::vector<double> >
         phaseSaturations(const Region&           reg,
                          const CellRange&        cells,
-                         const BlackoilPropertiesInterface& props,
+                         BlackoilPropertiesInterface& props,
-                         const std::vector< std::vector<double> >& phase_pressures);
+                         const std::vector<double> swat_init,
+                         std::vector< std::vector<double> >& phase_pressures);
 
 
 
@@ -254,7 +255,7 @@ namespace Opm
 
             class InitialStateComputer {
             public:
-                InitialStateComputer(const BlackoilPropertiesInterface& props,
+                InitialStateComputer(BlackoilPropertiesInterface& props,
                                      const Opm::DeckConstPtr            deck,
                                      const UnstructuredGrid&            G    ,
                                      const double                       grav = unit::gravity)
@@ -334,6 +335,19 @@ namespace Opm
                         }
                     }
                     
+
+                    // Check for presence of kw SWATINIT
+                    if (deck->hasKeyword("SWATINIT")) {
+                        const std::vector<double>& swat_init = deck->getKeyword("SWATINIT")->getSIDoubleData();
+                        swat_init_.resize(G.number_of_cells);
+                        const int* gc = G.global_cell;
+                        for (int c = 0; c < G.number_of_cells; ++c) {
+                            const int deck_pos = (gc == NULL) ? c : gc[c];
+                            swat_init_[c] = swat_init[deck_pos];
+                        }
+
+                    }
+
                     // Compute pressures, saturations, rs and rv factors.
                     calcPressSatRsRv(eqlmap, rec, props, G, grav);
 
@@ -360,15 +374,18 @@ namespace Opm
                 PVec sat_;
                 Vec rs_;
                 Vec rv_;
+                Vec swat_init_;
 
                 template <class RMap>
                 void
                 calcPressSatRsRv(const RMap&                       reg  ,
                                  const std::vector< EquilRecord >& rec  ,
-                                 const Opm::BlackoilPropertiesInterface& props,
+                                 Opm::BlackoilPropertiesInterface& props,
                                  const UnstructuredGrid&           G    ,
                                  const double grav)
                 {
+                    typedef Miscibility::NoMixing NoMix;
+
                     for (typename RMap::RegionId
                              r = 0, nr = reg.numRegions();
                          r < nr; ++r)
@@ -383,7 +400,7 @@ namespace Opm
                    
                         PVec press = phasePressures(G, eqreg, cells, grav);
 
-                        const PVec sat = phaseSaturations(eqreg, cells, props, press);
+                        const PVec sat = phaseSaturations(eqreg, cells, props, swat_init_, press);
 
                         const int np = props.numPhases();
                         for (int p = 0; p < np; ++p) {

opm/core/simulator/initStateEquil_impl.hpp

@@ -579,7 +579,8 @@ namespace Opm
         std::vector< std::vector<double> >
         phaseSaturations(const Region&           reg,
                          const CellRange&        cells,
-                         const BlackoilPropertiesInterface& props,
+                         BlackoilPropertiesInterface& props,
+                         const std::vector<double> swat_init,
                          std::vector< std::vector<double> >& phase_pressures)
         {
             const double z0   = reg.datum();
@@ -610,8 +611,14 @@ namespace Opm
                 double sw = 0.0;
                 if (water) {
                     const double pcov = phase_pressures[oilpos][local_index] - phase_pressures[waterpos][local_index];
+                    if (swat_init.empty()) { // Invert Pc to find sw
                         sw = satFromPc(props, waterpos, cell, pcov);
                         phase_saturations[waterpos][local_index] = sw;
+                    } else { // Scale Pc to reflect imposed sw
+                        sw = swat_init[cell];
+                        props.swatInitScaling(cell, pcov, sw);
+                        phase_saturations[waterpos][local_index] = sw;
+                    }
                 }
                 double sg = 0.0;
                 if (gas) {
@@ -736,7 +743,7 @@ namespace Opm
      * \param[in] gravity  Acceleration of gravity, assumed to be in Z direction.
      */
     void initStateEquil(const UnstructuredGrid& grid,
-                        const BlackoilPropertiesInterface& props,
+                        BlackoilPropertiesInterface& props,
                         const Opm::DeckConstPtr deck,
                         const double gravity,
                         BlackoilState& state)