Added free function initSaturation().

The state argument is of type SimulatorState& - and no longer a template
parameter.

opm/core/simulator/initState.hpp

@@ -22,6 +22,7 @@
 
 #include <opm/parser/eclipse/Deck/Deck.hpp>
 
+#include <opm/core/simulator/SimulatorState.hpp>
 struct UnstructuredGrid;
 
 namespace Opm
@@ -35,6 +36,18 @@ namespace Opm
     ///
     /// Functions for initializing a reservoir state.
 
+    /// Will initialize the first and second component of the
+    /// SATURATION field in all the cells in the set @cells. The
+    /// @props object will be queried, and depending on the value
+    /// @satType either the minimum or the maximum saturation is
+    /// applied to thee first component in the SATURATION field.
+    /// For the second component (1 - first_sat) is used.
+
+    enum ExtremalSat { MinSat, MaxSat };
+    template <class Props>
+    static void initSaturation(const std::vector<int>& cells , const Props& props , SimulatorState& state , ExtremalSat satType);
+
+
     /// Initialize a two-phase state from parameters.
     /// The following parameters are accepted (defaults):
     ///    - convection_testcase   (false) -- Water in the 'left' part of the grid.

opm/core/simulator/initState_impl.hpp

@@ -40,6 +40,38 @@
 namespace Opm
 {
 
+
+    template <class Props>
+    static void initSaturation(const std::vector<int>& cells , const Props& props , SimulatorState& state , ExtremalSat satType) {
+        const int num_phases = state.numPhases();
+        std::vector<double> min_sat(num_phases * cells.size());
+        std::vector<double> max_sat(num_phases * cells.size());
+        props.satRange(cells.size() ,cells.data() , min_sat.data() , max_sat.data());
+
+        {
+            std::vector<double> second_sat(cells.size());
+            std::vector<double> first_sat(cells.size());
+
+            for (size_t index = 0; index < cells.size(); index++) {
+                if (satType == MinSat) {
+                    first_sat[index] = min_sat[ num_phases * index];
+                    second_sat[index] = 1 - min_sat[ num_phases * index];
+                } else {
+                    first_sat[index] = max_sat[ num_phases * index];
+                    second_sat[index] = 1 - max_sat[ num_phases * index];
+                }
+            }
+
+            state.setCellDataComponent( "SATURATION" , 0 , cells , first_sat );
+            state.setCellDataComponent( "SATURATION" , 1 , cells , second_sat );
+        }
+    }
+
+
+        // Initialize saturations so that there is water below woc,
+        // and oil above.
+
+
     namespace
     {
 #ifdef __clang__
@@ -75,8 +107,8 @@ namespace Opm
 #pragma clang diagnostic pop
 #endif /* __clang__ */
 
+
         enum WaterInit { WaterBelow, WaterAbove };
-        enum ExtremalSat { MinSat, MaxSat };
 
         /// Will initialize the first and second component of the
         /// SATURATION field in all the cells in the set @cells. The
@@ -87,21 +119,26 @@ namespace Opm
 
         template <class Props, class State>
         static void initSaturation(const std::vector<int>& cells , const Props& props , State& state , ExtremalSat satType) {
-	    std::vector<double> min_sat(cells.size());
+	    std::vector<double> min_sat(state.numPhases() * cells.size());
-	    std::vector<double> max_sat(cells.size());
+	    std::vector<double> max_sat(state.numPhases() * cells.size());
-	    std::vector<double> second_sat(cells.size());
+            props.satRange(cells.size() ,cells.data() , min_sat.data() , max_sat.data());
-	    std::vector<double>* first_sat;
 
-	    props.satRange(cells.size() ,cells.data() , min_sat.data() , max_sat.data());
+            {
-	    if (satType == MinSat) {
+                std::vector<double> first_sat(cells.size());
-	      first_sat = &min_sat;
+                std::vector<double> second_sat(cells.size());
-	    } else {
-	      first_sat = &max_sat;
-	    }
 
-	    std::transform( first_sat->begin() , first_sat->end() , second_sat.begin() , [](double s) { return 1 - s; });
+                for (size_t index=0; index < cells.size(); index++) {
-	    state.setCellDataComponent( "SATURATION" , 0 , cells , *first_sat );
+                    if (satType == MinSat) {
-	    state.setCellDataComponent( "SATURATION" , 1 , cells , second_sat );
+                        first_sat[index]  =     min_sat[index * state.numPhases()];
+                        second_sat[index] = 1 - min_sat[index * state.numPhases()];
+                    } else {
+                        first_sat[index]  =     max_sat[index * state.numPhases()];
+                        second_sat[index] = 1 - max_sat[index * state.numPhases()];
+                    }
+                }
+                state.setCellDataComponent( "SATURATION" , 0 , cells , first_sat );
+                state.setCellDataComponent( "SATURATION" , 1 , cells , second_sat );
+            }
 	}
 
 

tutorials/tutorial3.cpp

@@ -37,6 +37,7 @@
 
 #include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
 
+#include <opm/core/simulator/initState.hpp>
 #include <opm/core/simulator/TwophaseState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 
@@ -267,18 +268,7 @@ try
     /// \internal [two-phase state]
     TwophaseState state;
     state.init(grid.number_of_cells , grid.number_of_faces, 2);
-    {
+    initSaturation( allcells , props , state , MinSat );
-        std::vector<double> min_sat(allcells.size());
-        std::vector<double> max_sat(allcells.size());
-        std::vector<double> second_sat(allcells.size());
-
-        props.satRange(allcells.size() ,allcells.data() , min_sat.data() , max_sat.data());
-
-        std::transform( min_sat.begin() , min_sat.end() , second_sat.begin() , [](double s) { return 1 - s; });
-        state.setCellDataComponent( "SATURATION" , 0 , allcells , min_sat );
-        state.setCellDataComponent( "SATURATION" , 1 , allcells , second_sat );
-
-    }
 
     /// \internal [two-phase state]
     /// \endinternal

tutorials/tutorial4.cpp

@@ -37,6 +37,7 @@
 
 #include <opm/core/transport/reorder/TransportSolverTwophaseReorder.hpp>
 
+#include <opm/core/simulator/initState.hpp>
 #include <opm/core/simulator/TwophaseState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 
@@ -214,18 +215,7 @@ try
     /// \internal[two-phase state]
     TwophaseState state;
     state.init(grid.number_of_cells , grid.number_of_faces, 2);
-    {
+    initSaturation( allcells , props , state , MinSat );
-        std::vector<double> min_sat(allcells.size());
-        std::vector<double> max_sat(allcells.size());
-        std::vector<double> second_sat(allcells.size());
-
-        props.satRange(allcells.size() ,allcells.data() , min_sat.data() , max_sat.data());
-
-        std::transform( min_sat.begin() , min_sat.end() , second_sat.begin() , [](double s) { return 1 - s; });
-        state.setCellDataComponent( "SATURATION" , 0 , allcells , min_sat );
-        state.setCellDataComponent( "SATURATION" , 1 , allcells , second_sat );
-
-    }
 
 
     /// \internal[two-phase state]