diff --git a/opm/simulators/wells/BlackoilWellModelGeneric.cpp b/opm/simulators/wells/BlackoilWellModelGeneric.cpp index 0c2496f8f..bfb480865 100644 --- a/opm/simulators/wells/BlackoilWellModelGeneric.cpp +++ b/opm/simulators/wells/BlackoilWellModelGeneric.cpp @@ -1239,7 +1239,7 @@ updateEclWells(const int timeStepIdx, } auto& ws = this->wellState().well(well_index); - this->wellState().updateStatus(well_index, well.getStatus()); + ws.updateStatus( well.getStatus() ); ws.reset_connection_factors(pd); this->prod_index_calc_[well_index].reInit(well); } diff --git a/opm/simulators/wells/PerfData.cpp b/opm/simulators/wells/PerfData.cpp index 00454d82f..955615ed7 100644 --- a/opm/simulators/wells/PerfData.cpp +++ b/opm/simulators/wells/PerfData.cpp @@ -24,8 +24,9 @@ namespace Opm { -PerfData::PerfData(std::size_t num_perf, bool injector_, std::size_t num_phases) +PerfData::PerfData(std::size_t num_perf, double pressure_first_connection_, bool injector_, std::size_t num_phases) : injector(injector_) + , pressure_first_connection(pressure_first_connection_) , pressure(num_perf) , rates(num_perf) , phase_rates(num_perf * num_phases) @@ -61,6 +62,7 @@ bool PerfData::try_assign(const PerfData& other) { if (this->injector != other.injector) return false; + this->pressure_first_connection = other.pressure_first_connection; this->pressure = other.pressure; this->rates = other.rates; this->phase_rates = other.phase_rates; diff --git a/opm/simulators/wells/PerfData.hpp b/opm/simulators/wells/PerfData.hpp index bd8920dd9..421593ce8 100644 --- a/opm/simulators/wells/PerfData.hpp +++ b/opm/simulators/wells/PerfData.hpp @@ -32,12 +32,13 @@ private: bool injector; public: - PerfData(std::size_t num_perf, bool injector_, std::size_t num_phases); + PerfData(std::size_t num_perf, double pressure_first_connection_, bool injector_, std::size_t num_phases); std::size_t size() const; bool empty() const; bool try_assign(const PerfData& other); + double pressure_first_connection; std::vector pressure; std::vector rates; std::vector phase_rates; diff --git a/opm/simulators/wells/SingleWellState.cpp b/opm/simulators/wells/SingleWellState.cpp index 437cd5f9e..b6b9278fc 100644 --- a/opm/simulators/wells/SingleWellState.cpp +++ b/opm/simulators/wells/SingleWellState.cpp @@ -22,16 +22,23 @@ namespace Opm { -SingleWellState::SingleWellState(const std::string& name_, const ParallelWellInfo& pinfo, bool is_producer, const std::vector& perf_input, std::size_t num_phases, double temp) +SingleWellState::SingleWellState(const std::string& name_, + const ParallelWellInfo& pinfo, + bool is_producer, + double pressure_first_connection, + const std::vector& perf_input, + const PhaseUsage& pu_, + double temp) : name(name_) , parallel_info(pinfo) , producer(is_producer) + , pu(pu_) , temperature(temp) - , well_potentials(num_phases) - , productivity_index(num_phases) - , surface_rates(num_phases) - , reservoir_rates(num_phases) - , perf_data(perf_input.size(), !is_producer, num_phases) + , well_potentials(pu_.num_phases) + , productivity_index(pu_.num_phases) + , surface_rates(pu_.num_phases) + , reservoir_rates(pu_.num_phases) + , perf_data(perf_input.size(), pressure_first_connection, !is_producer, pu_.num_phases) { for (std::size_t perf = 0; perf < perf_input.size(); perf++) { this->perf_data.cell_index[perf] = perf_input[perf].cell_index; @@ -143,7 +150,110 @@ double SingleWellState::sum_solvent_rates() const { return this->sum_connection_rates(this->perf_data.solvent_rates); } + +void SingleWellState::update_producer_targets(const Well& ecl_well, const SummaryState& st) { + const double bhp_safety_factor = 0.99; + const auto& prod_controls = ecl_well.productionControls(st); + if (prod_controls.hasControl(Well::ProducerCMode::THP)) + this->thp = prod_controls.thp_limit; + + + auto cmode_is_bhp = (prod_controls.cmode == Well::ProducerCMode::BHP); + auto bhp_limit = prod_controls.bhp_limit; + + if (ecl_well.getStatus() == Well::Status::STOP) { + if (cmode_is_bhp) + this->bhp = bhp_limit; + else + this->bhp = this->perf_data.pressure_first_connection; + + return; + } + + if (prod_controls.cmode == Well::ProducerCMode::GRUP) { + this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor; + return; + } + + switch (prod_controls.cmode) { + case Well::ProducerCMode::ORAT: + assert(this->pu.phase_used[BlackoilPhases::Liquid]); + this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate; + break; + case Well::ProducerCMode::WRAT: + assert(this->pu.phase_used[BlackoilPhases::Aqua]); + this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate; + break; + case Well::ProducerCMode::GRAT: + assert(this->pu.phase_used[BlackoilPhases::Vapour]); + this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate; + break; + default: + // Keep zero init. + break; + } + + if (cmode_is_bhp) + this->bhp = bhp_limit; + else + this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor; + +} + +void SingleWellState::update_injector_targets(const Well& ecl_well, const SummaryState& st) { + const double bhp_safety_factor = 1.01; + const auto& inj_controls = ecl_well.injectionControls(st); + + if (inj_controls.hasControl(Well::InjectorCMode::THP)) + this->thp = inj_controls.thp_limit; + + auto cmode_is_bhp = (inj_controls.cmode == Well::InjectorCMode::BHP); + auto bhp_limit = inj_controls.bhp_limit; + + if (ecl_well.getStatus() == Well::Status::STOP) { + if (cmode_is_bhp) + this->bhp = bhp_limit; + else + this->bhp = this->perf_data.pressure_first_connection; + + return; + } + + + if (inj_controls.cmode == Well::InjectorCMode::GRUP) { + this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor; + return; + } + + if (inj_controls.cmode == Well::InjectorCMode::RATE) { + auto inj_surf_rate = inj_controls.surface_rate; + switch (inj_controls.injector_type) { + case InjectorType::WATER: + assert(pu.phase_used[BlackoilPhases::Aqua]); + this->surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate; + break; + case InjectorType::GAS: + assert(pu.phase_used[BlackoilPhases::Vapour]); + this->surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate; + break; + case InjectorType::OIL: + assert(pu.phase_used[BlackoilPhases::Liquid]); + this->surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate; + break; + case InjectorType::MULTI: + // Not currently handled, keep zero init. + break; + } + } + + if (cmode_is_bhp) + this->bhp = bhp_limit; + else + this->bhp = this->perf_data.pressure_first_connection * bhp_safety_factor; +} } + + diff --git a/opm/simulators/wells/SingleWellState.hpp b/opm/simulators/wells/SingleWellState.hpp index 5655b0832..ddf10bd17 100644 --- a/opm/simulators/wells/SingleWellState.hpp +++ b/opm/simulators/wells/SingleWellState.hpp @@ -23,11 +23,13 @@ #include #include -#include #include #include + +#include #include #include +#include namespace Opm { @@ -35,13 +37,20 @@ struct PerforationData; class SingleWellState { public: - SingleWellState(const std::string& name, const ParallelWellInfo& pinfo, bool is_producer, const std::vector& perf_input, std::size_t num_phases, double temp); + SingleWellState(const std::string& name, + const ParallelWellInfo& pinfo, + bool is_producer, + double presssure_first_connection, + const std::vector& perf_input, + const PhaseUsage& pu, + double temp); std::string name; std::reference_wrapper parallel_info; Well::Status status{Well::Status::OPEN}; bool producer; + PhaseUsage pu; double bhp{0}; double thp{0}; double temperature{0}; @@ -65,6 +74,8 @@ public: /// PerforationData::connection_transmissibility_factor actually /// used (overwrites existing internal values). void reset_connection_factors(const std::vector& new_perf_data); + void update_producer_targets(const Well& ecl_well, const SummaryState& st); + void update_injector_targets(const Well& ecl_well, const SummaryState& st); void updateStatus(Well::Status status); void init_timestep(const SingleWellState& other); void shut(); diff --git a/opm/simulators/wells/WellState.cpp b/opm/simulators/wells/WellState.cpp index a5b420ba5..feffa5983 100644 --- a/opm/simulators/wells/WellState.cpp +++ b/opm/simulators/wells/WellState.cpp @@ -54,17 +54,15 @@ void WellState::base_init(const std::vector& cellPressures, } } -void WellState::initSingleProducer(const std::vector& cellPressures, - const Well& well, - const std::vector& well_perf_data, +void WellState::initSingleProducer(const Well& well, const ParallelWellInfo& well_info, + double pressure_first_connection, + const std::vector& well_perf_data, const SummaryState& summary_state) { - const double bhp_safety_factor = 0.99; const auto& pu = this->phase_usage_; - const int np = pu.num_phases; const double temp = 273.15 + 15.56; - auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, true, well_perf_data, static_cast(np), temp}); + auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, true, pressure_first_connection, well_perf_data, pu, temp}); if ( ws.perf_data.empty()) return; @@ -72,67 +70,21 @@ void WellState::initSingleProducer(const std::vector& cellPressures, ws.status = Well::Status::OPEN; } - const auto& prod_controls = well.isProducer() ? well.productionControls(summary_state) : Well::ProductionControls(0); - if (prod_controls.hasControl(Well::ProducerCMode::THP)) - ws.thp = prod_controls.thp_limit; - - - auto cmode_is_bhp = (prod_controls.cmode == Well::ProducerCMode::BHP); - auto bhp_limit = prod_controls.bhp_limit; - auto pressure_first_connection = well_info.broadcastFirstPerforationValue(cellPressures[ws.perf_data.cell_index[0]]); - - if (well.getStatus() == Well::Status::STOP) { - if (cmode_is_bhp) - ws.bhp = bhp_limit; - else - ws.bhp = pressure_first_connection; - - return; - } - - if (prod_controls.cmode == Well::ProducerCMode::GRUP) { - ws.bhp = pressure_first_connection * bhp_safety_factor; - return; - } - - switch (prod_controls.cmode) { - case Well::ProducerCMode::ORAT: - assert(pu.phase_used[BlackoilPhases::Liquid]); - ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate; - break; - case Well::ProducerCMode::WRAT: - assert(pu.phase_used[BlackoilPhases::Aqua]); - ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate; - break; - case Well::ProducerCMode::GRAT: - assert(pu.phase_used[BlackoilPhases::Vapour]); - ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate; - break; - default: - // Keep zero init. - break; - } - - if (cmode_is_bhp) - ws.bhp = bhp_limit; - else - ws.bhp = pressure_first_connection * bhp_safety_factor; + ws.update_producer_targets(well, summary_state); } -void WellState::initSingleInjector(const std::vector& cellPressures, - const Well& well, - const std::vector& well_perf_data, +void WellState::initSingleInjector(const Well& well, const ParallelWellInfo& well_info, + double pressure_first_connection, + const std::vector& well_perf_data, const SummaryState& summary_state) { - const double bhp_safety_factor = 1.01; const auto& pu = this->phase_usage_; - const int np = pu.num_phases; const auto& inj_controls = well.injectionControls(summary_state); const double temp = inj_controls.temperature; - auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, well.isProducer(), well_perf_data, static_cast(np), temp}); + auto& ws = this->wells_.add(well.name(), SingleWellState{well.name(), well_info, false, pressure_first_connection, well_perf_data, pu, temp}); if ( ws.perf_data.empty()) return; @@ -140,53 +92,7 @@ void WellState::initSingleInjector(const std::vector& cellPressures, ws.status = Well::Status::OPEN; } - if (inj_controls.hasControl(Well::InjectorCMode::THP)) - ws.thp = inj_controls.thp_limit; - - auto cmode_is_bhp = (inj_controls.cmode == Well::InjectorCMode::BHP); - auto bhp_limit = inj_controls.bhp_limit; - auto pressure_first_connection = well_info.broadcastFirstPerforationValue(cellPressures[ws.perf_data.cell_index[0]]); - - if (well.getStatus() == Well::Status::STOP) { - if (cmode_is_bhp) - ws.bhp = bhp_limit; - else - ws.bhp = pressure_first_connection; - - return; - } - - - if (inj_controls.cmode == Well::InjectorCMode::GRUP) { - ws.bhp = pressure_first_connection * bhp_safety_factor; - return; - } - - if (inj_controls.cmode == Well::InjectorCMode::RATE) { - auto inj_surf_rate = inj_controls.surface_rate; - switch (inj_controls.injector_type) { - case InjectorType::WATER: - assert(pu.phase_used[BlackoilPhases::Aqua]); - ws.surface_rates[pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate; - break; - case InjectorType::GAS: - assert(pu.phase_used[BlackoilPhases::Vapour]); - ws.surface_rates[pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate; - break; - case InjectorType::OIL: - assert(pu.phase_used[BlackoilPhases::Liquid]); - ws.surface_rates[pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate; - break; - case InjectorType::MULTI: - // Not currently handled, keep zero init. - break; - } - } - - if (cmode_is_bhp) - ws.bhp = bhp_limit; - else - ws.bhp = pressure_first_connection * bhp_safety_factor; + ws.update_injector_targets(well, summary_state); } @@ -200,10 +106,15 @@ void WellState::initSingleWell(const std::vector& cellPressures, if (well.isInjector() == well.isProducer()) throw std::logic_error(fmt::format("Well must be either producer or injector - logic error for well: {}", well.name())); + double pressure_first_connection = -1; + if (!well_perf_data.empty()) + pressure_first_connection = cellPressures[well_perf_data[0].cell_index]; + pressure_first_connection = well_info.broadcastFirstPerforationValue(pressure_first_connection); + if (well.isInjector()) - this->initSingleInjector(cellPressures, well, well_perf_data, well_info, summary_state); + this->initSingleInjector(well, well_info, pressure_first_connection, well_perf_data, summary_state); else - this->initSingleProducer(cellPressures, well, well_perf_data, well_info, summary_state); + this->initSingleProducer(well, well_info, pressure_first_connection, well_perf_data, summary_state); } diff --git a/opm/simulators/wells/WellState.hpp b/opm/simulators/wells/WellState.hpp index 688e9e224..8594de1b0 100644 --- a/opm/simulators/wells/WellState.hpp +++ b/opm/simulators/wells/WellState.hpp @@ -296,16 +296,16 @@ private: const ParallelWellInfo& well_info, const SummaryState& summary_state); - void initSingleProducer(const std::vector& cellPressures, - const Well& well, - const std::vector& well_perf_data, + void initSingleProducer(const Well& well, const ParallelWellInfo& well_info, + double pressure_first_connection, + const std::vector& well_perf_data, const SummaryState& summary_state); - void initSingleInjector(const std::vector& cellPressures, - const Well& well, - const std::vector& well_perf_data, + void initSingleInjector(const Well& well, const ParallelWellInfo& well_info, + double pressure_first_connection, + const std::vector& well_perf_data, const SummaryState& summary_state); }; diff --git a/tests/test_wellstate.cpp b/tests/test_wellstate.cpp index 38e70d25e..a31739ade 100644 --- a/tests/test_wellstate.cpp +++ b/tests/test_wellstate.cpp @@ -546,10 +546,10 @@ BOOST_AUTO_TEST_CASE(TESTSegmentState2) { BOOST_AUTO_TEST_CASE(TESTPerfData) { - Opm::PerfData pd1(3, true, 3); - Opm::PerfData pd2(3, true, 3); - Opm::PerfData pd3(2, true, 3); - Opm::PerfData pd4(3, false, 3); + Opm::PerfData pd1(3, 100, true, 3); + Opm::PerfData pd2(3, 100, true, 3); + Opm::PerfData pd3(2, 100, true, 3); + Opm::PerfData pd4(3, 100, false, 3); for (std::size_t i = 0; i < 3; i++) { @@ -575,9 +575,14 @@ BOOST_AUTO_TEST_CASE(TESTPerfData) { BOOST_AUTO_TEST_CASE(TestSingleWellState) { Opm::ParallelWellInfo pinfo; std::vector connections = {{0,1,1,0},{1,1,1,1},{2,1,1,2}}; - Opm::SingleWellState ws1("W1", pinfo, true, connections, 3, 1); - Opm::SingleWellState ws2("W2", pinfo, true, connections, 3, 2); - Opm::SingleWellState ws3("W3", pinfo, false, connections, 3, 3); + Opm::PhaseUsage pu; + + // This is totally bonkers, but the pu needs a complete deck to initialize properly + pu.num_phases = 3; + + Opm::SingleWellState ws1("W1", pinfo, true, 100, connections, pu, 1); + Opm::SingleWellState ws2("W2", pinfo, true, 100, connections, pu, 2); + Opm::SingleWellState ws3("W3", pinfo, false, 100, connections, pu, 3); ws1.bhp = 100; ws1.thp = 200;