Extract common parts of blackoil and incomp. state

Put the identical parts of the simulator state into a base class that
they can be referenced from when adressing the common fields.

CMakeLists_files.cmake

@@ -95,9 +95,11 @@ list (APPEND MAIN_SOURCE_FILES
 	opm/core/props/satfunc/SatFuncStone2.cpp
 	opm/core/props/satfunc/SaturationPropsBasic.cpp
 	opm/core/props/satfunc/SaturationPropsFromDeck.cpp
+	opm/core/simulator/BlackoilState.cpp
 	opm/core/simulator/SimulatorCompressibleTwophase.cpp
 	opm/core/simulator/SimulatorIncompTwophase.cpp
 	opm/core/simulator/SimulatorReport.cpp
+	opm/core/simulator/SimulatorState.cpp
 	opm/core/simulator/SimulatorTimer.cpp
 	opm/core/tof/DGBasis.cpp
 	opm/core/tof/TofReorder.cpp
@@ -304,8 +306,10 @@ list (APPEND PUBLIC_HEADER_FILES
 	opm/core/simulator/SimulatorCompressibleTwophase.hpp
 	opm/core/simulator/SimulatorIncompTwophase.hpp
 	opm/core/simulator/SimulatorReport.hpp
+	opm/core/simulator/SimulatorState.hpp
 	opm/core/simulator/SimulatorTimer.hpp
 	opm/core/simulator/TwophaseState.hpp
+	opm/core/simulator/TwophaseState_impl.hpp
 	opm/core/simulator/WellState.hpp
 	opm/core/simulator/initState.hpp
 	opm/core/simulator/initState_impl.hpp

opm/core/simulator/BlackoilState.cpp

@@ -0,0 +1,37 @@
+#include "BlackoilState.hpp"
+#include <opm/core/props/BlackoilPropertiesInterface.hpp>
+
+using namespace Opm;
+
+void
+BlackoilState::init(const UnstructuredGrid& g, int num_phases) {
+   SimulatorState::init(g, num_phases);
+   gor_.resize(g.number_of_cells, 0.) ;
+   // surfvol_ intentionally empty, left to initBlackoilSurfvol
+}
+
+/// Set the first saturation to either its min or max value in
+/// the indicated cells. The second saturation value s2 is set
+/// to (1.0 - s1) for each cell. Any further saturation values
+/// are unchanged.
+void
+BlackoilState::setFirstSat(const std::vector<int>& cells,
+                           const Opm::BlackoilPropertiesInterface& props,
+                           ExtremalSat es) {
+    SimulatorState::setFirstSat(cells, props, es);
+}
+
+bool
+BlackoilState::equals(const SimulatorState& other,
+                      double epsilon) const {
+    const BlackoilState* that = dynamic_cast <const BlackoilState*> (&other);
+    bool equal = that != 0;
+    equal = equal && SimulatorState::equals (other, epsilon);
+    equal = equal && SimulatorState::vectorApproxEqual(this->surfacevol(),
+                                                       that->surfacevol(),
+                                                       epsilon);
+    equal = equal && SimulatorState::vectorApproxEqual(this->gasoilratio(),
+                                                       that->gasoilratio(),
+                                                       epsilon);
+    return equal;
+}

opm/core/simulator/BlackoilState.hpp

@@ -22,37 +22,17 @@
 
 #include <opm/core/grid.h>
 #include <opm/core/props/BlackoilPropertiesInterface.hpp>
-#include <opm/core/utility/ErrorMacros.hpp>
+#include <opm/core/simulator/SimulatorState.hpp>
 #include <vector>
-#include <cmath>
 
 namespace Opm
 {
 
     /// Simulator state for a blackoil simulator.
-    class BlackoilState
+    class BlackoilState : public SimulatorState
     {
     public:
-
+        virtual void init(const UnstructuredGrid& g, int num_phases);
-        void init(const UnstructuredGrid& g, const int num_phases)
-        {
-            num_phases_ = num_phases;
-            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(num_phases * g.number_of_cells, 0.0);
-            for (int cell = 0; cell < g.number_of_cells; ++cell) {
-                // Defaulting the second saturation to 1.0.
-                // This will usually be oil in a water-oil case,
-                // gas in an oil-gas case.
-                // For proper initialization, one should not rely on this,
-                // but use available phase information instead.
-                sat_[num_phases*cell + 1] = 1.0;
-            }
-            gor_.resize(g.number_of_cells, 0.0);
-        }
-
-        enum ExtremalSat { MinSat, MaxSat };
 
         /// Set the first saturation to either its min or max value in
         /// the indicated cells. The second saturation value s2 is set
@@ -60,80 +40,20 @@ namespace Opm
         /// are unchanged.
         void setFirstSat(const std::vector<int>& cells,
                          const Opm::BlackoilPropertiesInterface& props,
-                         ExtremalSat es)
+                         ExtremalSat es);
-        {
-            if (cells.empty()) {
-                return;
-            }
-            const int n = cells.size();
-            assert(n > 0);
-            std::vector<double> smin(num_phases_*n);
-            std::vector<double> smax(num_phases_*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_[num_phases_*cell] = svals[num_phases_*ci];
-                sat_[num_phases_*cell + 1] = 1.0 - sat_[num_phases_*cell];
-            }
-        }
 
-        int numPhases() const
+        virtual bool equals(const SimulatorState& other,
-        {
+                            double epsilon = 1e-8) const;
-            return num_phases_;
-        }
 
-
-
-        bool equals(const BlackoilState& other, double epsilon = 1e-8) const {
-            bool equal = (num_phases_ == other.num_phases_);
-
-            equal = equal && (vectorApproxEqual( pressure() , other.pressure() , epsilon));
-            equal = equal && (vectorApproxEqual( facepressure() , other.facepressure() , epsilon));
-            equal = equal && (vectorApproxEqual( faceflux() , other.faceflux() , epsilon));
-            equal = equal && (vectorApproxEqual( surfacevol() , other.surfacevol() , epsilon));
-            equal = equal && (vectorApproxEqual( saturation() , other.saturation() , epsilon));
-            equal = equal && (vectorApproxEqual( gasoilratio() , other.gasoilratio() , epsilon));
-
-            return equal;
-        }
-
-
-        std::vector<double>& pressure    () { return press_ ; }
-        std::vector<double>& facepressure() { return fpress_; }
-        std::vector<double>& faceflux    () { return flux_  ; }
         std::vector<double>& surfacevol  () { return surfvol_; }
-        std::vector<double>& saturation  () { return sat_   ; }
         std::vector<double>& gasoilratio () { return gor_   ; }
 
-        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>& surfacevol  () const { return surfvol_; }
-        const std::vector<double>& saturation  () const { return sat_   ; }
         const std::vector<double>& gasoilratio () const { return gor_   ; }
 
     private:
-        int num_phases_;
-        std::vector<double> press_ ;
-        std::vector<double> fpress_;
-        std::vector<double> flux_  ;
         std::vector<double> surfvol_;
-        std::vector<double> sat_   ;
         std::vector<double> gor_   ;
-
-
-        static bool vectorApproxEqual(const std::vector<double>& v1, const std::vector<double>& v2 , double epsilon) {
-            if (v1.size() != v2.size())
-                return false;
-            
-            for (size_t i = 0; i < v1.size(); i++)
-                if (std::abs(v1[i] - v2[i]) > epsilon * (std::abs(v1[i]) + std::abs(v2[i])))
-                    return false;
-            
-            return true;
-        }
-
     };
 
 } // namespace Opm

opm/core/simulator/SimulatorState.cpp

@@ -0,0 +1,96 @@
+#include <opm/core/simulator/SimulatorState.hpp>
+#include <opm/core/grid.h>
+
+#include <cmath>
+
+using namespace Opm;
+
+bool
+SimulatorState::equals (const SimulatorState& other,
+                        double epsilon) const {
+    bool equal = (num_phases_ == other.num_phases_);
+
+    // if we use &=, then all the tests will be run regardless
+    equal = equal && vectorApproxEqual( pressure() , other.pressure() , epsilon);
+    equal = equal && vectorApproxEqual( facepressure() , other.facepressure() , epsilon);
+    equal = equal && vectorApproxEqual( faceflux() , other.faceflux() , epsilon);
+    equal = equal && vectorApproxEqual( saturation() , other.saturation() , epsilon);
+
+    return equal;
+}
+
+bool
+SimulatorState::vectorApproxEqual(const std::vector<double>& v1,
+                                  const std::vector<double>& v2,
+                                  double epsilon) {
+    if (v1.size() != v2.size()) {
+        return false;
+    }
+
+    for (size_t i = 0; i < v1.size(); i++) {
+        const double diff = std::abs(v1[i] - v2[i]);
+        const double scale = std::abs(v1[i]) + std::abs(v2[i]);
+        if (diff > epsilon * scale) {
+            return false;
+        }
+    }
+
+    return true;
+}
+
+void
+SimulatorState::init(const UnstructuredGrid& g, int num_phases)
+{
+    num_phases_ = num_phases;
+    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(num_phases_ * g.number_of_cells, 0.0);
+    for (int cell = 0; cell < g.number_of_cells; ++cell) {
+        // Defaulting the second saturation to 1.0.
+        // This will usually be oil in a water-oil case,
+        // gas in an oil-gas case.
+        // For proper initialization, one should not rely on this,
+        // but use available phase information instead.
+        sat_[num_phases_*cell + 1] = 1.0;
+    }
+}
+
+template <typename Props> void
+SimulatorState::setFirstSat(const std::vector<int>& cells,
+                            const Props& props,
+                            ExtremalSat es)
+{
+    if (cells.empty()) {
+        return;
+    }
+    int n = cells.size();
+    std::vector<double> smin(num_phases_*n);
+    std::vector<double> smax(num_phases_*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_[num_phases_*cell] = svals[num_phases_*ci];
+        sat_[num_phases_*cell + 1] = 1.0 - sat_[num_phases_*cell];
+    }
+}
+
+// template instantiations for all known (to this library) subclasses
+// of SimulatorState that will call this method. notice that there are
+// no empty angle brackets after "template" -- that would have been
+// specialization instead
+#include <opm/core/props/BlackoilPropertiesInterface.hpp>
+#include <opm/core/props/IncompPropertiesInterface.hpp>
+
+template void
+SimulatorState::setFirstSat <IncompPropertiesInterface> (
+        const std::vector<int> &cells,
+        const IncompPropertiesInterface &props,
+        ExtremalSat es);
+
+template void
+SimulatorState::setFirstSat <BlackoilPropertiesInterface> (
+        const std::vector<int> &cells,
+        const BlackoilPropertiesInterface &props,
+        ExtremalSat es);

opm/core/simulator/SimulatorState.hpp

@@ -0,0 +1,78 @@
+// Copyright (C) 2013 Uni Research AS
+// This file is licensed under the GNU General Public License v3.0
+
+#ifndef OPM_SIMULATORSTATE_HEADER_INCLUDED
+#define OPM_SIMULATORSTATE_HEADER_INCLUDED
+
+#include <vector>
+
+// forward declaration
+struct UnstructuredGrid;
+
+namespace Opm
+{
+    class SimulatorState
+    {
+    public:
+
+        virtual void init(const UnstructuredGrid& g, int num_phases);
+
+        enum ExtremalSat { MinSat, MaxSat };
+
+    protected:
+        /**
+         * Initialize the first saturation to maximum value. This method
+         * should be considered deprecated. Avoid to use it!
+         *
+         * \tparam Props Fluid and rock properties that pertain to this
+         *               kind of simulation. Currently, only Blackoil-
+         *               and IncompPropertiesInterface are supported.
+         */
+        template <typename Props>
+        void setFirstSat(const std::vector<int>& cells,
+                         const Props& props,
+                         ExtremalSat es);
+    public:
+        int numPhases() const { return num_phases_; }
+
+        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_   ; }
+
+        /**
+         * Compare this state with another, to see if they are different
+         * only within a small margin.
+         */
+        virtual bool equals(const SimulatorState& other,
+                            double epsilon = 1e-8) const;
+    private:
+        int num_phases_;
+        std::vector<double> press_ ;
+        std::vector<double> fpress_;
+        std::vector<double> flux_  ;
+        std::vector<double> sat_   ;
+
+    protected:
+        /**
+         * Check if two vectors are equal within a margin.
+         *
+         * @param epsilon Relative difference that is tolerated for the
+         *                vectors to still be considered equal.
+         *
+         * @return True if every element is within the margin, false if
+         *         there is at least one that is not.
+         */
+        static bool vectorApproxEqual(const std::vector<double>& v1,
+                                      const std::vector<double>& v2,
+                                      double epsilon);
+    };
+
+} // namespace Opm
+
+#endif // OPM_SIMULATORSTATE_HEADER_INCLUDED

opm/core/simulator/TwophaseState.hpp

@@ -20,77 +20,26 @@
 #ifndef OPM_TWOPHASESTATE_HEADER_INCLUDED
 #define OPM_TWOPHASESTATE_HEADER_INCLUDED
 
-#include <opm/core/grid.h>
 #include <opm/core/props/IncompPropertiesInterface.hpp>
-#include <vector>
+#include <opm/core/simulator/SimulatorState.hpp>
 
 namespace Opm
 {
 
     /// Simulator state for a two-phase simulator.
-    class TwophaseState
+    class TwophaseState : public SimulatorState
     {
     public:
-
-        void init(const UnstructuredGrid& g, int num_phases)
-        {
-            num_phases_ = num_phases;
-            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(num_phases_ * g.number_of_cells, 0.0);
-            for (int cell = 0; cell < g.number_of_cells; ++cell) {
-                // Defaulting the second saturation to 1.0.
-                // This will usually be oil in a water-oil case,
-                // gas in an oil-gas case.
-                // For proper initialization, one should not rely on this,
-                // but use available phase information instead.
-                sat_[num_phases_*cell + 1] = 1.0;
-            }
-        }
-
-        enum ExtremalSat { MinSat, MaxSat };
-
         void setFirstSat(const std::vector<int>& cells,
                          const Opm::IncompPropertiesInterface& props,
-                         ExtremalSat es)
+                         ExtremalSat es);
-        {
-            const int n = cells.size();
-            std::vector<double> smin(num_phases_*n);
-            std::vector<double> smax(num_phases_*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_[num_phases_*cell] = svals[num_phases_*ci];
-                sat_[num_phases_*cell + 1] = 1.0 - sat_[num_phases_*cell];
-            }
-        }
 
-        int numPhases() const
+        virtual bool equals (const SimulatorState& other,
-        {
+                             double epsilon = 1e-8) const;
-            return num_phases_;
-        }
-
-        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:
-        int num_phases_;
-        std::vector<double> press_ ;
-        std::vector<double> fpress_;
-        std::vector<double> flux_  ;
-        std::vector<double> sat_   ;
     };
 
 } // namespace Opm
 
+#include "TwophaseState_impl.hpp"
 
 #endif // OPM_TWOPHASESTATE_HEADER_INCLUDED

opm/core/simulator/TwophaseState_impl.hpp

@@ -0,0 +1,22 @@
+#ifndef OPM_TWOPHASESTATE_HEADER_INCLUDED
+#error Do not include this file directly!
+#endif
+
+namespace Opm {
+
+inline void
+TwophaseState::setFirstSat(const std::vector<int>& cells,
+                           const Opm::IncompPropertiesInterface& props,
+                           ExtremalSat es) {
+    SimulatorState::setFirstSat(cells, props, es);
+}
+
+inline bool
+TwophaseState::equals (const SimulatorState& other,
+                       double epsilon) const {
+    return dynamic_cast <const TwophaseState*>(&other)
+            ? SimulatorState::equals (other, epsilon)
+            : false;
+}
+
+} // namespace Opm