Merge remote-tracking branch 'upstream/master' into opm-parser-integrate

opm/core/props/BlackoilPhases.hpp

@@ -40,6 +40,40 @@ namespace Opm
         int phase_pos[MaxNumPhases];
     };
 
+    /// Check or assign presence of a formed, free phase.  Limited to
+    /// the 'BlackoilPhases'.
+    ///
+    /// Use a std::vector<PhasePresence> to represent the conditions
+    /// in an entire model.
+    class PhasePresence
+    {
+    public:
+        PhasePresence()
+            : present_(0)
+        {}
+
+        bool hasFreeWater() const { return present(BlackoilPhases::Aqua  ); }
+        bool hasFreeOil  () const { return present(BlackoilPhases::Liquid); }
+        bool hasFreeGas  () const { return present(BlackoilPhases::Vapour); }
+
+        void setFreeWater() { insert(BlackoilPhases::Aqua  ); }
+        void setFreeOil  () { insert(BlackoilPhases::Liquid); }
+        void setFreeGas  () { insert(BlackoilPhases::Vapour); }
+
+    private:
+        unsigned char present_;
+
+        bool present(const BlackoilPhases::PhaseIndex i) const
+        {
+            return present_ & (1 << i);
+        }
+
+        void insert(const BlackoilPhases::PhaseIndex i)
+        {
+            present_ |= (1 << i);
+        }
+    };
+
 } // namespace Opm
 
 #endif // OPM_BLACKOILPHASES_HEADER_INCLUDED

opm/core/simulator/initState_impl.hpp

@@ -30,6 +30,7 @@
 #include <opm/core/props/IncompPropertiesInterface.hpp>
 #include <opm/core/props/BlackoilPropertiesInterface.hpp>
 #include <opm/core/props/phaseUsageFromDeck.hpp>
+#include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 
 #include <iostream>
 #include <cmath>
@@ -606,14 +607,13 @@ namespace Opm
                 }
             }
         }
-    }
-
+    }   
     /// Initialize surface volume from pressure and saturation by z = As.
-    /// Here the RS factor is used to compute an intial z for the
+    /// Here the gas/oil ratio is used to compute an intial z for the
     /// computation of A.
     template <class Props, class State>
-    void initBlackoilSurfvol(const UnstructuredGrid& grid,
-                             const BlackoilPropertiesInterface& props,
+    void initBlackoilSurfvolUsingRS(const UnstructuredGrid& grid,
+                             const Props& props,
                              State& state)
     {
         if (props.numPhases() != 3) {
@@ -739,7 +739,8 @@ namespace Opm
                 int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c];
                 state.gasoilratio()[c] = rs_deck[c_deck];
             }
-            initBlackoilSurfvol(grid, props, state);
+            initBlackoilSurfvolUsingRS(grid, props, state);
+            computeSaturation(props,state);
         } else {
             OPM_THROW(std::runtime_error, "Temporarily, we require the RS field.");
         }

opm/core/utility/miscUtilitiesBlackoil.cpp

@@ -19,20 +19,26 @@
 
 
 #include "config.h"
+
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
-#include <opm/core/utility/Units.hpp>
+
 #include <opm/core/grid.h>
 #include <opm/core/wells.h>
+#include <opm/core/linalg/blas_lapack.h>
 #include <opm/core/props/BlackoilPropertiesInterface.hpp>
 #include <opm/core/simulator/BlackoilState.hpp>
 #include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
-#include <iostream>
+#include <opm/core/utility/Units.hpp>
+
 #include <algorithm>
-#include <functional>
 #include <cmath>
+#include <functional>
+#include <limits>
+#include <iostream>
 #include <iterator>
 
+
 namespace Opm
 {
 
@@ -277,6 +283,67 @@ namespace Opm
         }
     }
 
+    /// @brief Computes saturation from surface volume
+    void computeSaturation(const BlackoilPropertiesInterface& props,
+                          BlackoilState& state){
+
+        const int np = props.numPhases();
+        const int nc = props.numCells();
+        std::vector<double> allA(nc*np*np);
+        std::vector<int> allcells(nc);
+        for (int c = 0; c < nc; ++c) {
+            allcells[c] = c;
+        }
+
+        //std::vector<double> res_vol(np);
+        const std::vector<double>& z = state.surfacevol();
+
+        props.matrix(nc, &state.pressure()[0], &z[0], &allcells[0], &allA[0], 0);
+
+        // Linear solver.
+        MAT_SIZE_T n = np;
+        MAT_SIZE_T nrhs = 1;
+        MAT_SIZE_T lda = np;
+        std::vector<MAT_SIZE_T> piv(np);
+        MAT_SIZE_T ldb = np;
+        MAT_SIZE_T info = 0;
+
+
+        //double res_vol;
+        double tot_sat;
+        const double epsilon = std::sqrt(std::numeric_limits<double>::epsilon());
+
+        for (int c = 0; c < nc; ++c) {
+            double* A = &allA[c*np*np];
+            const double* z_loc = &z[c*np];
+            double* s = &state.saturation()[c*np];
+
+            for (int p = 0; p < np; ++p){
+                s[p] = z_loc[p];
+                }
+
+            dgesv_(&n, &nrhs, &A[0], &lda, &piv[0], &s[0], &ldb, &info);
+
+            tot_sat = 0;
+            for (int p = 0; p < np; ++p){
+                if (s[p] < epsilon) // saturation may be less then zero due to round of errors
+                    s[p] = 0;
+
+                tot_sat += s[p];
+            }
+
+            for (int p = 0; p < np; ++p){
+                s[p]  = s[p]/tot_sat;
+            }
+
+
+
+
+
+        }
+
+    }
+
 
     /// Compute two-phase transport source terms from well terms.
     /// Note: Unlike the incompressible version of this function,

opm/core/utility/miscUtilitiesBlackoil.hpp

@@ -136,6 +136,9 @@ namespace Opm
                            const double* saturation,
                            double* surfacevol);
 
+    /// Computes saturation from surface volume densities
+    void computeSaturation(const BlackoilPropertiesInterface& props,
+                          BlackoilState& state);
 
     /// Compute two-phase transport source terms from well terms.
     /// Note: Unlike the incompressible version of this function,