From cb4b698b515f64e2e588883153e5825e4eb5f66c Mon Sep 17 00:00:00 2001
From: Tor Harald Sandve <tor.harald.sandve@iris.no>
Date: Fri, 13 May 2016 09:10:13 +0200
Subject: [PATCH] Add method for calculating the initial hydroCarbonState

The following hydroCarbonState are used
    enum HydroCarbonState {
        GasOnly = 0,
        GasAndOil = 1,
        OilOnly = 2
    };
Cells almost filled with water are treated as a GasAndOil state
---
 CMakeLists_files.cmake                    |  1 +
 opm/core/simulator/BlackoilState.hpp      |  6 ++++
 opm/core/utility/initHydroCarbonState.hpp | 42 +++++++++++++++++++++++
 3 files changed, 49 insertions(+)
 create mode 100644 opm/core/utility/initHydroCarbonState.hpp

diff --git a/CMakeLists_files.cmake b/CMakeLists_files.cmake
index a5ee0621..d36438cb 100644
--- a/CMakeLists_files.cmake
+++ b/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
diff --git a/opm/core/simulator/BlackoilState.hpp b/opm/core/simulator/BlackoilState.hpp
index 51c6c21c..905b0943 100644
--- a/opm/core/simulator/BlackoilState.hpp
+++ b/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
     {
diff --git a/opm/core/utility/initHydroCarbonState.hpp b/opm/core/utility/initHydroCarbonState.hpp
new file mode 100644
index 00000000..c0be0f68
--- /dev/null
+++ b/opm/core/utility/initHydroCarbonState.hpp
@@ -0,0 +1,42 @@
+#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]);
+    if (!pu.phase_used[Gas]) {
+        return; // do nothing
+    }
+    std::vector<int>& hydroCarbonState = state.hydroCarbonState();
+    const int np = pu.num_phases;
+    hydroCarbonState.resize(num_cells);
+    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