From 88d1d9762dfca836c688f98f1c2e95f022e2a03c Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Halvor=20M=C3=B8ll=20Nilsen?= <Halvor.M.Nilsen@sintef.no>
Date: Mon, 27 Aug 2012 12:20:03 +0200
Subject: [PATCH] Addes support for initialisation of three phases using SWOF
 and SGOF

---
 opm/core/utility/initState_impl.hpp | 28 ++++++++++++++++++++--------
 1 file changed, 20 insertions(+), 8 deletions(-)

diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp
index c50db8a9..f6670434 100644
--- a/opm/core/utility/initState_impl.hpp
+++ b/opm/core/utility/initState_impl.hpp
@@ -512,11 +512,12 @@ namespace Opm
                            State& state)
     {
         const int num_phases = props.numPhases();
-        if (num_phases != 2) {
-            THROW("initStateFromDeck(): currently handling only two-phase scenarios.");
-        }
+        
         state.init(grid, num_phases);
         if (deck.hasField("EQUIL")) {
+            if (num_phases != 2) {
+                THROW("initStateFromDeck(): currently handling only two-phase scenarios.");
+            }
             // Set saturations depending on oil-water contact.
             const EQUIL& equil= deck.getEQUIL();
             if (equil.equil.size() != 1) {
@@ -535,11 +536,22 @@ namespace Opm
             const std::vector<double>& sw_deck = deck.getFloatingPointValue("SWAT");
             const std::vector<double>& p_deck = deck.getFloatingPointValue("PRESSURE");
             const int num_cells = grid.number_of_cells;
-            for (int c = 0; c < num_cells; ++c) {
-                int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
-                s[2*c] = sw_deck[c_deck];
-                s[2*c + 1] = 1.0 - s[2*c];
-                p[c] = p_deck[c_deck];
+            if(num_phases == 2){
+                for (int c = 0; c < num_cells; ++c) {
+                    int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                    s[2*c] = sw_deck[c_deck];
+                    s[2*c + 1] = 1.0 - s[2*c];
+                    p[c] = p_deck[c_deck];
+                }
+            }else{
+                const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
+                 for (int c = 0; c < num_cells; ++c) {
+                    int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
+                    s[2*c] = sw_deck[c_deck];
+                    s[2*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]);
+                    s[2*c + 2] = sg_deck[c_deck];
+                    p[c] = p_deck[c_deck];
+                 }
             }
         } else {
             THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE.");