From bdcf0291e07be19c0abe1dd908b741ad2667b929 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 09:52:13 +0200 Subject: [PATCH 01/10] Fix error message. --- opm/core/wells/WellsManager.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/wells/WellsManager.cpp b/opm/core/wells/WellsManager.cpp index 0215b053e..c31345fa4 100644 --- a/opm/core/wells/WellsManager.cpp +++ b/opm/core/wells/WellsManager.cpp @@ -544,7 +544,7 @@ namespace Opm cf[pu.phase_pos[BlackoilPhases::Liquid]] = 1.0; } else if (wci_line.injector_type_[0] == 'G') { if (!pu.phase_used[BlackoilPhases::Vapour]) { - THROW("Water phase not used, yet found water-injecting well."); + THROW("Gas phase not used, yet found gas-injecting well."); } cf[pu.phase_pos[BlackoilPhases::Vapour]] = 1.0; } From fa6b7729725a4d89d79699c7d0e32224b67367aa Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 09:53:11 +0200 Subject: [PATCH 02/10] Changed well initialization and property calculation. Bhp is now initialized to bhp target for bhp-controlled wells. Mobilities and pvt properties are now calculated from well perforation pressure and injection specifications for injectors, producers still use cell properties as before. --- opm/core/pressure/CompressibleTpfa.cpp | 39 +++++++++++++++++++++----- opm/core/simulator/WellState.hpp | 13 +++++++-- 2 files changed, 42 insertions(+), 10 deletions(-) diff --git a/opm/core/pressure/CompressibleTpfa.cpp b/opm/core/pressure/CompressibleTpfa.cpp index 954e07fa4..859585991 100644 --- a/opm/core/pressure/CompressibleTpfa.cpp +++ b/opm/core/pressure/CompressibleTpfa.cpp @@ -36,6 +36,7 @@ #include #include #include +#include namespace Opm { @@ -446,7 +447,7 @@ namespace Opm /// Compute per-iteration dynamic properties for wells. void CompressibleTpfa::computeWellDynamicData(const double /*dt*/, const BlackoilState& /*state*/, - const WellState& /*well_state*/) + const WellState& well_state) { // These are the variables that get computed by this function: // @@ -458,18 +459,42 @@ namespace Opm wellperf_A_.resize(nperf*np*np); wellperf_phasemob_.resize(nperf*np); // The A matrix is set equal to the perforation grid cells' - // matrix, for both injectors and producers. + // matrix for producers, computed from bhp and injection + // component fractions from // The mobilities are set equal to the perforation grid cells' - // mobilities, for both injectors and producers. + // mobilities for producers. + std::vector mu(np); for (int w = 0; w < nw; ++w) { + bool producer = (wells_->type[w] == PRODUCER); + const double* comp_frac = &wells_->comp_frac[np*w]; for (int j = wells_->well_connpos[w]; j < wells_->well_connpos[w+1]; ++j) { const int c = wells_->well_cells[j]; - const double* cA = &cell_A_[np*np*c]; double* wpA = &wellperf_A_[np*np*j]; - std::copy(cA, cA + np*np, wpA); - const double* cM = &cell_phasemob_[np*c]; double* wpM = &wellperf_phasemob_[np*j]; - std::copy(cM, cM + np, wpM); + if (producer) { + const double* cA = &cell_A_[np*np*c]; + std::copy(cA, cA + np*np, wpA); + const double* cM = &cell_phasemob_[np*c]; + std::copy(cM, cM + np, wpM); + } else { + 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]; + } + // Hack warning: comp_frac is used as a component + // surface-volume variable in calls to matrix() and + // viscosity(), but as a saturation in the call to + // relperm(). This is probably ok as long as injectors + // only inject pure fluids. + props_.matrix(1, &perf_p, comp_frac, &c, wpA, NULL); + props_.viscosity(1, &perf_p, comp_frac, &c, &mu[0], NULL); + ASSERT(std::fabs(std::accumulate(comp_frac, comp_frac + np, 0.0) - 1.0) < 1e-6); + props_.relperm (1, comp_frac, &c, wpM , NULL); + for (int phase = 0; phase < np; ++phase) { + wpM[phase] /= mu[phase]; + } + } } } } diff --git a/opm/core/simulator/WellState.hpp b/opm/core/simulator/WellState.hpp index 59d9f0602..7ddc886e2 100644 --- a/opm/core/simulator/WellState.hpp +++ b/opm/core/simulator/WellState.hpp @@ -37,10 +37,17 @@ namespace Opm if (wells) { const int nw = wells->number_of_wells; bhp_.resize(nw); - // Initialize bhp to be pressure in first perforation cell. + // Initialize bhp to be target pressure + // if bhp-controlled well, otherwise set + // to pressure in first perforation cell. for (int w = 0; w < nw; ++w) { - const int cell = wells->well_cells[wells->well_connpos[w]]; - bhp_[w] = state.pressure()[cell]; + const WellControls* ctrl = wells->ctrls[w]; + if (ctrl->type[ctrl->current] == BHP) { + bhp_[w] = ctrl->target[ctrl->current]; + } else { + const int cell = wells->well_cells[wells->well_connpos[w]]; + bhp_[w] = state.pressure()[cell]; + } } perfrates_.resize(wells->well_connpos[nw]); } From 67b5f007fde0ae8b5351e5e4cb4e9dacf4cda32d Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 10:40:36 +0200 Subject: [PATCH 03/10] Made initialization from SWAT/SGAS etc. more robust and general. --- opm/core/utility/initState_impl.hpp | 66 ++++++++++++++++++++++------- 1 file changed, 51 insertions(+), 15 deletions(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index 7d01ade0e..c139bf0e0 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(): Uuser 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. From 65447604ae680821eb960dbeb068d5dc972406de Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Fri, 14 Sep 2012 20:56:08 +0200 Subject: [PATCH 04/10] Typo fix. --- opm/core/utility/initState_impl.hpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/opm/core/utility/initState_impl.hpp b/opm/core/utility/initState_impl.hpp index c139bf0e0..575e18a41 100644 --- a/opm/core/utility/initState_impl.hpp +++ b/opm/core/utility/initState_impl.hpp @@ -515,7 +515,7 @@ namespace Opm const int num_phases = props.numPhases(); const PhaseUsage pu = phaseUsageFromDeck(deck); if (num_phases != pu.num_phases) { - THROW("initStateFromDeck(): Uuser specified property object with " << num_phases << " phases, " + THROW("initStateFromDeck(): user specified property object with " << num_phases << " phases, " "found " << pu.num_phases << " phases in deck."); } state.init(grid, num_phases); From 2017481a58d650f1b2b182a543848879eda60f32 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Sep 2012 14:33:57 +0200 Subject: [PATCH 05/10] Improve diagnostic output if crossflow is detected. --- opm/core/utility/miscUtilities.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) 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); From 50a23c0f5d7e0d2a254585b84b55f72ad2ac400f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Sep 2012 14:34:33 +0200 Subject: [PATCH 06/10] Support shut wells in compressible tpfa solver. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 60 ++++++++++++++++------ 1 file changed, 44 insertions(+), 16 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index eb8dc8f34..8eecbfed6 100644 --- a/opm/core/pressure/tpfa/cfs_tpfa_residual.c +++ b/opm/core/pressure/tpfa/cfs_tpfa_residual.c @@ -741,6 +741,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 +838,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 +879,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 +901,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 +914,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, From e3388575d620f05946405e4b5a4152cff7d9c7bc Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Atgeirr=20Fl=C3=B8=20Rasmussen?= Date: Thu, 20 Sep 2012 14:35:03 +0200 Subject: [PATCH 07/10] Fix treatment of WELOPEN keyword. Now you can actually shut and open wells with WELOPEN. The following caveats apply: - this may interact improperly with group controls, - dynamic usage of WCONINJE/WCONPROD should not be mixed with WELOPEN. --- opm/core/wells/WellsManager.cpp | 51 ++++++++++++++++----------------- 1 file changed, 25 insertions(+), 26 deletions(-) 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. From bdcd5236bd9d592410d3caa480236de168226bbf Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 20 Sep 2012 15:03:38 +0200 Subject: [PATCH 08/10] Don't crash on models without wells. The user will legitimately want to run models that do not specify wells (e.g., using boundary conditions). While we do not yet fully support that configuration (no wells), we absolutely must not crash by dereferencing null pointers or generating pointers into ::empty() std::vector<>s. This commit installs the required guards needed to avoid said failure mode. --- opm/core/pressure/CompressibleTpfa.cpp | 25 ++++++++++++++----------- 1 file changed, 14 insertions(+), 11 deletions(-) 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], From b137eb0b65dfaecf01db6bef9f2ba5cbc84b28b4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?B=C3=A5rd=20Skaflestad?= Date: Thu, 20 Sep 2012 15:48:48 +0200 Subject: [PATCH 09/10] Interpret ``wells != 0 && wells->W == 0'' as ``no wells''. The CompressibleTpfa class always passes a non-null `forces->wells' object to the constructor, assembly, and reconstruction routines but uses ``forces->wells->W == 0'' to signify a simulation model without wells. This is, arguably, an error in the CompressibleTpfa class but one that does not require a lot of work to support in the cfs_tpfa_residual module. Insert the extra tests in an effort to honour the ``liberal in what you accept, strict in what you produce'' principle. --- opm/core/pressure/tpfa/cfs_tpfa_residual.c | 24 +++++++++++++--------- 1 file changed, 14 insertions(+), 10 deletions(-) diff --git a/opm/core/pressure/tpfa/cfs_tpfa_residual.c b/opm/core/pressure/tpfa/cfs_tpfa_residual.c index eb8dc8f34..023e3c2d1 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; @@ -1127,8 +1129,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 +1195,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 +1300,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); } From ab21d44c9ab40a25db5041b892cf0c61de7197e0 Mon Sep 17 00:00:00 2001 From: Roland Kaufmann Date: Fri, 21 Sep 2012 15:06:47 +0200 Subject: [PATCH 10/10] Disable warning for using DUNE's FieldVector::size In DUNE 2.2 FieldVector::size changed from being a member to being a method. A compatibility warning is issued if you include the relevant headers. This warning can be silenced for DUNE modules by using passing the option --enable-fieldvector-size-is-method to ./configure. This patch effectively does the same, but through a macro definition. --- opm/core/linalg/LinearSolverIstl.cpp | 4 ++++ 1 file changed, 4 insertions(+) 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