Multiple changes dealing with initialization. Work in progress.

- Moved simulator state class into its own file.
 - Using new initState...() methods in spu_2p.cpp
 - No longer controlled by 'scenario' parameter.

opm/core/TwophaseState.hpp

@@ -0,0 +1,85 @@
+/*
+  Copyright 2012 SINTEF ICT, Applied Mathematics.
+
+  This file is part of the Open Porous Media project (OPM).
+
+  OPM is free software: you can redistribute it and/or modify
+  it under the terms of the GNU General Public License as published by
+  the Free Software Foundation, either version 3 of the License, or
+  (at your option) any later version.
+
+  OPM is distributed in the hope that it will be useful,
+  but WITHOUT ANY WARRANTY; without even the implied warranty of
+  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+  GNU General Public License for more details.
+
+  You should have received a copy of the GNU General Public License
+  along with OPM.  If not, see <http://www.gnu.org/licenses/>.
+*/
+
+#ifndef OPM_TWOPHASESTATE_HEADER_INCLUDED
+#define OPM_TWOPHASESTATE_HEADER_INCLUDED
+
+#include <opm/core/grid.h>
+#include <opm/core/fluid/IncompPropertiesInterface.hpp>
+#include <vector>
+
+namespace Opm
+{
+
+    class TwophaseState
+    {
+    public:
+
+        void init(const UnstructuredGrid& g)
+        {
+            press_.resize(g.number_of_cells, 0.0);
+            fpress_.resize(g.number_of_faces, 0.0);
+            flux_.resize(g.number_of_faces, 0.0);
+            sat_.resize(2 * g.number_of_cells, 0.0);
+        }
+
+        enum ExtremalSat { MinSat, MaxSat };
+
+        void setWaterSat(const std::vector<int>& cells,
+                         const Opm::IncompPropertiesInterface& props,
+                         ExtremalSat es)
+        {
+            const int n = cells.size();
+            std::vector<double> smin(2*n);
+            std::vector<double> smax(2*n);
+            props.satRange(n, &cells[0], &smin[0], &smax[0]);
+            const double* svals = (es == MinSat) ? &smin[0] : &smax[0];
+            for (int ci = 0; ci < n; ++ci) {
+                const int cell = cells[ci];
+                sat_[2*cell] = svals[2*ci];
+                sat_[2*cell + 1] = 1.0 - sat_[2*cell];
+            }
+        }
+
+        int numPhases() const
+        {
+            return 2;
+        }
+
+        std::vector<double>& pressure    () { return press_ ; }
+        std::vector<double>& facepressure() { return fpress_; }
+        std::vector<double>& faceflux    () { return flux_  ; }
+        std::vector<double>& saturation  () { return sat_   ; }
+
+        const std::vector<double>& pressure    () const { return press_ ; }
+        const std::vector<double>& facepressure() const { return fpress_; }
+        const std::vector<double>& faceflux    () const { return flux_  ; }
+        const std::vector<double>& saturation  () const { return sat_   ; }
+
+    private:
+        std::vector<double> press_ ;
+        std::vector<double> fpress_;
+        std::vector<double> flux_  ;
+        std::vector<double> sat_   ;
+    };
+
+} // namespace Opm
+
+
+#endif // OPM_TWOPHASESTATE_HEADER_INCLUDED

opm/core/utility/initState_impl.hpp

@@ -159,6 +159,7 @@ namespace Opm
         if (state.numPhases() != 2) {
             THROW("initStateTwophaseFromDeck(): state must have two phases.");
         }
+        state.init(grid);
         const int num_cells = props.numCells();
         const bool convection_testcase = param.getDefault("convection_testcase", false);
         const bool segregation_testcase = param.getDefault("segregation_testcase", false);
@@ -179,6 +180,13 @@ namespace Opm
             const double init_p = param.getDefault("ref_pressure", 100)*unit::barsa;
             std::fill(state.pressure().begin(), state.pressure().end(), init_p);
         } else if (segregation_testcase) {
+            // Warn against error-prone usage.
+            if (gravity == 0.0) {
+                std::cout << "**** Warning: running gravity segregation scenario, but gravity is zero." << std::endl;
+            }
+            if (grid.cartdims[2] <= 1) {
+                std::cout << "**** Warning: running gravity segregation scenario, which expects nz > 1." << std::endl;
+            }
             // Initialise water saturation to max *above* water-oil contact.
             const double woc = param.get<double>("water_oil_contact");
             initWaterOilContact(grid, props, woc, WaterAbove, state);
@@ -187,6 +195,13 @@ namespace Opm
             double dens[2] = { props.density()[1], props.density()[0] };
             initHydrostaticPressure(grid, dens, woc, gravity, woc, ref_p, state);
         } else if (param.has("water_oil_contact")) {
+            // Warn against error-prone usage.
+            if (gravity == 0.0) {
+                std::cout << "**** Warning: running gravity convection scenario, but gravity is zero." << std::endl;
+            }
+            if (grid.cartdims[2] <= 1) {
+                std::cout << "**** Warning: running gravity convection scenario, which expects nz > 1." << std::endl;
+            }
             // Initialise water saturation to max below water-oil contact.
             const double woc = param.get<double>("water_oil_contact");
             initWaterOilContact(grid, props, woc, WaterBelow, state);
@@ -240,6 +255,7 @@ namespace Opm
         if (state.numPhases() != 2) {
             THROW("initStateTwophaseFromDeck(): state must have two phases.");
         }
+        state.init(grid);
         if (deck.hasField("EQUIL")) {
             // Set saturations depending on oil-water contact.
             const EQUIL& equil= deck.getEQUIL();