From 5deeed1dac302fe450c273c26308cb40d85a5941 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= <atgeirr@sintef.no>
Date: Mon, 21 May 2012 14:03:56 +0200
Subject: [PATCH] Added initBlackoilSurfvol() function.

---
 examples/sim_wateroil.cpp           |  1 +
 opm/core/utility/initState.hpp      | 15 ++++++++++++
 opm/core/utility/initState_impl.hpp | 37 +++++++++++++++++++++++++++++
 3 files changed, 53 insertions(+)

diff --git a/examples/sim_wateroil.cpp b/examples/sim_wateroil.cpp
index d4417175..1e8df834 100644
--- a/examples/sim_wateroil.cpp
+++ b/examples/sim_wateroil.cpp
@@ -201,6 +201,7 @@ main(int argc, char** argv)
         } else {
             initStateFromDeck(*grid->c_grid(), *props, deck, gravity[2], state);
         }
+        initBlackoilSurfvol(*grid->c_grid(), *props, state);
     } else {
         // Grid init.
         const int nx = param.getDefault("nx", 100);
diff --git a/opm/core/utility/initState.hpp b/opm/core/utility/initState.hpp
index e49a692c..2d22769e 100644
--- a/opm/core/utility/initState.hpp
+++ b/opm/core/utility/initState.hpp
@@ -96,6 +96,21 @@ namespace Opm
                            const double gravity,
                            State& state);
 
+    /// Initialize a two-phase water-oil blackoil state from input deck.
+    /// If EQUIL is present:
+    ///   - saturation is set according to the water-oil contact,
+    ///   - pressure is set to hydrostatic equilibrium.
+    /// Otherwise:
+    ///   - saturation is set according to SWAT,
+    ///   - pressure is set according to PRESSURE.
+    /// In addition, this function sets surfacevol.
+    template <class Props, class State>
+    void initBlackoilStateFromDeck(const UnstructuredGrid& grid,
+                                   const Props& props,
+                                   const EclipseGridParser& deck,
+                                   const double gravity,
+                                   State& state);
+
 
 } // namespace Opm
 
diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp
index 389d5d40..f49db2a5 100644
--- a/opm/core/utility/initState_impl.hpp
+++ b/opm/core/utility/initState_impl.hpp
@@ -550,6 +550,43 @@ namespace Opm
     }
 
 
+
+
+
+    /// Initialize surface volume from pressure and saturation by z = As.
+    template <class Props, class State>
+    void initBlackoilSurfvol(const UnstructuredGrid& grid,
+                             const Props& props,
+                             State& state)
+    {
+        state.surfacevol() = state.saturation();
+        const int np = props.numPhases();
+        const int nc = grid.number_of_cells;
+        std::vector<double> allA(nc*np*np);
+        std::vector<int> allcells(nc);
+        for (int c = 0; c < nc; ++c) {
+            allcells[c] = c;
+        }
+        // Assuming that using the saturation as z argument here does not change
+        // the outcome. This is not guaranteed unless we have only a single phase
+        // per cell.
+        props.matrix(nc, &state.pressure()[0], &state.surfacevol()[0], &allcells[0], &allA[0], 0);
+        for (int c = 0; c < nc; ++c) {
+            // Using z = As
+            double* z = &state.surfacevol()[c*np];
+            const double* s = &state.saturation()[c*np];
+            const double* A = &allA[c*np*np];
+            for (int row = 0; row < np; ++row) {
+                z[row] = 0.0;
+                for (int col = 0; col < np; ++col) {
+                    // Recall: A has column-major ordering.
+                    z[row] += A[row + np*col]*s[col];
+                }
+            }
+        }
+    }
+
+
 } // namespace Opm