Merge remote-tracking branch 'upstream/master'

opm/core/grid/cpgpreprocess/preprocess.h

@@ -56,8 +56,8 @@ extern "C" {
     struct grdecl {
         int           dims[3]; /**< Cartesian box dimensions. */
         const double *coord;   /**< Pillar end-points. */
-        const double *zcorn;   /**< Explicit "active" map.  May be NULL.*/
+        const double *zcorn;   /**< Corner-point depths. */
-        const int    *actnum;  /**< Corner-point depths. */
+        const int    *actnum;  /**< Explicit "active" map.  May be NULL.*/
     };
 
     /**

opm/core/pressure/CompressibleTpfa.cpp

@@ -632,6 +632,23 @@ namespace Opm
                             &state.pressure()[0],
                             &state.faceflux()[0],
                             &state.facepressure()[0]);
+
+        // Compute well perforation pressures (not done by the C code).
+        if (wells_ != 0) {
+            const int nw = wells_->number_of_wells;
+            const int np = props_.numPhases();
+            for (int w = 0; w < nw; ++w) {
+                const double* comp_frac = &wells_->comp_frac[np*w];
+                for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) {
+                    const double bhp = well_state.bhp()[w];
+                    double perf_p = bhp;
+                    for (int phase = 0; phase < np; ++phase) {
+                        perf_p += wellperf_gpot_[np*j + phase]*comp_frac[phase];
+                    }
+                    well_state.perfPress()[j] = perf_p;
+                }
+            }
+        }
     }
 
 } // namespace Opm

opm/core/simulator/SimulatorCompressibleTwophase.cpp

@@ -196,6 +196,7 @@ namespace Opm
             if (!file) {
                 THROW("Failed to open " << fname.str());
             }
+            file.precision(15);
             const std::vector<double>& d = *(it->second);
             std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
         }
@@ -434,9 +435,8 @@ namespace Opm
                 computePorevolume(grid_, props_.porosity(), *rock_comp_, state.pressure(), porevol);
             }
 
-            // Process transport sources (to include bdy terms and well flows).
+            // Process transport sources from well flows.
-            Opm::computeTransportSource(grid_, src_, state.faceflux(), 1.0,
+            Opm::computeTransportSource(props_, wells_, well_state, transport_src);
-                                        wells_, well_state.perfRates(), transport_src);
 
             // Solve transport.
             transport_timer.start();
@@ -449,9 +449,7 @@ namespace Opm
                 tsolver_.solve(&state.faceflux()[0], &state.pressure()[0],
                                &initial_porevol[0], &porevol[0], &transport_src[0], stepsize,
                                state.saturation(), state.surfacevol());
-                Opm::computeInjectedProduced(props_,
+                Opm::computeInjectedProduced(props_, state, transport_src, stepsize, injected, produced);
-                                             state.pressure(), state.surfacevol(), state.saturation(),
-                                             transport_src, stepsize, injected, produced);
                 if (gravity_ != 0 && use_segregation_split_) {
                     tsolver_.solveGravity(columns_, stepsize, state.saturation(), state.surfacevol());
                 }

opm/core/simulator/SimulatorIncompTwophase.cpp

@@ -244,6 +244,7 @@ namespace Opm
             if (!file) {
                 THROW("Failed to open " << fname.str());
             }
+            file.precision(15);
             const std::vector<double>& d = *(it->second);
             std::copy(d.begin(), d.end(), std::ostream_iterator<double>(file, "\n"));
         }

opm/core/simulator/WellState.hpp

@@ -49,7 +49,8 @@ namespace Opm
                         bhp_[w] = state.pressure()[cell];
                     }
                 }
-                perfrates_.resize(wells->well_connpos[nw]);
+                perfrates_.resize(wells->well_connpos[nw], 0.0);
+                perfpress_.resize(wells->well_connpos[nw], -1e100);
             }
         }
 
@@ -61,9 +62,14 @@ namespace Opm
         std::vector<double>& perfRates() { return perfrates_; }
         const std::vector<double>& perfRates() const { return perfrates_; }
 
+        /// One pressure per well connection.
+        std::vector<double>& perfPress() { return perfpress_; }
+        const std::vector<double>& perfPress() const { return perfpress_; }
+
     private:
         std::vector<double> bhp_;
         std::vector<double> perfrates_;
+        std::vector<double> perfpress_;
     };
 
 } // namespace Opm

opm/core/transport/reorder/TransportModelCompressibleTwophase.cpp

@@ -152,7 +152,7 @@ namespace Opm
             B_cell = 1.0/tm.A_[np*np*cell + 0];
             double src_flux       = -tm.source_[cell];
             bool src_is_inflow = src_flux < 0.0;
-            influx  =  src_is_inflow ? src_flux : 0.0;
+            influx  =  src_is_inflow ? B_cell* src_flux : 0.0;
             outflux = !src_is_inflow ? src_flux : 0.0;
             comp_term = (tm.porevolume_[cell] - tm.porevolume0_[cell])/tm.porevolume0_[cell];
             dtpv    = tm.dt_/tm.porevolume0_[cell];

opm/core/utility/miscUtilitiesBlackoil.cpp

@@ -21,7 +21,10 @@
 #include <opm/core/utility/miscUtilitiesBlackoil.hpp>
 #include <opm/core/utility/Units.hpp>
 #include <opm/core/grid.h>
+#include <opm/core/newwells.h>
 #include <opm/core/fluid/BlackoilPropertiesInterface.hpp>
+#include <opm/core/simulator/BlackoilState.hpp>
+#include <opm/core/simulator/WellState.hpp>
 #include <opm/core/utility/ErrorMacros.hpp>
 #include <algorithm>
 #include <functional>
@@ -31,53 +34,70 @@
 namespace Opm
 {
 
-    /// @brief Computes injected and produced volumes of all phases.
+    /// @brief Computes injected and produced surface volumes of all phases.
     /// Note 1: assumes that only the first phase is injected.
     /// Note 2: assumes that transport has been done with an
     ///         implicit method, i.e. that the current state
     ///         gives the mobilities used for the preceding timestep.
-    /// @param[in]  props     fluid and rock properties.
+    /// Note 3: Gives surface volume values, not reservoir volumes
-    /// @param[in]  p         pressure (one value per cell)
+    ///         (as the incompressible version of the function does).
-    /// @param[in]  z         surface-volume values (for all P phases)
+    ///         Also, assumes that transport_src is given in surface volumes
-    /// @param[in]  s         saturation values (for all P phases)
+    ///         for injector terms!
-    /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  props           fluid and rock properties.
-    /// @param[in]  dt        timestep used
+    /// @param[in]  state           state variables (pressure, sat, surfvol)
-    /// @param[out] injected  must point to a valid array with P elements,
+    /// @param[in]  transport_src   if < 0: total resv outflow, if > 0: first phase surfv inflow
-    ///                       where P = s.size()/src.size().
+    /// @param[in]  dt              timestep used
-    /// @param[out] produced  must also point to a valid array with P elements.
+    /// @param[out] injected        must point to a valid array with P elements,
+    ///                             where P = s.size()/src.size().
+    /// @param[out] produced        must also point to a valid array with P elements.
     void computeInjectedProduced(const BlackoilPropertiesInterface& props,
-                                 const std::vector<double>& press,
+                                 const BlackoilState& state,
-                                 const std::vector<double>& z,
+                                 const std::vector<double>& transport_src,
-                                 const std::vector<double>& s,
-                                 const std::vector<double>& src,
                                  const double dt,
                                  double* injected,
                                  double* produced)
     {
-        const int num_cells = src.size();
+        const int num_cells = transport_src.size();
-        const int np = s.size()/src.size();
+        if (props.numCells() != num_cells) {
-        if (int(s.size()) != num_cells*np) {
+            THROW("Size of transport_src vector does not match number of cells in props.");
-            THROW("Sizes of s and src vectors do not match.");
         }
+        const int np = props.numPhases();
+        if (int(state.saturation().size()) != num_cells*np) {
+            THROW("Sizes of state vectors do not match number of cells.");
+        }
+        const std::vector<double>& press = state.pressure();
+        const std::vector<double>& s = state.saturation();
+        const std::vector<double>& z = state.surfacevol();
         std::fill(injected, injected + np, 0.0);
         std::fill(produced, produced + np, 0.0);
         std::vector<double> visc(np);
         std::vector<double> mob(np);
+        std::vector<double> A(np*np);
+        std::vector<double> prod_resv_phase(np);
+        std::vector<double> prod_surfvol(np);
         for (int c = 0; c < num_cells; ++c) {
-            if (src[c] > 0.0) {
+            if (transport_src[c] > 0.0) {
-                injected[0] += src[c]*dt;
+                injected[0] += transport_src[c]*dt;
-            } else if (src[c] < 0.0) {
+            } else if (transport_src[c] < 0.0) {
-                const double flux = -src[c]*dt;
+                const double flux = -transport_src[c]*dt;
                 const double* sat = &s[np*c];
                 props.relperm(1, sat, &c, &mob[0], 0);
                 props.viscosity(1, &press[c], &z[np*c], &c, &visc[0], 0);
+                props.matrix(1, &press[c], &z[np*c], &c, &A[0], 0);
                 double totmob = 0.0;
                 for (int p = 0; p < np; ++p) {
                     mob[p] /= visc[p];
                     totmob += mob[p];
                 }
+                std::fill(prod_surfvol.begin(), prod_surfvol.end(), 0.0);
                 for (int p = 0; p < np; ++p) {
-                    produced[p] += (mob[p]/totmob)*flux;
+                    prod_resv_phase[p] = (mob[p]/totmob)*flux;
+                    for (int q = 0; q < np; ++q) {
+                        prod_surfvol[q] += prod_resv_phase[p]*A[q + np*p];
+                    }
+                }
+                for (int p = 0; p < np; ++p) {
+                    produced[p] += prod_surfvol[p];
                 }
             }
         }
@@ -251,4 +271,58 @@ namespace Opm
         }
     }
 
+
+    /// Compute two-phase transport source terms from well terms.
+    /// Note: Unlike the incompressible version of this function,
+    ///       this version computes surface volume injection rates,
+    ///       production rates are still total reservoir volumes.
+    /// \param[in]  props         Fluid and rock properties.
+    /// \param[in]  wells         Wells data structure.
+    /// \param[in]  well_state    Well pressures and fluxes.
+    /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
+    ///                           (+) positive  inflow of first (water) phase (surface volume),
+    ///                           (-) negative  total outflow of both phases (reservoir volume).
+    void computeTransportSource(const BlackoilPropertiesInterface& props,
+                                const Wells* wells,
+                                const WellState& well_state,
+                                std::vector<double>& transport_src)
+    {
+        int nc = props.numCells();
+        transport_src.clear();
+        transport_src.resize(nc, 0.0);
+        // Well contributions.
+        if (wells) {
+            const int nw = wells->number_of_wells;
+            const int np = wells->number_of_phases;
+            if (np != 2) {
+                THROW("computeTransportSource() requires a 2 phase case.");
+            }
+            std::vector<double> A(np*np);
+            for (int w = 0; w < nw; ++w) {
+                const double* comp_frac = wells->comp_frac + np*w;
+                for (int perf = wells->well_connpos[w]; perf < wells->well_connpos[w + 1]; ++perf) {
+                    const int perf_cell = wells->well_cells[perf];
+                    double perf_rate = well_state.perfRates()[perf];
+                    if (perf_rate > 0.0) {
+                        // perf_rate is a total inflow reservoir rate, we want a surface water rate.
+                        if (wells->type[w] != INJECTOR) {
+                            std::cout << "**** Warning: crossflow in well "
+                                      << w << " perf " << perf - wells->well_connpos[w]
+                                      << " ignored. Reservoir rate was "
+                                      << perf_rate/Opm::unit::day << " m^3/day." << std::endl;
+                            perf_rate = 0.0;
+                        } else {
+                            ASSERT(std::fabs(comp_frac[0] + comp_frac[1] - 1.0) < 1e-6);
+                            perf_rate *= comp_frac[0]; // Water reservoir volume rate.
+                            props.matrix(1, &well_state.perfPress()[perf], comp_frac, &perf_cell, &A[0], 0);
+                            perf_rate *= A[0];         // Water surface volume rate.
+                        }
+                    }
+                    transport_src[perf_cell] += perf_rate;
+                }
+            }
+        }
+    }
+
+
 } // namespace Opm

opm/core/utility/miscUtilitiesBlackoil.hpp

@@ -22,36 +22,40 @@
 
 #include <vector>
 
-struct UnstructuredGrid;
+struct Wells;
 
 namespace Opm
 {
 
     class BlackoilPropertiesInterface;
+    class BlackoilState;
+    class WellState;
 
-    /// @brief Computes injected and produced volumes of all phases.
+
+    /// @brief Computes injected and produced surface volumes of all phases.
     /// Note 1: assumes that only the first phase is injected.
     /// Note 2: assumes that transport has been done with an
     ///         implicit method, i.e. that the current state
     ///         gives the mobilities used for the preceding timestep.
-    /// @param[in]  props     fluid and rock properties.
+    /// Note 3: Gives surface volume values, not reservoir volumes
-    /// @param[in]  p         pressure (one value per cell)
+    ///         (as the incompressible version of the function does).
-    /// @param[in]  z         surface-volume values (for all P phases)
+    ///         Also, assumes that transport_src is given in surface volumes
-    /// @param[in]  s         saturation values (for all P phases)
+    ///         for injector terms!
-    /// @param[in]  src       if < 0: total outflow, if > 0: first phase inflow.
+    /// @param[in]  props           fluid and rock properties.
-    /// @param[in]  dt        timestep used
+    /// @param[in]  state           state variables (pressure, sat, surfvol)
-    /// @param[out] injected  must point to a valid array with P elements,
+    /// @param[in]  transport_src   if < 0: total resv outflow, if > 0: first phase surfv inflow
-    ///                       where P = s.size()/src.size().
+    /// @param[in]  dt              timestep used
-    /// @param[out] produced  must also point to a valid array with P elements.
+    /// @param[out] injected        must point to a valid array with P elements,
+    ///                             where P = s.size()/src.size().
+    /// @param[out] produced        must also point to a valid array with P elements.
     void computeInjectedProduced(const BlackoilPropertiesInterface& props,
-                                 const std::vector<double>& p,
+                                 const BlackoilState& state,
-                                 const std::vector<double>& z,
+                                 const std::vector<double>& transport_src,
-                                 const std::vector<double>& s,
-                                 const std::vector<double>& src,
                                  const double dt,
                                  double* injected,
                                  double* produced);
 
+
     /// @brief Computes total mobility for a set of saturation values.
     /// @param[in]  props     rock and fluid properties
     /// @param[in]  cells     cells with which the saturation values are associated
@@ -66,6 +70,7 @@ namespace Opm
                               const std::vector<double>& s,
                               std::vector<double>& totmob);
 
+
     /// @brief Computes total mobility and omega for a set of saturation values.
     /// @param[in]  props     rock and fluid properties
     /// @param[in]  cells     cells with which the saturation values are associated
@@ -131,6 +136,22 @@ namespace Opm
                            const double* saturation,
                            double* surfacevol);
 
+
+    /// Compute two-phase transport source terms from well terms.
+    /// Note: Unlike the incompressible version of this function,
+    ///       this version computes surface volume injection rates,
+    ///       production rates are still total reservoir volumes.
+    /// \param[in]  props         Fluid and rock properties.
+    /// \param[in]  wells         Wells data structure.
+    /// \param[in]  well_state    Well pressures and fluxes.
+    /// \param[out] transport_src The transport source terms. They are to be interpreted depending on sign:
+    ///                           (+) positive  inflow of first (water) phase (surface volume),
+    ///                           (-) negative  total outflow of both phases (reservoir volume).
+    void computeTransportSource(const BlackoilPropertiesInterface& props,
+                                const Wells* wells,
+                                const WellState& well_state,
+                                std::vector<double>& transport_src);
+
 } // namespace Opm
 
 #endif // OPM_MISCUTILITIESBLACKOIL_HEADER_INCLUDED