From 913054c47369279b1d074d71f08757fa48ee9536 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Mon, 3 Sep 2012 13:54:50 +0200
Subject: [PATCH] Added more checks in 3-phase init code.

---
 opm/core/utility/initState_impl.hpp | 14 +++++++++-----
 1 file changed, 9 insertions(+), 5 deletions(-)

diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp
index f66704348..9a635e93f 100644
--- a/opm/core/utility/initState_impl.hpp
+++ b/opm/core/utility/initState_impl.hpp
@@ -512,11 +512,10 @@ namespace Opm
                            State& state)
     {
         const int num_phases = props.numPhases();
-        
         state.init(grid, num_phases);
         if (deck.hasField("EQUIL")) {
             if (num_phases != 2) {
-                THROW("initStateFromDeck(): currently handling only two-phase scenarios.");
+                THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios.");
             }
             // Set saturations depending on oil-water contact.
             const EQUIL& equil= deck.getEQUIL();
@@ -536,22 +535,27 @@ 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;
-            if(num_phases == 2){
+            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{
+            } else if (num_phases == 3) {
+                if (!deck.hasField("SGAS")) {
+                    THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found).");
+                }
                 const std::vector<double>& sg_deck = deck.getFloatingPointValue("SGAS");
-                 for (int c = 0; c < num_cells; ++c) {
+                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(): init with SWAT etc. only available with 2 or 3 phases.");
             }
         } else {
             THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE.");