Add well potentials to SingleWellState

This commit is contained in:
Joakim Hove 2021-08-05 10:57:15 +02:00
parent e84eaa3179
commit 2f504536f4
10 changed files with 32 additions and 42 deletions

View File

@ -1326,8 +1326,9 @@ namespace Opm {
// Store it in the well state
// potentials is resized and set to zero in the beginning of well->ComputeWellPotentials
// and updated only if sucessfull. i.e. the potentials are zero for exceptions
auto& ws = this->wellState().well(well->indexOfWell());
for (int p = 0; p < np; ++p) {
this->wellState().wellPotentials(well->indexOfWell())[p] = std::abs(potentials[p]);
ws.well_potentials[p] = std::abs(potentials[p]);
}
}

View File

@ -367,8 +367,7 @@ namespace Opm
const int np = number_of_phases_;
const double sign = well_copy.well_ecl_.isInjector() ? 1.0 : -1.0;
for (int phase = 0; phase < np; ++phase){
well_state_copy.wellRates(well_copy.index_of_well_)[phase]
= sign * well_state_copy.wellPotentials(well_copy.index_of_well_)[phase];
well_state_copy.wellRates(well_copy.index_of_well_)[phase] = sign * ws.well_potentials[phase];
}
well_copy.scaleSegmentRatesWithWellRates(well_state_copy);

View File

@ -21,9 +21,10 @@
namespace Opm {
SingleWellState::SingleWellState(bool is_producer, double temp)
SingleWellState::SingleWellState(bool is_producer, std::size_t num_phases, double temp)
: producer(is_producer)
, temperature(temp)
, well_potentials(num_phases)
{}

View File

@ -20,6 +20,8 @@
#ifndef OPM_SINGLE_WELL_STATE_HEADER_INCLUDED
#define OPM_SINGLE_WELL_STATE_HEADER_INCLUDED
#include <vector>
#include <opm/parser/eclipse/EclipseState/Schedule/Well/Well.hpp>
#include <opm/parser/eclipse/EclipseState/Schedule/Events.hpp>
@ -27,7 +29,7 @@ namespace Opm {
class SingleWellState {
public:
SingleWellState(bool is_producer, double temp);
SingleWellState(bool is_producer, std::size_t num_phases, double temp);
Well::Status status{Well::Status::OPEN};
bool producer;
@ -36,6 +38,7 @@ public:
double temperature{};
double dissolved_gas_rate{0};
double vaporized_oil_rate{0};
std::vector<double> well_potentials;
Events events;
Well::InjectorCMode injection_cmode{Well::InjectorCMode::CMODE_UNDEFINED};
Well::ProducerCMode production_cmode{Well::ProducerCMode::CMODE_UNDEFINED};

View File

@ -1383,9 +1383,10 @@ namespace WellGroupHelpers
continue;
}
const auto& ws = wellState.well(well_index);
// add contribution from wells unconditionally
for (int phase = 0; phase < np; phase++) {
pot[phase] += wefac * wellState.wellPotentials(well_index)[phase];
pot[phase] += wefac * ws.well_potentials[phase];
}
}
@ -1431,8 +1432,8 @@ namespace WellGroupHelpers
{
// the well is found and owned
int well_index = it->second[0];
const auto wpot = wellState.wellPotentials(well_index);
const auto& ws = wellState.well(well_index);
const auto& wpot = ws.well_potentials;
if (pu.phase_used[BlackoilPhases::Liquid] > 0)
oilpot = wpot[pu.phase_pos[BlackoilPhases::Liquid]];

View File

@ -693,9 +693,10 @@ updateWellTestStateEconomic(const WellState& well_state,
bool rate_limit_violated = false;
const auto& quantity_limit = econ_production_limits.quantityLimit();
const auto& ws = well_state.well(this->index_of_well_);
if (econ_production_limits.onAnyRateLimit()) {
if (quantity_limit == WellEconProductionLimits::QuantityLimit::POTN)
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellPotentials(index_of_well_).data(), deferred_logger);
rate_limit_violated = checkRateEconLimits(econ_production_limits, ws.well_potentials.data(), deferred_logger);
else {
rate_limit_violated = checkRateEconLimits(econ_production_limits, well_state.wellRates(index_of_well_).data(), deferred_logger);
}

View File

@ -246,6 +246,7 @@ namespace Opm
deferred_logger.info(" well " + this->name() + " is being tested for economic limits");
WellState well_state_copy = well_state;
auto& ws = well_state_copy.well(this->indexOfWell());
updateWellStateWithTarget(simulator, group_state, well_state_copy, deferred_logger);
calculateExplicitQuantities(simulator, well_state_copy, deferred_logger);
@ -270,7 +271,7 @@ namespace Opm
}
const int np = well_state_copy.numPhases();
for (int p = 0; p < np; ++p) {
well_state_copy.wellPotentials(this->indexOfWell())[p] = std::abs(potentials[p]);
ws.well_potentials[p] = std::abs(potentials[p]);
}
this->updateWellTestState(well_state_copy, simulation_time, /*writeMessageToOPMLog=*/ false, welltest_state_temp, deferred_logger);
this->closeCompletions(welltest_state_temp);
@ -470,6 +471,7 @@ namespace Opm
// create a copy of the well_state to use. If the operability checking is sucessful, we use this one
// to replace the original one
WellState well_state_copy = well_state;
auto& ws = well_state_copy.well(this->indexOfWell());
// TODO: well state for this well is kind of all zero status
// we should be able to provide a better initialization
@ -510,7 +512,7 @@ namespace Opm
}
const int np = well_state_copy.numPhases();
for (int p = 0; p < np; ++p) {
well_state_copy.wellPotentials(this->indexOfWell())[p] = std::abs(potentials[p]);
ws.well_potentials[p] = std::abs(potentials[p]);
}
well_state = well_state_copy;
} else {
@ -673,7 +675,7 @@ namespace Opm
double total_rate = std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = well_state.wellPotentials(well_index)[p];
well_state.wellRates(well_index)[p] = ws.well_potentials[p];
}
}
break;
@ -690,7 +692,7 @@ namespace Opm
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = well_state.wellPotentials(well_index)[p];
well_state.wellRates(well_index)[p] = ws.well_potentials[p];
}
}
break;
@ -876,7 +878,7 @@ namespace Opm
// using the well potentials
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = -well_state.wellPotentials(well_index)[p];
well_state.wellRates(well_index)[p] = -ws.well_potentials[p];
}
}
break;
@ -896,7 +898,7 @@ namespace Opm
double total_rate = -std::accumulate(rates.begin(), rates.end(), 0.0);
if (total_rate <= 0.0){
for (int p = 0; p<np; ++p) {
well_state.wellRates(well_index)[p] = -well_state.wellPotentials(well_index)[p];
well_state.wellRates(well_index)[p] = -ws.well_potentials[p];
}
}
break;
@ -939,14 +941,15 @@ namespace Opm
{
const int np = this->number_of_phases_;
std::vector<double> scaling_factor(np);
const auto& ws = well_state.well(this->index_of_well_);
double total_potentials = 0.0;
for (int p = 0; p<np; ++p) {
total_potentials += well_state.wellPotentials(this->index_of_well_)[p];
total_potentials += ws.well_potentials[p];
}
if (total_potentials > 0) {
for (int p = 0; p<np; ++p) {
scaling_factor[p] = well_state.wellPotentials(this->index_of_well_)[p] / total_potentials;
scaling_factor[p] = ws.well_potentials[p] / total_potentials;
}
return scaling_factor;
}

View File

@ -44,7 +44,6 @@ void WellState::base_init(const std::vector<double>& cellPressures,
this->parallel_well_info_.clear();
this->wellrates_.clear();
this->segment_state.clear();
this->well_potentials_.clear();
this->productivity_index_.clear();
this->wells_.clear();
{
@ -90,10 +89,9 @@ void WellState::initSingleWell(const std::vector<double>& cellPressures,
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.isProducer(), temp});
auto& ws = this->wells_.add(well.name(), SingleWellState{well.isProducer(), static_cast<std::size_t>(np), temp});
this->parallel_well_info_.add(well.name(), well_info);
this->wellrates_.add(well.name(), std::vector<double>(np, 0));
this->well_potentials_.add(well.name(), std::vector<double>(np, 0));
const int num_perf_this_well = well_info->communication().sum(well_perf_data.size());
this->segment_state.add(well.name(), SegmentState{});
this->perfdata.add(well.name(), PerfData{static_cast<std::size_t>(num_perf_this_well), well.isInjector(), this->phase_usage_});
@ -334,7 +332,6 @@ void WellState::init(const std::vector<double>& cellPressures,
auto it = prevState->wellMap().find(well.name());
if (it != end)
{
const int newIndex = w;
const int oldIndex = it->second[ 0 ];
const auto& prev_well = prevState->well(oldIndex);
@ -358,11 +355,7 @@ void WellState::init(const std::vector<double>& cellPressures,
wellRates(w) = prevState->wellRates(oldIndex);
wellReservoirRates(w) = prevState->wellReservoirRates(oldIndex);
// Well potentials
for (int p=0; p < np; p++) {
this->wellPotentials(newIndex)[p] = prevState->wellPotentials(oldIndex)[p];
}
new_well.well_potentials = prev_well.well_potentials;
// perfPhaseRates
const int num_perf_old_well = (*it).second[ 2 ];
@ -487,7 +480,7 @@ WellState::report(const int* globalCellIdxMap,
}
const auto& reservoir_rates = this->well_reservoir_rates_[well_index];
const auto& well_potentials = this->well_potentials_[well_index];
const auto& well_potentials = ws.well_potentials;
const auto& wpi = this->productivity_index_[well_index];
const auto& wv = this->wellRates(well_index);

View File

@ -186,14 +186,6 @@ public:
return this->productivity_index_[well_index];
}
std::vector<double>& wellPotentials(std::size_t well_index) {
return this->well_potentials_[well_index];
}
const std::vector<double>& wellPotentials(std::size_t well_index) const {
return this->well_potentials_[well_index];
}
template<class Comm>
void communicateGroupRates(const Comm& comm);
@ -360,10 +352,6 @@ private:
// Productivity Index
WellContainer<std::vector<double>> productivity_index_;
// Well potentials
WellContainer<std::vector<double>> well_potentials_;
data::Segment
reportSegmentResults(const PhaseUsage& pu,
const int well_id,

View File

@ -577,9 +577,9 @@ GAS
BOOST_AUTO_TEST_CASE(TestSingleWellState) {
Opm::SingleWellState ws1(true, 1);
Opm::SingleWellState ws2(true, 2);
Opm::SingleWellState ws3(false, 3);
Opm::SingleWellState ws1(true, 3, 1);
Opm::SingleWellState ws2(true, 3, 2);
Opm::SingleWellState ws3(false, 3, 3);
ws1.bhp = 100;
ws1.thp = 200;