Merge pull request #1087 from babrodtk/somax_in_restart

Added somax to input/output to restart file

opm/autodiff/BlackoilModelBase.hpp

@@ -103,6 +103,8 @@ namespace Opm {
             ADB rsSat; // Saturated gas-oil ratio
             ADB rvSat; // Saturated oil-gas ratio
 
+            std::vector<double> soMax; // Maximum oil saturation
+
             std::array<V, fipValues> fip;
         };
 

opm/autodiff/BlackoilModelBase_impl.hpp

@@ -463,6 +463,7 @@ typedef Eigen::Array<double,
         : rq(num_phases)
         , rsSat(ADB::null())
         , rvSat(ADB::null())
+        , soMax()
         , fip()
     {
     }
@@ -662,6 +663,7 @@ typedef Eigen::Array<double,
                     } else {
                         state.rv = sd_.rvSat;
                     }
+                    sd_.soMax = fluid_.satOilMax();
                 }
             }
             else {

opm/autodiff/BlackoilModelEbos.hpp

@@ -1247,6 +1247,9 @@ namespace Opm {
             VectorType& Pb = simData.getCellData( "PBUB" );
             VectorType& Pd = simData.getCellData( "PDEW" );
 
+            simData.registerCellData( "SOMAX", 1 );
+            VectorType& somax = simData.getCellData( "SOMAX" );
+
             std::vector<int> failed_cells_pb;
             std::vector<int> failed_cells_pd;
             const auto& gridView = ebosSimulator().gridView();
@@ -1271,6 +1274,8 @@ namespace Opm {
 
                 temperature[cellIdx] = fs.temperature(FluidSystem::oilPhaseIdx).value();
 
+                somax[cellIdx] = ebosSimulator().model().maxOilSaturation(cellIdx);
+
                 if (aqua_active) {
                     saturation[ satIdx + aqua_pos ] = fs.saturation(FluidSystem::waterPhaseIdx).value();
                     bWater[cellIdx] = fs.invB(FluidSystem::waterPhaseIdx).value();

opm/autodiff/BlackoilPropsAdFromDeck.cpp

@@ -168,10 +168,10 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
         if (deck.hasKeyword("VAPPARS") && deck.hasKeyword("VAPOIL") && deck.hasKeyword("DISGAS")) {
             vap1_ = deck.getKeyword("VAPPARS").getRecord(0).getItem(0).get< double >(0);
             vap2_ = deck.getKeyword("VAPPARS").getRecord(0).getItem(1).get< double >(0);
-            satOilMax_.resize(number_of_cells, 0.0);
         } else if (deck.hasKeyword("VAPPARS")) {
             OPM_THROW(std::runtime_error, "Input has VAPPARS, but missing VAPOIL and/or DISGAS\n");
         }
+        satOilMax_.resize(number_of_cells, 0.0);
 
         SaturationPropsFromDeck* ptr
             = new SaturationPropsFromDeck();
@@ -879,19 +879,22 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
     /// Update for max oil saturation.
     void BlackoilPropsAdFromDeck::updateSatOilMax(const std::vector<double>& saturation)
     {
-        if (!satOilMax_.empty()) {
+        const int n = satOilMax_.size();
-            const int n = satOilMax_.size();
+        const int np = phase_usage_.num_phases;
-            const int np = phase_usage_.num_phases;
+        const int posOil = phase_usage_.phase_pos[Oil];
-            const int posOil = phase_usage_.phase_pos[Oil];
+        const double* s = saturation.data();
-            const double* s = saturation.data();
+        for (int i=0; i<n; ++i) {
-            for (int i=0; i<n; ++i) {
+            if (satOilMax_[i] < s[np*i+posOil]) {
-                if (satOilMax_[i] < s[np*i+posOil]) {
+                satOilMax_[i] = s[np*i+posOil];
-                    satOilMax_[i] = s[np*i+posOil];
-                }
             }
         }
     }
 
+    const std::vector<double>& BlackoilPropsAdFromDeck::satOilMax() const
+    {
+        return satOilMax_;
+    }
+
     /// 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.
@@ -918,7 +921,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
                                            const std::vector<int>& cells,
                                            const double vap) const
     {
-        if (!satOilMax_.empty() && vap > 0.0) {
+        if (vap > 0.0) {
             const int n = cells.size();
             V factor = V::Ones(n, 1);
             const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());
@@ -943,7 +946,7 @@ BlackoilPropsAdFromDeck::BlackoilPropsAdFromDeck(const BlackoilPropsAdFromDeck&
                                            const std::vector<int>& cells,
                                            const double vap) const
     {
-        if (!satOilMax_.empty() && vap > 0.0) {
+        if (vap > 0.0) {
             const int n = cells.size();
             V factor = V::Ones(n, 1);
             const double eps_sqrt = std::sqrt(std::numeric_limits<double>::epsilon());

opm/autodiff/BlackoilPropsAdFromDeck.hpp

@@ -362,6 +362,9 @@ namespace Opm
         /// Update for max oil saturation.                  
         void updateSatOilMax(const std::vector<double>& saturation);
 
+        /// Returns the max oil saturation
+        const std::vector<double>& satOilMax() const;
+
         /// 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/Compat.cpp

@@ -174,6 +174,10 @@ void solutionToSim( const data::Solution& sol,
         state.getCellData("SSOL") = sol.data("SSOL");
     }
 
+    if ( sol.has( "SOMAX" ) ) {
+        state.registerCellData("SOMAX", 1);
+        state.getCellData("SOMAX") = sol.data("SOMAX");
+    }
 }
 
 

opm/autodiff/SimulatorFullyImplicitBlackoilEbos.hpp

@@ -149,6 +149,17 @@ public:
             // 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_);
+            {
+                const int num_cells = Opm::UgGridHelpers::numCells(grid());
+
+                typedef std::vector<double> VectorType;
+
+                const VectorType& somax = state.getCellData( "SOMAX" );
+
+                for (int cellIdx = 0; cellIdx < num_cells; ++cellIdx) {
+                    ebosSimulator_.model().setMaxOilSaturation(somax[cellIdx], cellIdx);
+                }
+            }
         }
 
         // Create timers and file for writing timing info.

opm/autodiff/SimulatorFullyImplicitBlackoilOutput.hpp

@@ -420,11 +420,12 @@ namespace Opm
                          ExtraData& extra )
     {
         std::map<std::string, UnitSystem::measure> solution_keys {{"PRESSURE" , UnitSystem::measure::pressure},
-                                                                  {"SWAT"     , UnitSystem::measure::identity},
+                                                                  {"SWAT" , UnitSystem::measure::identity},
-                                                                  {"SGAS"     , UnitSystem::measure::identity},
+                                                                  {"SGAS" , UnitSystem::measure::identity},
-                                                                  {"TEMP"     , UnitSystem::measure::temperature},
+                                                                  {"TEMP" , UnitSystem::measure::temperature},
-                                                                  {"RS"       , UnitSystem::measure::gas_oil_ratio},
+                                                                  {"RS" , UnitSystem::measure::gas_oil_ratio},
-                                                                  {"RV"       , UnitSystem::measure::oil_gas_ratio}};
+                                                                  {"RV" , UnitSystem::measure::oil_gas_ratio},
+                                                                  {"SOMAX", UnitSystem::measure::identity}};
         std::map<std::string, bool> extra_keys {{"OPMEXTRA" , false}};
 
         if (restart_double_si_) {
@@ -550,6 +551,8 @@ namespace Opm
             addToSimData( simData, "RSSAT", sd.rsSat );
             addToSimData( simData, "RVSAT", sd.rvSat );
 
+            addToSimData( simData, "SOMAX", sd.soMax );
+
             return simData;
         }
 
@@ -774,6 +777,13 @@ namespace Opm
                         data::TargetType::RESTART_AUXILIARY);
             }
 
+            if (sd.hasCellData("SOMAX")) {
+                output.insert("SOMAX",
+                        Opm::UnitSystem::measure::identity,
+                        std::move( sd.getCellData("SOMAX") ),
+                        data::TargetType::RESTART_AUXILIARY);
+            }
+
             //Warn for any unhandled keyword
             if (log) {
                 for (auto& keyValue : rstKeywords) {

opm/polymer/fullyimplicit/FullyImplicitCompressiblePolymerSolver.hpp

@@ -82,6 +82,7 @@ namespace Opm {
                 : rq(num_phases)
                 , rsSat(ADB::null())
                 , rvSat(ADB::null())
+                , soMax()
                 , fip()
             {
             }
@@ -92,6 +93,7 @@ namespace Opm {
             std::vector<ReservoirResidualQuant> rq;
             ADB rsSat;
             ADB rvSat;
+            std::vector<double> soMax;
             std::array<V, fipValues> fip;
         };