Merge pull request #3659 from joakim-hove/refactor-wellstate-init

Refactor wellstate init
This commit is contained in:
Bård Skaflestad 2021-11-03 13:16:34 +01:00 committed by GitHub
commit 2971a3f051
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
8 changed files with 284 additions and 215 deletions

View File

@ -1472,9 +1472,10 @@ updateEclWells(const int timeStepIdx,
++pdIter;
}
}
auto& ws = this->wellState().well(well_index);
this->wellState().updateStatus(well_index, well.getStatus());
this->wellState().resetConnectionTransFactors(well_index, pd);
ws.updateStatus( well.getStatus() );
ws.reset_connection_factors(pd);
this->prod_index_calc_[well_index].reInit(well);
}
}
@ -2329,7 +2330,8 @@ runWellPIScaling(const int timeStepIdx,
++pdIter;
}
}
this->wellState().resetConnectionTransFactors(well_index, pd);
auto& ws = this->wellState().well(well_index);
ws.reset_connection_factors(pd);
this->prod_index_calc_[well_index].reInit(well);
};

View File

@ -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;

View File

@ -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<double> pressure;
std::vector<double> rates;
std::vector<double> phase_rates;

View File

@ -18,19 +18,35 @@
*/
#include <opm/simulators/wells/SingleWellState.hpp>
#include <opm/simulators/wells/PerforationData.hpp>
namespace Opm {
SingleWellState::SingleWellState(const ParallelWellInfo& pinfo, bool is_producer, std::size_t num_perf, std::size_t num_phases, double temp)
: parallel_info(pinfo)
SingleWellState::SingleWellState(const std::string& name_,
const ParallelWellInfo& pinfo,
bool is_producer,
double pressure_first_connection,
const std::vector<PerforationData>& 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(num_perf, !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;
this->perf_data.connection_transmissibility_factor[perf] = perf_input[perf].connection_transmissibility_factor;
this->perf_data.satnum_id[perf] = perf_input[perf].satnum_id;
this->perf_data.ecl_index[perf] = perf_input[perf].ecl_index;
}
}
void SingleWellState::init_timestep(const SingleWellState& other) {
@ -56,6 +72,9 @@ void SingleWellState::shut() {
std::fill(this->surface_rates.begin(), this->surface_rates.end(), 0);
std::fill(this->reservoir_rates.begin(), this->reservoir_rates.end(), 0);
std::fill(this->productivity_index.begin(), this->productivity_index.end(), 0);
auto& connpi = this->perf_data.prod_index;
connpi.assign(connpi.size(), 0);
}
void SingleWellState::stop() {
@ -67,6 +86,51 @@ void SingleWellState::open() {
this->status = Well::Status::OPEN;
}
void SingleWellState::updateStatus(Well::Status new_status) {
switch (new_status) {
case Well::Status::OPEN:
this->open();
break;
case Well::Status::SHUT:
this->shut();
break;
case Well::Status::STOP:
this->stop();
break;
default:
throw std::logic_error("Invalid well status");
}
}
void SingleWellState::reset_connection_factors(const std::vector<PerforationData>& new_perf_data) {
if (this->perf_data.size() != new_perf_data.size()) {
throw std::invalid_argument {
"Size mismatch for perforation data in well " + this->name
};
}
for (std::size_t conn_index = 0; conn_index < new_perf_data.size(); conn_index++) {
if (this->perf_data.cell_index[conn_index] != static_cast<std::size_t>(new_perf_data[conn_index].cell_index)) {
throw std::invalid_argument {
"Cell index mismatch in connection "
+ std::to_string(conn_index)
+ " of well "
+ this->name
};
}
if (this->perf_data.satnum_id[conn_index] != new_perf_data[conn_index].satnum_id) {
throw std::invalid_argument {
"Saturation function table mismatch in connection "
+ std::to_string(conn_index)
+ " of well "
+ this->name
};
}
this->perf_data.connection_transmissibility_factor[conn_index] = new_perf_data[conn_index].connection_transmissibility_factor;
}
}
double SingleWellState::sum_connection_rates(const std::vector<double>& connection_rates) const {
return this->parallel_info.get().sumPerfValues(connection_rates.begin(), connection_rates.end());
@ -86,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;
}
}

View File

@ -23,22 +23,34 @@
#include <functional>
#include <vector>
#include <opm/simulators/wells/SegmentState.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
#include <opm/simulators/wells/SegmentState.hpp>
#include <opm/simulators/wells/PerfData.hpp>
#include <opm/simulators/wells/ParallelWellInfo.hpp>
#include <opm/core/props/BlackoilPhases.hpp>
namespace Opm {
struct PerforationData;
class SingleWellState {
public:
SingleWellState(const ParallelWellInfo& pinfo, bool is_producer, std::size_t num_perf, 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<PerforationData>& perf_input,
const PhaseUsage& pu,
double temp);
std::string name;
std::reference_wrapper<const ParallelWellInfo> parallel_info;
Well::Status status{Well::Status::OPEN};
bool producer;
PhaseUsage pu;
double bhp{0};
double thp{0};
double temperature{0};
@ -55,6 +67,16 @@ public:
Well::ProducerCMode production_cmode{Well::ProducerCMode::CMODE_UNDEFINED};
/// Special purpose method to support dynamically rescaling a well's
/// CTFs through WELPI.
///
/// \param[in] new_perf_data New perforation data. Only
/// PerforationData::connection_transmissibility_factor actually
/// used (overwrites existing internal values).
void reset_connection_factors(const std::vector<PerforationData>& 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();
void stop();

View File

@ -18,6 +18,7 @@
along with OPM. If not, see <http://www.gnu.org/licenses/>.
*/
#include <fmt/format.h>
#include <config.h>
#include <opm/simulators/wells/WellState.hpp>
@ -53,7 +54,46 @@ void WellState::base_init(const std::vector<double>& cellPressures,
}
}
void WellState::initSingleProducer(const Well& well,
const ParallelWellInfo& well_info,
double pressure_first_connection,
const std::vector<PerforationData>& well_perf_data,
const SummaryState& summary_state) {
const auto& pu = this->phase_usage_;
const double temp = 273.15 + 15.56;
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;
if (well.getStatus() == Well::Status::OPEN) {
ws.status = Well::Status::OPEN;
}
ws.update_producer_targets(well, summary_state);
}
void WellState::initSingleInjector(const Well& well,
const ParallelWellInfo& well_info,
double pressure_first_connection,
const std::vector<PerforationData>& well_perf_data,
const SummaryState& summary_state) {
const auto& pu = this->phase_usage_;
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, false, pressure_first_connection, well_perf_data, pu, temp});
if ( ws.perf_data.empty())
return;
if (well.getStatus() == Well::Status::OPEN) {
ws.status = Well::Status::OPEN;
}
ws.update_injector_targets(well, summary_state);
}
@ -63,134 +103,18 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
const ParallelWellInfo& well_info,
const SummaryState& summary_state)
{
assert(well.isInjector() || well.isProducer());
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);
// Set default zero initial well rates.
// May be overwritten below.
const auto& pu = this->phase_usage_;
const int np = pu.num_phases;
double temp = well.isInjector() ? well.injectionControls(summary_state).temperature : 273.15 + 15.56;
auto& ws = this->wells_.add(well.name(), SingleWellState{well_info, well.isProducer(), well_perf_data.size(), static_cast<std::size_t>(np), temp});
if ( ws.perf_data.empty())
return;
const auto inj_controls = well.isInjector() ? well.injectionControls(summary_state) : Well::InjectionControls(0);
const auto prod_controls = well.isProducer() ? well.productionControls(summary_state) : Well::ProductionControls(0);
const bool is_bhp = well.isInjector() ? (inj_controls.cmode == Well::InjectorCMode::BHP)
: (prod_controls.cmode == Well::ProducerCMode::BHP);
const double bhp_limit = well.isInjector() ? inj_controls.bhp_limit : prod_controls.bhp_limit;
const bool is_grup = well.isInjector() ? (inj_controls.cmode == Well::InjectorCMode::GRUP)
: (prod_controls.cmode == Well::ProducerCMode::GRUP);
const double inj_surf_rate = well.isInjector() ? inj_controls.surface_rate : 0.0; // To avoid a "maybe-uninitialized" warning.
const double local_pressure = well_perf_data.empty() ?
0 : cellPressures[well_perf_data[0].cell_index];
const double global_pressure = well_info.broadcastFirstPerforationValue(local_pressure);
if (well.getStatus() == Well::Status::OPEN) {
ws.status = Well::Status::OPEN;
}
if (well.getStatus() == Well::Status::STOP) {
// Stopped well:
// 1. Rates: zero well rates.
// 2. Bhp: assign bhp equal to bhp control, if
// applicable, otherwise assign equal to
// first perforation cell pressure.
if (is_bhp) {
ws.bhp = bhp_limit;
} else {
ws.bhp = global_pressure;
}
} else if (is_grup) {
// Well under group control.
// 1. Rates: zero well rates.
// 2. Bhp: initialize bhp to be a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
const double safety_factor = well.isInjector() ? 1.01 : 0.99;
ws.bhp = safety_factor * global_pressure;
} else {
// Open well, under own control:
// 1. Rates: initialize well rates to match
// controls if type is ORAT/GRAT/WRAT
// (producer) or RATE (injector).
// Otherwise, we cannot set the correct
// value here and initialize to zero rate.
auto& rates = ws.surface_rates;
if (well.isInjector()) {
if (inj_controls.cmode == Well::InjectorCMode::RATE) {
switch (inj_controls.injector_type) {
case InjectorType::WATER:
assert(pu.phase_used[BlackoilPhases::Aqua]);
rates[pu.phase_pos[BlackoilPhases::Aqua]] = inj_surf_rate;
break;
case InjectorType::GAS:
assert(pu.phase_used[BlackoilPhases::Vapour]);
rates[pu.phase_pos[BlackoilPhases::Vapour]] = inj_surf_rate;
break;
case InjectorType::OIL:
assert(pu.phase_used[BlackoilPhases::Liquid]);
rates[pu.phase_pos[BlackoilPhases::Liquid]] = inj_surf_rate;
break;
case InjectorType::MULTI:
// Not currently handled, keep zero init.
break;
}
} else {
// Keep zero init.
}
} else {
assert(well.isProducer());
// Note negative rates for producing wells.
switch (prod_controls.cmode) {
case Well::ProducerCMode::ORAT:
assert(pu.phase_used[BlackoilPhases::Liquid]);
rates[pu.phase_pos[BlackoilPhases::Liquid]] = -prod_controls.oil_rate;
break;
case Well::ProducerCMode::WRAT:
assert(pu.phase_used[BlackoilPhases::Aqua]);
rates[pu.phase_pos[BlackoilPhases::Aqua]] = -prod_controls.water_rate;
break;
case Well::ProducerCMode::GRAT:
assert(pu.phase_used[BlackoilPhases::Vapour]);
rates[pu.phase_pos[BlackoilPhases::Vapour]] = -prod_controls.gas_rate;
break;
default:
// Keep zero init.
break;
}
}
// 2. Bhp: initialize bhp to be target pressure if
// bhp-controlled well, otherwise set to a
// little above or below (depending on if
// the well is an injector or producer)
// pressure in first perforation cell.
if (is_bhp) {
ws.bhp = bhp_limit;
} else {
const double safety_factor = well.isInjector() ? 1.01 : 0.99;
ws.bhp = safety_factor * global_pressure;
}
}
// 3. Thp: assign thp equal to thp target/limit, if such a limit exists,
// otherwise keep it zero.
const bool has_thp = well.isInjector() ? inj_controls.hasControl(Well::InjectorCMode::THP)
: prod_controls.hasControl(Well::ProducerCMode::THP);
const double thp_limit = well.isInjector() ? inj_controls.thp_limit : prod_controls.thp_limit;
if (has_thp) {
ws.thp = thp_limit;
}
if (well.isInjector())
this->initSingleInjector(well, well_info, pressure_first_connection, well_perf_data, summary_state);
else
this->initSingleProducer(well, well_info, pressure_first_connection, well_perf_data, summary_state);
}
void WellState::init(const std::vector<double>& cellPressures,
const Schedule& schedule,
const std::vector<Well>& wells_ecl,
@ -242,14 +166,8 @@ void WellState::init(const std::vector<double>& cellPressures,
auto& perf_data = ws.perf_data;
const int num_perf_this_well = perf_data.size();
const int global_num_perf_this_well = ecl_well.getConnections().num_open();
const auto& perf_input = well_perf_data[w];
for (int perf = 0; perf < num_perf_this_well; ++perf) {
perf_data.cell_index[perf] = perf_input[perf].cell_index;
perf_data.connection_transmissibility_factor[perf] = perf_input[perf].connection_transmissibility_factor;
perf_data.satnum_id[perf] = perf_input[perf].satnum_id;
perf_data.ecl_index[perf] = perf_input[perf].ecl_index;
if (wells_ecl[w].getStatus() == Well::Status::OPEN) {
for (int p = 0; p < this->numPhases(); ++p) {
perf_data.phase_rates[this->numPhases()*perf + p] = ws.surface_rates[p] / double(global_num_perf_this_well);
@ -751,27 +669,12 @@ void WellState::shutWell(int well_index)
{
auto& ws = this->well(well_index);
ws.shut();
auto& perf_data = ws.perf_data;
auto& connpi = perf_data.prod_index;
connpi.assign(connpi.size(), 0);
}
void WellState::updateStatus(int well_index, Well::Status status)
{
switch (status) {
case Well::Status::OPEN:
this->openWell(well_index);
break;
case Well::Status::SHUT:
this->shutWell(well_index);
break;
case Well::Status::STOP:
this->stopWell(well_index);
break;
default:
throw std::logic_error("Invalid well status");
}
auto& ws = this->well(well_index);
ws.updateStatus(status);
}
@ -911,40 +814,6 @@ void WellState::updateWellsDefaultALQ( const std::vector<Well>& wells_ecl )
}
}
void WellState::resetConnectionTransFactors(const int well_index,
const std::vector<PerforationData>& new_perf_data)
{
auto& ws = this->well(well_index);
auto& perf_data = ws.perf_data;
if (perf_data.size() != new_perf_data.size()) {
throw std::invalid_argument {
"Size mismatch for perforation data in well "
+ std::to_string(well_index)
};
}
for (std::size_t conn_index = 0; conn_index < new_perf_data.size(); conn_index++) {
if (perf_data.cell_index[conn_index] != static_cast<std::size_t>(new_perf_data[conn_index].cell_index)) {
throw std::invalid_argument {
"Cell index mismatch in connection "
+ std::to_string(conn_index)
+ " of well "
+ std::to_string(well_index)
};
}
if (perf_data.satnum_id[conn_index] != new_perf_data[conn_index].satnum_id) {
throw std::invalid_argument {
"Saturation function table mismatch in connection "
+ std::to_string(conn_index)
+ " of well "
+ std::to_string(well_index)
};
}
perf_data.connection_transmissibility_factor[conn_index] = new_perf_data[conn_index].connection_transmissibility_factor;
}
}
const ParallelWellInfo&
WellState::parallelWellInfo(std::size_t well_index) const
@ -956,3 +825,5 @@ WellState::parallelWellInfo(std::size_t well_index) const
template void WellState::updateGlobalIsGrup<ParallelWellInfo::Communication>(const ParallelWellInfo::Communication& comm);
template void WellState::communicateGroupRates<ParallelWellInfo::Communication>(const ParallelWellInfo::Communication& comm);
} // namespace Opm

View File

@ -197,18 +197,6 @@ public:
bool wellIsOwned(const std::string& wellName) const;
/// Special purpose method to support dynamically rescaling a well's
/// CTFs through WELPI.
///
/// \param[in] well_index Process-local linear index of single well.
/// Must be in the range 0..numWells()-1.
///
/// \param[in] well_perf_data New perforation data. Only
/// PerforationData::connection_transmissibility_factor actually
/// used (overwrites existing internal values).
void resetConnectionTransFactors(const int well_index,
const std::vector<PerforationData>& well_perf_data);
void updateStatus(int well_index, Well::Status status);
void openWell(int well_index);
@ -324,7 +312,17 @@ private:
const ParallelWellInfo& well_info,
const SummaryState& summary_state);
void initSingleProducer(const Well& well,
const ParallelWellInfo& well_info,
double pressure_first_connection,
const std::vector<PerforationData>& well_perf_data,
const SummaryState& summary_state);
void initSingleInjector(const Well& well,
const ParallelWellInfo& well_info,
double pressure_first_connection,
const std::vector<PerforationData>& well_perf_data,
const SummaryState& summary_state);
};

View File

@ -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++) {
@ -574,9 +574,15 @@ BOOST_AUTO_TEST_CASE(TESTPerfData) {
BOOST_AUTO_TEST_CASE(TestSingleWellState) {
Opm::ParallelWellInfo pinfo;
Opm::SingleWellState ws1(pinfo, true, 10, 3, 1);
Opm::SingleWellState ws2(pinfo, true, 10, 3, 2);
Opm::SingleWellState ws3(pinfo, false, 10, 3, 3);
std::vector<Opm::PerforationData> connections = {{0,1,1,0},{1,1,1,1},{2,1,1,2}};
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;