diff --git a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp index e2457f39b..271f93951 100644 --- a/opm/core/fluid/BlackoilPropertiesFromDeck.cpp +++ b/opm/core/fluid/BlackoilPropertiesFromDeck.cpp @@ -28,8 +28,8 @@ namespace Opm { rock_.init(deck, grid); pvt_.init(deck, 200); - SaturationPropsFromDeck* ptr - = new SaturationPropsFromDeck(); + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); 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("threephase_model", "simple"); - bool use_stone2 = (threephase_model == "stone2"); if (sat_samples > 1) { - if (use_stone2) { + if (threephase_model == "stone2") { SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); - } else { + } else if (threephase_model == "simple") { SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + 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* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); - } else { + } else if (threephase_model == "simple") { SaturationPropsFromDeck* ptr = new SaturationPropsFromDeck(); satprops_.reset(ptr); ptr->init(deck, grid, sat_samples); + } else if (threephase_model == "gwseg") { + SaturationPropsFromDeck* ptr + = new SaturationPropsFromDeck(); + satprops_.reset(ptr); + ptr->init(deck, grid, sat_samples); + } else { + THROW("Unknown threephase_model: " << threephase_model); } } diff --git a/opm/core/fluid/SaturationPropsFromDeck.hpp b/opm/core/fluid/SaturationPropsFromDeck.hpp index 0a0599ec8..ea2acb342 100644 --- a/opm/core/fluid/SaturationPropsFromDeck.hpp +++ b/opm/core/fluid/SaturationPropsFromDeck.hpp @@ -26,6 +26,7 @@ #include #include #include +#include #include struct UnstructuredGrid; diff --git a/opm/core/linalg/LinearSolverIstl.cpp b/opm/core/linalg/LinearSolverIstl.cpp index 3c4435e8e..a3f9583ee 100644 --- a/opm/core/linalg/LinearSolverIstl.cpp +++ b/opm/core/linalg/LinearSolverIstl.cpp @@ -26,6 +26,10 @@ #include +// 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 #include diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 859585991..42b3fba75 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/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 wellperf_A_; // std::vector wellperf_phasemob_; const int np = props_.numPhases(); - const int nw = wells_->number_of_wells; - const int nperf = wells_->well_connpos[nw]; + const int nw = (wells_ != 0) ? wells_->number_of_wells : 0; + 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(&grid_); CompletionData completion_data; - completion_data.gpot = &wellperf_gpot_[0]; - completion_data.A = &wellperf_A_[0]; - completion_data.phasemob = &wellperf_phasemob_[0]; + completion_data.gpot = ! wellperf_gpot_.empty() ? &wellperf_gpot_[0] : 0; + completion_data.A = ! wellperf_A_.empty() ? &wellperf_A_[0] : 0; + completion_data.phasemob = ! wellperf_phasemob_.empty() ? &wellperf_phasemob_[0] : 0; cfs_tpfa_res_wells wells_tmp; wells_tmp.W = const_cast(wells_); wells_tmp.data = &completion_data; @@ -599,9 +599,9 @@ namespace Opm { UnstructuredGrid* gg = const_cast(&grid_); CompletionData completion_data; - completion_data.gpot = const_cast(&wellperf_gpot_[0]); - completion_data.A = const_cast(&wellperf_A_[0]); - completion_data.phasemob = const_cast(&wellperf_phasemob_[0]); + completion_data.gpot = ! wellperf_gpot_.empty() ? const_cast(&wellperf_gpot_[0]) : 0; + completion_data.A = ! wellperf_A_.empty() ? const_cast(&wellperf_A_[0]) : 0; + completion_data.phasemob = ! wellperf_phasemob_.empty() ? const_cast(&wellperf_phasemob_[0]) : 0; cfs_tpfa_res_wells wells_tmp; wells_tmp.W = const_cast(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], diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index eb8dc8f34..0e9b6f7d2 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/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,23 +840,25 @@ 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 ], + h, &res, &w2c, &w2w); + break; - switch (ctrl->type[ ctrl->current ]) { - case BHP : - welleq_coeff_bhp(np, pw - ctrl->target[ ctrl->current ], - h, &res, &w2c, &w2w); - break; + case RESERVOIR_RATE: + assert (W->number_of_phases == np); + welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w); + break; - case RESERVOIR_RATE: - assert (W->number_of_phases == np); - welleq_coeff_resv(np, h, ctrl, &res, &w2c, &w2w); - break; - - case SURFACE_RATE: - assert (0); /* Surface rate currently not supported */ - break; + case SURFACE_RATE: + assert (0); /* Surface rate currently not supported */ + break; + } } /* Assemble completion contributions */ @@ -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); - assemble_completion_to_cell(c, nc + w, np, dt, h); + 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) { - assert (wells->W != NULL); + if ((wells != NULL) && (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) && - (wpress != NULL) && (wflux != NULL)) { + if ((forces != NULL) && + (forces->wells != NULL) && + (forces->wells->W != NULL) && (wpress != NULL) && (wflux != NULL)) { compute_wflux(np, forces->wells, pmobc, cpress, wpress, wflux); } diff --git a/opm/core/utility/miscUtilities.cpp b/opm/core/utility/miscUtilities.cpp index aeafbfca7..53e997731 100644 --- a/opm/core/utility/miscUtilities.cpp +++ b/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); diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index c31345fa4..95cb0369a 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -716,32 +716,31 @@ namespace Opm #endif 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_; - std::map::const_iterator it = well_names_to_index.find(wellname); - 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; - if (line.openshutflag_ == "SHUT") { - // We currently don't care if the well is open or not. - /// \TODO Should this perhaps be allowed? I.e. should it be if(well_shut) { shutwell(); } else { /* do nothing*/ }? - //ASSERT(w_->ctrls[index]->current < 0); - } else if (line.openshutflag_ == "OPEN") { - //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; - - - } + 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_; + std::map::const_iterator it = well_names_to_index.find(wellname); + if (it == well_names_to_index.end()) { + THROW("Trying to open/shut well with name: \"" << wellname<<"\" but it's not registered under WELSPECS."); + } + const int index = it->second; + if (line.openshutflag_ == "SHUT") { + int& cur_ctrl = w_->ctrls[index]->current; + if (cur_ctrl >= 0) { + cur_ctrl = ~cur_ctrl; + } + ASSERT(w_->ctrls[index]->current < 0); + } else if (line.openshutflag_ == "OPEN") { + 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."); + } + } } // Build the well_collection_ well group hierarchy.