Compute initial saturations from surface volumes

Add new function is added that computes saturation from surface volumes
solving z = As for each cell. This function is used to compute an
intial guess to the saturations in initState_impl.hpp.

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,8 +607,7 @@ 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
     /// computation of A.
@@ -740,6 +740,7 @@ namespace Opm
                 state.gasoilratio()[c] = rs_deck[c_deck];
             }
             initBlackoilSurfvol(grid, props, state);
+            computeSaturation(props,state);
         } else {
             OPM_THROW(std::runtime_error, "Temporarily, we require the RS field.");
         }

opm/core/utility/miscUtilitiesBlackoil.cpp

@@ -32,6 +32,8 @@
 #include <functional>
 #include <cmath>
 #include <iterator>
+#include <opm/core/linalg/blas_lapack.h>
+
 
 namespace Opm
 {
@@ -277,6 +279,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,