Preparations for Recording Dynamic Well State

Coalesce blocks with same conditions, split long lines, and apply
'const' where appropriate.  While here, also tighten the "rate = 0"
criterion to include denormalised numbers.
This commit is contained in:
Bård Skaflestad 2021-03-01 15:17:52 +01:00
parent cd7c4d9927
commit 65edfb702c
4 changed files with 163 additions and 114 deletions

View File

@ -268,11 +268,12 @@ namespace Opm {
/// Returns true if the well was actually found and shut.
bool forceShutWellByNameIfPredictionMode(const std::string& wellname, const double simulation_time);
void updateEclWell(int timeStepIdx, int well_index);
void updateEclWell(int timeStepIdx, const std::string& wname);
void updateEclWell(const int timeStepIdx, const int well_index);
void updateEclWell(const int timeStepIdx, const std::string& wname);
bool hasWell(const std::string& wname);
double wellPI(int well_index) const;
double wellPI(const int well_index) const;
double wellPI(const std::string& well_name) const;
protected:
Simulator& ebosSimulator_;
@ -345,7 +346,7 @@ namespace Opm {
WellTestState wellTestState_;
std::unique_ptr<GuideRate> guideRate_;
std::map<std::string, double> node_pressures_; // Storing network pressures for output.
std::map<std::string, double> node_pressures_{}; // Storing network pressures for output.
// used to better efficiency of calcuation
mutable BVector scaleAddRes_;

View File

@ -298,8 +298,8 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
beginTimeStep() {
beginTimeStep()
{
updatePerforationIntensiveQuantities();
Opm::DeferredLogger local_deferredLogger;
@ -425,7 +425,10 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx, const double simulationTime, Opm::DeferredLogger& deferred_logger) {
BlackoilWellModel<TypeTag>::wellTesting(const int timeStepIdx,
const double simulationTime,
Opm::DeferredLogger& deferred_logger)
{
const auto& wtest_config = schedule()[timeStepIdx].wtest_config();
if (wtest_config.size() != 0) { // there is a WTEST request
@ -435,7 +438,9 @@ namespace Opm {
std::vector< Scalar > B_avg(numComponents(), Scalar());
computeAverageFormationFactor(B_avg);
const auto& wellsForTesting = wellTestState_.updateWells(wtest_config, wells_ecl_, simulationTime);
const auto wellsForTesting = wellTestState_
.updateWells(wtest_config, wells_ecl_, simulationTime);
for (const auto& testWell : wellsForTesting) {
const std::string& well_name = testWell.first;
@ -444,9 +449,12 @@ namespace Opm {
// some preparation before the well can be used
well->init(&phase_usage_, depth_, gravity_, local_num_cells_, B_avg);
const Well& wellEcl = schedule().getWell(well_name, timeStepIdx);
double well_efficiency_factor = wellEcl.getEfficiencyFactor();
WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx), schedule(), timeStepIdx, well_efficiency_factor);
WellGroupHelpers::accumulateGroupEfficiencyFactor(schedule().getGroup(wellEcl.groupName(), timeStepIdx),
schedule(), timeStepIdx, well_efficiency_factor);
well->setWellEfficiencyFactor(well_efficiency_factor);
well->setVFPProperties(vfp_properties_.get());
well->setGuideRate(guideRate_.get());
@ -459,11 +467,16 @@ namespace Opm {
}
}
// called at the end of a report step
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
endReportStep() {
endReportStep()
{
// Clear the communication data structures for above values.
for (auto&& pinfo : local_parallel_well_info_)
{
@ -471,17 +484,26 @@ namespace Opm {
}
}
// called at the end of a report step
template<typename TypeTag>
const SimulatorReportSingle&
BlackoilWellModel<TypeTag>::
lastReport() const {return last_report_; }
// called at the end of a time step
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
timeStepSucceeded(const double& simulationTime, const double dt) {
timeStepSucceeded(const double& simulationTime, const double dt)
{
// time step is finished and we are not any more at the beginning of an report step
report_step_starts_ = false;
@ -735,12 +757,13 @@ namespace Opm {
if (nw > 0) {
well_container.reserve(nw);
for (int w = 0; w < nw; ++w) {
const Well& well_ecl = wells_ecl_[w];
const std::string& well_name = well_ecl.name();
// A new WCON keywords can re-open a well that was closed/shut due to Physical limit
if ( wellTestState_.hasWellClosed(well_name)) {
if (this->wellTestState_.hasWellClosed(well_name)) {
// TODO: more checking here, to make sure this standard more specific and complete
// maybe there is some WCON keywords will not open the well
if (well_state_.effectiveEventsOccurred(w)) {
@ -760,7 +783,8 @@ namespace Opm {
// something like wellTestState_.hasWell(well_name)?
bool wellIsStopped = false;
if (wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::ECONOMIC) ||
wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::PHYSICAL) ) {
wellTestState_.hasWellClosed(well_name, WellTestConfig::Reason::PHYSICAL))
{
if (well_ecl.getAutomaticShutIn()) {
// shut wells are not added to the well container
well_state_.shutWell(w);
@ -783,28 +807,41 @@ namespace Opm {
// (prediction type) rate control is zero, then it is effectively shut.
if (!well_ecl.getAllowCrossFlow() && well_ecl.isProducer() && well_ecl.predictionMode()) {
const auto& summaryState = ebosSimulator_.vanguard().summaryState();
auto prod_controls = well_ecl.productionControls(summaryState);
const auto prod_controls = well_ecl.productionControls(summaryState);
auto is_zero = [](const double x)
{
return std::isfinite(x) && !std::isnormal(x);
};
bool zero_rate_control = false;
switch (prod_controls.cmode) {
case Well::ProducerCMode::ORAT:
zero_rate_control = (prod_controls.oil_rate == 0.0);
zero_rate_control = is_zero(prod_controls.oil_rate);
break;
case Well::ProducerCMode::WRAT:
zero_rate_control = (prod_controls.water_rate == 0.0);
zero_rate_control = is_zero(prod_controls.water_rate);
break;
case Well::ProducerCMode::GRAT:
zero_rate_control = (prod_controls.gas_rate == 0.0);
zero_rate_control = is_zero(prod_controls.gas_rate);
break;
case Well::ProducerCMode::LRAT:
zero_rate_control = (prod_controls.liquid_rate == 0.0);
zero_rate_control = is_zero(prod_controls.liquid_rate);
break;
case Well::ProducerCMode::RESV:
zero_rate_control = (prod_controls.resv_rate == 0.0);
zero_rate_control = is_zero(prod_controls.resv_rate);
break;
default:
// Might still have zero rate controls, but is pressure controlled.
zero_rate_control = false;
break;
}
if (zero_rate_control) {
// Treat as shut, do not add to container.
local_deferredLogger.info(" Well shut due to zero rate control and disallowing crossflow: " + well_ecl.name());
@ -867,10 +904,11 @@ namespace Opm {
const auto& perf_data = this->well_perf_data_[wellID];
// Cater for case where local part might have no perforations.
const int pvtreg = perf_data.empty() ?
0 : pvt_region_idx_[perf_data.front().cell_index];
const auto pvtreg = perf_data.empty()
? 0 : pvt_region_idx_[perf_data.front().cell_index];
const auto& parallel_well_info = *local_parallel_well_info_[wellID];
auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
const auto global_pvtreg = parallel_well_info.broadcastFirstPerforationValue(pvtreg);
return std::make_unique<WellType>(this->wells_ecl_[wellID],
parallel_well_info,
@ -1426,11 +1464,13 @@ namespace Opm {
}
}
}
logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown, "computeWellPotentials() failed.", terminal_output_);
logAndCheckForExceptionsAndThrow(deferred_logger, exception_thrown,
"computeWellPotentials() failed.",
terminal_output_);
// Store it in the well state
well_state_.wellPotentials() = well_potentials;
}
@ -2536,7 +2576,7 @@ namespace Opm {
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateEclWell(int timeStepIdx, int well_index)
updateEclWell(const int timeStepIdx, const int well_index)
{
const auto& schedule = this->ebosSimulator_.vanguard().schedule();
const auto& wname = this->wells_ecl_[well_index].name();
@ -2556,23 +2596,36 @@ namespace Opm {
this->prod_index_calc_[well_index].reInit(well);
}
template<typename TypeTag>
void
BlackoilWellModel<TypeTag>::
updateEclWell(int timeStepIdx, const std::string& wname) {
auto well_iter = std::find_if( this->wells_ecl_.begin(), this->wells_ecl_.end(), [wname] (const auto& well) -> bool { return well.name() == wname;});
if (well_iter == this->wells_ecl_.end())
throw std::logic_error("Could not find well: " + wname);
updateEclWell(const int timeStepIdx, const std::string& wname)
{
auto well_iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(),
[&wname](const auto& well) -> bool
{
return well.name() == wname;
});
if (well_iter == this->wells_ecl_.end()) {
throw std::logic_error { "Could not find well: " + wname };
}
auto well_index = std::distance(this->wells_ecl_.begin(), well_iter);
this->updateEclWell(timeStepIdx, well_index);
}
template<typename TypeTag>
double
BlackoilWellModel<TypeTag>::
wellPI(int well_index) const
wellPI(const int well_index) const
{
const auto& pu = this->phase_usage_;
const auto np = this->numPhases();
@ -2603,14 +2656,24 @@ namespace Opm {
}
}
template<typename TypeTag>
double
BlackoilWellModel<TypeTag>::
wellPI(const std::string& well_name) const
{
auto well_iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(), [&well_name](const Well& well) { return well.name() == well_name; });
if (well_iter == this->wells_ecl_.end())
throw std::logic_error("Could not find well: " + well_name);
auto well_iter = std::find_if(this->wells_ecl_.begin(), this->wells_ecl_.end(),
[&well_name](const Well& well)
{
return well.name() == well_name;
});
if (well_iter == this->wells_ecl_.end()) {
throw std::logic_error { "Could not find well: " + well_name };
}
auto well_index = std::distance(this->wells_ecl_.begin(), well_iter);
return this->wellPI(well_index);
@ -2630,20 +2693,20 @@ namespace Opm {
auto hasWellPIEvent = [this, timeStepIdx](const int well_index) -> bool
{
return this->schedule()[timeStepIdx]
.wellgroup_events().hasEvent(this->wells_ecl_[well_index].name(),
return this->schedule()[timeStepIdx].wellgroup_events()
.hasEvent(this->wells_ecl_[well_index].name(),
ScheduleEvents::Events::WELL_PRODUCTIVITY_INDEX);
};
auto rescaleWellPI =
[this, timeStepIdx](const int well_index,
const double newWellPI) -> void
{
{
const auto& wname = this->wells_ecl_[well_index].name();
auto& schedule = this->ebosSimulator_.vanguard().schedule(); // Mutable
schedule.applyWellProdIndexScaling(wname, timeStepIdx, newWellPI);
}
this->updateEclWell(timeStepIdx, well_index);
};
@ -2672,8 +2735,7 @@ namespace Opm {
const auto nw = this->numLocalWells();
for (auto wellID = 0*nw; wellID < nw; ++wellID) {
if (hasWellPIEvent(wellID)) {
const auto newWellPI = this->wellPI(wellID);
rescaleWellPI(wellID, newWellPI);
rescaleWellPI(wellID, this->wellPI(wellID));
}
}

View File

@ -2295,11 +2295,11 @@ namespace Opm
// Sum with communication in case of distributed well.
const auto& comm = this->parallel_well_info_.communication();
if (comm.size() > 1)
{
if (comm.size() > 1) {
comm.sum(wellPI, np);
}
assert (static_cast<int>(subsetPerfID) == this->number_of_perforations_ &&
assert ((static_cast<int>(subsetPerfID) == this->number_of_perforations_) &&
"Internal logic error in processing connections for PI/II");
}

View File

@ -553,9 +553,11 @@ namespace Opm
data::Wells res = WellState::report(pu, globalCellIdxMap);
const int nw = this->numWells();
if( nw == 0 ) return res;
const int np = pu.num_phases;
if (nw == 0) {
return res;
}
const int np = pu.num_phases;
using rt = data::Rates::opt;
std::vector<rt> phs(np);
@ -571,58 +573,42 @@ namespace Opm
phs.at( pu.phase_pos[Gas] ) = rt::gas;
}
/* this is a reference or example on **how** to convert from
* WellState to something understood by opm-output. it is intended
* to be properly implemented and maintained as a part of
* simulators, as it relies on simulator internals, details and
* representations.
*/
// This is a reference or example on **how** to convert from
// WellState to something understood by opm-common's output
// layer. It is intended to be properly implemented and
// maintained as a part of simulators, as it relies on simulator
// internals, details and representations.
for (const auto& wt : this->wellMap()) {
const auto w = wt.second[ 0 ];
const auto& pwinfo = *parallel_well_info_[w];
if ((this->status_[w] == Well::Status::SHUT) || !pwinfo.isOwner())
{
continue;
}
auto& well = res.at(wt.first);
//well.injectionControl = static_cast<int>(this->currentInjectionControls()[ w ]);
//well.productionControl = static_cast<int>(this->currentProductionControls()[ w ]);
const int well_rate_index = w * pu.num_phases;
if (pu.phase_used[Water]) {
well.rates.set( rt::reservoir_water, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Water]] );
const auto i = well_rate_index + pu.phase_pos[Water];
well.rates.set(rt::reservoir_water, this->well_reservoir_rates_[i]);
well.rates.set(rt::productivity_index_water, this->productivity_index_[i]);
well.rates.set(rt::well_potential_water, this->well_potentials_[i]);
}
if (pu.phase_used[Oil]) {
well.rates.set( rt::reservoir_oil, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Oil]] );
const auto i = well_rate_index + pu.phase_pos[Oil];
well.rates.set(rt::reservoir_oil, this->well_reservoir_rates_[i]);
well.rates.set(rt::productivity_index_oil, this->productivity_index_[i]);
well.rates.set(rt::well_potential_oil, this->well_potentials_[i]);
}
if (pu.phase_used[Gas]) {
well.rates.set( rt::reservoir_gas, this->well_reservoir_rates_[well_rate_index + pu.phase_pos[Gas]] );
}
if ( pu.phase_used[Water] ) {
well.rates.set( rt::productivity_index_water, this->productivity_index_[well_rate_index + pu.phase_pos[Water]] );
}
if ( pu.phase_used[Oil] ) {
well.rates.set( rt::productivity_index_oil, this->productivity_index_[well_rate_index + pu.phase_pos[Oil]] );
}
if ( pu.phase_used[Gas] ) {
well.rates.set( rt::productivity_index_gas, this->productivity_index_[well_rate_index + pu.phase_pos[Gas]] );
}
if ( pu.phase_used[Water] ) {
well.rates.set( rt::well_potential_water, this->well_potentials_[well_rate_index + pu.phase_pos[Water]] );
}
if ( pu.phase_used[Oil] ) {
well.rates.set( rt::well_potential_oil, this->well_potentials_[well_rate_index + pu.phase_pos[Oil]] );
}
if ( pu.phase_used[Gas] ) {
well.rates.set( rt::well_potential_gas, this->well_potentials_[well_rate_index + pu.phase_pos[Gas]] );
const auto i = well_rate_index + pu.phase_pos[Gas];
well.rates.set(rt::reservoir_gas, this->well_reservoir_rates_[i]);
well.rates.set(rt::productivity_index_gas, this->productivity_index_[i]);
well.rates.set(rt::well_potential_gas, this->well_potentials_[i]);
}
if (pu.has_solvent || pu.has_zFraction) {