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/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index 7d01ade0e..575e18a41 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -29,6 +29,7 @@ #include #include #include +#include #include namespace Opm @@ -512,15 +513,23 @@ namespace Opm State& state) { const int num_phases = props.numPhases(); + const PhaseUsage pu = phaseUsageFromDeck(deck); + if (num_phases != pu.num_phases) { + THROW("initStateFromDeck(): user specified property object with " << num_phases << " phases, " + "found " << pu.num_phases << " phases in deck."); + } state.init(grid, num_phases); if (deck.hasField("EQUIL")) { if (num_phases != 2) { THROW("initStateFromDeck(): EQUIL-based init currently handling only two-phase scenarios."); } + if (pu.phase_used[BlackoilPhases::Vapour]) { + THROW("initStateFromDeck(): EQUIL-based init currently handling only oil-water scenario (no gas)."); + } // Set saturations depending on oil-water contact. const EQUIL& equil= deck.getEQUIL(); if (equil.equil.size() != 1) { - THROW("No region support yet."); + THROW("initStateFromDeck(): No region support yet."); } const double woc = equil.equil[0].water_oil_contact_depth_; initWaterOilContact(grid, props, woc, WaterBelow, state); @@ -528,37 +537,64 @@ namespace Opm const double datum_z = equil.equil[0].datum_depth_; const double datum_p = equil.equil[0].datum_depth_pressure_; initHydrostaticPressure(grid, props, woc, gravity, datum_z, datum_p, state); - } else if (deck.hasField("SWAT") && deck.hasField("PRESSURE")) { - // Set saturations from SWAT, pressure from PRESSURE. + } else if (deck.hasField("PRESSURE")) { + // Set saturations from SWAT/SGAS, pressure from PRESSURE. std::vector& s = state.saturation(); std::vector& p = state.pressure(); - const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& p_deck = deck.getFloatingPointValue("PRESSURE"); const int num_cells = grid.number_of_cells; if (num_phases == 2) { - for (int c = 0; c < num_cells; ++c) { - int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; - s[2*c] = sw_deck[c_deck]; - s[2*c + 1] = 1.0 - s[2*c]; - p[c] = p_deck[c_deck]; + if (!pu.phase_used[BlackoilPhases::Aqua]) { + // oil-gas: we require SGAS + if (!deck.hasField("SGAS")) { + THROW("initStateFromDeck(): missing SGAS keyword in 2-phase init"); + } + const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + gpos] = sg_deck[c_deck]; + s[2*c + opos] = 1.0 - sg_deck[c_deck]; + p[c] = p_deck[c_deck]; + } + } else { + // water-oil or water-gas: we require SWAT + if (!deck.hasField("SWAT")) { + THROW("initStateFromDeck(): missing SWAT keyword in 2-phase init"); + } + const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int nwpos = (wpos + 1) % 2; + for (int c = 0; c < num_cells; ++c) { + int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; + s[2*c + wpos] = sw_deck[c_deck]; + s[2*c + nwpos] = 1.0 - sw_deck[c_deck]; + p[c] = p_deck[c_deck]; + } } } else if (num_phases == 3) { - if (!deck.hasField("SGAS")) { - THROW("initStateFromDeck(): missing SGAS keyword in 3-phase init (only SWAT and PRESSURE found)."); + const bool has_swat_sgas = deck.hasField("SWAT") && deck.hasField("SGAS"); + if (!has_swat_sgas) { + THROW("initStateFromDeck(): missing SGAS or SWAT keyword in 3-phase init."); } + const int wpos = pu.phase_pos[BlackoilPhases::Aqua]; + const int gpos = pu.phase_pos[BlackoilPhases::Vapour]; + const int opos = pu.phase_pos[BlackoilPhases::Liquid]; + const std::vector& sw_deck = deck.getFloatingPointValue("SWAT"); const std::vector& sg_deck = deck.getFloatingPointValue("SGAS"); for (int c = 0; c < num_cells; ++c) { int c_deck = (grid.global_cell == NULL) ? c : grid.global_cell[c]; - s[3*c] = sw_deck[c_deck]; - s[3*c + 1] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); - s[3*c + 2] = sg_deck[c_deck]; + s[3*c + wpos] = sw_deck[c_deck]; + s[3*c + opos] = 1.0 - (sw_deck[c_deck] + sg_deck[c_deck]); + s[3*c + gpos] = sg_deck[c_deck]; p[c] = p_deck[c_deck]; } } else { THROW("initStateFromDeck(): init with SWAT etc. only available with 2 or 3 phases."); } } else { - THROW("initStateFromDeck(): we must either have EQUIL, or both SWAT and PRESSURE."); + THROW("initStateFromDeck(): we must either have EQUIL, or PRESSURE and SWAT/SOIL/SGAS."); } // Finally, init face pressures. 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.