Added initBlackoilSurfvol() function.

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);

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
 

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