Merge pull request #1015 from totto82/hydrocarbonState

Add member hydroCarbonState to the BlackoilState

CMakeLists_files.cmake

@@ -377,6 +377,7 @@ list (APPEND PUBLIC_HEADER_FILES
         opm/core/utility/Event.hpp
         opm/core/utility/Event_impl.hpp
         opm/core/utility/Factory.hpp
+        opm/core/utility/initHydroCarbonState.hpp
         opm/core/utility/MonotCubicInterpolator.hpp
         opm/core/utility/NonuniformTableLinear.hpp
         opm/core/utility/NullStream.hpp

opm/core/simulator/BlackoilState.cpp

@@ -21,15 +21,18 @@ BlackoilState::BlackoilState( size_t num_cells , size_t num_faces , size_t num_p
 }
 
 BlackoilState::BlackoilState( const BlackoilState& other )
-    : SimulationDataContainer(other)
+    : SimulationDataContainer(other),
+      hydrocarbonstate_(other.hydroCarbonState())
 {
     setBlackoilStateReferencePointers();
+
 }
 
 BlackoilState& BlackoilState::operator=( const BlackoilState& other )
 {
     SimulationDataContainer::operator=(other);
     setBlackoilStateReferencePointers();
+    hydrocarbonstate_ = other.hydroCarbonState();
     return *this;
 }
 

opm/core/simulator/BlackoilState.hpp

@@ -30,6 +30,12 @@
 namespace Opm
 {
 
+    enum HydroCarbonState {
+        GasOnly = 0,
+        GasAndOil = 1,
+        OilOnly = 2
+    };
+
     /// Simulator state for a blackoil simulator.
     class BlackoilState : public SimulationDataContainer
     {
@@ -62,10 +68,12 @@ namespace Opm
         std::vector<double>& surfacevol  () { return *surfacevol_ref_;  }
         std::vector<double>& gasoilratio () { return *gasoilratio_ref_; }
         std::vector<double>& rv ()          { return *rv_ref_;          }
+        std::vector<HydroCarbonState>& hydroCarbonState() { return hydrocarbonstate_;  }
 
         const std::vector<double>& surfacevol  () const { return *surfacevol_ref_;  }
         const std::vector<double>& gasoilratio () const { return *gasoilratio_ref_; }
         const std::vector<double>& rv ()          const { return *rv_ref_;          }
+        const std::vector<HydroCarbonState>& hydroCarbonState() const { return hydrocarbonstate_;  }
 
     private:
         void setBlackoilStateReferencePointers();
@@ -73,6 +81,9 @@ namespace Opm
         std::vector<double>* gasoilratio_ref_;
         std::vector<double>* rv_ref_;
 
+        // A vector storing the hydro carbon state.
+        std::vector<HydroCarbonState> hydrocarbonstate_;
+
 
     };
 } // namespace Opm

opm/core/utility/initHydroCarbonState.hpp

@@ -0,0 +1,44 @@
+#ifndef INITHYDROCARBONSTATE_HPP
+#define INITHYDROCARBONSTATE_HPP
+
+#include "opm/core/simulator/BlackoilState.hpp"
+
+namespace Opm
+{
+
+void initHydroCarbonState(BlackoilState& state, const PhaseUsage& pu, const int num_cells) {
+    enum { Oil = BlackoilPhases::Liquid, Gas = BlackoilPhases::Vapour, Water = BlackoilPhases::Aqua };
+    // hydrocarbonstate is only used when gas and oil is present
+    assert(pu.phase_used[Oil]);
+    std::vector<HydroCarbonState>& hydroCarbonState = state.hydroCarbonState();
+    hydroCarbonState.resize(num_cells);
+    if (!pu.phase_used[Gas]) {
+        // hydroCarbonState should only be used when oil and gas is present. Return OilOnly to avoid potential trouble.
+        std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::OilOnly);
+        return;
+    }
+    const int np = pu.num_phases;
+    std::fill(hydroCarbonState.begin(), hydroCarbonState.end(), HydroCarbonState::GasAndOil);
+
+    // set hydrocarbon state
+    const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
+    const std::vector<double>& saturation = state.saturation();
+    for (int c = 0; c < num_cells; ++c) {
+        if (pu.phase_used[Water]) {
+            if ( saturation[c*np + pu.phase_pos[ Water ]] > (1.0 - epsilon)) {
+                continue; // cases (almost) filled with water is treated as GasAndOil case;
+            }
+        }
+        if ( saturation[c*np + pu.phase_pos[ Gas ]] == 0.0) {
+            hydroCarbonState[c] = HydroCarbonState::OilOnly;
+            continue;
+        }
+        if ( saturation[c*np + pu.phase_pos[ Oil ]] == 0.0) {
+            hydroCarbonState[c] = HydroCarbonState::GasOnly;
+        }
+    }
+}
+
+
+} // namespace Opm
+#endif // INITHYDROCARBONSTATE_HPP