Merge branch 'master' into ert

Conflicts:
	Makefile.am

opm/core/fluid/BlackoilPropertiesFromDeck.cpp

@@ -28,8 +28,8 @@ namespace Opm
     {
         rock_.init(deck, grid);
         pvt_.init(deck, 200);
-        SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
+        SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
-            = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
+            = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
         satprops_.reset(ptr);
         ptr->init(deck, grid, 200);
 
@@ -50,30 +50,43 @@ namespace Opm
         // Unfortunate lack of pointer smartness here...
         const int sat_samples = param.getDefault("sat_tab_size", 200);
         std::string threephase_model = param.getDefault<std::string>("threephase_model", "simple");
-        bool use_stone2 = (threephase_model == "stone2");
         if (sat_samples > 1) {
-            if (use_stone2) {
+            if (threephase_model == "stone2") {
                 SaturationPropsFromDeck<SatFuncStone2Uniform>* ptr
                     = new SaturationPropsFromDeck<SatFuncStone2Uniform>();
                 satprops_.reset(ptr);
                 ptr->init(deck, grid, sat_samples);
-            } else {
+            } else if (threephase_model == "simple") {
                 SaturationPropsFromDeck<SatFuncSimpleUniform>* ptr
                     = new SaturationPropsFromDeck<SatFuncSimpleUniform>();
                 satprops_.reset(ptr);
                 ptr->init(deck, grid, sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegUniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegUniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, grid, sat_samples);
+            } else {
+                THROW("Unknown threephase_model: " << threephase_model);
             }
         } else {
-            if (use_stone2) {
+            if (threephase_model == "stone2") {
                 SaturationPropsFromDeck<SatFuncStone2Nonuniform>* ptr
                     = new SaturationPropsFromDeck<SatFuncStone2Nonuniform>();
                 satprops_.reset(ptr);
                 ptr->init(deck, grid, sat_samples);
-            } else {
+            } else if (threephase_model == "simple") {
                 SaturationPropsFromDeck<SatFuncSimpleNonuniform>* ptr
                     = new SaturationPropsFromDeck<SatFuncSimpleNonuniform>();
                 satprops_.reset(ptr);
                 ptr->init(deck, grid, sat_samples);
+            } else if (threephase_model == "gwseg") {
+                SaturationPropsFromDeck<SatFuncGwsegNonuniform>* ptr
+                    = new SaturationPropsFromDeck<SatFuncGwsegNonuniform>();
+                satprops_.reset(ptr);
+                ptr->init(deck, grid, sat_samples);
+            } else {
+                THROW("Unknown threephase_model: " << threephase_model);
             }
         }
 

opm/core/fluid/SaturationPropsFromDeck.hpp

@@ -26,6 +26,7 @@
 #include <opm/core/fluid/blackoil/BlackoilPhases.hpp>
 #include <opm/core/fluid/SatFuncStone2.hpp>
 #include <opm/core/fluid/SatFuncSimple.hpp>
+#include <opm/core/fluid/SatFuncGwseg.hpp>
 #include <vector>
 
 struct UnstructuredGrid;

opm/core/linalg/LinearSolverIstl.cpp

@@ -26,6 +26,10 @@
 
 #include <opm/core/utility/have_boost_redef.hpp>
 
+// Silence compatibility warning from DUNE headers since we don't use
+// the deprecated member anyway (in this compilation unit)
+#define DUNE_COMMON_FIELDVECTOR_SIZE_IS_METHOD 1
+
 // TODO: clean up includes.
 #include <dune/common/deprecated.hh>
 #include <dune/istl/bvector.hh>

opm/core/pressure/CompressibleTpfa.cpp

@@ -123,7 +123,7 @@ namespace Opm
                                  WellState& well_state)
     {
         const int nc = grid_.number_of_cells;
-        const int nw = wells_->number_of_wells;
+        const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
 
         // Set up dynamic data.
         computePerSolveDynamicData(dt, state, well_state);
@@ -454,8 +454,8 @@ namespace Opm
         // std::vector<double> wellperf_A_;
         // std::vector<double> wellperf_phasemob_;
         const int np = props_.numPhases();
-        const int nw = wells_->number_of_wells;
+        const int nw = (wells_ != 0) ? wells_->number_of_wells : 0;
-        const int nperf = wells_->well_connpos[nw];
+        const int nperf = (wells_ != 0) ? wells_->well_connpos[nw] : 0;
         wellperf_A_.resize(nperf*np*np);
         wellperf_phasemob_.resize(nperf*np);
         // The A matrix is set equal to the perforation grid cells'
@@ -512,9 +512,9 @@ namespace Opm
         const double* z = &state.surfacevol()[0];
         UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
         CompletionData completion_data;
-        completion_data.gpot = &wellperf_gpot_[0];
+        completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0;
-        completion_data.A = &wellperf_A_[0];
+        completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0;
-        completion_data.phasemob = &wellperf_phasemob_[0];
+        completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0;
         cfs_tpfa_res_wells wells_tmp;
         wells_tmp.W = const_cast<Wells*>(wells_);
         wells_tmp.data = &completion_data;
@@ -599,9 +599,9 @@ namespace Opm
     {
         UnstructuredGrid* gg = const_cast<UnstructuredGrid*>(&grid_);
         CompletionData completion_data;
-        completion_data.gpot = const_cast<double*>(&wellperf_gpot_[0]);
+        completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast<double*>(&wellperf_gpot_[0]) : 0;
-        completion_data.A = const_cast<double*>(&wellperf_A_[0]);
+        completion_data.A = ! wellperf_A_.empty() ? const_cast<double*>(&wellperf_A_[0]) : 0;
-        completion_data.phasemob = const_cast<double*>(&wellperf_phasemob_[0]);
+        completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast<double*>(&wellperf_phasemob_[0]) : 0;
         cfs_tpfa_res_wells wells_tmp;
         wells_tmp.W = const_cast<Wells*>(wells_);
         wells_tmp.data = &completion_data;
@@ -609,6 +609,9 @@ namespace Opm
         forces.wells = &wells_tmp;
         forces.src = NULL;
 
+        double* wpress = ! well_state.bhp      ().empty() ? & well_state.bhp      ()[0] : 0;
+        double* wflux  = ! well_state.perfRates().empty() ? & well_state.perfRates()[0] : 0;
+
         cfs_tpfa_res_flux(gg,
                           &forces,
                           props_.numPhases(),
@@ -617,9 +620,9 @@ namespace Opm
                           &face_phasemob_[0],
                           &face_gravcap_[0],
                           &state.pressure()[0],
-                          &well_state.bhp()[0],
+                          wpress,
                           &state.faceflux()[0],
-                          &well_state.perfRates()[0]);
+                          wflux);
         cfs_tpfa_res_fpress(gg,
                             props_.numPhases(),
                             &htrans_[0],

opm/core/pressure/tpfa/cfs_tpfa_residual.c

@@ -142,7 +142,7 @@ impl_allocate(struct UnstructuredGrid   *G       ,
     nnu    = G->number_of_cells;
     nwperf = 0;
 
-    if (wells != NULL) {                    assert (wells->W != NULL);
+    if ((wells != NULL) && (wells->W != NULL)) {
         nnu    += wells->W->number_of_wells;
         nwperf  = wells->W->well_connpos[ wells->W->number_of_wells ];
     }
@@ -185,13 +185,15 @@ construct_matrix(struct UnstructuredGrid   *G    ,
 /* ---------------------------------------------------------------------- */
 {
     int    f, c1, c2, w, i, nc, nnu;
+    int    wells_present;
     size_t nnz;
 
     struct CSRMatrix *A;
 
     nc = nnu = G->number_of_cells;
-    if (wells != NULL) {
+
-        assert (wells->W != NULL);
+    wells_present = (wells != NULL) && (wells->W != NULL);
+    if (wells_present) {
         nnu += wells->W->number_of_wells;
     }
 
@@ -214,7 +216,7 @@ construct_matrix(struct UnstructuredGrid   *G    ,
             }
         }
 
-        if (wells != NULL) {
+        if (wells_present) {
             /* Well <-> cell connections */
             struct Wells *W = wells->W;
 
@@ -252,7 +254,7 @@ construct_matrix(struct UnstructuredGrid   *G    ,
             }
         }
 
-        if (wells != NULL) {
+        if (wells_present) {
             /* Fill well <-> cell connections */
             struct Wells *W = wells->W;
 
@@ -741,6 +743,29 @@ assemble_completion_to_cell(int c, int wdof, int np, double dt,
 }
 
 
+/* ---------------------------------------------------------------------- */
+static void
+welleq_coeff_shut(int np, struct cfs_tpfa_res_data *h,
+                  double *res, double *w2c, double *w2w)
+/* ---------------------------------------------------------------------- */
+{
+    int           p;
+    double        fwi;
+    const double *dpflux_w;
+
+    /* Sum reservoir phase flux derivatives set by
+     * compute_darcyflux_and_deriv(). */
+    dpflux_w = h->pimpl->flux_work + (1 * np);
+    for (p = 0, fwi = 0.0; p < np; p++) {
+        fwi += dpflux_w[ p ];
+    }
+
+    *res = 0.0;
+    *w2c = 0.0;
+    *w2w = fwi;
+}
+
+
 /* ---------------------------------------------------------------------- */
 static void
 welleq_coeff_bhp(int np, double dp, struct cfs_tpfa_res_data *h,
@@ -815,9 +840,10 @@ assemble_completion_to_well(int w, int c, int nc, int np,
     ctrl = W->ctrls[ w ];
 
     if (ctrl->current < 0) {
-        assert (0);             /* Shut wells currently not supported */
+        /* Interpreting a negative current control index to mean a shut well */
+        welleq_coeff_shut(np, h, &res, &w2c, &w2w);
     }
-
+    else {
         switch (ctrl->type[ ctrl->current ]) {
         case BHP :
             welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ],
@@ -833,6 +859,7 @@ assemble_completion_to_well(int w, int c, int nc, int np,
             assert (0);             /* Surface rate currently not supported */
             break;
         }
+    }
 
     /* Assemble completion contributions */
     wdof = nc + w;
@@ -854,7 +881,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells   *wells ,
                       struct cfs_tpfa_res_data    *h     )
 {
     int           w, i, c, np, np2, nc;
-    int           is_neumann;
+    int           is_neumann, is_open;
     double        pw, dp;
     double       *WI, *gpot, *pmobp;
     const double *Ac, *dAc;
@@ -876,6 +903,7 @@ assemble_well_contrib(struct cfs_tpfa_res_wells   *wells ,
 
     for (w = i = 0; w < W->number_of_wells; w++) {
         pw = wpress[ w ];
+        is_open = W->ctrls[w]->current >= 0;
 
         for (; i < W->well_connpos[w + 1];
              i++, gpot += np, pmobp += np) {
@@ -888,7 +916,9 @@ assemble_well_contrib(struct cfs_tpfa_res_wells   *wells ,
 
             init_completion_contrib(i, np, Ac, dAc, h->pimpl);
 
+            if (is_open) {
                 assemble_completion_to_cell(c, nc + w, np, dt, h);
+            }
 
             /* Prepare for RESV controls */
             compute_darcyflux_and_deriv(np, WI[i], dp, pmobp, gpot,
@@ -1127,8 +1157,7 @@ cfs_tpfa_res_construct(struct UnstructuredGrid   *G      ,
         nf     = G->number_of_faces;
         nwperf = 0;
 
-        if (wells != NULL) {
+        if ((wells != NULL) && (wells->W != NULL)) {
-            assert (wells->W != NULL);
             nwperf = wells->W->well_connpos[ wells->W->number_of_wells ];
         }
 
@@ -1194,7 +1223,9 @@ cfs_tpfa_res_assemble(struct UnstructuredGrid     *G        ,
         assemble_cell_contrib(G, c, h);
     }
 
-    if ((forces != NULL) && (forces->wells != NULL)) {
+    if ((forces           != NULL) &&
+        (forces->wells    != NULL) &&
+        (forces->wells->W != NULL)) {
         compute_well_compflux_and_deriv(forces->wells, cq->nphases,
                                         cpress, wpress, h->pimpl);
 
@@ -1297,8 +1328,9 @@ cfs_tpfa_res_flux(struct UnstructuredGrid    *G        ,
 {
     compute_flux(G, np, trans, pmobf, gravcap_f, cpress, fflux);
 
-    if ((forces != NULL) && (forces->wells != NULL) &&
+    if ((forces           != NULL) &&
-        (wpress != NULL) && (wflux         != NULL)) {
+        (forces->wells    != NULL) &&
+        (forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) {
 
         compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux);
     }

opm/core/utility/miscUtilities.cpp

@@ -374,7 +374,10 @@ namespace Opm
                     if (perf_rate > 0.0) {
                         // perf_rate is a total inflow rate, we want a water rate.
                         if (wells->type[w] != INJECTOR) {
-                            std::cout << "**** Warning: crossflow in well with index " << w << " ignored." << std::endl;
+                            std::cout << "**** Warning: crossflow in well "
+                                      << w << " perf " << perf - wells->well_connpos[w]
+                                      << " ignored. 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);

opm/core/wells/WellsManager.cpp

@@ -717,7 +717,6 @@ namespace Opm
 
         if (deck.hasField("WELOPEN")) {
             const WELOPEN& welopen = deck.getWELOPEN();
-
             for (size_t i = 0; i < welopen.welopen.size(); ++i) {
                 WelopenLine line = welopen.welopen[i];
                 std::string wellname = line.well_;
@@ -725,22 +724,22 @@ namespace Opm
                 if (it == well_names_to_index.end()) {
                     THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS.");
                 }
-        		int index = it->second;
+                const int index = it->second;
                 if (line.openshutflag_ == "SHUT") {
-        			// We currently don't care if the well is open or not.
+                    int& cur_ctrl = w_->ctrls[index]->current;
-        			/// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }?
+                    if (cur_ctrl >= 0) {
-        			//ASSERT(w_->ctrls[index]->current < 0);
+                        cur_ctrl = ~cur_ctrl;
+                    }
+                    ASSERT(w_->ctrls[index]->current < 0);
                 } else if (line.openshutflag_ == "OPEN") {
-        			//ASSERT(w_->ctrls[index]->current >= 0);
+                    int& cur_ctrl = w_->ctrls[index]->current;
+                    if (cur_ctrl < 0) {
+                        cur_ctrl = ~cur_ctrl;
+                    }
+                    ASSERT(w_->ctrls[index]->current >= 0);
                 } else {
                     THROW("Unknown Open/close keyword: \"" << line.openshutflag_<< "\". Allowed values: OPEN, SHUT.");
                 }
-
-        		// We revert back to it's original control.
-        		// Note that this is OK as ~~ = id.
-        		w_->ctrls[index]->current = ~w_->ctrls[index]->current;
-
-
             }
         }